ロジスティック回帰４

さて、今回は、多重ロジスティック回帰を勉強します
多重とつくと、オドロオドロしいのですが、まあ、重回帰分析みたいなものです。つまり、複数の独立変数それぞれの、目的変数に対する効果を、計算するというだけのことです。

こういうと、何だか単純な話ですね。まあ、単純な話なのですが
ただし、理論まで単純かどうかは、僕の力では分かりません。

そんなわけで、早速、分析してみましょう
今まで、カブトガニのメスにツガイのオスがいるかどうかを、カブトガニのメスの甲羅の幅から推測してきました。その結果、甲羅の幅が大きいほど、ペアとなるオス、つまり、彼氏がいる確率が高そうだ、と言うことが分かりました。

さて、ツガイのオスがいる確率は、メスの甲羅の幅だけに影響されているのでしょうか。
一般に、カブトガニのメスは、甲羅の色が明るいほど、若いそうです。
そこで、甲羅の色が、ツガイとなるオスがいる確率に依拠するかを確認してみます。人間の女の子は、若いほど価値があるような気がしますが、カブトガニでもそうなのか、確認してみましょう。

==========================================================

もう一度、データ行列を確認します

No. color spine width satell weight
1 3 3 28.3 8 3050
2 4 3 22.5 0 1550
3 2 1 26 9 2300

そうこんな感じでした。
さて、今回は、colorを、説明変数に取り込みます
ただし、colorを質的変数として扱うので、colorからダミー行列を
作成します

ダミー行列としては、こんな感じになります

No. bright mid mid.dark
1 0 1 0
2 0 0 1
3 1 0 0
4 0 0 1
5 0 0 1

colorは、2-5の値をとる変数でした
また、2が一番明るくて、数が大きくなるほど、色が暗くなります
そこで、この4つの水準を、bright,mid,mid.darkの4つの変数として再表現します。つまり、あるレコード、つまり、カブトガニの個体のcolorが、2ならは、brightが1になると言う感じで、質的変数に変換します

なお、colorの5番目の水準がないと思うかもしれません
がしかし、それは、質的変数が全てゼロの場合に相当しますので
大丈夫です

ちなみに、colorを質的変数に変換するのは、こんな感じで、処理しました。まず、ダミー行列xを使って

for (i in 1 : 4){
for(j in 1 : 173){
if(crab.data$color[j]== ( i + 1 ) ){
x[j,i] <- 1
}

else{
x[j,i] <- 0
}
}
}


color変数を質的変数に変換しておきます
その上で、ダミー行列xをこんな感じで、もとのデータ行列に足しておきましょう

crab.data <- data.frame(crab.data,x)

さて、これで準備は整いましたね
さっそく、多重ロジスティック回帰を行います

glm(formula = satell.flag ~ width + bright + mid + mid.dark,
family = binomial)

とすれば、計算できます
いつものように、結果を確認します

> summary(m.glm.result)

Call:
glm(formula = satell.flag ~ width + bright + mid + mid.dark,
family = binomial)

Deviance Residuals:
Min 1Q Median 3Q Max
-2.1124 -0.9848 0.5243 0.8513 2.1413

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -12.7151 2.7617 -4.604 4.14e-06 ***
width 0.4680 0.1055 4.434 9.26e-06 ***
bright 1.3299 0.8525 1.560 0.1188
mid 1.4023 0.5484 2.557 0.0106 *
mid.dark 1.1061 0.5921 1.868 0.0617 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 225.76 on 172 degrees of freedom
Residual deviance: 187.46 on 168 degrees of freedom
AIC: 197.46

Number of Fisher Scoring iterations: 4

なるほどって感じです
式で書けば、

logit(π) =-12.7151+ 0.4680*width+1.3299*bright+mid*1.4023+mid.dark*1.1061

です。
ちなみに

logit(π)=ln[π/(1-π)]

のように、定義されています。

要するに、甲羅の色が明るくなると、それなりに、ツガイのオスが増えるようです。ただし、明るい色の間だけでは、甲羅の色の効果は、それ程強くないようです

続いて、尤度比検定の結果を確認します

> anova(m.glm.result,test="Chisq")
Analysis of Deviance Table

Model: binomial, link: logit

Response: satell.flag

Terms added sequentially (first to last)


