71. 回帰分析と重回帰分析


回帰分析を行なうために以下の関数が用意されている.

ここで対象となるモデルは以下のような線形モデルである.

上式をベクトル表記すると y = Xb + e となる.このときの y は応答ベクトル,X は説明変数のベクトル(モデル行列)で,x0 は切片項(要素が全て 1 である列ベクトル)となっている.

回帰分析と重回帰分析

関数 lm() により線形モデルの当てはめを行うことが出来る.この関数により,回帰分析や分散分析,そして共分散分析を行うことが出来る.

詳しい解説は『工学のためのデータサイエンス入門』(間瀬・神保・鎌倉・金藤 共著,数理工学社) を参照のこと.分散分析や非線形回帰についても詳しい解説が載っている.

関数 lm() の書式と引数

書式は以下の通りである.重要なのはモデルを指定する引数 formula と,データ名を指定する引数 data である.

 
 lm(formula, data, subset, weights, na.action,
    method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
    singular.ok = TRUE, contrasts = NULL, offset = NULL, ...)

引数 formula と回帰のモデル式の対応表を挙げる.

formula

回帰におけるモデル式

y ~ x

モデル式 y = a + bx + ε( ε は誤差項)について,目的変数 y と説明変数 x をベクトルで指定する.

y ~ x1 + x2

モデル式 y = a + b1x1 + b2x2 + ε( ε は誤差項)について,目的変数 y と説明変数 x1,x2 をベクトルで指定する.

y ~ x1 * x2

交互作用項を含んだモデル式( x1:x2 でもよい) y = a + b1x1 + b2x2 + b3x1x2 + ε( ε は誤差項)について,目的変数 y と説明変数 x1,x2 をベクトルで指定する.

y ~ x1 + x2 + x1*x2

上と同じ交互作用モデル.

y ~ (x1 + x2)^2

上と同じ交互作用モデル.

y ~ x - 1

切片(定数)項を除外する( + 0 でも可).

y ~ 1 + x + I(x^2)

多項式回帰 : y = b0 + b1x1 + b2x2 + ε( ε は誤差項).I(x^2) は poly(x,2) でも可.

y ~ x | z

z で条件付けしたときの y の x への単回帰.

y ~ . , data = データ名

あるデータに目的変数 y と説明変数 x1,・・・ が入っている場合で,モデル式がy = a + b1x1 + ・・・ + ε( ε は誤差項)である場合は,まず,目的変数 y をベクトルで指定し,右辺は「 y 以外」という意味で . (ピリオド)を指定することも出来る.

モデル情報を取り出す関数

lm() の返り値は当てはめられたモデルのオブジェクトとなっている.このオブジェクト obj を以下の関数に与えることで,新たな統計処理を行うことが出来る.この中からよく使うものを説明する.

関数

機能

AIC(obj)

AIC(赤池情報量規準)を求める.

anova(obj1 , obj2)

モデルを比較して分散分析表を生成する.

coefficients(obj)

回帰係数 (行列) を抽出.coef(obj) と省略できる.

deviance(obj)

重みつけられた残差平方和.

formula(obj)

モデル式を抽出.

logLik(obj)

対数尤度を求める.

plot(obj)

残差,当てはめ値などの 4 種類のプロットを生成.

predict(obj, newdata=data.frame)

提供されるデータフレームが元のものと同じラベルを持つことを強制する.値は data.frame 中の非ランダムな変量に対する予測値のベクトルまたは行列となる.

print(obj)

オブジェクトの簡略版を表示.

residuals(obj)

適当に重みつけられた残差 (の行列) を抽出する.resid(obj) と省略できる.df.residuals(obj) も参照されたい.

step(obj)

階層を保ちながら,項を加えたり減らしたりして適当なモデルを選ぶ.この探索で見つかった最大の AIC (赤池情報量規準)  を持つモデルが返される.

summary(obj)

回帰分析の完全な要約が表示される.

回帰分析の例

単回帰モデルの当てはめを行う.

 
 result <- lm(Volume ~ Girth, data = trees)                   # 回帰分析を行う
 summary(result)                                              # 分析結果の要約

関数 lm() で回帰分析を行った要約結果を変数 result に代入した場合,result$fstatistic は実際の F 検定結果と自由度は返すが p 値は返さない.自力で F 統計量の p 値を計算する場合は以下のようにすればよい.

 
 result <- summary( result <- lm(Volume ~ Girth, data = trees) )
 f.stat <- result$fstatistic
 p.value <- 1-pf(f.stat["value"],f.stat["numdf"],f.stat["dendf"])
 p.value

value 
    0

両対数を取ってから線形回帰モデルの当てはめを行うこともできる.

 
 result1 <- lm(log(Volume) ~ log(Girth), data = trees)
 summary(result1)

続いて,関数 anova() で分散分析表を書くことが出来る.

 
 anova(result1)

Analysis of Variance Table
Response: log(Volume)
           Df Sum Sq Mean Sq F value    Pr(>F)    
log(Girth)  1 7.9254  7.9254  599.72 < 2.2e-16 ***
Residuals  29 0.3832  0.0132                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

 aov(result1)

Call:
   aov(formula = result1)
