<< 目次へ
最終更新日 ： 2009.5.12
重回帰モデル
重回帰モデルは単回帰モデルを拡張したモデルで、単回帰モデルにおいて独立変数は1つだけであったが、重回帰モデルでは2つ以上の独立変数を含むようなモデル式になる。例えば、3つの連続型である独立変数を用いる場合は次式のように、かなり複雑なモデルが考えられる。
Y = X1 + X2 + X3 + X1*X2 + X1*X3 + X2*X3 + X1*X2*X3 ・・・(1)
これは2次の交互作用項（X1*X2, X1*X3, X2*X3）に加えて、3次の交互作用項（X1*X2*X3）が含まれているモデルである。一般線形モデルにおける重回帰モデルは基本的に線形関係が仮定できることが前提となっているので、非線形モデル、例えば次のようなモデルは扱うことができない。
Y = X1 + X2 + X3 + X1^2 + X2^2 + X3^3 + X1*X2 + X1*X3 + X2*X3 + X1*X2*X3 ・・・(2)
これは(1)のモデルから、さらに2次項（X1^2, X2^2, X3^2）が追加されている。重回帰分析を行うにあたって、まず最初に注意すべきことは「3つの独立変数を用いる場合でもこれだけ複雑なモデルを考えなければならない」ということである。
モデル選択の基本は最も複雑なモデルから最も単純なモデルへ、というものであるから、最初に最も複雑なモデルを解析することになる。例えば、4つの独立変数を用いる場合は
Y = X1 + X2 + X3 + X4 + X1^2 + X2^2 + X3^2 + X4^2 + X1*X2 + X1*X3 + X1*X4 + X2*X3 + X3*X4 ・・・(3) X1*X2*X3 + X1*X2*X4 + X2*X3*X4 X1*X2*X3*X4
というかなり複雑なモデルも考えられるわけである。もっともRには自動的に変数を選択して最適なモデルを選んでくれるようなstep()という関数が用意されているので、手間は掛からないのだがむやみに多くの変数を投入することはあまりおススメできない。
また、これに関連して覚えておくべきことがある。それは重回帰分析ではn/3個以上の母数を推定することはやるべきではない、ということである。nは得られているデータ数（サンプルサイズ）のことで、例えば50人分のデータが得られている場合は50/3=16.6となるので、16個の母数しか推定すべきではないのである。
式(3)のモデルでは18個の母数が存在するので、このようなモデルを扱うのであれば、あらかじめn/3=18となるようなn個のデータを取っておく必要があるというわけだ。ちなみに、これは最低でも54人以上のデータが必要となる。
重回帰の練習のため、表1のデータを用いて解析をしてみよう。これはオックスフォード大学出版会のいちページ（http://www.oup.com/grafenhails/）から自由に観覧できるようになっている。表1のデータはChapter 07のワークシートにあるTreesというシートより得ることができる。
ちなみに、このデータは高さと直径を使って木の体積を予測しようと取られたものである。高さはHEIGHT、直径はDIAMETER、そして体積がVOLUMEである。
| DIAMETER | HEIGHT | VOLUME |
| 8.3 | 70 | 10.3 |
| 8.6 | 65 | 10.3 |
| 8.8 | 63 | 10.2 |
| 10.5 | 72 | 16.4 |
| ・ ・ ・ （以下省略） |
・ ・ ・ （以下省略） |
・ ・ ・ （以下省略） |
まずはstep()を使わずに手作業で変数選択を行っていく方法を紹介しよう。重回帰分析を行う際には必ず各変数間の相関関係を散布図によって確認しておく必要がある。そしてこれは関数pairs()によって容易に実行することができる。
---
# データセットの準備
# ワークシートの任意のセル範囲を適当なテキストエディタにコピペする
# Trees.txtとしてDドライブに保存する
dat <- read.table("D://Trees.txt", header=TRUE)
attach(dat)
pairs(dat, panel=panel.smooth)
---
基本的に一番下にあるVOLUMEの行だけを見ればよい。この2つは縦軸にVOLUMEを、横軸にそれぞれDIAMETERとHEIGHTをとった場合の散布図となっている。また、pairs()の引数にpanel=panel.smoothを指定することによって、図中の赤い曲線が描かれるようになる。これを特にスムージング曲線などといったりする。
VOLUMEとDIAMETERは明らかに正の相関関係が見られるようである。VOLUMEとHEIGHTも同様に正の相関関係が見られるといえるだろう。ただし、ばらつきの程度は少し大きいようである。
続いてmgcvライブラリに含まれているgam()を使って、変数間に曲線関係が存在するかどうかを確かめてみるという方法がある。gam()は一般化加法モデルを扱うための関数であり、s()はノンパラメトリック平滑化を行うための関数である。これらを理解するとなると、またそれはそれで時間がかかるので、ここでは「グラフをみて曲線関係があるかどうか確認できるのだ」ということだけ気に留めておけばよいだろう。
--- > library(mgcv) #パッケージmgcvの呼び出し > model <- gam(VOLUME ~ s(DIAMETER) + s(HEIGHT)) > par(mfrow=c(1,2)) #プロットウィンドウを1行2列で表示させる > plot(model) #平滑化曲線のプロット Press return for next page.... > par(mfrow=c(1,1)) #プロットウィンドウを元に戻す（1行1列での表示） ---
DIAMETERの方に描かれているのは、曲線になっているがこの程度では曲線関係があるとはいえない。この程度とはどの程度かという疑問もあるだろうが、これは経験的な判断となってしまう部分も多い。ただ、この時点ではあくまでも事前確認に過ぎないので、それほど厳密な判断を求める必要はない。
さて、一方でHEIGHTの方はどうだろうか。これは明らかに曲線関係があるとはいえないだろう。結局、どちらも曲線関係がみられないようなので、線形モデルとして扱って問題ないといえる。もし仮に、ここで明らかな曲線関係が見られるようであると、非線形モデルとして扱う必要が出てくるかもしれない（モデルに2次項DIAMETER^2とHEIGHT^2を含めることになる）。
従属変数VOLUMEと独立変数DIAMETER、HEIGHTそれぞれの間に曲線関係がないと分かったところで、次に行うべきことは複雑な交互作用が存在するかどうかを確かめることである。これはモデルに多くの独立変数を含める場合には極めて重要な作業となる。
これを行うためにtreeパッケージに含まれるtree()という樹木モデルに当てはめるための関数を用いる。なお、treeパッケージは標準で実装されていないため、ユーザ自らでダウンロードする必要がある。
---
#treeパッケージのインストール
> install.packages("tree")
#treeパッケージの呼び出し
> library(tree)
#樹木モデルの指定
> model2 <- tree(VOLUME ~ ., data = dat)
#樹木図を描く
> plot(model2)
> text(model2)
---
 |
