統計ER

統計ソフトRの使い方を中心に、統計解析方法の解説をするブログ。ありそうでなかなか見つからないサンプルサイズ計算などニッチな方法について紹介しています。

統計ソフトRで重回帰分析において簡単に変数選択する方法は?

ブログランキングに参加しています。
まずはぽちぽちっとお願いします。
↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓
にほんブログ村 科学ブログ 数学へ

重回帰分析で変数選択は苦労するもの。 統計ソフトRで重回帰分析をする場合に、 いいプログラムを見つけた。 その名もbestglm。

AIC, BIC(デフォルト), BICqなどの Information Criterion 情報規準を使って ベストの変数の組み合わせを見つけてくれるすぐれもの。 bestglmパッケージのznuclearデータで試してみよう。

bestglmの使い方

事前にインストールしておいたあと bestglmを呼び出す。

install.packages("bestglm")

 

library(bestglm)

  znuclearデータフレームは、原発のデータ。

Data on 32 nuclear power plants. The response variable is cost and there are ten covariates.

  bestglmの場合、data.frame内で、 カラムの順番に決まりがある。 まず独立変数を並べて、 一番最後のカラムに従属変数を置く。   このdata.frameではcostが従属変数だ。 以下のような順番に並んでいる必要がある。 右端がcost。

> round(znuclear, 2)
    date    T1    T2 capacity PR NE CT BW     N PT  cost
1  68.58  0.19 -1.72    -0.62  0  1  0  0  0.89  0  0.17
2  67.33 -1.17  1.01     1.14  0  0  1  0 -1.80  0  0.13
3  67.33 -1.17  1.91     1.14  1  0  1  0 -1.80  0  0.07
4  68.00 -0.79  0.50     1.14  0  1  1  0  0.73  0  1.09
5  68.00 -0.79  1.40     1.14  1  1  1  0  0.73  0  1.05
6  67.92 -0.11 -1.11    -1.79  0  1  1  0 -0.68  0 -0.59
7  68.17 -0.43 -1.23     0.10  0  0  0  0 -0.16  0 -1.22
8  68.42  0.19 -0.25    -2.26  0  0  0  0 -1.80  0 -0.81
9  68.42  0.47 -0.66     0.10  1  0  0  0 -0.16  0  0.15
10 68.33 -0.43  0.85    -0.05  0  1  1  1 -1.09  0  1.24
11 68.58 -0.43  0.23    -1.45  0  0  0  0 -0.68  0 -0.55
12 68.75 -0.11 -1.59    -0.06  0  1  0  0  0.02  0 -0.18
13 68.42  0.47  0.05    -1.67  0  0  1  0 -1.09  0 -0.12
14 68.92  0.98 -1.00     1.08  0  0  0  0  0.18  0  0.37
15 68.92 -0.11  0.32     0.23  0  0  0  1  1.02  0 -0.24
16 68.42 -0.79  0.50    -0.12  0  0  0  0 -0.68  0 -0.05
17 69.50  1.21 -0.15     0.21  0  1  0  0  1.08  0  1.33
18 68.42  0.47  1.25    -1.67  1  0  1  0 -1.09  0 -1.05
19 69.17  0.47  0.50     1.23  0  0  0  0 -1.80  0  1.89
20 68.92  0.73 -0.25     1.08  1  0  0  0  0.32  0  0.34
21 68.75 -0.79  0.76     0.52  0  0  1  1  0.96  0  0.73
22 70.92  2.02 -0.45     0.13  1  1  0  0  1.25  0  1.15
23 69.67  0.73 -0.25    -0.08  0  0  1  0  1.14  0  0.97
24 70.08  1.43 -0.35     0.09  1  0  0  0 -0.68  0  0.91
25 70.42  1.43 -1.98    -1.61  0  0  1  0  1.20  0  0.25
26 71.08  1.64 -0.45     1.38  0  0  1  0  1.30  0  1.27
27 67.25 -0.11  0.14    -0.30  0  0  0  0  0.32  1 -1.94
28 67.17 -1.60 -1.47     0.09  0  0  1  0  0.18  1 -1.07
29 67.83 -0.43  0.14     0.40  0  0  0  1  0.64  1 -1.10
30 67.83 -0.43  0.85     0.40  1  0  0  1  0.64  1 -1.14
31 67.25 -0.11  0.93    -0.30  1  0  0  0  0.32  1 -1.81
32 67.83 -2.62  1.55     0.40  1  0  0  1  0.64  1 -1.23

  bestglm()の結果を見てみる。 BICで検討したBest Modelが表示される。

