用量反応曲線を描いているのだが、複数の曲線からそれぞれEC50を得たとき、比較したい。
よくある(?)手法として、EC50のSDを求めて、それが曲線間でかぶる範囲があるかないかで差を判定する、という手法があるらしいのだが、これが頭の中で物議を醸していた（のだが放置していた）。
用量反応曲線をやり始めてから、この問題に直面したときに放置しておいたのだが、用量反応曲線について聞かれることがあったのでちょっと調べた。
多くの研究者はGraphPad Prismというソフトウェアを使って描くらしいのだが、いかんせん有料。
すると、Analyzing Dose-Response DataというPDF説明書を見つけたので、これを使ってみる。
実測データを0 ~ 100に変換するまでは省略。
変換データはこちら。
各濃度で3回実験を行なっている。欠損値はNAとしている。
これを地味に入力した。
conc <- c(10, 8, 7.523, 7, 6.523, 6, 5.523, 5, 4.523, 4, 3.523) * -1 #濃度は既に対数スケール resp <- rbind( c(0, 0, 0, 0, 0, 0), c(2.831, 0, 9.524, NA, NA, NA), c(28.571, 23.81, 35.714, 2.804, 0, 5.807), c(45.238, 42.857, 52.381, 14.019, 8.411, 16.822), c(66.667, 71.429, 66.667, 22.43, 28.037, 33.645), c(76.19, 83.333, 83.333, 47.664, 56.075, 64.486), c(83.333, 97.619, 90.476, 64.486, 78.505, 81.308), c(92.857, 102.381, 104.762, 75.701, 86.916, 95.327), c(NA, NA, NA, 84.112, 95.327, 106.542), c(92.857, 102.381, 104.762, 89.720, 103.738, 106.542), c(NA, NA, NA, 89.72, 103.738, NA) ) resp_a <- c(resp[,1:3]) #阻害剤が入っていない実験 resp_b <- c(resp[,4:6]) #阻害剤が入っている実験
いつものように4パラメータのロジスティック推定を行う。
濃度は既に対数変換しているので
fct = L.4(fixed = c(NA, 0, 100, NA))
関数で、最大・最小漸近線はそれぞれ100と0に固定する。
実を言うと、実験を複数回していても、回帰するときに引数formulaを
resp ~ conc
で渡せば、
robust
a character string specifying the rho function for robust estimation. Default is non-robust least squares estimation ("mean"). Available robust methods are: median estimation ("median"), least median of squares ("lms"), least trimmed squares ("lts"), metric trimming ("trimmed"), metric winsorizing ("winsor") and Tukey's biweight ("tukey").
で勝手に計算してくれたらしい。
conc <- rep(conc, 3) #respに対応する濃度ベクトルにする。 res <- list(0, 0) #推定を格納する。 res[[1]] <- drm(resp_a ~ conc, na.action=na.omit, fct=L.4(fixed=c(NA, 0, 100,NA), names=c("hill", "", "", "EC50"))) res[[2]] <- drm(resp_b ~ conc, na.action=na.omit, fct=L.4(fixed=c(NA, 0, 100,NA), names=c("hill", "", "", "EC50"))) EC50 <- matrix(0, nr=2, nc=2) #EC50とSDを格納する。 for(i in 1:2){ EC50[, i] <- summary(res[[i]])$coefficients[2, 1:2] } dimnames(EC50) <- list(c("EC50", "std"), c(colnames(resp))) xlim <- range(conc) ylim <- range(resp, na.rm=TRUE) #typeを指定しなければ、robustで指定したデータだけがプロットされる。 #曲線だけ描きたければ、"none" plot(res[[1]], log="", xlim=xlim, ylim=ylim, col=1, pch=1, type="all") par(new=TRUE) plot(res[[2]], log="", xlim=xlim, ylim=ylim, col=2, pch=2, type="all") #SD幅を描く。 arrows(EC50[1, 1] - EC50[2, 1], 50, EC50[1, 1] + EC50[2, 1], 50, length=0, col=1) arrows(EC50[1, 2] - EC50[2, 2], 50, EC50[1, 2] + EC50[2, 2], 50, length=0, col=2)
EC50
[,1] [,2]
EC50 -6.91263828 -6.06909176
std 0.04251684 0.04780971
プリズムが返す推定値
一致しそうだ。
PDFによると
Users frequently object to the fact that the standard error for EC50 is not reported. These standard error values express the uncertainty inherent in estimating the parameters by curve fitting. For the equation fitted here, the midpoint of the curve is log EC50, not EC50. Don’t try to express the standard error for EC50 obtained in this way using a single number―in particular, it is not the antilog of the standard error of log EC50. We suggest that you use the 95-percent confidence interval for EC50, which Prism does report.
普通、この手の解析でSDは返って来ないのだけれど、Prismなら返してくれるんだぜ(ﾄﾞﾔｧ 的なことが書いてある。
最後の一文は重要で、比較したいのならPrismが返してくれる95%信頼区間を使え、とある。
エラーバーがかぶるかかぶらないかの統計学的有意差判定は、20111212MIKUセミナーでも少し触れた。そもそも、SDがかぶらないから有意という判定は正直なところ、誰が言い出したのかは不明である。
SDで比較ができるのもよくわからない。
さて、"drc"の結果からp値なるものを取り出そうとしたことは多々あるが、できていないし、2群の比較となるとよくわからない。
しかし、EC50（言うなれば平均）、SD、実験回数（サンプル数）がわかっていれば、ゴリ押しで信頼区間は求められる。
n <- 3 #サンプル数 ciw <- qt(0.975, n) * EC50["std", ] / sqrt(n) #信頼区間の幅 EC50["EC50",] - ciw #信頼区間下限 EC50["EC50",] + ciw #信頼区間上限
[1] -6.990758 -6.156937 [1] -6.834518 -5.981247
となり、Prismの結果と一致する。
ここで、
conc <- 10 ^ conc fct = LL.4(fixed = c(NA, 0, 100, NA))
とし、ログーロジスティック推定をしたとすると、
EC50
[,1] [,2]
EC50 1.222823e-07 8.528973e-07
std 1.197101e-08 9.389035e-08
となり、SDの値が異なる。
ただし、プロットおよびSD幅は変わらない。
このまま、95%信頼区間を出しにかかれば、
[1] 1.002869e-07 6.803844e-07 [1] 1.442777e-07 1.025410e-06
となり、対数を取っていない信頼区間を求められる。
これはPrismの値と一致してそうだが、おそらく、計算の都合か何かで若干異なっている。