樹木図というのは、簡単にいえばある基準より大きいか小さいかを判断して枝分かれさせていくような手法である。上の樹木図を見てみると、1番上の「DIAMETER < 16.15」と表示されている部分がある。右の方へ伸びている枝は、DIAMETERが16.15よりも大きい場合のVOLUMEの平均値を表している。
# DIAMETER>16.15の場合のVOLUMEのデータ > VOLUME[DIAMETER > 16.15] [1] 42.6 55.4 55.7 58.3 51.5 51.0 77.0 # これの平均値 樹木図の55.93という値と同じことがわかる > mean(VOLUME[DIAMETER > 16.15]) [1] 55.92857
逆に左側へ伸びている枝はDIAMETERが16.15以下のVOLUMEのデータの構造を表している。まず「DIAMETER<12.45」という基準からさらに枝分かれしているが、右側の30.53という値はDIAMETERが16.15以下という条件の下で、さらにDIAMETER < 12.45の場合のVOLUMEの平均値である。
# DIAMETERが16.15以下、かつDIAMETERが12.45より大きい場合のVOLUMEのデータ > VOLUME[DIAMETER <= 16.15 & DIAMETER > 12.45] [1] 22.2 33.8 27.4 25.7 24.9 34.5 31.7 36.3 38.3 # これの平均値 DIAMETER<12.45の右枝の30.53と同じである > mean(VOLUME[DIAMETER <= 16.15 & DIAMETER > 12.45]) [1] 30.53333
DIAMETER < 10.75という枝についても同様である。実をいうと、今回の例では用いている説明変数の数が2つだけなので、わざわざ樹木図を描いてみるまでもない。しかし、説明変数の数が3つ以上になると複雑な交互作用が存在する可能性が出てくるので、この樹木図を描いてみることが重要となる。なぜなら、樹木図は容易に交互作用を発見することができるからである。
それでは交互作用の発見のために、例としてRに最初から用意されているswissのデータセットを用いて同じように樹木図を描いてみよう。swissは6個の変数をもったデータセットで、これは：
> data(swiss)
とすることによって使用することができるようになる。また、このデータセットにある6個の変数は以下のようなものである。
このデータセットについて樹木図を描いてみる：
> result <- tree(Fertility ~ ., data=swiss) > plot(result) > text(result)
 |