Df Deviance Resid. Df Resid. Dev P(>|Chi|)
NULL 172 225.759
width 1 31.306 171 194.453 2.204e-08
bright 1 0.084 170 194.368 0.771
mid 1 3.259 169 191.109 0.071
mid.dark 1 3.652 168 187.457 0.056


たしかに、甲羅の幅は、よく効いています
また、甲羅の色が明るいということ自体、つまり、bright+mid+mid.darkは、widthだけに比べても、少ないながらにも効果が ありそうです(6.995 P<0.07)。しかし、明るい色の間で効果があるかといわれれば、何だかビミョーってのが、正直な所でしょう

では、この結果を、グラフで確認します
前と同じように、予測確率を、プロットします
ただし、今回は、甲羅の色ごとに予測確率をプロットします

まず、予測確率を、こんな感じで計算しておきます

>exp(-12.715+0.4680*crab.data$width+1.3299*crab.data$bright)/(1+exp(-12.715+0.4680*crab.data$width+1.3299*crab.data$bright))

この結果を、元のCrab.dataあたりにでも、記憶させましょう
もちろん、この計算を、甲羅の色ごとにおこないますよ

さて、計算が終わったら、グラフ化します

> plot(width[width.order],fit.bright[width.order],xlim=c(21,34),ylim=c(0,1),xlab="",ylab="",type="l")

> par(new=T)

> plot(width[width.order],fit.mid[width.order],xlim=c(21,34),ylim=c(0,1),xlab="",ylab="",type="l",col="red")

> par(new=T)

> plot(width[width.order],fit.mid.dark[width.order],xlim=c(21,34),ylim=c(0,1),xlab="",ylab="",type="l",col="skyblue")

> par(new=T)
> plot(width[width.order],fit.dark[width.order],xlim=c(21,34),ylim=c(0,1),xlab="Crab.Width",ylab="Probability",type="l",col="navy")

たぶんこんな感じです
で、その結果が、上のグラフです。
たしかに、甲羅の色が明るいと、区別が付かないように見えますね

====================================================
カテゴリカルデータ解析入門
Alan Agresti
サイエンティスト社

p.171あたりに挑戦

次回は、モデル選択に挑戦します

=====================================================

さて、今回は、上に加えて、回帰木にも挑戦しました
回帰木とは、SPSSで言う所のAnswer Treeってヤツです
一般的には、決定木(Desicion Tree)と言うほうが解り易いかもしれません

その結果が、上の画像です
crab.cateとありますが、これは、甲羅の幅を、小さい方からカテゴリー変数にしたものです

あんまり、よく理解していないのですが、

>par(xpd=T)
> plot( rpart(satell.flag~width.cate+color,crab.data),unifor=T,branch=0.7,margin=0.05)
> text( rpart(satell.flag~width.cate+color,crab.data),use.n=T,all.leaves=F)
>

な感じで計算しました
画像に出ている、a,b,c,d,e,f,g,hは、甲羅の幅のカテゴリーで、

"23.24""23.25-24.25" "24.25-25.25" "25.25-26.25" "26.25-27.25" "27.25-28.25" "28.25-29.25" "29.25"

と言う風に分類されています
まあ、何となくうまくいった感じですね

結果をざっくり解釈すると、
ツガイとなるオスがいる確率が高いのは、甲羅の色が明るい＝若いメスで、かつ、甲羅の幅が広い個体。甲羅の幅が広くても、甲羅の色が暗い＝年寄りのメスは、ツガイとなるオスが少なくなる
年寄りのメスにオスがいる確率は、甲羅の幅が非常に小さいメスと同程度

こういうと怒られそうですが、
なんだか、カブトガニの世界も、人間の世界と同じですね
女の子は、若い方が価値があるのは、進化の結果なのかも知れません


===========================================
Rでまなぶデータマイニング
熊谷悦生ら
九天社(2007)

P.106あたりを挑戦

この本は、Rに慣れるにはいいけど、
計算方法学ぶには、不適切
計算原理とか、計算の応用方法は、完全無視
その意味で、計算原理を知ってる人向けだね

他の本を読め、ってことか

あと、さっそく、某データを使って、分析してたら
スゲーーー当たり前の結果が出てしまった、しらけてしまった
まあ、使い手の技量不足を露呈しましたってことですね

また、そのうち、がんばります


===========================================