No.16864　Rを用いた９０％信頼区間の算出方法　　【赤岳】　2012/05/08(Tue) 11:17

いつもここで勉強させていただいています。
さて，Rを用いて混合効果モデルを行い，以下のような結果を出しましたが，3つの固定効果それぞれの推定値，SE等は出ていますが，90％信頼区間を出すコマンドの書き方がわかりません。
どなたか，ご教授いただけませんでしょうか。よろしくお願いします。
なお，「推定値±t*SE」による計算式ではなく，Rの関数を使った方法を知りたいです。

> model.Cmax<-lme(Cmax.log10~Treatment+Group+Dmpoint, random=~1|Volunteer, method="ML",data=set)
> summary(model.Cmax)
Linear mixed-effects model fit by maximum likelihood
Data: set
AIC BIC logLik
-34.2166 -28.87437 23.1083

Random effects:
Formula: ~1 | Volunteer
(Intercept) Residual
StdDev: 0.04816989 0.05230218

Fixed effects: Cmax.log10 ~ Treatment + Group + Dmpoint
Value Std.Error DF t-value p-value
(Intercept) 0.8001070 0.03605658 7 22.190317 0.0000
Treatmentpht 0.0786912 0.02813088 7 2.797324 0.0266
GroupB -0.0740434 0.04619337 7 -1.602902 0.1530
Dmpointp2 0.0045478 0.02813088 7 0.161667 0.8761
（以下，省略）
No.16865　Re: Rを用いた90％信頼区間の算出方法　　【青木繁伸】　2012/05/08(Tue) 13:50

lme R confidence interval で検索すると interval 関数のヘルプが見つかります。
fm1 <- lme(distance ~ age, Orthodont, random = ~ age | Subject)
intervals(fm1)
No.16866　Re: Rを用いた90％信頼区間の算出方法　　【赤岳】　2012/05/08(Tue) 14:49

青木先生，ご教示ありがとうございます。コマンドの設定はわかりました。
Treatmentにおける90％CIの結果は，下記のようになりました。
しかし，検算のため，式：Mean±t(0.05,n=9)*SEで90％CI求めると，一致しませんでした。
この試験は，2x2のクロスオーバー試験で，9例に2つの薬剤をクロスで投与したデザインです。この検算ではt分布のコマンドをqt(0.05,7)としましたが，これが問題なのでしょうか。
不一致の原因はなにか，9例のクロスの場合，t分布の自由度はいくらになるのか，ご教示いただきたくお願いいたします。

> intervals(model.Cmax,level=0.9)
Approximate 90% confidence intervals
Fixed effects:
lower est. upper
(Intercept) 0.73986147 0.800107023 0.860352574
Treatmentpht 0.03168838 0.078691176 0.125693972
（以下，省略）

検算
> 0.078691176-qt(0.05,7)*0.02813088
[1] 0.1319873
> 0.078691176+qt(0.05,7)*0.02813088
[1] 0.02539501
No.16867　Re: Rを用いた90％信頼区間の算出方法　　【青木繁伸】　2012/05/08(Tue) 15:48

テストデータ（あなたの本当のデータでなくて良い）を含めて，追試が出来る（コンソールにコピーするだけでちゃんと動く）形で質問してください。
No.16868　Re: Rを用いた90％信頼区間の算出方法　　【青木繁伸】　2012/05/08(Tue) 15:55

プログラム nlme:::intervals.lme をみると，
-qt((1 - level)/2, objectfixDFX) * sqrt(diag(objectvarFix))のような部分があるので，lmeが返すobjectのfixDFX 要素を見ればそれが自由度でしょう。
No.16869　Re: Rを用いた90％信頼区間の算出方法　　【赤岳】　2012/05/08(Tue) 16:47

青木先生，
失礼しました。とりあえず検討データと実施したコマンド（混合効果モデル，90%CI，90%CIの手計算）をお示しします。

> set<-read.csv("DATA.csv")
> set
Subject Treatment Group Seq Cmax
1 a Ref A p2 0.7857568
2 a New A p1 0.7727616
3 b Ref A p2 0.7764106
4 b New A p1 0.8024316
5 c Ref A p2 0.9704863
6 c New A p1 0.7907073
7 d Ref B p1 0.9389698
8 d New B p2 0.7498136
9 e Ref B p1 0.7665616
10 e New B p2 0.7079107
11 f Ref A p2 0.9895388
12 f New A p1 0.8409838
13 g Ref A p2 0.8945376
14 g New A p1 0.7936507
15 h Ref B p1 0.8358807
16 h New B p2 0.8065191
17 i Ref B p1 0.6776070
18 i New B p2 0.6582023

> model.Cmax<-lme(Cmax~Treatment+Group+Seq, random=~1|Subject, method="ML",data=set)
> summary(model.Cmax)
Linear mixed-effects model fit by maximum likelihood
Data: set
AIC BIC logLik
-34.2166 -28.87437 23.1083

Random effects:
Formula: ~1 | Subject
(Intercept) Residual
StdDev: 0.04816989 0.05230218

