と査読で言われた。
状況設定としてはこんな感じである。いま、とある検査項目 （連続量）について、予後を予測するマーカーになるか検討したい。ここで、デザインとしてはとりあえず研究を立ち上げてみました、というような後ろ向き観察研究で、予後転帰が確定している症例について、診断確定時の検査値を取ってきた、という状況である。ここで、サンプル数は予後が良かった群も悪かった群も20 くらいである。
以外にも、40個くらいの項目がとってある。というのも、これらは「取れるから取った」だけの値であり、これらによって予後が予測できるかというと、3個くらいは、確実に予後に関連していると言われているものがあるが、ほかはあっ（察し）という感じである。
注目している検査項目 が、ROCAUC 0.75 くらいだった、としよう。他の40個近くの検査項目は、患者背景のばらつきを検定してください、と言われているので、何も考えずtでもwilcox でもカイ自乗でもいいので検定するのだが、無作為抽出・割付した研究ではないので、当然、いくつかは有意差がある。
査読者が言うには、これらの有意な検査項目と、新規マーカー を比較して、
が有意にマーカーなり得るか検討しろ、ということらしい。
いまいち腑に落ちなかったのでシミュレーションで考えてみる。
いま、査読者が言っていることで疑問に思っているのが、ある検査項目 に着目した時、1番目の群の検査結果
と対照となる群
について
・2群間で分布が異なっていること（患者背景で有意差があるかどうか、について言っている）
・分布が異なるとき、2群の判別性能をROC で評価するとどうなるのか
ということでいいと思われる。
まず、 と
には差がない（と神様だけが知っている）とき、帰無仮説検定では有意水準
だけ、有意な結果を生じる。
どちらも平均0、分散1 の正規分布から、2群分適当にランダムサンプリングでデータを作り、この標本について、t検定とROC 解析をしてp値とAUC の値を記録する。
差がないときのシミュレーション結果は、図の緑点である。p=0.05 のところにタテ線を引いた。
2群には差がないとわかっているにもかかわらず、偶然にもp<0.05 になり有意になってしまった結果については、ROCAUC は0.7 相当ある。
有意にならなかった場合でも、ROCAUC は0.6くらいは出てしまう。
さて、査読者はご丁寧にも、「ROC 解析のAUC は、0.7 程度で有用な判別性能があると言われています」というようなご高説を賜ってくださったので、これから鑑みるに、2群の間に差がないデータでも、有用な判別性能を持つ場合があるので、患者背景で大量に検定した結果についてROC 解析をやるべき、と言っているのであろう()
次に、 と
には差がある（と神様だけが知っている）ときを考える。ここで、2群には1 の差がある（両群とも分散1）と設定したとき、検出力は0.84 程度になることが計算でわかる。
このときの結果が図の青点である。差がない場合よりも、小さいp 値が出やすいので青点は左に偏っている。
しかし、AUC については0.7 程度が最も密度が高い。しかも、差がない緑に比べてAUC が高い傾向にあるのかというと、そうでもないのでAUC に対して効果量は影響しなさそう。
となると、帰無仮説検定を絨毯爆撃して有意な項目についてROC 解析をすると、AUC 0.7 のような有用() な結果になってしまいがちな感じがするのだが、査読者がそうのたまう()のでやった()
library(pROC) N <- c(20, 18) alpha <- 0.05 niter <- 10000 # 差がない res <- matrix(0, niter, 2) for(i in seq(niter)){ g <- unlist(mapply(rep, 1:2, N)) d <- unlist(mapply(rnorm, n=N, m=0)) res[i, 1] <- t.test(d ~ g)$p.value res[i, 2] <- roc(g, d)$auc } # 差がある res1 <- matrix(0, niter, 2) M <- c(1, 0) for(i in seq(niter)){ g <- unlist(mapply(rep, 1:2, N)) d <- unlist(mapply(rnorm, n=N, m=M)) res1[i, 1] <- t.test(d ~ g)$p.value res1[i, 2] <- roc(g, d)$auc } mean(res1[,1] < alpha) r <- rbind(res, res1) cols <- c(rgb(0, 1, 0, 0.3), rgb(0, 0, 1, 0.3)) xl <- c(0, 1) yl <- range(r[,2], 1) plot(r, col=rep(cols, each=niter), pch=16, cex=0.7, xlab="p value", ylab="AUC", xlim=xl, ylim=yl) abline(v=alpha, h=0.5, lty=3) legend("topright", legend=c(expression(Delta==0), expression(Delta==1)), col=cols, pch=16, cex=2)
