統計ER

統計ソフトRの使い方を中心に、統計解析方法の解説をするブログ。ありそうでなかなか見つからないサンプルサイズ計算などニッチな方法について紹介しています。

真の要因か見かけの要因かを見分ける回帰直線の差の検定

ブログランキングに参加しています。
まずはぽちぽちっとお願いします。
↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓
にほんブログ村 科学ブログ 数学へ

ある健康指標が地域差がありそうだ。こんなとき健康指標の平均値の差の検定がしたくなるもの。しかし、健康指標の多くは年齢に関係している。そして地域も年齢構成が異なることが多い。

本当に健康指標に差があるのか、それとも単に年齢構成の違いだけか。それを確認するのに回帰直線の差の検定が使える。手法としては、共分散分析を使うことになる。

統計ソフトRに組み込まれている関数を使った場合と使わない場合、両方でやってみた。

引用書籍

医学への統計学 13. 交絡因子の調整 13.1 共分散分析

最新第3版はこちら。

例題のデータ

ag は年齢、bp は血圧、.a はA地区、.b はB地区を表している。

ag.a <- c(44,36,57,39,32,41,46,41,40,46,41,62,44,30,44,52,49,48,40,48,33,46,44,35,38,40,52,43,47,55,48,25)
bp.a <- c(190,140,270,126,98,212,161,181,115,135,198,217,150,64,151,189,219,181,203,229,139,133,159,227,182,165,164,199,161,210,213,152)
ag.b <- c(45,58,66,62,64,50,54,66,50,35,61,57,43,68,48,66,44,59,59,51,52)
bp.b <- c(221,197,125,235,253,167,260,234,172,168,220,191,127,329,141,284,225,215,197,127,207)

血圧を比較すると

平均値の差の検定(t検定)で血圧を比較すると、統計学的有意に異なるという結果。等分散を仮定しても、仮定しなくても(Welchの方法)でも同じ結果。ちなみに標準偏差は10くらい違う。

> t.test(bp.a,bp.b,var=T)

        Two Sample t-test

data:  bp.a and bp.b
t = -2.3816, df = 51, p-value = 0.02101
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -58.269538  -4.965581
sample estimates:
mean of x mean of y 
 172.9062  204.5238 

> t.test(bp.a,bp.b)

        Welch Two Sample t-test

data:  bp.a and bp.b
t = -2.2786, df = 36.531, p-value = 0.02863
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -59.74434  -3.49078
sample estimates:
mean of x mean of y 
 172.9062  204.5238 

> c(sd(bp.a),sd(bp.b))
[1] 43.03364 53.17764

年齢と血圧の散布図

血圧など健康指標に関する調査は、必ず年齢との関係を見たほうがいい。一般的には加齢とともに健康状態が悪くなる。散布図を描いてみる。

plot(ag.a, bp.a, las=1, xlim=c(0,100), ylim=c(50,350),
     xlab="Age", ylab="Blood Pressure", pch=2)
points(ag.b, bp.b)
title("True or not? Difference between two areas")

A地区が三角、B地区が丸で示してある。ほとんど違いはないように見えるがどうだろうか?

f:id:toukeier:20190107191537p:plain

二つの地区の年齢と血圧の回帰直線を求めてみる

linear.reg <- function(x,y){
 xbar <- mean(x)
 ybar <- mean(y)
 n <- length(x)
 b <- sum((x-xbar)*(y-ybar))/sum((x-xbar)^2)
 a <- ybar-b*xbar
 yhat <- a+b*x
 R2 <- sum((yhat-ybar)^2)/sum((y-ybar)^2)
 VE <- 1/(n-2)*sum((y-(a+b*x))^2)
 rse <- sqrt(VE)
 SEb <- sqrt(VE/sum((x-xbar)^2))
 Tb <- b/SEb
 pb <- 2*pt(abs(Tb), n-2, lower.tail=F)
 bl <- b-qt(1-0.05/2, n-2)*SEb
 bu <- b+qt(1-0.05/2, n-2)*SEb
 SEa <- sqrt(VE*(1/n+xbar^2/sum((x-xbar)^2)))
 Ta <- a/SEa
 pa <- 2*pt(abs(Ta), n-2, lower.tail=F)
 al <- a-qt(1-0.05/2, n-2)*SEa
 au <- a+qt(1-0.05/2, n-2)*SEa
 list(
      c(Slope=b,"t_Slope"=Tb, "p_Slope"=pb, "LL_Slope"=bl, "UL_Slope"=bu,
      Intercept=a,"t_Intcp"=Ta, "p_Intcp"=pa, "LL_Intcp"=al, "UL_Intcp"=au,
      "Error Var."=VE,"Resid. SE"=rse,"R^2"=R2))
}