Fixed effects: Cmax ~ Treatment + Group + Seq
Value Std.Error DF t-value p-value
(Intercept) 0.8001070 0.03605658 7 22.190317 0.0000
TreatmentRef 0.0786912 0.02813088 7 2.797324 0.0266
GroupB -0.0740434 0.04619337 7 -1.602902 0.1530
Seqp2 0.0045478 0.02813088 7 0.161667 0.8761
（以下，省略）

> intervals(model.Cmax, level=0.9)
Approximate 90% confidence intervals

Fixed effects:
lower est. upper
(Intercept) 0.73986147 0.800107023 0.860352574
TreatmentRef 0.03168838 0.078691176 0.125693972
GroupB -0.15122615 -0.074043440 0.003139275
Seqp2 -0.04245497 0.004547826 0.051550623
attr(,"label")
[1] "Fixed effects:"
（以下，省略）

> 0.0786912+qt(0.05,7)*0.02813088 ＃手計算で上記Treatmentの解からValueとSEを持ってきています。やはり，信頼区間が手計算と関数での計算値が異なっています。
[1] 0.02539504
> 0.0786912-qt(0.05,7)*0.02813088 ＃手計算で上記Treatmentの解からValueとSEを持ってきています。
[1] 0.1319874
No.16870　Re: Rを用いた90％信頼区間の算出方法　　【青木繁伸】　2012/05/08(Tue) 17:30