Terms:
                log(Girth) Residuals
Sum of Squares    7.925446  0.383244
Deg. of Freedom          1        29
Residual standard error: 0.1149578 
Estimated effects may be unbalanced

以下はモデルの情報を取り出す一つの例である.以下は,result1 で使ったモデルに,さらに説明変数 log(Height) を取り入れて新たなモデルの当てはめを行っている.説明項に多少の変更を加えた場合でモデルの当てはめをやり直す場合には,この関数 update() を用いるのが便利である.

 
 result1 <- lm(log(Volume) ~ log(Girth), data = trees)
 result2 <- update(result1, ~ . + log(Height), data = trees)
 summary(result2)

Call:
lm(formula = log(Volume) ~ log(Girth) + log(Height), data = trees)
Residuals:
      Min        1Q    Median        3Q       Max 
-0.168561 -0.048488  0.002431  0.063637  0.129223 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -6.63162    0.79979  -8.292 5.06e-09 ***
log(Girth)   1.98265    0.07501  26.432  < 2e-16 ***
log(Height)  1.11712    0.20444   5.464 7.81e-06 ***

関数 step() を用いて,result2 のモデルについて,モデルから変数を 1 つ削った場合と変数を削らない場合の解析を全通り行い,各結果の AIC を比較している.

 
 result3 <- step(result2)

Start:  AIC= -152.69 
 log(Volume) ~ log(Girth) + log(Height) 
 
              Df Sum of Sq      RSS      AIC
<none>                        0.185 -152.685
- log(Height)  1     0.198    0.383 -132.185
- log(Girth)   1     4.628    4.813  -53.743

 summary(result3)

Call:
lm(formula = log(Volume) ~ log(Girth) + log(Height), data = trees)
Residuals:
      Min        1Q    Median        3Q       Max 
-0.168561 -0.048488  0.002431  0.063637  0.129223 
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -6.63162    0.79979  -8.292 5.06e-09 ***
log(Girth)   1.98265    0.07501  26.432  < 2e-16 ***
log(Height)  1.11712    0.20444   5.464 7.81e-06 ***

説明変数に非説明変数以外の全てを指定する場合は .(ピリオド)を指定し,切片項を取り除く場合は$-1$を指定する.

 
 result <- lm(log(Volume) ~ .-1, data = trees)
 summary(result)

Call:
lm(formula = log(Volume) ~ . - 1, data = trees)
Residuals:
      Min        1Q    Median        3Q       Max 
-0.188133 -0.072707 -0.001769  0.062948  0.174312 
Coefficients:
       Estimate Std. Error t value Pr(>|t|)    
Girth  0.144695   0.006381   22.68  < 2e-16 ***
Height 0.017830   0.001138   15.66 1.09e-15 ***

パッケージ wle を入れることで,関数 mle.cp(),mle.aic(),mle.cv() を用いることによりそれぞれ Cp 統計量,AIC,クロス・バリデーション規準により重回帰モデルの選択を行う.

 
 result <- mle.cp(log(Volume) ~ log(Girth) + log(Height), data = trees)  # cp 統計量
 result <- mle.cv(log(Volume) ~ log(Girth) + log(Height), data = trees)  # CV
 result <- mle.aic(log(Volume) ~ log(Girth) + log(Height), data = trees) # AIC
 summary(result)

Call:
mle.aic(formula = log(Volume) ~ log(Girth) + log(Height), data = trees)
Akaike Information Criterion (AIC):
     (Intercept) log(Girth) log(Height)      aic
[1,]           1          1           1  -64.560
[2,]           1          1           0  -36.700
[3,]           0          1           1    2.196
[4,]           0          1           0  169.000
[5,]           1          0           1  632.100
[6,]           0          0           1  976.300
[7,]           1          0           0 1158.000
Printed the first  7  best models 

パッケージ wle を入れることで,関数 wle.cp(),wle.aic(),wle.cv() を用いることによりそれぞれ 重み付けされた Cp 統計量,AIC,クロス・バリデーション規準により重回帰モデルの選択を行うことも出来る. 

 


分析結果の取り出し

検定結果を変数などに代入したい場合は,検定結果を適当な変数に代入した後,関数 names() で中身を調べればよい.

 
 out1 <- t.test(y1, y2,var.equal=T)
 names(out1)

[1] "statistic"   "parameter"   "p.value"     "conf.int"    "estimate"   
[6] "null.value"  "alternative" "method"      "data.name"
  
 out1$p.value     # p 値の取り出し

[1] 0.005408422

ただし,回帰分析結果を取り出す場合は,関数 names() を使うだけでは取り出せない場合があり,試行錯誤が必要となる場面がしばしば出てくる・・・.

 
 dummy1 <- data.frame(group=c(rep(1,50), rep(2,50)), y=c(y1,y2) )
 out2 <- summary( lm(y ~ group, data=dummy1) )
 out2$coefficients[ ,"Pr(>|t|)"]               # ラベル名は names(out2) で調べまくる

(Intercept)        group 
0.073337026  0.005408422

 out2$coefficients[2,"Pr(>|t|)"]               # 何とか p 値を取り出す・・・

[1] 0.005408422