「統計学関連なんでもあり」の過去ログ--- 046

No.20696　多変量解析でのAUCやCIのブートストラップ　　【おおはし】　2013/12/25(Wed) 18:20

Rで多変量ロジスティック回帰分析を行い，AUCとオッズ比についてブートストラップにより信頼区間を求めることを考えておりますが，以下のようなコードで正しいでしょうか。
オッズ比については，最後の段階でexpを通すのが正しいと思われたのですが，ご教授いただけますと幸いです。warningsが出る分については無視してもかまわないでしょうか。

library(pROC)
library(ROCR)
library(boot)
library(epicalc)
library(snow)

auc_calc <- function(data, indices, formula, response) {
  d=data[indices,]
  afit <- glm(formula, data=d, family=binomial)
  afit.predict <- ROCR::prediction(fitted(afit), d[response])
  aperf <- performance(afit.predict, "auc")
  return(aperf@y.values[[1]])
}

m_oddsratio_calc <- function(data, indices, formula) {
  d=data[indices,]
  afit <- glm(formula, data=d, family=binomial)
  return(coef(afit))
}

analysis <- function(data, formula, response) {
  dpres <- glm(formula=formula, family="binomial", data=data)
  pred  <- ROCR::prediction(fitted(dpres), data[response])
  perf  <- ROCR::performance(pred,"tpr","fpr")
  M <- max(perf@y.values[[1]]-perf@x.values[[1]])
  d_cutoff <- which(perf@y.values[[1]]-perf@x.values[[1]]==M)
  cat("** Cutoff: ")
  print(perf@alpha.values[[1]][d_cutoff])
  cat("** Sensitivity: ")
  print(perf@y.values[[1]][d_cutoff])
  cat("** Specificity: ")
  print(1-perf@x.values[[1]][d_cutoff])
  perf1 <- ROCR::performance(pred,"auc")
  cat("** AUC: ")
  print(perf1@y.values[[1]])
  cat("** Bootstrapped AUC CI:\n")
  set.seed(3141592)
  boot.results <- boot(data=data, statistic=auc_calc, R=1000, formula=formula, response=response)
  print(boot.ci(boot.results,type = c("perc")))

  cat("** Odds ratio\n")
  print(logistic.display(dpres, crude.p.value=TRUE, decimal=4))

  cat("** Bootstrapped odds ratio\n")
  boot.results <- boot(data=data, statistic=m_oddsratio_calc, R=1000, formula=formula)
  or_results <- c()
  for(it in 1:length(boot.results$t0)) # n=vn variates
  {
    btor <- boot.ci(boot.results, type = c("perc"), index=it, h=exp)
    or_results <- rbind(or_results, c(btor$percent[4:5]))
  }
  or_result <- data.frame(t(or_results))
  colnames(or_results) <- c("(percentile L,", "percentile H)")
  rownames(or_results) <- names(boot.results$t0)
  print(or_results)
}

library(ISwR)
data(juul)
juul$menarche <- factor(juul$menarche, labels = c("No","Yes"))
juul.girl <- subset(juul, complete.cases(age) & complete.cases(menarche) & complete.cases(igf1))
analysis(juul.girl, menarche~age+igf1, "menarche")

No.20697　Re: 多変量解析でのAUCやCIのブートストラップ　　【青木繁伸】　2013/12/25(Wed) 18:44

warnings は
glm.fit: fitted probabilities numerically 0 or 1 occurred
ですね。
「完全分離 (complete separation) 」が起きているのでしょう。
そのような場合を含めるか含めないか，どうすればよいかは私には分かりません。

このプログラムが正しいかどうかは見ていません。
正確な答えが分かっているデータに対してこのプログラムを動かし，ほぼ同じ解が得られるかどうかを確かめればよいでしょう。

忌憚なく意見を申し上げれば，あるプログラムが正しいかどうかを検証するのは，「かったるい」。

ノウハウをお教えしましょう。
コメントのないプログラムというのは読む気がしないということ。
コメントといっても，プログラムのオウム返しのコメントは何の意味もないと言うこと。
上でいうコメントは，単なる注釈というのではなく，ここでは何を求めるためにどういうアルゴリズムで計算を進めているんですというような，プログラムをトレースする人の役に立つ注釈ということです。特に，パフォーマンスを得るために凝ったことをしているような箇所では，本質的にはどういうことをしていますというようなことの説明が必要ということです。

ああだこうだいっても，正しい結論が出ていればプログラムは正しい，そうでなければそうでない。当たり前のことです。

正しい答えが出ないとき，プログラムのどこが間違いかを指摘するより，正しいプログラムを新たに書いた方が効率が良い（モティペーションがあれば）。
プログラムが正しい答えを出すなら，より効率的なプログラムにするための指摘をするのは労多くして功少なし（それで良いじゃない）。
ということで，結論は，「答えの分かっているデータで，正しい答えが出るかどうか」に尽きるわけです。
それを判定するのは，プログラムを書いた本人。

プログラムを書けるような人ならば，人に頼らないこと。

No.20718　Re: 多変量解析でのAUCやCIのブートストラップ　　【おおはし】　2013/12/31(Tue) 15:33

青木先生

ご連絡ありがとうございます。
入力と答えの分かっているペアが入手できず，見掛け上はおかしな値ではなさそうなのですが，確証が得られず心配な事案です。
もう少し頑張ってみようと思います。

今後とも宜しくお願いいたします。