数学いらずの医科統計学PART7 CHAPTER34のモデルの紹介の話。
モデルをサンプルデータに適合させると、モデルのそれぞれのパラメータに対して推定値を信頼区間が得られる。
とあって、昔苦労した記憶があるのだが、drc パッケージ中にこの面倒をやってくれる
confint
という関数があったのでやってみる。
library(drc) ryegrass.m1 <- drm(rootl ~ conc, data = ryegrass, fct = LL.4()) plot(ryegrass.m1, type="all") conf0 <- confint(ryegrass.m1) confint(ryegrass.m1, "e") # "e" のラベルをつけたパラメータだけ計算する。
conf0 2.5 % 97.5 % b:(Intercept) 2.01211606 3.9523221 c:(Intercept) 0.03878752 0.9240389 d:(Intercept) 7.39961398 8.1863026 e:(Intercept) 2.67052621 3.4453837
EC50の95%信頼区間もプロットすれば
plot(ryegrass.m1, type="all") segments(conf0[4, 1], mean(ryegrass.m1$parm[2:3, ]), conf0[4, 2], lwd=4, col=2)
というわけで昔苦労した記憶に挑戦する。
res の list オブジェクトまで作ったら
lapply(res, confint) [[1]] 2.5 % 97.5 % hill:(Intercept) -2.100733 -1.567433 EC50:(Intercept) -7.000203 -6.825073 [[2]] 2.5 % 97.5 % hill:(Intercept) -2.339949 -1.642450 EC50:(Intercept) -6.167189 -5.970994
となり、Prismが計算したEC50の95%信頼区間に一致する。
しかしHill係数についてはよくわからない。たぶんパラメータの置き方によるのだと思うが、薬理学的な比較で重要なのはEC50なのでこれでよしとしておこう。
今回の信頼区間推定で重要なのは、各濃度点での観測数がばらばらでも、曲線全体の当てはめとして推定してくれるっぽいことだと思う。
というわけで気になったのでやってみるが、昔酸素解離曲線(1), (2)というのをやった。
これらでは、酸素飽和度SpO2の一組に対して、動脈酸素分圧PaO2の一組しかデータがない。けれども、曲線全体の当てはめを考えて信頼区間を推定してくれる。
sig, sig2 のオブジェクトを作成したら、
sig A 'drc' model. Call: drm(formula = SpO2 ~ PaO2, fct = L.4()) Coefficients: b:(Intercept) c:(Intercept) d:(Intercept) e:(Intercept) -0.07582 -18.87027 96.79567 21.62124
confint(sig) 2.5 % 97.5 % b:(Intercept) -0.09006235 -0.06158466 c:(Intercept) -31.14474790 -6.59578972 d:(Intercept) 94.54421684 99.04712703 e:(Intercept) 17.76837267 25.47411610
sig2 A 'drc' model. Call: drm(formula = SpO2 ~ PaO2, fct = LL.4(fixed = c(NA, NA, 100, NA))) Coefficients: b:(Intercept) c:(Intercept) e:(Intercept) -2.77761 -0.02002 26.96873
confint(sig2) 2.5 % 97.5 % b:(Intercept) -2.8483718 -2.7068436 c:(Intercept) -0.8729649 0.8329262 e:(Intercept) 26.6181168 27.3193404