linreg.a <- unlist(linear.reg(x=ag.a, y=bp.a))
linreg.b <- unlist(linear.reg(x=ag.b, y=bp.b))

linreg.a
linreg.b

結果はこちら。Slopeが傾き、Interceptが切片。A地区の傾きが3.119、B地区の傾きが3.038、A地区の切片が37.8、B地区の切片が37.0となった。ほどんど同じと言ってもいいような気がする。

> linreg.a
        Slope       t_Slope       p_Slope      LL_Slope      UL_Slope 
 3.119131e+00  3.811036e+00  6.400747e-04  1.447639e+00  4.790623e+00 
    Intercept       t_Intcp       p_Intcp      LL_Intcp      UL_Intcp 
 3.780888e+01  1.049872e+00  3.021611e-01 -3.573914e+01  1.113569e+02 
   Error Var.     Resid. SE           R^2 
 1.289388e+03  3.590805e+01  3.262060e-01 
> linreg.b
        Slope       t_Slope       p_Slope      LL_Slope      UL_Slope 
   3.03780189    2.62659370    0.01661368    0.61710301    5.45850077 
    Intercept       t_Intcp       p_Intcp      LL_Intcp      UL_Intcp 
  37.01073387    0.57304659    0.57333336  -98.16912384  172.19059158 
   Error Var.     Resid. SE           R^2 
2183.76192669   46.73073856    0.26638079 

散布図に二つの回帰直線を描き入れてみる

abline(a=linreg.a[6], b=linreg.a[1],lty=2)
abline(a=linreg.b[6], b=linreg.b[1])

点線がA地区、実線がB地区。ほとんど重なっているのがわかる。

f:id:toukeier:20190107192009p:plain

要するに年齢の違いなのでは?

では、先ほどの血圧が地区間統計学的有意に異なっていたのはなぜか?年齢に違いはないのか?平均値の差の検定をやってみると、統計学的有意に異なるという結果だった。等分散を仮定しても、仮定しなくても結果は同様。

> t.test(ag.a, ag.b, var=T)

        Two Sample t-test

data:  ag.a and ag.b
t = -5.0438, df = 51, p-value = 6.139e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -16.539166  -7.121548
sample estimates:
mean of x mean of y 
 43.31250  55.14286 

> t.test(ag.a, ag.b)

        Welch Two Sample t-test

data:  ag.a and ag.b
t = -4.9007, df = 38.726, p-value = 1.738e-05
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -16.714272  -6.946443
sample estimates:
mean of x mean of y 
 43.31250  55.14286 

> c(sd(ag.a),sd(ag.b))
[1] 7.879895 9.034853

年齢と血圧の二地区の母回帰直線の傾きが同じかどうか、切片が同じかどうか

回帰直線の傾きが同じと言えるかどうか?同じと言えるなら、さらに切片はどうか?切片が同じであれば、平行でもなく、同一の回帰直線でよい。つまり、二つの回帰直線には差がなく、血圧の差は、地域差ではなく、年齢分布差ということになる。

