最終更新日： 2009.5.12
swissのデータセットを用いた重回帰分析の例
■ まずは変数同士の関係を散布図で確認する。
> pairs(swiss, panel=panel.smooth)
■ つづいて曲線関係がないかを加法モデルを当てはめて確認する。
# gam()を用いて、一般加法モデルの当てはめを行う > gam.result <- gam(Fertility ~ s(Agriculture) + s(Examination) + s(Education) + + s(Catholic) + s(Infant.Mortality)) # 新しいグラフウィンドウ > windows() # グラフウィンドウを2行3列に分割 > par(mfrow=c(2,3)) # スムージング曲線をプロット > plot(gam.result)
これを見ると、どうやらExaminationとEducationとに曲線関係があるようである。したがって、モデルにExamination^2とEducation^2という2次項を組み入れて考える必要がある。
■ 樹木図は描いてもよいが、、、 描くなら以下のようにする。
> tree.result <- tree(Fertility ~ ., data=swiss) > windows() > plot(tree.result) > text(tree.result)
それでは実際に解析の手順（の1例）を示す。
まず曲線関係が疑われるので、交互作用を含まないモデルにExamination^2とEducation^2の2つの2次項を入れて分析してみる。出力結果は係数表の部分のみで、あとの部分は省略してある。これによれば、どちらの2次項も有意ではないがEducation^2についてはp=0.14、Examination^2についてはp=0.06と、これらの項を取り除くには少々の不安が残る（gam()によって得られたグラフをみても、曲線関係がありそうなのは明らかである）。
> lm.res1 <- lm(Fertility ~ . + I(Education^2) + I(Examination^2), data=swiss)
> summary(lm.res1)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 65.59858 10.72032 6.119 3.54e-07 ***
Agriculture -0.12237 0.07346 -1.666 0.10375
Examination -1.52247 0.70188 -2.169 0.03623 *
Education -0.37282 0.41586 -0.896 0.37549
Catholic 0.07795 0.03724 2.093 0.04289 *
Infant.Mortality 1.37971 0.40653 3.394 0.00159 **
I(Education^2) -0.01241 0.00828 -1.498 0.14208
I(Examination^2) 0.03542 0.01854 1.910 0.06347 .
---
確実に曲線関係が認められないわけでもなさそうなので、2次項を含めて、交互作用項を含めたモデルを解析してみる。モデル式の記述で ~ . という部分があるが、これは全ての変数を用いるという意味である。ピリオド（.）は全ての変数という意味である。また、それに ^2 と付け加えられているのは2次の交互作用を指定するという意味で、2次項（累乗項）を意味しているのではないことに注意すること。
まず明らかに有意でない交互作用項を取り除くことからはじめる。不要な交互作用項を取り除くことによって単一の変数が有意になることもあるので、ここでは交互作用項のみを取り除くべきである。ここではAgriculture:Education、Examination:Catholic、Education:Catholicは明らかに有意ではない。
> lm.result1 <- lm(Fertility ~ .^2 + I(Education^2) + I(Examination^2), data=swiss)
> summary(lm.result1)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 259.679228 68.381854 3.797 0.000691 ***
Agriculture -1.966543 0.708890 -2.774 0.009581 **
Examination -8.943802 3.799544 -2.354 0.025565 *
Education 0.568031 3.944942 0.144 0.886505
Catholic -0.380349 0.431479 -0.882 0.385294
Infant.Mortality -5.904220 3.283151 -1.798 0.082541 .
I(Education^2) -0.039096 0.029672 -1.318 0.197956
I(Examination^2) 0.010867 0.057384 0.189 0.851122
Agriculture:Examination 0.032384 0.017709 1.829 0.077747 .
Agriculture:Education 0.003455 0.019102 0.181 0.857714
Agriculture:Catholic 0.003454 0.002995 1.153 0.258170
Agriculture:Infant.Mortality 0.053439 0.031096 1.719 0.096365 .
Examination:Education 0.106770 0.063367 1.685 0.102733
Examination:Catholic -0.004868 0.015071 -0.323 0.749038
Examination:Infant.Mortality 0.284293 0.160229 1.774 0.086515 .
Education:Catholic 0.001836 0.013028 0.141 0.888911
Education:Infant.Mortality -0.123877 0.168556 -0.735 0.468286
Catholic:Infant.Mortality 0.015813 0.017787 0.889 0.381290
---
モデルを更新する場合にはupdate()を用いる。なお、update()を用いずに、度々lm()で新たなモデルを解析しなおしてはいけない。すると今度はEducation:Infant.MortalityとCatholic:Infant.Mortalityが不要なのでこれらを取り除くことにする。
> lm.result2 <- update(lm.result1, ~. - Agriculture:Education - Examination:Catholic -
+ Education:Catholic)
> summary(lm.result2)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 252.084973 59.495468 4.237 0.000179 ***
Agriculture -1.889901 0.641679 -2.945 0.005971 **
Examination -9.165186 3.191009 -2.872 0.007178 **
Education 0.898174 2.982577 0.301 0.765256
Catholic -0.369139 0.338110 -1.092 0.283086
Infant.Mortality -5.342534 2.612026 -2.045 0.049115 *
I(Education^2) -0.039263 0.020120 -1.951 0.059808 .
I(Examination^2) 0.027501 0.036875 0.746 0.461229
Agriculture:Examination 0.034448 0.014430 2.387 0.023054 *
Agriculture:Catholic 0.003538 0.002372 1.491 0.145714
Agriculture:Infant.Mortality 0.049577 0.027941 1.774 0.085527 .
Examination:Education 0.092610 0.040554 2.284 0.029174 *
Examination:Infant.Mortality 0.262715 0.119639 2.196 0.035465 *
Education:Infant.Mortality -0.117195 0.147425 -0.795 0.432502
Catholic:Infant.Mortality 0.012764 0.012953 0.985 0.331834
---
それではEducation:Infant.MortalityとCatholic:Infant.Mortalityの2つの交互作用項を取り除いたモデルを見てみよう。どうやら、さっきからExamination^2の項があることによって、Catholicの項が有意にならないような雰囲気である。具体的にいうと、Examination^2をモデルに組み込むことによって、Chatholicの説明率を下げているのである。
> lm.result3 <- update(lm.result2, ~ . - Education:Infant.Mortality - Catholic:Infant.Mortality)
> summary(lm.result3)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 236.690821 53.193661 4.450 8.78e-05 ***
Agriculture -2.047557 0.597358 -3.428 0.00161 **
Examination -6.731892 2.072904 -3.248 0.00262 **
Education -1.474926 0.524988 -2.809 0.00817 **
Catholic -0.070114 0.141614 -0.495 0.62371
Infant.Mortality -4.833092 2.284621 -2.115 0.04179 *
I(Education^2) -0.032127 0.015133 -2.123 0.04112 *
I(Examination^2) 0.022293 0.036147 0.617 0.54152
Agriculture:Examination 0.032234 0.013913 2.317 0.02667 *
Agriculture:Catholic 0.003053 0.002285 1.336 0.19043
Agriculture:Infant.Mortality 0.060197 0.024142 2.493 0.01768 *
Examination:Education 0.083535 0.035719 2.339 0.02537 *
Examination:Infant.Mortality 0.160592 0.064505 2.490 0.01784 *
---
そこでExamnination^2の項を取り除いて、再度、モデルを解析してみよう。どうやら、Examination^2の項がCatholicの説明率を下げているようではないようである。すると、Chatolicの説明率を下げているのはAgriculture:Catholicなのだろうか、と疑えるのでExamination^2はそのままにして、Agriculture:Catholicの交互作用項のみを取り除いてみる。
> lm.result4 <- update(lm.result3, ~ . - I(Examination^2))
> summary(lm.result4)
Call:
lm(formula = Fertility ~ Agriculture + Examination + Education +
Catholic + Infant.Mortality + I(Education^2) + Agriculture:Examination +
Agriculture:Catholic + Agriculture:Infant.Mortality + Examination:Education +
Examination:Infant.Mortality, data = swiss)
Residuals:
Min 1Q Median 3Q Max
-11.7636 -2.8725 0.2969 2.4269 14.0546
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 228.385099 51.003026 4.478 7.68e-05 ***
Agriculture -1.918621 0.554603 -3.459 0.001442 **
Examination -5.902973 1.564075 -3.774 0.000597 ***
Education -1.642260 0.445441 -3.687 0.000764 ***
Catholic -0.031804 0.126131 -0.252 0.802400
Infant.Mortality -4.814805 2.264116 -2.127 0.040582 *
I(Education^2) -0.037335 0.012446 -3.000 0.004953 **
Agriculture:Examination 0.027386 0.011378 2.407 0.021499 *
Agriculture:Catholic 0.002491 0.002077 1.199 0.238453
Agriculture:Infant.Mortality 0.058042 0.023675 2.452 0.019355 *
Examination:Education 0.098905 0.025360 3.900 0.000416 ***
Examination:Infant.Mortality 0.160501 0.063931 2.511 0.016827 *
---
やはり、Catholicの説明率を下げていたのはAgriculture:Catholicだったようである。しかし、それでも、結局Examination^2の項は不要なので取り除いてみよう。
> lm.result4 <- update(lm.result3, ~ . - Agriculture:Catholic)
> summary(lm.result4)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 226.449866 53.225390 4.255 0.000149 ***
Agriculture -1.786269 0.570728 -3.130 0.003519 **
Examination -5.808849 1.976168 -2.939 0.005789 **
Education -1.425913 0.529545 -2.693 0.010802 *
Catholic 0.113057 0.035844 3.154 0.003299 **
Infant.Mortality -5.381088 2.272558 -2.368 0.023547 *
I(Education^2) -0.040106 0.014060 -2.853 0.007233 **
I(Examination^2) 0.003042 0.033521 0.091 0.928215
Agriculture:Examination 0.018650 0.009603 1.942 0.060213 .
Agriculture:Infant.Mortality 0.063326 0.024296 2.606 0.013351 *
Examination:Education 0.090751 0.035702 2.542 0.015607 *
Examination:Infant.Mortality 0.176270 0.064135 2.748 0.009405 **
---
すると、全ての項が有意になったのでこれで最終的なモデルが出来上がったことになる。今回の場合、もっと初期の段階でExamination^2の2次項を取り除いてもよかったかもしれない。しかし、それをしなかったのは最初のgam()によって得られた図より、曲線関係が強く疑われたためである。ちなみに、step()関数を用いてモデル選択をすることもでき、その方法は間瀬ほか「工学のためのデータサイエンス入門」数理工学社（2004）p131-136に詳しい手順が書かれている。
step()はAICをモデル選択の判断基準にしているため、最終的に選択されるモデルはこのページで示したモデルとは異なる。ここで得られたモデルのAICは314であるのに対して、step()で得られるモデルのAICは317となっている（小数点以下は切り捨てた）。AICはその値が小さいほどよいモデルであることを示しているので、step()を用いてモデル選択をした参考書の例よりも、ここで示した方法（変数増減法）によって得られたモデルのほうがよいといえる。
> lm.result5 <- update(lm.result4, ~ . - Examination^2)
> summary(lm.result5)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 225.408202 51.251845 4.398 9.3e-05 ***
Agriculture -1.773153 0.544463 -3.257 0.002460 **
Examination -5.701934 1.564507 -3.645 0.000838 ***
Education -1.454521 0.419562 -3.467 0.001381 **
Catholic 0.113804 0.034403 3.308 0.002139 **
Infant.Mortality -5.361763 2.231173 -2.403 0.021534 *
I(Education^2) -0.040713 0.012197 -3.338 0.001971 **
Agriculture:Examination 0.018269 0.008518 2.145 0.038783 *
Agriculture:Infant.Mortality 0.062883 0.023470 2.679 0.011054 *
Examination:Education 0.093029 0.025033 3.716 0.000684 ***
Examination:Infant.Mortality 0.175787 0.063028 2.789 0.008395 **
---
