Rで一般化線形モデル

今回は一般化線形モデルについて書いてみようと思います。
ただ、自分には完全に理解できる頭はないので、要点だけ掻い摘んで紹介したいと思います。
(なら書かなきゃいいという考え方もありますが、自分の頭の中を整理するために書かせていただきます。)

一般化線形モデルとは?

以前紹介したことのある単回帰分析重回帰分析は、
残差と呼ばれるデータのばらつき(実際は回帰線(予測値)と実際値の差)が正規分布であることが前提でした。
(単回帰分析や重回帰分析は、一般線形モデルといいます。)
しかし実際のデータは正規分布以外の分布(二項分布やポアソン分布など)に従う場合もあります。
そういった正規分布以外の分布でも扱うことができるのが「一般化線形モデル」(Generalized Linear Model、略式:GLM)です。
ちなみに、一般化線形モデルについては、こちらのサイトの情報を参考にさせていただきました。

Rで一般化線形モデルへの適用を行うには?

Rで一般化線形モデルを適用するには、glmという関数を利用します。
使い方は下記のとおりです。

glm(モデル式, family = 目的変数の分布, data = データフレーム)



モデル式には「目的変数~説明変数1+説明変数2+…」を、データフレームは対象のデータを指定します。
familyには目的変数の分布を指定します、目的変数の分布は下記のものがあります。

family 説明
binomial 二項分布。ロジティクス回帰分析などで利用。
poisson ポアソン分布
Gamma ガンマ分布
gaussian 正規分布。重回帰分析と同じ


それでは、一般線形モデル(lm関数)と一般化線形モデル(glm関数)を比較しながらRで実行していきましょう。
まずは一般線形モデルです。

#データを取得。不必要な列は除外して変数に格納
data(airquality)
airq2<-airquality[,1:4]
airq2
Ozone Solar.R Wind Temp
1 41 190 7.4 67
2 36 118 8.0 72
3 12 149 12.6 74
4 18 313 11.5 62
5 NA NA 14.3 56
6 28 NA 14.9 66
7 23 299 8.6 65
8 19 99 13.8 59
9 8 19 20.1 61
10 NA 194 8.6 69
11 7 NA 6.9 74
12 16 256 9.7 69
13 11 290 9.2 66
14 14 274 10.9 68
15 18 65 13.2 58
16 14 334 11.5 64
17 34 307 12.0 66
18 6 78 18.4 57
19 30 322 11.5 68
20 11 44 9.7 62
21 1 8 9.7 59
22 11 320 16.6 73
23 4 25 9.7 61
24 32 92 12.0 61
25 NA 66 16.6 57
26 NA 266 14.9 58
27 NA NA 8.0 57
28 23 13 12.0 67
29 45 252 14.9 81
30 115 223 5.7 79
31 37 279 7.4 76
32 NA 286 8.6 78
33 NA 287 9.7 74
34 NA 242 16.1 67
35 NA 186 9.2 84
36 NA 220 8.6 85
37 NA 264 14.3 79
38 29 127 9.7 82
39 NA 273 6.9 87
40 71 291 13.8 90
41 39 323 11.5 87
42 NA 259 10.9 93
43 NA 250 9.2 92
44 23 148 8.0 82
45 NA 332 13.8 80
46 NA 322 11.5 79
47 21 191 14.9 77
48 37 284 20.7 72
49 20 37 9.2 65
50 12 120 11.5 73
51 13 137 10.3 76
52 NA 150 6.3 77
53 NA 59 1.7 76
54 NA 91 4.6 76
55 NA 250 6.3 76
56 NA 135 8.0 75
57 NA 127 8.0 78
58 NA 47 10.3 73
59 NA 98 11.5 80
60 NA 31 14.9 77
61 NA 138 8.0 83
62 135 269 4.1 84
63 49 248 9.2 85
64 32 236 9.2 81
65 NA 101 10.9 84
66 64 175 4.6 83
67 40 314 10.9 83
68 77 276 5.1 88
69 97 267 6.3 92
70 97 272 5.7 92
71 85 175 7.4 89
72 NA 139 8.6 82
73 10 264 14.3 73
74 27 175 14.9 81
75 NA 291 14.9 91
76 7 48 14.3 80
77 48 260 6.9 81
78 35 274 10.3 82
79 61 285 6.3 84
80 79 187 5.1 87
81 63 220 11.5 85
82 16 7 6.9 74
83 NA 258 9.7 81
84 NA 295 11.5 82
85 80 294 8.6 86
86 108 223 8.0 85
87 20 81 8.6 82
88 52 82 12.0 86
89 82 213 7.4 88
90 50 275 7.4 86
91 64 253 7.4 83
92 59 254 9.2 81
93 39 83 6.9 81
94 9 24 13.8 81
95 16 77 7.4 82
96 78 NA 6.9 86
97 35 NA 7.4 85
98 66 NA 4.6 87
99 122 255 4.0 89
100 89 229 10.3 90
101 110 207 8.0 90
102 NA 222 8.6 92
103 NA 137 11.5 86
104 44 192 11.5 86
105 28 273 11.5 82
106 65 157 9.7 80
107 NA 64 11.5 79
108 22 71 10.3 77
109 59 51 6.3 79
110 23 115 7.4 76
111 31 244 10.9 78
112 44 190 10.3 78
113 21 259 15.5 77
114 9 36 14.3 72
115 NA 255 12.6 75
116 45 212 9.7 79
117 168 238 3.4 81
118 73 215 8.0 86
119 NA 153 5.7 88
120 76 203 9.7 97
121 118 225 2.3 94
122 84 237 6.3 96
123 85 188 6.3 94
124 96 167 6.9 91
125 78 197 5.1 92
126 73 183 2.8 93
127 91 189 4.6 93
128 47 95 7.4 87
129 32 92 15.5 84
130 20 252 10.9 80
131 23 220 10.3 78
132 21 230 10.9 75
133 24 259 9.7 73
134 44 236 14.9 81
135 21 259 15.5 76
136 28 238 6.3 77
137 9 24 10.9 71
138 13 112 11.5 71
139 46 237 6.9 78
140 18 224 13.8 67
141 13 27 10.3 76
142 24 238 10.3 68
143 16 201 8.0 82
144 13 238 12.6 64
145 23 14 9.2 71
146 36 139 10.3 81
147 7 49 10.3 69
148 14 20 16.6 63
149 30 193 6.9 70
150 NA 145 13.2 77
151 14 191 14.3 75
152 18 131 8.0 76
153 20 223 11.5 68