プログラムを見ると，あなたの計算方法とは違うようですね。
intervals <- function (object, level = 0.95, which = c("all", "var-cov", "fixed"), 
    ...) 
{
    which <- match.arg(which)
    val <- list()
    if (which != "var-cov") {
        est <- fixef(object)
        len <- -qt((1 - level)/2, object$fixDF$X) * sqrt(diag(object$varFix))
        vfix <- array(c(est - len, est, est + len), c(length(est), 
            3), list(names(est), c("lower", "est.", "upper")))
        attr(vfix, "label") <- "Fixed effects:"
        val <- list(fixed = vfix)
    }
    :
となっています。逐次実行してみると以下のようになるでしょう。
> model.Cmax<-lme(Cmax~Treatment+Group+Seq, random=~1|Subject, method="ML",data=set)
> object <- model.Cmax # intervals 関数内では object なので，簡単のために
> level <- 0.9
> (est <- fixef(object)) # 推定値
 (Intercept) TreatmentRef       GroupB        Seqp2 
 0.800107000  0.078691185 -0.074043410  0.004547835 
> (len <- -qt((1 - level)/2, object$fixDF$X) * sqrt(diag(object$varFix))) # t * SE
 (Intercept) TreatmentRef       GroupB        Seqp2 
  0.06024554   0.04700278   0.07718270   0.04700278 
> (1 - level)/2 # α
[1] 0.05
> object$fixDF$X # 自由度
 (Intercept) TreatmentRef       GroupB        Seqp2 
           7            7            7            7 
> -qt((1 - level)/2, object$fixDF$X) # level=0.9 に対するクオンタイル値
 (Intercept) TreatmentRef       GroupB        Seqp2 
    1.894579     1.894579     1.894579     1.894579 
> object$varFix # この対角成分の平方根が SE
               (Intercept)  TreatmentRef        GroupB         Seqp2
(Intercept)   0.0010111708 -2.735517e-04 -7.376191e-04 -2.735517e-04
TreatmentRef -0.0002735517  6.154913e-04 -1.020711e-19 -6.838792e-05
GroupB       -0.0007376191 -1.020711e-19  1.659643e-03 -2.539673e-19
Seqp2        -0.0002735517 -6.838792e-05 -2.539673e-19  6.154913e-04
> sqrt(diag(object$varFix)) # 対角成分の平方根
 (Intercept) TreatmentRef       GroupB        Seqp2 
  0.03179891   0.02480910   0.04073872   0.02480910 
> ( vfix <- array(c(est - len, est, est + len), c(length(est), 
+             3), list(names(est), c("lower", "est.", "upper")))) # Est ± t * SE
                   lower         est.       upper
(Intercept)   0.73986146  0.800107000 0.860352538
TreatmentRef  0.03168840  0.078691185 0.125693968
GroupB       -0.15122611 -0.074043410 0.003139291
Seqp2        -0.04245495  0.004547835 0.051550618
> attr(vfix, "label") <- "Fixed effects:"
> (val <- list(fixed = vfix))r # これが最終結果
$fixed
                   lower         est.       uppe
(Intercept)   0.73986146  0.800107000 0.860352538
TreatmentRef  0.03168840  0.078691185 0.125693968
GroupB       -0.15122611 -0.074043410 0.003139291
Seqp2        -0.04245495  0.004547835 0.051550618
attr(,"label")
[1] "Fixed effects:"

> intervals(model.Cmax, level=0.9) # intervals 関数で返されるものを念のために
Approximate 90% confidence intervals

 Fixed effects:
                   lower         est.       upper
(Intercept)   0.73986146  0.800107000 0.860352538
TreatmentRef  0.03168840  0.078691185 0.125693968
GroupB       -0.15122611 -0.074043410 0.003139291
Seqp2        -0.04245495  0.004547835 0.051550618
attr(,"label")
[1] "Fixed effects:"

 Random Effects:
  Level: Subject 
                     lower       est.      upper
sd((Intercept)) 0.02496539 0.04816988 0.09294217

 Within-group standard error:
     lower       est.      upper 
0.03549324 0.05230217 0.07707149
No.16871　Re: Rを用いた90％信頼区間の算出方法　　【赤岳】　2012/05/08(Tue) 17:48

青木先生，
fixDFXを見ると，やはり自由度は7のようです。ご指摘の計算方法が違うということですが，sqrt(diag(objectvarFix)が解で表示されているSEとは違うということでしょうか。
すみません。ご指摘について行けておりません。
No.16872　Re: Rを用いた90％信頼区間の算出方法　　【青木繁伸】　2012/05/08(Tue) 17:54

この上のコメントに実行結果を付けました。

> sqrt(diag(object$varFix)が解で表示されているSEとは違うということでしょうか。

解で表示されているSEとは違いますね。
No.16873　Re: Rを用いた90％信頼区間の算出方法　　【赤岳】　2012/05/08(Tue) 18:18

青木先生，
ご検討いただきありがとうございます。ほんとに感謝いたします。
結局，sqrt(diag(object$varFix)は0.02480910で，これを用いると確かに手計算と一致しました。
しかし，Fixed effectの解とかなり違います。
どう理解したらよろしいのでしょうか。どちらの結果を用いるべきなのでしょうか。
ちなみにSASのOutput（私自身使ったことはありません）では，Fixed effectの解に表示されるSEが用いられて計算されているようです。
No.16874　Re: Rを用いた90％信頼区間の算出方法　　【青木繁伸】　2012/05/08(Tue) 18:32

summary で表示される SE が何なのかは，nlme:::summary.lme を見れば分かります。
    fixed <- fixef(object)
    stdFixed <- sqrt(diag(as.matrix(object$varFix)))
    object$corFixed <- array(t(object$varFix/stdFixed)/stdFixed, 
        dim(object$varFix), list(names(fixed), names(fixed)))
    if (adjustSigma && object$method == "ML") 
        stdFixed <- stdFixed * sqrt(object$dims$N/(object$dims$N - 
            length(stdFixed)))
の部分で，adjustSigma && object$method == "ML"が TRUE なので，adjust されているのです。
つまり，sqrt(objectdimsN/(objectdimsN - length(stdFixed)))が，掛けられている。今の場合だと sqrt(18/14 ) 倍されている。

> どちらの結果を用いるべきなのでしょうか。
> SASのOutput（私自身使ったことはありません）では，Fixed effectの解に表示されるSEが用いられて計算されているようです。

どちらが良いか，同じかなどは，その分析方法について，お勉強しないとわからないのでしょうね。
No.16875　Re: Rを用いた90％信頼区間の算出方法　　【赤岳】　2012/05/08(Tue) 19:03

青木先生
端的なヒントをありがとうございます。
モデル式の最尤推定法のmethodを”ML”から”REML”に変更すると，手計算の値と一致しました。
MLとREMLの違いはよくわかりませんが，methodは，REMLを指定した方が安心できることがわかりました。

> model.Cmax<-lme(Cmax~Treatment+Group+Seq, random=~1|Subject, method="REML",data=set)
> summary(model.Cmax)
（省略）
Fixed effects: Cmax ~ Treatment + Group + Seq
Value Std.Error DF t-value p-value
(Intercept) 0.8001070 0.03605659 7 22.190315 0.0000
TreatmentRef 0.0786912 0.02813088 7 2.797324 0.0266
（省略）

> intervals(model.Cmax,level=0.9)
Approximate 90% confidence intervals
Fixed effects:
lower est. upper
(Intercept) 0.73179499 0.800107023 0.86841906
TreatmentRef 0.02539502 0.078691176 0.13198733
（省略）

＃Fixed effectの解のSEを用いた手計算の結果
> 0.0786912-qt(0.05,7)*0.02813088
[1] 0.1319874
> 0.0786912+qt(0.05,7)*0.02813088
[1] 0.02539504
No.16876　Re: Rを用いた90％信頼区間の算出方法　　【青木繁伸】　2012/05/09(Wed) 07:35

> モデル式の最尤推定法のmethodを”ML”から”REML”に変更すると，手計算の値と一致しました。
> MLとREMLの違いはよくわかりませんが，methodは，REMLを指定した方が安心できることがわかりました。

method は，If "REML" the model is fit by maximizing the restricted log-likelihood. If "ML" the log-likelihood is maximized. Defaults to "REML".
No.16880　Re: Rを用いた90％信頼区間の算出方法　　【赤岳】　2012/05/09(Wed) 10:56

青木先生，

制限付き最尤法と最尤法の違いはよくわかりません。勉強してみます。
REMLがDefaultなんですね。ご教示ありがとうございました。
また，詳細に数々のご教授をいただきほんとにありがとうございました。
今後ともよろしくご指導お願いいたします。

● 「統計学関連なんでもあり」の過去ログ--- 045 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る