「統計学関連なんでもあり」の過去ログ--- 046

No.20931　lmeのグループごとのestimate　　【赤岳】　2014/03/09(Sun) 17:40

いつもお世話になっております。
混合効果モデルにより，treatment0と1の経時的変化（Slope=time*treatment)を比較しました。
さらに，treatmentごとのSlopeの推定値±SEを抽出したいのですが，コマンドが判らず，困っています。
treatment0のSlopeとSEは，下記のtimeに当たる -1.10331±0.296558と思います。
treatment1のSlopeは，-1.10331＋0.48221= -0.62111だと思いますが，SEは出てきません。自分で計算しないといけないのでしょうか。

どなたかご教示いただけませんか。
念のため，コマンドとOutput，sample dataを示します。（画面が長くなりすみません）

> Dat<-read.csv("data.sample.csv")
> model <- lme(effect ~ time + treatment + time*treatment, random= ~1|subject, method="REML", data=Dat)
> summary(model)
Linear mixed-effects model fit by REML
Data: Dat
AIC BIC logLik
407.9816 419.3326 -197.9908

Random effects:
Formula: ~1 | subject
(Intercept) Residual
StdDev: 42.29267 6.873126

Fixed effects: effect ~ time + treatment + time * treatment
Value Std.Error DF t-value p-value
(Intercept) 141.55022 13.208265 39 10.716791 0.0000
time -1.10331 0.296558 39 -3.720388 0.0006
treatment -1.56119 8.598637 39 -0.181563 0.8569
time:treatment 0.48221 0.467510 39 1.031447 0.3087

> Dat
subject time effect treatment
1 4446 0 172.20 1
2 4446 3 163.80 1
3 4446 6 166.95 1
4 4446 9 172.20 1
5 4446 12 163.80 1
6 4420 0 131.25 0
7 4420 3 124.95 0
8 4420 6 121.80 0
9 4420 9 116.55 0
10 4420 12 116.55 0
11 4404 0 102.90 1
12 4404 3 113.40 1
13 4404 6 78.75 1
14 4404 9 93.45 1
15 4404 12 99.75 1
16 4391 0 76.65 0
17 4391 6 71.40 0
18 4391 9 61.95 0
19 4390 0 100.80 1
20 4390 3 119.70 1
21 4390 6 123.90 1
22 4390 9 105.00 1
23 4390 12 98.70 1
24 4378 0 122.85 0
25 4378 3 129.15 0
26 4378 6 137.55 0
27 4378 9 120.75 0
28 4378 12 118.65 0
29 4348 0 202.65 0
30 4348 3 198.45 0
31 4348 6 199.50 0
32 4348 9 199.50 0
33 4348 12 187.95 0
34 4347 0 183.75 0
35 4347 3 185.85 0
36 4347 6 174.30 0
37 4347 9 181.65 0
38 4347 12 176.40 1
39 4337 0 186.90 1
40 4337 3 196.35 1
41 4337 6 194.25 1
42 4337 9 191.10 1
43 4337 12 179.55 1
44 4334 0 143.85 0
45 4334 3 134.40 0
46 4334 6 128.10 0
47 4334 9 131.25 0
48 4334 12 119.70 0
49 4328 0 111.30 0
50 4328 3 96.60 0
51 4328 6 95.55 0
52 4328 9 97.65 0
53 4328 12 94.50 0

No.20938　Re: lmeのグループごとのestimate　　【赤岳】　2014/03/10(Mon) 18:43

Groupごとのestimateを求めるコードがないとすると，自力でtreatment1の分散を求まなければなりません。
先の内容を整理しました。
今回のモデルを，ランダム効果を除いたモデル式で書くと
Y(time,treatment) = μ + β1*time + β2*treatment + β3*time*treatment
となるかと思います。
treatmentは，0（プラセボ）又は1（実薬）なので，
Y(time,treatment=0) = μ + β1*time
Y(time,treatment=1) = μ + β1*time + β2 + β3*time
後者を整理すれば
Y(time,treatment=1) = (μ+β2) + (β1+β3)*time
すなわち，
treatmen=0のときは，切片μ，Slope=β1の直線
treatmen=1のときは，切片μ＋β2，Slope=β1+β3の直線，
となります。
これらの値は先に記載したようにデフォルトの結果を使って求まりますが，（β1＋β3）の分散が厄介です。
Var(β1+β3)=Var(β1)+Var(β3)+2*Cov(β1,β3)の関係がありますので，
Covが分かれば，Varは分かります。
また，
Var(β1+β3)=Var(β1)+Var(β3)+2*Cov(β1,β3)の関係がありますので，
これらから，こつこつと求めていくしかありません。
もっと，コマンドを使って求めることはできないものでしょうか。
どなたか，ご教示ください。