#各列の関係性を確認
pairs(airq2,panel=panel.smooth,lwd=2)
一般化線形モデル_pairs

#重回帰分析を実施
airq2.lm<-lm(Ozone~.,data=airq2) #結果を出力(Q-Qプロットで、分布に従っているか確認)
#Q-QプロットはX軸上に観測した累積パーセント、Y軸上に期待累積パーセントを持つグラフを表示します。
#すべての散布データが参照線に近ければ、そのデータセットは与えられた分布に従うということができます。
#qqnorm(x):xの正規確率プロット。横軸に標準正規分布のパーセンタイル、縦軸に観測値となる散布図を出力
# 正規分布に従っていれば直線に並ぶはず。
#qqline(x):xが完全に正規分布に従っている場合の予想線
qqnorm(resid(airq2.lm))# 残差の正規確率プロット(残差:目的変数の実測値と予測値の差。回帰式からどれだけ外れているかってこと)
qqline(resid(airq2.lm),col=”red”)# 残差の正規確率の参照線
一般化線形モデル_QQプロット


さて、出力されたQ-Qプロットのグラフを確認してみると、残差が正規分布に従っていないことが分かりますね。
仮に正規分布に従っているのであれば、Q-Qプロットのグラフの赤い参照線上に点が集まるはずですが、外れていますしね。


というわけで、同じデータに対して一般化線形モデルの適用を行っていきたいと思います。
一般化線形モデルの当てはまりの良さはAICという指標を利用します。
AICに関する説明はWikipediaを参照してください。

#重回帰分析の結果のAICを出力する。
#AIC:赤池情報量規準。統計モデルの良さを評価するための指標。
#値が小さいほどあてはまりが良いとされますが、
#相対的な評価として用いてるため、「○以下であることが望ましい」というような基準はありません。
AIC(airq2.lm)
[1] 998.7171

#他の分布で一般化線形モデルで回帰分析行う
#glmを利用すると、対象データの確率分布を指定して回帰分析を行える。(lmは正規分布のみ)
#binomial:二項分布。ロジティクス回帰分析などで利用。
#poisson:ポアソン分布。
#Gamma:ガンマ分布
#gaussian:正規分布
#使い方は「result<-glm(モデル式,data=対象データ,family=確率分布,link=リンク関数(省略化)」となる
#下の例では、結果を別変数に入れるのが面倒なので、無視。AICで当てはまりの良さを評価する
AIC(glm(Ozone~Solar.R+Wind+Temp,data=airq2,family=binomial))#2項分布のデータではないのでエラーになる
Error in eval(expr, envir, enclos) :
y の値は [0, 1] 中でなければなりません
AIC(glm(Ozone~Solar.R+Wind+Temp,data=airq2,family=poisson))
[1] 1344.839
AIC(glm(Ozone~Solar.R+Wind+Temp,data=airq2,family=Gamma))
[1] 939.8778
AIC(glm(Ozone~Solar.R+Wind+Temp,data=airq2,family=gaussian))#重回帰分析の結果と同じ値になる
[1] 998.7171

#当てはまりのよいモデルの結果を変数に格納
result <- glm(Ozone~Solar.R+Wind+Temp,data=airq2,family=Gamma)
summary(result)#適用後の詳細情報を参照

Call:
glm(formula = Ozone ~ Solar.R + Wind + Temp, family = Gamma,
data = airq2)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.91385 -0.44403 -0.09604 0.29322 1.25132

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.061e-01 1.530e-02 6.934 3.23e-10 ***
Solar.R -6.825e-05 1.779e-05 -3.836 0.000212 ***
Wind 1.442e-03 3.471e-04 4.156 6.55e-05 ***
Temp -9.627e-04 1.569e-04 -6.137 1.45e-08 ***

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 0.2610159)

Null deviance: 71.950 on 110 degrees of freedom
Residual deviance: 29.177 on 107 degrees of freedom
(42 observations deleted due to missingness)
AIC: 939.88

Number of Fisher Scoring iterations: 6

plot(result)#様々な情報をグラフ化してくれるので参照


一般化線形モデルの適用結果の確認

さて、一般化線形モデルの適用が終わったら、結果を確認します。
上記のようにsummary関数を利用して、適用結果の詳細を確認します。
注目すべきは「Coefficients」です。
どの項目もP値が非常に小さいため、有意な項目といえます。
「Estimate」が係数になるので、どの項目がどれだけ影響を与えているかわかりますね。
これで一般化線形モデルに当てはめたモデル式を取得できます。

最後に

一般化線形モデルは少しとっつきにくかったですが、使ってみると一般線形モデル(lm関数)よりも使い勝手がいい気がします。
こちらの方が汎用的に使用できそうですしね。
他にも便利なものがありそうですし、これからも臆せずいろいろ触ってみようと思います。


コメントを残す

メールアドレスが公開されることはありません。

CAPTCHA