> out <- bestglm(znuclear)
> out
BIC
BICq equivalent for q in (0.349204366418933, 0.716418902103362)
Best Model:
               Estimate Std. Error   t value     Pr(>|t|)
(Intercept) -38.7480703 7.91826983 -4.893502 4.910313e-05
date          0.5620284 0.11445901  4.910303 4.701224e-05
capacity      0.4759804 0.07818015  6.088252 2.310934e-06
NE            0.6588957 0.19616044  3.358963 2.510375e-03
CT            0.3714664 0.15987847  2.323430 2.858187e-02
N            -0.2277672 0.10786682 -2.111560 4.489115e-02
PT           -0.5982476 0.30044058 -1.991235 5.748951e-02
 
summary()でBest modelのモデルの有意性が確認できる。
> summary(out)
Fitting algorithm:  BIC-leaps
Best Model:
           df  deviance
Null Model 25  4.436699
Full Model 31 31.000000

        likelihood-ratio test - GLM

data:  H0: Null Model vs. H1: Best Fit BIC-leaps
X = 26.563, df = 6, p-value = 0.0001748

  out$Subsetsで、検討したモデル一覧が表示される。 *がBest model。

> out$Subsets
    (Intercept)  date    T1    T2 capacity    PR    NE    CT    BW     N    PT logLikelihood        BIC
0          TRUE FALSE FALSE FALSE    FALSE FALSE FALSE FALSE FALSE FALSE FALSE     0.5079792  -1.015958
1          TRUE FALSE FALSE FALSE    FALSE FALSE FALSE FALSE FALSE FALSE  TRUE    10.2059259 -16.946116
2          TRUE FALSE FALSE FALSE     TRUE FALSE FALSE FALSE FALSE FALSE  TRUE    17.8241085 -28.716745
3          TRUE  TRUE FALSE FALSE     TRUE FALSE FALSE FALSE FALSE FALSE  TRUE    23.3113617 -36.225516
4          TRUE  TRUE FALSE FALSE     TRUE FALSE  TRUE FALSE FALSE FALSE  TRUE    26.6826218 -39.502300
5          TRUE  TRUE FALSE FALSE     TRUE FALSE  TRUE  TRUE FALSE  TRUE FALSE    29.2577991 -41.186919
6*         TRUE  TRUE FALSE FALSE     TRUE FALSE  TRUE  TRUE FALSE  TRUE  TRUE    31.6132054 -42.431995
7          TRUE  TRUE FALSE FALSE     TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE    32.1063164 -39.952482
8          TRUE  TRUE FALSE  TRUE     TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE    33.2254075 -38.724928
9          TRUE  TRUE  TRUE  TRUE     TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE    33.2836564 -35.375690
10         TRUE  TRUE  TRUE  TRUE     TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE    33.3647536 -32.072148

  out$BestModelsで、Bestから5番目までが表示される。

> out$BestModels
  date    T1    T2 capacity    PR   NE   CT    BW     N    PT Criterion
1 TRUE FALSE FALSE     TRUE FALSE TRUE TRUE FALSE  TRUE  TRUE -42.43200
2 TRUE FALSE FALSE     TRUE FALSE TRUE TRUE FALSE  TRUE FALSE -41.18692
3 TRUE FALSE FALSE     TRUE FALSE TRUE TRUE FALSE FALSE  TRUE -40.64612
4 TRUE FALSE FALSE     TRUE FALSE TRUE TRUE  TRUE  TRUE  TRUE -39.95248
5 TRUE FALSE FALSE     TRUE  TRUE TRUE TRUE FALSE  TRUE  TRUE -39.84694

  out$BestModelで、Best modelのcoefficientsだけが表示される。

> out$BestModel

Call:
lm(formula = y ~ ., data = data.frame(Xy[, c(bestset[-1], FALSE), 
    drop = FALSE], y = y))

Coefficients:
(Intercept)         date     capacity           NE           CT            N           PT  
   -38.7481       0.5620       0.4760       0.6589       0.3715      -0.2278      -0.5982  

  Best modelのresidualを用いてQQ plotを描くことができる。

qqnorm(resid(out$BestModel), ylab="residuals, best model")

f:id:toukeier:20201005213955p:plain

まとめ

重回帰分析において情報規準によってベストの変数選択を教えてくれる bestglm()を紹介した。 もちろん計算上もっとも当てはまりのよいモデルを教えてくれているだけで、 変数の選択には変数の意味合いは全く加味されていない。 研究者はこの結果を参考に、 自分で最終解析変数セットを決めなければならない。