2009.4.11
2009.3.9より
２本の回帰直線の比較を、（ｘと２つのサンプルのちがいを表す名義変数のあいだの）交互作用を検討することによって行っている場合をときどき見かける。そこで仮定されていることをもう少し詳しく見てみた。。
まずは例に使ったデータを以下に示す。第１のサンプルのxが以下のxx1である、
> xx1
[1] 1.901628 1.189416 2.497983 2.970987 3.958775 4.822747 4.565434 1.689338
[9] 8.571783 1.564368 7.491317 2.822647
第１のサンプルのyが以下のyy1である、
> yy1
[1] 2.248822 1.979994 3.750488 3.545448 5.246537 5.256942 5.351267 2.810869
[9] 9.463105 2.548663 8.418361 2.930062
第２のサンプルのxが以下のxx2である、
> xx2
[1] 9.980802 8.324157 6.803240 6.054728 8.228597 8.778725 7.613337 1.163588
[9] 2.651272 7.439294 3.565445 7.699280
第２のサンプルのyが以下のyy2である、
> yy2
[1] 16.8921915 14.6140562 11.4044017 10.2141847 15.3765631 18.5378851
[7] 16.0148672 -0.0994805 2.2876513 12.9308352 4.9445873 15.0245021
両方のサンプルをあわせて１つにしたものの、xが以下のxx10、yがyy10、両サンプルを区別するためのダミー変数がff01である（データ自体は上記と同じ）、
> yy10
[1] 2.2488221 1.9799937 3.7504875 3.5454482 5.2465367 5.2569419
[7] 5.3512673 2.8108693 9.4631050 2.5486629 8.4183612 2.9300625
[13] 16.8921915 14.6140562 11.4044017 10.2141847 15.3765631 18.5378851
[19] 16.0148672 -0.0994805 2.2876513 12.9308352 4.9445873 15.0245021
> xx10
[1] 1.901628 1.189416 2.497983 2.970987 3.958775 4.822747 4.565434 1.689338
[9] 8.571783 1.564368 7.491317 2.822647 9.980802 8.324157 6.803240 6.054728
[17] 8.228597 8.778725 7.613337 1.163588 2.651272 7.439294 3.565445 7.699280
> ff01
[1] 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1
まず、それぞれのサンプルへの回帰直線あてはめである。以下では関数lm()を使っているが、関数glmにおいて、family=gaussianかつlink=identityとしても同じである。
第１のサンプル
> resc1<-lm(yy1~xx1)
> summary(resc1)
Call:
lm(formula = yy1 ~ xx1)
Residuals:
Min 1Q Median 3Q Max
-0.6755 -0.2501 0.0360 0.2488 0.4927
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.75276 0.21113 3.565 0.00514 **
xx1 1.01069 0.04901 20.624 1.59e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3829 on 10 degrees of freedom
Multiple R-squared: 0.977, Adjusted R-squared: 0.9747
F-statistic: 425.4 on 1 and 10 DF, p-value: 1.589e-09
第２のサンプル、
> resc2<-lm(yy2~xx2)
> summary(resc2)
Call:
lm(formula = yy2 ~ xx2)
Residuals:
Min 1Q Median 3Q Max
-2.2354 -0.6950 -0.1526 0.3850 2.1049
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.8689 1.0185 -2.817 0.0183 *
xx2 2.2039 0.1452 15.176 3.13e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.294 on 10 degrees of freedom
Multiple R-squared: 0.9584, Adjusted R-squared: 0.9542
F-statistic: 230.3 on 1 and 10 DF, p-value: 3.125e-08
続いて、両方を一緒に分析したものである。まず、サンプルの区別なくプールしたもの、つまり、両サンプルをプールして１本の直線をあてはめたものである、
> resc10<-lm(yy10~xx10)
> summary(resc10)
Call:
lm(formula = yy10 ~ xx10)
Residuals:
Min 1Q Median 3Q Max
-5.0815 -0.6858 0.2465 1.2380 3.6026
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.6356 0.8539 -1.916 0.0685 .
xx10 1.8876 0.1467 12.866 1.03e-11 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.019 on 22 degrees of freedom
Multiple R-squared: 0.8827, Adjusted R-squared: 0.8774
F-statistic: 165.5 on 1 and 22 DF, p-value: 1.030e-11
次にサンプルのちがいを表すダミー変数も説明変数に入れたもの。共分散分析での主役になる、２本の平行な線をあてはめた場合にあたる。
> resc101<-lm(yy10~xx10+ff01)
> summary(resc101)
Call:
lm(formula = yy10 ~ xx10 + ff01)
Residuals:
Min 1Q Median 3Q Max
-3.25833 -0.75121 -0.01644 1.33139 3.22872
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.7225 0.7582 -2.272 0.0337 *
xx10 1.6851 0.1511 11.151 2.78e-10 ***
ff01 2.2390 0.8488 2.638 0.0154 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.791 on 21 degrees of freedom
Multiple R-squared: 0.9119, Adjusted R-squared: 0.9035
F-statistic: 108.7 on 2 and 21 DF, p-value: 8.371e-12
そして、交互作用もあるものである。すぐ上の例では、２本の直線の傾きは共通だが、こちらでは傾きも切片もちがってかまわない。
> resc10i1<-lm(yy10~xx10*ff01)
> summary(resc10i1)
Call:
lm(formula = yy10 ~ xx10 * ff01)
Residuals:
Min 1Q Median 3Q Max
-2.235364 -0.468314 0.004677 0.248784 2.104908
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.7528 0.5260 1.431 0.167851
xx10 1.0107 0.1221 8.278 6.85e-08 ***
ff01 -3.6216 0.9170 -3.950 0.000792 ***
xx10:ff01 1.1932 0.1624 7.347 4.23e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.9541 on 20 degrees of freedom
Multiple R-squared: 0.9762, Adjusted R-squared: 0.9726
F-statistic: 273.2 on 3 and 20 DF, p-value: < 2.2e-16
最後の交互作用もある場合で、ff01が0の場合を計算してみると、切片0.7528、回帰係数（傾き）1.0107となる。
サンプルの差を表すのに0と1ではなく、以下のような、1と2のちがいを使ってみる
> ff12
[1] 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
交互作用がある場合は、
> lm(yy10~xx10*ff12)
Call:
lm(formula = yy10 ~ xx10 * ff12)
Coefficients:
(Intercept) xx10 ff12 xx10:ff12
4.3744 -0.1825 -3.6216 1.1932
これも同様に、ff12が1の場合と2の場合の切片や回帰係数（傾き）を計算してみると、上記の、第1のサンプル単独な場合、第2のサンプル単独な場合とそれぞれ同じである。
２本の、傾きも切片もちがってかまわない直線をあてはめるのだから、交互作用（サンプルのちがいを表す２値的な変数との）が入っている場合は、それぞれ別々に回帰直線を当てはめた場合と同じかというと、そうではない。対数尤度を見てみる。
> logLik(resc1)
'log Lik.' -4.414788 (df=3)
> logLik(resc2)
'log Lik.' -19.02389 (df=3)
> logLik(resc10i1)
'log Lik.' -30.73781 (df=5)
となり、個別に当てはめた方がだいぶ大きい（よくあてはまっている）。
少しまわり道をそれて、これらの対数尤度がどのように計算されているか確認してみよう。まず、第１のサンプル単独の場合である。
残差分散を計算する。12は標本数である。
> v1<-sum((fitted(resc1)-yy1)^2)/12
> esd1<-sqrt(v1)
正規分布母密度関数の関数dnormを使えば
> sum(dnorm((fitted(resc1)-yy1),mean=0,sd=esd1,log=TRUE))
[1] -4.414788
と対数尤度が計算でき、これは関数lmが計算したものと同じである。
同様に、第２のサンプル単独の場合
> v2<-sum((fitted(resc2)-yy2)^2)/12
> esd2<-sqrt(v2)
> sum(dnorm((fitted(resc2)-yy2),mean=0,sd=esd2,log=TRUE))
[1] -19.02389
こちらも、関数lmが計算したものと同じである。
最後に交互作用がある場合である、
> v10i1<-sum((fitted(resc10i1)-yy10)^2)/24
> esd10i1<-sqrt(v10i1)
> sum(dnorm((fitted(resc10i1)-yy10),mean=0,sd=esd10i1,log=TRUE))
[1] -30.73781
これも関数lmが計算したものと同じである。すると、この場合には、どちらのサンプルにも同じ残差分散の値を使っていることが分かる。
残差分散の値は
> v1
[1] 0.1222021
> v2
[1] 1.394828
> v10i1
[1] 0.7585152
である。
（サンプルのちがいを表す２値的な変数との）交互作用が入っている場合には残差分散は１つの値を使っているのに対して、それぞれ別々に回帰直線を当てはめた場合ではサンプルによりちがった値を使っているのだった。確率モデルとしてみると、正規線形モデルでの回帰直線当てはめは、切片、傾き、残差分散の３つのパラメーターを持つモデル使っていることになるが、交互作用が入っている場合には切片２つ、傾き２つ、残差分散１つで合計の５つのパラメーター、それぞれ別々に回帰直線を当てはめる場合では切片、傾き、残差分散がそれぞれ２つで合計の６つのパラメーターということになる。ざきほどの対数尤度の差は、この残差分散の差である。つまりサンプルにより直線の周りのばらつきの程度には差があるのに、共通した値を使っているため、確率モデルとして見るとあてはまりが悪く対数尤度が低くなっていたのであった。