モデルのパラメータの信頼区間推定

MikuHatsune2013-03-27

数学いらずの医科統計学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