本書第 6 章「GLMの応用範囲をひろげる」をもとにした。
第 3 章では観測データのあてはめ（確率分布の指定）にポアソン分布を用いたが、今回は二項分布で、二項分布を使った GLM の一つであるロジスティック回帰分析を行う。
観測データと分析
本章の観測データは以下のとおり。
> df<-read.csv("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/binomial/data4a.csv")
> summary(df)
N y x f
Min. :8 Min. :0.00 Min. : 7.660 C:50
1st Qu.:8 1st Qu.:3.00 1st Qu.: 9.338 T:50
Median :8 Median :6.00 Median : 9.965
Mean :8 Mean :5.08 Mean : 9.967
3rd Qu.:8 3rd Qu.:8.00 3rd Qu.:10.770
Max. :8 Max. :8.00 Max. :12.440
以下は本書の図 6.2 に相当する、観測データのプロット。
> theme_set(theme_gray(base_size=12, base_family="HiraKakuProN-W3"))
> p1<-ggplot(df, aes(x=x,y=y,shape=f,colour=f))+geom_point()+xlab("体サイズ x")+ylab("生存種子数 y")
> p2<-ggplot(df,aes(x=f,y=y))+geom_boxplot()
> library(grid)
> grid.newpage()
> pushViewport(viewport(layout = grid.layout(1,2)))
> print(p1,vp=viewport(layout.pos.row=1,layout.pos.col = 1))
> print(p2,vp=viewport(layout.pos.row=1,layout.pos.col = 2))
このプロット図からわかるのは
- 体サイズが大きいほど、生存種子数が多い
- 肥料を与えた方が、生存種子数が多い
この分析では「生存種子数が多い」を「生存確率が高い」とする。
ここでは、本書の解説とは逆に分析結果から示す。以下は図 6.7に相当、左図は肥料なし（赤）、右図は肥料あり（青）。描画方法は末尾の「Rコード」参照。
曲線は分析結果のもの、x軸の体サイズ 10 で比較しているように、肥料ありの方が種子数が多い（= 生存確率が高い）のは明確。
リンク関数等の詳細は次回取り上げることにして、ここでは分析結果から示す。GLM による推定値は以下とおり。
> fit<-glm(cbind(y,N-y)~x + f, data=d, family = binomial)
> fit
Coefficients:
(Intercept) x fT
-19.536 1.952 2.022
オッヅとリスク
ロジスティック関数の逆関数のロジット関数 logit(qi) より
qi が今回の「生存確率」で、qi / (1 - qi) はオッヅ（odds）で「（生存する確率）／（生存しない確率）」。
上記最後の ∝ 関係式は、切片の定数 exp(β1) を省略したもの。この式から、体サイズ xi の大きさが「1 単位増加」した影響を調べると
odds ∝ exp(1.95(xi + 1))exp(2.02fi)
∝ exp(1.95xi)exp(1.95)exp(2.02fi)
> exp(1.95)
[1] 7.028688
7.5 倍の増加が予測される。
このように「GLMの目的と特徴」で示したように、かけ算として要因の効果の解釈しやすくなっている。
「オッヅとリスク」の関係は、一般にも耳にする例は
モデル選択
各モデルを AIC のみで評価すると
> library(MASS)
> stepAIC(fit)
Df Deviance AIC
<none> 123.03 272.21
- f 1 217.17 364.35
- x 1 490.58 637.76
モデル：AIC
一定：644.4
f ：637.76
x ：364.35
x+f ：272.21
なお「一定モデル」は以下のように、モデル指定して AIC を得た。
> glm(cbind(y,N-y)~1, data=d, family = binomial)
「x + f モデル」が AIC の観点から最良となった。
「二項分布を使った GLM」に続く。
R コード
y 値の算出で、predict の結果が「生存確率」なので 8 を掛けている。
> fit<-glm(cbind(y,N-y)~x + f, data=d, family = binomial)
> df.C<-subset(d,subset = (f=="C"))
> df.T<-subset(d,subset = (f=="T"))
> xx<-seq(min(d$x),max(d$x),length=100)
> yy<- predict(fit, data.frame(x=xx,f=rep("C",100)),type="response")*8
> pA<-ggplot()+geom_point(data=df.C,aes(x=x,y=y),color="red")+ geom_line(data=data.frame(x=xx,y=yy),aes(x=x,y=y))+xlab("体サイズ x")+ylab("生存種子数 y")+geom_vline(xintercept = 10,linetype="longdash")
> yy<- predict(fit, data.frame(x=xx,f=rep("T",100)),type="response")*8
> pB<-ggplot()+geom_point(data=df.T,aes(x=x,y=y),color="blue")+ geom_line(data=data.frame(x=xx,y=yy),aes(x=x,y=y))+xlab("体サイズ x")+ylab("生存種子数 y")+geom_vline(xintercept = 10,linetype="longdash")
> grid.newpage()
> pushViewport(viewport(layout = grid.layout(1,2)))
> print(pA,vp=viewport(layout.pos.row=1,layout.pos.col = 1))
> print(pB,vp=viewport(layout.pos.row=1,layout.pos.col = 2))
リンク関数等の詳細は次回取り上げることにして、ここでは分析結果から示す。GLM による推定値は以下とおり。
> fit<-glm(cbind(y,N-y)~x + f, data=d, family = binomial)
> fit
Coefficients:
(Intercept) x fT
-19.536 1.952 2.022
オッヅとリスク
上記最後の ∝ 関係式は、切片の定数 exp(β1) を省略したもの。この式から、体サイズ xi の大きさが「1 単位増加」した影響を調べると
odds ∝ exp(1.95(xi + 1))exp(2.02fi)
∝ exp(1.95xi)exp(1.95)exp(2.02fi)
> exp(1.95)
[1] 7.028688
つまり 7 倍「生存確率」が増加する。
また「肥料の有無」では
また「肥料の有無」では
> exp(2.02)
[1] 7.538325
このように「GLMの目的と特徴」で示したように、かけ算として要因の効果の解釈しやすくなっている。
「オッヅとリスク」の関係は、一般にも耳にする例は
個人 i の生活習慣 X の効果をあらわす係数が βs であるとして、発病確率をロジスティック回帰で調べたら、その最尤推定値 βs = 1.95 が得られたとしましょう。この場合は
となるので、病気になるオッヅ比は（近似的には「リスク」）は exp 1.95 倍（7倍）ぐらいと見つもられます。P.125
モデル選択
各モデルを AIC のみで評価すると
> library(MASS)
> stepAIC(fit)
Df Deviance AIC
<none> 123.03 272.21
- f 1 217.17 364.35
- x 1 490.58 637.76
モデル：AIC
一定：644.4
f ：637.76
x ：364.35
x+f ：272.21
なお「一定モデル」は以下のように、モデル指定して AIC を得た。
> glm(cbind(y,N-y)~1, data=d, family = binomial)
「二項分布を使った GLM」に続く。
R コード
y 値の算出で、predict の結果が「生存確率」なので 8 を掛けている。
> fit<-glm(cbind(y,N-y)~x + f, data=d, family = binomial)
> df.C<-subset(d,subset = (f=="C"))
> df.T<-subset(d,subset = (f=="T"))
> xx<-seq(min(d$x),max(d$x),length=100)
> yy<- predict(fit, data.frame(x=xx,f=rep("C",100)),type="response")*8
> pA<-ggplot()+geom_point(data=df.C,aes(x=x,y=y),color="red")+ geom_line(data=data.frame(x=xx,y=yy),aes(x=x,y=y))+xlab("体サイズ x")+ylab("生存種子数 y")+geom_vline(xintercept = 10,linetype="longdash")
> yy<- predict(fit, data.frame(x=xx,f=rep("T",100)),type="response")*8
> pB<-ggplot()+geom_point(data=df.T,aes(x=x,y=y),color="blue")+ geom_line(data=data.frame(x=xx,y=yy),aes(x=x,y=y))+xlab("体サイズ x")+ylab("生存種子数 y")+geom_vline(xintercept = 10,linetype="longdash")
> grid.newpage()
> pushViewport(viewport(layout = grid.layout(1,2)))
> print(pA,vp=viewport(layout.pos.row=1,layout.pos.col = 1))
> print(pB,vp=viewport(layout.pos.row=1,layout.pos.col = 2))
/*Tex
*/
0 件のコメント:
コメントを投稿