端的にいうと分かれて伸びている枝の長さが長い変数は重要な変数である。だから、パッと見てEducationとCatholicという2つの変数はFertilityを説明するのに大きく関与している変数だと分かる。それはともかくとして、この樹木図からどうしたら交互作用があるかどうかを、ある程度でも予測することができるだろうか。これも端的にいうことしかできないのだが、視覚的には枝分かれの数が複雑（多い）場合は最初から、交互作用が複雑であろうな、と注意深く分析していかなければならない。最初に描いた樹木図に比べると、今回描いた樹木図は複雑であるから、交互作用にも注意しなければならないであろう（ということしか実は分からない・・・）。
樹木図は大まかなデータの構造を見せてくれるが、どの変数同士の交互作用が認められるかは実際にモデル解析をしてみないと分からない。まぁ、こういうデータの見方（表現方法）を知っておくのも重要かと(^_^;)
【swissのデータセットを用いた解析は別ページを参照のこと】
ここまでの内容を一旦、まとめておこう。
これは重回帰分析を行う前にやっておくべきことである。これをすることによって、あらかじめ曲線関係の有無と交互作用の有無が把握できるからである。しかし、実際にそれが認められるかどうかはモデル解析によって確かめる必要があることはいうまでもない。
今回用いた表1のデータについて、明らかに曲線関係はなさそうなので、モデルに2次項を含める必要はないだろう。これが分かっているだけでも推定すべきパラメータを節約することができたわけである。もっとも、2つの変数しかもたないので最も複雑なモデルから考えていったとしても、それほど手間にはならないだろう。しかしながら、扱う変数が多くなるほど考えられる全てのモデルについて解析を行うことは難しいものとなる。
試しに最も複雑なモデル、すなわち2次項を含んだモデルを解析してみよう。
---
#lm()を用いて最も複雑なモデルを解析する
> model3 <- lm(VOLUME ~ DIAMETER + HEIGHT + I(DIAMETER^2) + I(HEIGHT^2) + DIAMETER*HEIGHT)
> summary(model3) #係数表の出力
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.60706 62.90855 0.105 0.9172
DIAMETER -5.12160 2.46674 -2.076 0.0483 *
HEIGHT 0.29491 1.77852 0.166 0.8696
I(DIAMETER^2) 0.16393 0.10089 1.625 0.1167
I(HEIGHT^2) -0.00494 0.01312 -0.376 0.7097
DIAMETER:HEIGHT 0.06628 0.05671 1.169 0.2535
---
案の定、2次項（DIAMETER^2, HEIGHT^2）は有意ではない。したがって、今度は2次項を削除して、改めてモデル解析を行う。こういう場合、update()という関数が便利である。
---
# 2次項を削除してモデルを更新する
# それを新たにmodel4という変数に代入
> model4 <- update(model3, ~ . - I(DIAMETER^2) - I(HEIGHT^2))
> summary(model4) # 係数表の出力
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 69.39632 23.83575 2.911 0.00713 **
DIAMETER -5.85585 1.92134 -3.048 0.00511 **
HEIGHT -1.29708 0.30984 -4.186 0.00027 ***
DIAMETER:HEIGHT 0.13465 0.02438 5.524 7.48e-06 ***
---
今度は全ての項が有意になっているので、最も当てはまりの良いモデルだといえる。
連続変数間の交互作用
2要因の分散分析ではカテゴリカル型の変数間の交互作用について説明したが、連続変数間の交互作用の解釈はそれに比べてかなり簡単である。なぜなら、連続型の変数の場合は演算子*は意味どおり「かける」という意味だからである。だから、交互作用項を含んだモデル
VOLUME = DIAMETER + HEIGHT + DIAMETER*HEIGHT
は、DIAMETER*HEIGHTというベクトル同士の積を表しているのであある。例えば、DEという変数にDIAMETERとHEDIGTの積を代入してそれをモデルに含めて解析しても同じ結果が得られるのである。
---
# DEにDIAMETERとHEIGHTの積を代入
> DE <- DIAMETER*HEIGHT
> DE
[1] 581.0 559.0 554.4 756.0 866.7 896.4 726.0 825.0 888.0 840.0 892.7
[12] 866.4 866.4 807.3 900.0 954.6 1096.5 1143.8 972.7 883.2 1092.0 1136.0
[23] 1073.0 1152.0 1255.1 1401.3 1435.0 1432.0 1440.0 1440.0 1792.2
# DIAMETER*HEIGHRの代わりにDEを追加
> model5 <- lm(VOLUME ~ DIAMETER + HEIGHT + DE))
> summary(model5) # 係数表の出力
Call:
lm(formula = VOLUME ~ DIAMETER + HEIGHT + DE)
Residuals:
Min 1Q Median 3Q Max
-6.5821 -1.0673 0.3026 1.5641 4.6649
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 69.39632 23.83575 2.911 0.00713 **
DIAMETER -5.85585 1.92134 -3.048 0.00511 **
HEIGHT -1.29708 0.30984 -4.186 0.00027 ***
DE 0.13465 0.02438 5.524 7.48e-06 ***
---
結局、樹木の直径（DIAMETER）が大きくなるほど樹木の高さ（HEIGHT）のもつ体積への影響の度合いが大きくなるということを意味している。
回帰診断
回帰分析を行ったあとにやるべきことがある。それは、
を確認することである。
分散の均一性については、標準化された残差をプロットする（下図の左下の図）ことによって確認することができる。図の見方としては、データポイントが全体に満遍なく散布していることが望ましく、適合値（Fitted values）が大きくなるにつれて分散も大きくなるような傾向を示すようだとまずいわけである。
続いて誤差の正規性についてはQ-Qプロットとよばれるものを用いればよい（右上の図）。これは描かれたデータポイントが直線にぴったりと乗っている状態が最も好ましい。あまりにもデータポイントが直線から離れている（全体的に曲線的に見える）ようだと正規性に疑いがあると判断できる。
ちなみに、左上の図は残差プロットである。これは実測値と予測値の差（残差）を示している。そして右下の図ははずれ値を見つけるためのものだと思えばよいだろう。この図において、特別に突起した線があるようだと、その番号のデータははずれ値であり、結果に良くない影響を及ぼしているかもしれないと判断するわけである。だから、今回のケースだと18番目のデータを削除して分析してみた方がよいかもしれない（ただ、ここではそれを確認しない）。
--- > par(mfrow=c(2,2)) # プロットウィンドウを4分割にする > plot(model4, which=1:4) # 回帰診断図のプロット > par(mfrow=c(1,1)) # プロットウィンドウの分割区分を元に戻す ---
参考文献
