2008-12-14

ロジスティック回帰1

ロジスティック回帰(Logistic Regression)とは、各レコードが「成功」・「失敗」のような二値応答変数に対して、そのオッズ(Odds)の対数を独立変数から最尤法を使って パラメータを計算するモデルです。オッズの対数のことを、ロジット(Logit)と呼ぶので、Logisticとよびます、たぶん。

Rのような統計計算ソフトに、よい形をしたデータを放り込むと、回帰係数と定数にあたるパラメータが出力されます。そして、そのパラメータを用い て、成功確率を逆算することで、独立変数が二値応答変数に与える効果を、定量的に調べることができます。成功確率は、Logitを、独立変数の関数に変形 することで求められます。

(具体例)
メスのカブトガニが、パートナーとなるオスがいる確率を、メスのカブトガニの甲羅の幅から推定してみる。

(データ)
No. color spine width satell weight satell.flag
1 3 3 28.3 8 3050 1
2 4 3 22.5 0 1550 0
3 2 1 26.0 9 2300 1
4 4 3 24.8 0 2100 0
5 4 3 26.0 4 2600 1
6 3 3 23.8 0 2100 0


No. :レコード番号
color :甲羅の色(2-5) 2;明るい
spine :縁棘の状態(1-5) 1;正常
width :甲羅の幅(cm)
weight:kg
Satell:サテライト数(あるメスにくっついているオスの数)
Satell.flg:サテライトの有無 1;あり

(手順)

(1)まず、Rにcsvファイルを読み込む
crab.data <- read.csv("ファイルの場所") (2)どんな感じで取り込んだかを見る head(crab.data) (3)早速ロジスティック回帰 crab.result <- glm(crab.data$satell.flag~crab.data$width,family = binomial) glm(二値応答~独立変数,family=binomial) (メモ) glmは、一般化線型モデルの関数を表す familyを指定すれば、ポアソン回帰とかも出来たはず、よく分からんが (4)結果 Summary(crab.result) ==================================================== Call: glm(formula = crab.data$satell.flag ~ crab.data$width, family = binomial) Deviance Residuals: Min 1Q Median 3Q Max -2.0281 -1.0458 0.5480 0.9066 1.6941 Coefficients: Estimate Std. Error z value Pr(>|z|)
(Intercept) -12.3508 2.6287 -4.698 2.62e-06 ***
crab.data$width 0.4972 0.1017 4.887 1.02e-06 ***
---
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: 194.45 on 171 degrees of freedom
AIC: 198.45

Number of Fisher Scoring iterations: 4

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

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -12.3508 2.6287 -4.698 2.62e-06 ***
crab.data$width 0.4972 0.1017 4.887 1.02e-06 ***

Interceptsとcrab.data$widthのEstimeteが定数項(α)と関数の傾き(β)ですね


(5)当てはめ値の計算

このαとβを用いて、当てはめ値(甲羅の幅ごとの、連れ合いのオスがいる確率の理論値)を計算する。変数satel.fittedvalに代入しておく


satel.fittedval <- exp(a + b*crab.data$width)/(1 + exp(a + b*crab.data$width))

計算結果を、最初のデータと合体

crab.data <- cbind(crab.data,satel.fittedval)

(6)列crab.data$widthとcrab.data$satel.fittedvalで散布図を描く

plot(crab.data$width,crab.data$satel.fittedval)

すると、画像のグラフを得る

ちゃんとS字型になっている
当てはめ値だから、当たり前だけど

とりあえず、初歩的なところは、処理できました
ばんざーーい

(参考文献)
カテゴリカルデータ解析入門
(An Introduction to Categorical Data Analysis)
Alan Agresti
サイエンティスト社(2003)

のp.143あたりを挑戦

(感想)
Agrestiて、カテゴリカルデータ分析業界では、凄いえらい先生らしい
あと、Rも、色々できるみたいだけど、まだよく分からん
特にグラフ
よい本を探そうかな

0 件のコメント:

コメントを投稿