reg.line.diff <- function(X1, Y1, X2, Y2){
 X <- c(X1,X2)
 Y <- c(Y1,Y2)
 X1bar <- mean(X1)
 X2bar <- mean(X2)
 Xbar  <- mean(X)
 Y1bar <- mean(Y1)
 Y2bar <- mean(Y2)
 Ybar  <- mean(Y)
 SSXY1 <- sum((X1-X1bar)*(Y1-Y1bar))
 SSXY2 <- sum((X2-X2bar)*(Y2-Y2bar))
 SSXY  <- sum((X-Xbar)*(Y-Ybar))
 SSX1  <- sum((X1-X1bar)^2)
 SSX2  <- sum((X2-X2bar)^2)
 SSX   <- sum((X-Xbar)^2)
 SSY1  <- sum((Y1-Y1bar)^2)
 SSY2  <- sum((Y2-Y2bar)^2)
 SSY   <- sum((Y-Ybar)^2)

 Delta1 <- SSY1-(SSXY1)^2/SSX1 + SSY2-(SSXY2)^2/SSX2
 Delta2 <- SSY1+SSY2 - (SSXY1+SSXY2)^2/(SSX1+SSX2)
 n  <- length(X1)
 m  <- length(X2)
 N  <- n+m
 Fb <- (Delta2-Delta1)/(Delta1/(N-4))
 pb <- pf(Fb, 1, N-4, lower.tail=F)
 b1 <- SSXY1/SSX1
 b2 <- SSXY2/SSX2
 a1 <- Y1bar-b1*X1bar
 a2 <- Y2bar-b2*X2bar
 bE <- (SSXY1+SSXY2)/(SSX1+SSX2)
 
 Delta3 <- SSY-(SSXY)^2/SSX
 Fa  <- (Delta3-Delta2)/(Delta2/(N-3))
 pa  <- pf(Fa, 1, N-3, lower.tail=F)
 aE1 <- Y1bar-bE*X1bar
 aE2 <- Y2bar-bE*X2bar
 bT  <- SSXY/SSX
 aT  <- Ybar-bT*Xbar

 a.diff <- aE2-aE1
 a.diff.l <- a.diff-qt(1-0.05/2, N-3)*sqrt((1/n+1/m+(X1bar-X2bar)^2/(SSX1+SSX2))*(Delta2/(N-3)))
 a.diff.u <- a.diff+qt(1-0.05/2, N-3)*sqrt((1/n+1/m+(X1bar-X2bar)^2/(SSX1+SSX2))*(Delta2/(N-3)))
 
 list(c(Delta1=Delta1, Delta2=Delta2, F_Slopes=Fb, p_Slopes=pb, 
        Delta3=Delta3, F_Intcps=Fa, p_Intcps=pa,
        Slope1=b1, Slope2=b2, Intcp1=a1, Intcp2=a2,
        Slope_Common=bE, Intcp1_Slope_Common=aE1, Intcp2_Slope_Common=aE2, 
        Diff_Intcps=a.diff, Diff_Intcps_LL=a.diff.l, Diff_Intcps_UL=a.diff.u,
        Slop_Total=bT, Intcp_Total=aT))
}
 
reg.line.diff1 <- unlist(reg.line.diff(X1=ag.a, Y1=bp.a, X2=ag.b, Y2=bp.b))

reg.line.diff1

結果はこちら。二つのSlopeが同じという帰無仮説のp値は0.953で統計学的有意でなく、二つのInterceptが同じという帰無仮説もp値が0.727で統計学的有意でなかった。ということは、傾きも切片も共通の値を使えることになる。全データの傾きは2.946、切片は44.04であった。

> reg.line.diff1
             Delta1              Delta2            F_Slopes 
       8.017313e+04        8.017897e+04        3.571058e-03 
           p_Slopes              Delta3            F_Intcps 
       9.525912e-01        8.037725e+04        1.236460e-01 
           p_Intcps              Slope1              Slope2 
       7.265909e-01        3.119131e+00        3.037802e+00 
             Intcp1              Intcp2        Slope_Common 
       3.780888e+01        3.701073e+01        3.081808e+00 
Intcp1_Slope_Common Intcp2_Slope_Common         Diff_Intcps 
       3.942545e+01        3.458412e+01       -4.841327e+00 
     Diff_Intcps_LL      Diff_Intcps_UL          Slop_Total 
      -3.249540e+01        2.281274e+01        2.945611e+00 
        Intcp_Total 
       4.404461e+01 

共通の傾き・切片の回帰直線を描き入れてみる

abline(a=reg.line.diff1[19], b=reg.line.diff1[18], col=2)

散布図に共通の傾き・切片の回帰直線を赤で描き入れた。ほとんど同一の直線になった。結局、両地区とも年齢と血圧の関係は同じで、もし両地区の年齢構成が同じであれば、平均値の差の検定でも両地区には差がないという結果だっただろうと想像できる。つまり、血圧の差を生んだのは、年齢構成の差だったということだ。

f:id:toukeier:20190107193108p:plain

lm()を使ってみるとどうなるか

それでは、組み込み関数 lm() を使ってみるとどうなるか、やってみる。Anova() を使うために、car パッケージを使えるようにする。インストールしていない場合はインストールする。

データの準備

ar.a <- rep("A",length(ag.a))
ar.b <- rep("B",length(ag.b))

dat <- data.frame(ag=c(ag.a,ag.b),bp=c(bp.a,bp.b),ar=c(ar.a,ar.b))

library(car)

toukeier.hatenablog.com

傾きを等しいと考えてもいいか

年齢と地区の交互作用項を入れたモデルで共分散分析を行う。

