二つのデータセットがあって、二つの回帰直線が描けたとき、そのあとどうすればいいか？
そのあとは、傾きが同じと言えるか？さらには切片が同じと言えるか？と進んでいく。
二つの回帰直線の差を検定してみる。
- 回帰直線の差の検定のためのサンプルデータ
- 回帰直線の差の検定の前にまず散布図を描いてみる
- 回帰直線を求めてみる
- 回帰直線を描き入れてみる
- 回帰直線の差の検定（傾きが同じかどうか）
- 回帰直線の傾きが共通で切片が異なるという結果を回帰直線にして描き入れてみる
- 回帰直線の差の検定を lm() で実施してみるとどうなるか
- まとめ
- 引用書籍
回帰直線の差の検定のためのサンプルデータ
データはこちら。C地区とO地区のNO2濃度（X）とせき・たん有症割合（Y）のデータ。
NO2.c <- c(0.025,0.030,0.016,0.028,0.021,0.041,0.012,0.015,0.030,0.018,0.023) pre.c <- c(4.1,6.8,2.6,11.5,2.8,10.7,2.7,2.8,3.5,5.0,3.6) are.c <- rep("c",length(NO2.c)) NO2.o <- c(0.024,0.033,0.030,0.017,0.022,0.027,0.016,0.017,0.020,0.022) pre.o <- c(6.8,8.7,11.4,4.5,5.6,10.5,5.2,6.8,6.0,8.0) are.o <- rep("o",length(NO2.o))
回帰直線の差の検定の前にまず散布図を描いてみる
散布図を描くにはplot()を使う。
O地区は、points()を使って、三角（pch=2）でプロットしている。
plot(NO2.c, pre.c, las=1, xlim=c(0,0.05), ylim=c(0,12), xlab="Annual average of NO2 concentration", ylab="Prevalence of cough and/or sputum (%)") title(" Difference between two regression lines ") points(NO2.o, pre.o, pch=2)
回帰直線を求めてみる
function(){} で自作した関数で回帰直線の傾きと切片を求める。
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)) } linear.reg(x=NO2.c, y=pre.c) linear.reg(x=NO2.o, y=pre.o) linreg.c <- unlist(linear.reg(x=NO2.c, y=pre.c)) linreg.o <- unlist(linear.reg(x=NO2.o, y=pre.o))
結果はこちら。
> linear.reg(x=NO2.c, y=pre.c) [[1]] Slope t_Slope p_Slope LL_Slope UL_Slope 279.85418266 3.23205342 0.01028875 83.98051896 475.72784635 Intercept t_Intcp p_Intcp LL_Intcp UL_Intcp -1.48929394 -0.69132430 0.50680073 -6.36257401 3.38398613 Error Var. Resid. SE R^2 5.32855590 2.30836650 0.53718391 > linear.reg(x=NO2.o, y=pre.o) [[1]] Slope t_Slope p_Slope LL_Slope UL_Slope 317.540322581 3.780409062 0.005385004 123.844538638 511.236106524 Intercept t_Intcp p_Intcp LL_Intcp UL_Intcp 0.110080645 0.055902022 0.956790802 -4.430835999 4.650997290 Error Var. Resid. SE R^2 2.099679940 1.449027239 0.641118694
結果の読み方は、こちらも参照のこと。
https://toukeier.hatenablog.com/entry/2019/01/05/132758toukeier.hatenablog.com
回帰直線を描き入れてみる
先ほどの散布図に二つの回帰直線を描き入れてみる。
Ｃ地区が実線で、Ｏ地区が点線（lty=2）だ。
abline(a=linreg.c[6], b=linreg.c[1]) abline(a=linreg.o[6], b=linreg.o[1], lty=2)
描き入れた図はこちら。
回帰直線の差の検定（傾きが同じかどうか）
上の図を見ると、二つの回帰直線の傾きはとても近いが、同じにして話を進めてもいいかどうかの検定。
Rのスクリプトは以下の通り。
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))) conf <- 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=NO2.c, Y1=pre.c, X2=NO2.o, Y2=pre.o)) reg.line.diff1
スクリプト内で、ナントカbarは平均、SSナントカは平方和（平均との差の二乗和）や積和（Xの平均との差とYの平均との差の積の和）。
それらを使って、Delta1, Delta2, Delta3を計算している。
Delta1は、傾きが異なるという仮説（対立仮説）のもとでの残差平方和。
Delta2は、傾きは等しいという仮説（帰無仮説）のもとでの残差平方和。
Delta3は、傾きも切片も等しいという仮説の下での残差平方和。
傾きの検定はDelta1とDelta2を使ってF分布で検定する。
F値の分子の自由度は地区の数（2）マイナス 1=1、分母の自由度は、全体のサンプルサイズ（N）マイナス パラメータの数（4つ）。
4つのパラメータとは、Intercept、Slope、地区、Slopeと地区の交互作用（地区ごとにSlopeが異なるという仮説）の4つ。
切片の検定は、Delta2とDelta3を使ってF分布で検定。
F値の分子の自由度は先ほどと同じく地区の数（2）マイナス 1=1、分母の自由度は、N マイナス パラメータの数（3）。
3つのパラメータは、Intercept、Slope、地区の3つ。
結果は、以下の通り。
> reg.line.diff1 Delta1 Delta2 F_Slopes 64.754442586 65.052361246 0.078212660 p_Slopes Delta3 F_Intcps 0.783108529 96.837793841 8.795034888 p_Intcps Slope1 Slope2 0.008280926 279.854182655 317.540322581 Intcp1 Intcp2 Slope_Common -1.489293937 0.110080645 290.976955534 Intcp1_Slope_Common Intcp2_Slope_Common Diff_Intcps -1.751184680 0.715725414 2.466910094 Diff_Intcps_LL Diff_Intcps_UL Slop_Total 0.719300276 4.214519912 281.451309098 Intcp_Total -0.355561311
傾きが等しいという帰無仮説は棄却されなかった（p_Slopes ≒ 0.78）ので、同じ傾きにして話を進めてもよいことになる。
共通の傾きは Slope_Common ≒ 290.98 である。
また、傾きが等しい上で、切片が等しいという帰無仮説は棄却された（p_Intcps ≒ 0.0083）ので、切片は異なって、C地区とO地区の回帰直線は並行だが別々がよいということになった。
切片の差は Diff_Intcps ≒ 2.47 で約2.5%だ。
95%信頼区間は0.72％から4.2％（Diff_Intcps_LL ≒ 0.72, Diff_Intcps_UL ≒ 4.21）。
回帰直線の傾きが共通で切片が異なるという結果を回帰直線にして描き入れてみる
結果を図に描き入れてみる。
共通の傾きの線は赤（col=2）で表示している。
C地区は実線で、O地区は点線だ。別々の回帰直線（黒）とほとんど違わないことがわかる。
それならばシンプルなほうがいい。
つまり、NO2濃度と有症割合の関係は地区によって異ならないが、各地区での有症割合レベルが高いか低いかの違いがあるというシンプルなほうがいい。
シンプルなほうが、ほかの地区にも応用が利く可能性が高い。
abline(a=reg.line.diff1[13], b=reg.line.diff1[12], col=2) abline(a=reg.line.diff1[14], b=reg.line.diff1[12], col=2, lty=2)
回帰直線の差の検定を lm() で実施してみるとどうなるか
統計ソフトRに組み込まれている関数 lm() を使うとどんなふうにできるか見てみる。
まずデータを少し加工する
dat <- data.frame(NO2=c(NO2.c,NO2.o),pre=c(pre.c,pre.o),are=c(are.c,are.o))
交互作用項を入れたモデル（傾きが別々のモデル）
上記でDelta1を計算したのと同等なのが、交互作用項を入れたモデルになる。
lm1 <- lm(pre~NO2*are, data=dat) library(car) Anova(lm1)
Anova()を使うのでcar パッケージをインストールして、使えるようにしておく。
https://toukeier.hatenablog.com/entry/%3Fp%3D1064toukeier.hatenablog.com
> Anova(lm1) Anova Table (Type II tests) Response: pre Sum Sq Df F value Pr(>F) NO2 85.373 1 22.4129 0.0001918 *** are 31.785 1 8.3446 0.0102038 * NO2:are 0.298 1 0.0782 0.7831085 Residuals 64.754 17 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
NO2:are が交互作用項。
p値は0.7831085である。
上記の傾きが等しい帰無仮説の検定（p_Slopes=0.783108529）と結果が同じことが確認できる。
ResidualsはDelta1に一致する。
交互作用が有意でないので交互作用項を入れないモデルにする
このモデルが二つの回帰直線の傾きが同じと考えた場合に一致する。
スクリプトは交互作用を外したシンプルなものだ。
lm2 <- lm(pre~NO2+are, data=dat) summary(lm2) confint(lm2) Anova(lm2)
結果は以下の通り。
> summary(lm2) Call: lm(formula = pre ~ NO2 + are, data = dat) Residuals: Min 1Q Median 3Q Max -3.4781 -1.3413 -0.1781 0.9595 5.1038 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1.7512 1.5217 -1.151 0.264870 NO2 290.9770 59.8680 4.860 0.000126 *** areo 2.4669 0.8318 2.966 0.008281 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.901 on 18 degrees of freedom Multiple R-squared: 0.6324, Adjusted R-squared: 0.5915 F-statistic: 15.48 on 2 and 18 DF, p-value: 0.0001227
NO2のEstimateが、上記で共通の傾きと称した値（Slope_Common=290.976955534）に一致している。
また、areoのEstimateが上記で言っていた切片の差（Diff_Intcps=2.466910094）に一致している。
InterceptのEstimateと足し合わせると、O地区の切片になる。
InterceptのEstimateはC地区の切片（Intcp1_Slope_Common=-1.751184680）と一致する。
95％信頼区間の結果は以下の通り。
> confint(lm2) 2.5 % 97.5 % (Intercept) -4.9481583 1.445789 NO2 165.1990206 416.754891 areo 0.7193003 4.214520
areoの95%信頼区間は、切片の差の信頼区間（Diff_Intcps_LL=0.719300276, Diff_Intcps_UL=4.214519912）に一致する。
このモデルの残差平方和がDelta2=65.052361246と一致しているのがわかる。
> Anova(lm2) Anova Table (Type II tests) Response: pre Sum Sq Df F value Pr(>F) NO2 85.373 1 23.623 0.0001257 *** are 31.785 1 8.795 0.0082809 ** Residuals 65.052 18 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
もし地域差がなければ地域の要因も取り除いて切片も同じと考えてもよい
地域で有意に異なることがわかったのでここで終了だが、もし地域差もない、つまり切片も異ならないのであれば、完全に一つの回帰直線にしてしまってもよい。
lm3 <- lm(pre~NO2, data=dat) Anova(lm3)
この時の残差平方和がDelta3=96.837793841なのである。
> Anova(lm3) Anova Table (Type II tests) Response: pre Sum Sq Df F value Pr(>F) NO2 80.105 1 15.717 0.000831 *** Residuals 96.838 19 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
まとめ
二つの回帰直線の差の検定を引用書籍に沿って確認してみた。
傾きが異なるかどうかは交互作用項を lm() 内に投入すればよいが、lm()を使わない場合はどうやるか、どんな残差平方和を使って検定するのかが詳細に確認できた。
二つのデータセットがあって、二つの回帰直線が引けたが、さてこの後どうするか？と悩んでしまっている人の次の一手になると思う。
引用書籍
ネタ本は「医学への統計学」6.2.3 二つの母回帰直線の差の検定
最新版はこちら。
