ロジスティック回帰から計算される確率の信頼区間

最近忙しくて更新ができないー。

でも、ついさっきデスマコロシアムには挑戦した。
正解なら今回もRで最短のはず。ワクワク
→12/13　0：30現在最短だ～
12/15盛大に抜かれた。どうやっても勝てる気がしない。

CodeIQ | 『第8回デスマコロシアム』問題 集計報告 @tbpgr #CodeIQ - tbpg’s programming memo

論文にロジスティック回帰から計算される確率の信頼区間が載っていた。
lm()ならば、predict関数で簡単に予測値の95%CIが出せる。
しかしglm()でのロジスティック回帰などだと単純には出力されない。

今回、この論文の表から完全に元データが再現できるではないか。

こういう、”答えがある”データは有難い。ということで、再現＋内的validationの追加まで試みる
かなりやっつけ仕事。グラフのタイトルなんかもつけてないのであしからず。

General thoracic and cardiovascular surgery 10/2014; DOI: 10.1007/s11748-014-0487-6
A simple risk scoring system for predicting acute exacerbation of interstitial pneumonia after... - Abstract - Europe PubMed Central

うーむ。院外だと課金がいるのか。

この論文はSASを使用したようだが、ROC曲線とかC-indexとか、予測値の９５％CIは
一致したのでよしとする。

sco <- 0:22  #元データの再現
pt <- c(2,0,0,24,21,11,8,226,64,50,33,375,64,83,37,5,8,5,2,3,0,0,1)
inc <- c(0,0,0,0,1,0,0,5,1,6,3,43,8,13,11,2,3,0,2,1,0,0,1)
score <- NULL
ae <- NULL
for (i in 1:23){
      score <- c(score, rep(sco[i],pt[i]))
      ae <- c(ae, rep(1,inc[i]), rep(0, (pt[i] - inc[i]))) 
}

dat <- data.frame(score=score, ae=ae)

###fit
fit1 <- glm(ae ~ score, data=dat, family=binomial)
summary(fit1)

###c index
library(pROC)
ROC <- roc(ae ~ score, data=dat, ci=TRUE, direction="auto")
plot(ROC, print.thres="best", print.thres.best.method="youden", 
     print.thres.best.weights=c(1, 0.5), grid=TRUE)
coords(ROC, "all")
print(ROC)

##予測値の信頼区間　lmだとpredictで行けるのにglmだと無理
new1 <- data.frame(score=seq(0,23,by=0.1))
pred1 <- predict(fit1, newdata=new1, type="link", dispersion=NULL, se=TRUE)

plot(seq(0,23,by=0.1),plogis(pred1$fit), type="l", xlim=c(0,23), ylim=c(0,1))
lines(seq(0,23,by=0.1), plogis(pred1$fit + pred1$se.fit*1.96), type="l", lty=3, xlim=c(0,23), ylim=c(0,1))
lines(seq(0,23,by=0.1), plogis(pred1$fit - pred1$se.fit*1.96), type="l", lty=3, xlim=c(0,23), ylim=c(0,1))
abline(h=seq(0,1,by=0.1),lty=3, v=22)

####rms　HarrellのBootstrapを利用したモデルの内的妥当性のチェック
library(rms)
fit2 <- lrm(ae ~ score, data=dat, x=TRUE,y=TRUE) #logistic reg. for rms
fit2

cali <- calibrate(fit2, method="boot",  B=300, predy=seq(.01,.99,length=100)) # usually B=200 or 300
windows()
plot(cali)
summary(cali)