こにしき(言葉・日本社会・教育)

関西学院大学(2016.04~)の寺沢拓敬のブログです(専門:言語社会学)。

独立変数にダミー変数が含まれるロジスティック回帰モデルから平均生起率を取り出すRスクリプト

スクリプトを走らせた結果

> data(infert) #自然・人工流産後の不妊症に関するデータ
> 
> #まず、不妊症の有無の予測モデルを作成
> res <- glm(case ~ age + parity + education + spontaneous + induced,
+             data=infert, family="binomial")
> summary(res)

Call:
glm(formula = case ~ age + parity + education + spontaneous + 
    induced, family = "binomial", data = infert)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.7603  -0.8162  -0.4956   0.8349   2.6536  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -1.14924    1.41220  -0.814   0.4158    
age               0.03958    0.03120   1.269   0.2046    
parity           -0.82828    0.19649  -4.215 2.49e-05 ***
education6-11yrs -1.04424    0.79255  -1.318   0.1876    
education12+ yrs -1.40321    0.83416  -1.682   0.0925 .  
spontaneous       2.04591    0.31016   6.596 4.21e-11 ***
induced           1.28876    0.30146   4.275 1.91e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 316.17  on 247  degrees of freedom
Residual deviance: 257.80  on 241  degrees of freedom
AIC: 271.8

Number of Fisher Scoring iterations: 4

> 
> #以下で何度も使うものを取り出す
> b <- res$coef  #係数
> M <- res$model #分析に使用されたデータフレーム(欠損値除外済み)
> 
> 
> ### ダミー変数化された名義尺度から平均値をとりだす ###
> 
> #まず取り出したい変数の名前を入れる。ここでは education
> Var <- "education"
> x <- cbind( 
+     table(M[,Var])/nrow(M)      #元データの分布(シェア)
+    ,c(0,b[grep(Var,names(b))])  #係数
+        )
> #各行ごとに積を出して、全部足す
> b1 <- sum(apply(x,1,prod),na.rm=T)
> 
> 
> #たとえば、他を一定としたときの age による変動
> 
> b.mean <- b1 +              # 上で算出したダミー変数の平均値           
+      b["(Intercept)"]     + #係数×平均値。以下同じ
+      b["parity"]*mean(M$parity) +
+      b["spontaneous"]*mean(M$spontaneous) +
+      b["induced"]*mean(M$induced)
>  
> x.age <- c(min(age):max(age))
> b.predi <- b.mean + b["age"]*x.age
> 
> p.predi <- 1/(1 + exp(-1*b.predi))
> 
> ### 描画
> plot(x.age,100*p.predi,type="p")
>  grid()