> 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()