lm1 <- lm(bp ~ ag*ar, data=dat)
summary(lm1)
Anova(lm1)

結果は、以下の通り。交互作用項 ag:arB統計学的有意ではないので、傾きを共通としても問題ない。Anova()のResidualsが上記のDelta1(傾きが異なるときの残差平方和)に一致しているのがわかる。

> summary(lm1)

Call:
lm(formula = bp ~ ag * ar, data = dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-112.506  -23.408   -2.317   27.068   85.419 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) 37.80888   40.56783   0.932  0.35591   
ag           3.11913    0.92197   3.383  0.00142 **
arB         -0.79815   69.07337  -0.012  0.99083   
ag:arB      -0.08133    1.36097  -0.060  0.95259   
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 40.45 on 49 degrees of freedom
Multiple R-squared:  0.3669,    Adjusted R-squared:  0.3282 
F-statistic: 9.467 on 3 and 49 DF,  p-value: 4.848e-05

> Anova(lm1)
Anova Table (Type II tests)

Response: bp
          Sum Sq Df F value    Pr(>F)    
ag         33787  1 20.6498 3.619e-05 ***
ar           198  1  0.1212    0.7292    
ag:ar          6  1  0.0036    0.9526    
Residuals  80173 49                      
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

傾きは等しいと考えてよく、共通の傾きを使う

交互作用項を外して、年齢と地区の線形結合モデルとする。

lm2 <- lm(bp ~ ag+ar, data=dat)
summary(lm2)
Anova(lm2)

結果はこちら。B地区の切片は、A地区と比べると-4.84と推定されたが、統計学的に有意ではなく、ゼロと考えてもよい。Anova()のResidualsはDelta2(傾き共通で切片が異なるときの残差平方和)と一致している。上記 lm() を使わない方法では、このDelta2と前のDelta1を使って、傾きを同じと考えてよいかの検定をしている。

> summary(lm2)

Call:
lm(formula = bp ~ ag + ar, data = dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-112.983  -23.270   -2.574   27.057   84.853 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  39.4255    29.9289   1.317    0.194    
ag            3.0818     0.6714   4.590 3.01e-05 ***
arB          -4.8413    13.7681  -0.352    0.727    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 40.04 on 50 degrees of freedom
Multiple R-squared:  0.3669,    Adjusted R-squared:  0.3416 
F-statistic: 14.49 on 2 and 50 DF,  p-value: 1.089e-05

> Anova(lm2)
Anova Table (Type II tests)

Response: bp
          Sum Sq Df F value    Pr(>F)    
ag         33787  1 21.0697 3.009e-05 ***
ar           198  1  0.1236    0.7266    
Residuals  80179 50                      
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

切片も共通でかまわないという結果の下、回帰直線は一つでよい

切片も共通で構わないとなれば、モデルはもっとも単純になる。単回帰モデルだ。

lm3 <- lm(bp ~ ag, data=dat)
summary(lm3)
Anova(lm3)

結果はこちら。ag の Estimate が共通の傾きで2.946、切片はInterceptの44.04となる。Delta3(傾きも切片も共通の時の残差平方和)が Anova(lm3)の結果中のResidualsと一致している。切片が同じでよいかどうかを、lm()を使わないやり方では、このDelta3と前のDelta2を使って検定している。

> summary(lm3)

Call:
lm(formula = bp ~ ag, data = dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-113.455  -22.652   -3.727   28.294   84.654 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  44.0446    26.6599   1.652    0.105    
ag            2.9456     0.5437   5.418 1.65e-06 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 39.7 on 51 degrees of freedom
Multiple R-squared:  0.3653,    Adjusted R-squared:  0.3529 
F-statistic: 29.35 on 1 and 51 DF,  p-value: 1.647e-06

> Anova(lm3)
Anova Table (Type II tests)

Response: bp
          Sum Sq Df F value    Pr(>F)    
ag         46264  1  29.355 1.647e-06 ***
Residuals  80377 51                      
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

回帰直線の差の検定の結果の見方はこちらも参照

toukeier.hatenablog.com

まとめ

二つの回帰直線の差の検定をつかって、真の要因か見かけの要因かを見分ける方法を実行した。線形モデルの関数 lm() を使って、共分散分析をすればよいだけなのだが、今回の場合のように順を追ってみていくと、共分散分析モデルの意味合いがとてもクリアに理解できる。交絡因子の調整の重要さを改めて思い知った。

ある健康指標について、年齢調整をしても地域差があるのかどうか確認したい、というような場合にとても重要な方法だ。