Lesson7 Generalized Liner Modeling
Kosugitti
2016年12月02日
毎回やる儀式
- インターネットに接続しておいてください（このサイトが見えていればOK）。
- Rstudioを起動し，この授業のプロジェクトを開いてください。(File>Open Project)
- File>New File>R scriptで新しいコードファイルを作ってください。
野球のデータを読み込んで，バッターだけのデータを取り出す，というのをやっておきましょう。
baseball <- read.csv("baseball2016.csv")
# 投手かどうか
baseball$Pos2 <- ifelse(baseball$position=="投手",1,2)
baseball$Pos2 <- factor(baseball$Pos2,labels=c("投手","野手"))
b.data <- subset(baseball,baseball$Pos2=="野手")ただの回帰じゃ，ダメなのよ
今までの話で，潜在変数を含んだパスモデルを書くことで様々な表現ができることがわかりましたが，これにはとある制限があります。それは，従属変数が連続的で正規分布していること，というものです。
library(ggplot2)
p <- ggplot(data=data.frame(X=c(-4,4)), aes(x=X))
p <- p + stat_function(fun=dnorm)
p <- p + stat_function(fun=dnorm, args=list(mean=2, sd=.5))
p
正規分布は，平均と標準偏差で様々な形が表現できますが，データが全て正規分布しているわけではありません。 例えばバッターのデータでも，こんな分布はどうでしょうか。
g <- ggplot(data=b.data,aes(x=HR)) + geom_histogram(binwidth=0.8)
g
mean(b.data$HR)## [1] 10.90123sd(b.data$HR)## [1] 10.81042ホームランは平均10，標準偏差9ぐらいですが，マイナス本数を打つことがありえないので正規分布とは言えません。 ですので，ホームラン本数を従属変数にする回帰分析は不適切な仮定を置いていることになります。
HR.reg1 <- lm(HR~height+weight,data=b.data)
summary(HR.reg1)##
## Call:
## lm(formula = HR ~ height + weight, data = b.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.045 -4.868 -1.488 2.655 32.906
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -75.1901 42.1345 -1.785 0.078227 .
## height 0.2383 0.2785 0.856 0.394690
## weight 0.4919 0.1371 3.588 0.000579 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.882 on 78 degrees of freedom
## Multiple R-squared: 0.3419, Adjusted R-squared: 0.325
## F-statistic: 20.26 on 2 and 78 DF, p-value: 8.203e-08g2 <- ggplot(data=b.data,aes(x=HR.reg1$model$HR,y=HR.reg1$fitted.values))
g2 <- g2 + geom_point() + geom_smooth(method="lm")
g2
こうした時，分布を様々な形に変えられるように線形モデルを一般化することができます。これが 一般化線形モデリング です。
心理学においては，心理変数の多くが正規分布するため，正規分布の例をあまり考えてきませんでした。なので，明らかに正規分布していない場合でも正規分布を仮定したモデルを当てはめ，適合度が低いと嘆いていることすらありました。少し補正を考えたとしても，友人の数や反応時間のように歪んだ分布の場合は，従属変数の値を対数変換して使うとか,割合のデータのように0から1の間に限定されているものは，従属変数の値を角変換して使う，といった手法が一般的でした。しかしこれからは，分布の形を見極めて，適した分布のモデルを使うのがスタンダードになってきます。皆さんも自分の予測したいデータの「分布の形」をしっかり考えて（そのためには散布図を書くことが必須），モデルを当てはめてください。
様々な分布
分布には様々なものがありますが，以下のようなものを覚えておけばいいでしょう。familyが分布を，linkはそれを線形モデルに書き換えるためのつなげ方を表現しています。
| Family | Link | 分布 | 使用例 |
|---|---|---|---|
| gaussian | idenitity | 正規分布 | 身長や成績など生態学的データ |
| binomial | logit | 二項分布 | 成功率(失敗・成功のどちらかなど) |
| poisson | logit | ポアソン分布 | 整数で数えられるデータ |
| Gamma | log | ガンマ分布 | スペルミスが出てくるまでの語数など |
分布については第三回広島ベイズ塾の資料 などを参考にしてください。
一般化線形モデリングをやってみる
ロジスティック回帰分析
結果が0/1などの二値変数であれば，二項分布を使った回帰分析をしますが，これは特に ロジスティック回帰分析 と呼ばれています。 ロジスティック関数は次のような形をしています。
x <- seq(-4,4,0.01)
plot(x,exp(x)/(1+exp(x)),type="l")
縦軸が0から1になっていることに注目してください。この関数は0から1までの滑らかなカーブを描いているわけです。だから二値変数に当てはめやすいのですね。
ロジスティック回帰分析はこのように，従属変数が名義尺度水準であるときに使われる回帰分析であると言えます。 実際の例を見てみましょう。MASSパッケージに入っているbirthwtデータを使います。
## low age lwt race
## Min. :0.0000 Min. :14.00 Min. : 80.0 Min. :1.000
## 1st Qu.:0.0000 1st Qu.:19.00 1st Qu.:110.0 1st Qu.:1.000
## Median :0.0000 Median :23.00 Median :121.0 Median :1.000
## Mean :0.3122 Mean :23.24 Mean :129.8 Mean :1.847
## 3rd Qu.:1.0000 3rd Qu.:26.00 3rd Qu.:140.0 3rd Qu.:3.000
## Max. :1.0000 Max. :45.00 Max. :250.0 Max. :3.000
## smoke ptl ht ui
## Min. :0.0000 Min. :0.0000 Min. :0.00000 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:0.0000
## Median :0.0000 Median :0.0000 Median :0.00000 Median :0.0000
## Mean :0.3915 Mean :0.1958 Mean :0.06349 Mean :0.1481
## 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:0.00000 3rd Qu.:0.0000
## Max. :1.0000 Max. :3.0000 Max. :1.00000 Max. :1.0000
## ftv bwt
## Min. :0.0000 Min. : 709
## 1st Qu.:0.0000 1st Qu.:2414
## Median :0.0000 Median :2977
## Mean :0.7937 Mean :2945
## 3rd Qu.:1.0000 3rd Qu.:3487
## Max. :6.0000 Max. :4990このデータは、1986年マサチューセッツ州で得られた，低体重児が出産される原因を調査するデータです。 変数lowが1になっているのがTRUE=低体重児，0が健常児です。説明に使われる他の変数は次の通り。
- 母親の年齢(age)
- 母親の体重(lwt)
- 人種(race,1=white,2=black,3=other)
- 喫煙するかどうか(smoke)
- 早産の回数(plt)
- 高血圧かどうか(ht)
- 子宮過敏症の有無(ui)
- 妊娠初期に検査した回数(ftv)
- 生まれた時の体重(bwt)
です。
ロジスティック回帰の実際
低体重かどうかが0/1のデータですので，これを従属変数にし，母親の体重を独立変数にした回帰を考えます。 通常の回帰は次のようになります。
result.lm <- lm(low~lwt,data=birthwt)
summary(result.lm)##
## Call:
## lm(formula = low ~ lwt, data = birthwt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.4277 -0.3375 -0.2859 0.6136 0.8687
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.646733 0.146012 4.429 1.61e-05 ***
## lwt -0.002577 0.001095 -2.354 0.0196 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4591 on 187 degrees of freedom
## Multiple R-squared: 0.02877, Adjusted R-squared: 0.02358
## F-statistic: 5.54 on 1 and 187 DF, p-value: 0.01962plot(birthwt$lwt,birthwt$low)
par(new=T)
plot(birthwt$lwt,result.lm$fitted.value,col="blue")
従属変数が二値ですから，線形回帰が当てはまるはずもありません。では次にロジスティック回帰分析をしてみます。
result.glm <- glm(low~lwt,data=birthwt,family=binomial)
summary(result.glm)##
## Call:
## glm(formula = low ~ lwt, family = binomial, data = birthwt)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0951 -0.9022 -0.8018 1.3609 1.9821
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.99831 0.78529 1.271 0.2036
## lwt -0.01406 0.00617 -2.279 0.0227 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 228.69 on 187 degrees of freedom
## AIC: 232.69
##
## Number of Fisher Scoring iterations: 4plot(birthwt$lwt,birthwt$low)
par(new=T)
plot(birthwt$lwt,result.glm$fitted.value,col="red")
幾分，滑らかな当てはまりになりました。
出力結果が0/1の間にありますので，予測値は「確率」のように解釈することができます。予測値の一部を見てみますと，それぞれが低体重児を出産する確率であると考えられます。
head(result.glm$fitted.value)## 85 86 87 88 89 91
## 0.1736052 0.2349235 0.3827710 0.3728574 0.3761505 0.3219314また，回帰係数はそのまま解釈するのではなく，ロジスティック関数に入れて考える必要があります。次のように入力しましょう。
exp(result.glm$coefficients)## (Intercept) lwt
## 2.7137035 0.9860401ここから，体重が一単位(1ポンド=454g=0.45kg)増えると，低体重児が生まれる確率が0.98倍になる，と解釈できます(こうした確率の比率のことを特に オッズ比 と言います)。
使い方は回帰分析と同じ
説明変数が複数になっても，基本的な考え方は回帰分析と同じです。
result.glm2 <- glm(low~age+lwt+smoke+ptl+ht+ui,data=birthwt,family=binomial)
summary(result.glm2)##
## Call:
## glm(formula = low ~ age + lwt + smoke + ptl + ht + ui, family = binomial,
## data = birthwt)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0765 -0.8086 -0.6246 1.0280 2.0222
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.381863 1.088914 1.269 0.20443
## age -0.042226 0.034585 -1.221 0.22212
## lwt -0.014318 0.006654 -2.152 0.03141 *
## smoke 0.550765 0.343648 1.603 0.10900
## ptl 0.593158 0.348433 1.702 0.08869 .
## ht 1.863640 0.686376 2.715 0.00662 **
## ui 0.736751 0.456508 1.614 0.10655
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 208.77 on 182 degrees of freedom
## AIC: 222.77
##
## Number of Fisher Scoring iterations: 4exp(result.glm2$coefficients)## (Intercept) age lwt smoke ptl ht
## 3.9823150 0.9586532 0.9857836 1.7345794 1.8096941 6.4471597
## ui
## 2.0891364result.glm3 <- step(result.glm2)## Start: AIC=222.77
## low ~ age + lwt + smoke + ptl + ht + ui
##
## Df Deviance AIC
## - age 1 210.31 222.31
## <none> 208.77 222.77
## - ui 1 211.33 223.33
## - smoke 1 211.33 223.33
## - ptl 1 211.78 223.78
## - lwt 1 213.97 225.97
## - ht 1 216.54 228.54
##
## Step: AIC=222.31
## low ~ lwt + smoke + ptl + ht + ui
##
## Df Deviance AIC
## <none> 210.31 222.31
## - ptl 1 212.83 222.83
## - smoke 1 213.01 223.01
## - ui 1 213.15 223.15
## - lwt 1 216.63 226.63
## - ht 1 218.45 228.45summary(result.glm3)##
## Call:
## glm(formula = low ~ lwt + smoke + ptl + ht + ui, family = binomial,
## data = birthwt)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0738 -0.7877 -0.6416 1.0657 1.9836
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.562843 0.859815 0.655 0.51272
## lwt -0.015493 0.006597 -2.348 0.01886 *
## smoke 0.563972 0.342368 1.647 0.09950 .
## ptl 0.533933 0.341417 1.564 0.11785
## ht 1.905592 0.685990 2.778 0.00547 **
## ui 0.769617 0.452910 1.699 0.08927 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 210.31 on 183 degrees of freedom
## AIC: 222.31
##
## Number of Fisher Scoring iterations: 4統計的に有意な変数は何か，それぞれの係数がどの程度オッズを変えるのか，モデルの適合度はどうか(AIC)などをみて考えることができます。変数選択をすることもできます。
ポアソン分布
もう一つ，別の分布の例をあげましょう。ポアソン分布を取り上げます。
ポアソン分布は整数で数えられるデータ，いわゆるカウントデータに使えます。バッターのデータの場合，ホームランの本数などがそうです。
HR.reg2 <- glm(HR ~ height + weight,data=b.data,family=poisson)
HR.reg2##
## Call: glm(formula = HR ~ height + weight, family = poisson, data = b.data)
##
## Coefficients:
## (Intercept) height weight
## -2.387608 0.004253 0.044250
##
## Degrees of Freedom: 80 Total (i.e. Null); 78 Residual
## Null Deviance: 827.4
## Residual Deviance: 565.4 AIC: 859.5g3 <- ggplot(data=b.data,aes(x=HR.reg1$model$HR,y=HR.reg2$fitted.values))
g3 <- g2 + geom_point() + geom_smooth(method="lm")
g3
ここでも，係数はlogit変換されていますから，解釈するためにはexp関数で元に戻す必要があります。
exp(HR.reg2$coefficients)## (Intercept) height weight
## 0.09184911 1.00426173 1.04524385体重が1単位増えるとHRが1本増える，というところでしょうか。
落穂拾い
- 0がとても多いカウントデータの場合は「0過剰ポアソン分布」を，0がないカウントデータの場合は「0切断ポアソン分布」があります。関数はパッケージを使うので説明は割愛します
- ロジスティック回帰は0/1のデータですが，複数の名義尺度水準変数に対しては 多項ロジスティック を，順序カテゴリカルな変数に対しては 順序ロジスティック回帰分析 を当てはめます。そういうものがあります。
- 今回は回帰分析の話だけでしたが，SEMのモデルでも従属変数がどのような分布をするか，指定することができます。
- 大事なのは 分布を意識して分析する ということです。
本日の課題
- ピッチャーのデータを作り，ピッチャーが打たれたHRの数を従属変数にしたGLMをやってみてください。独立変数はなんでも良いいですが，データにふさわしい分布族を選ぶこと。