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

No.20765　混合効果モデル　　【赤岳】　2014/01/22(Wed) 18:21

いつもご教示ありがとうございます。
さて，困っています。時系列の繰り返し測定値データを混合効果モデルで行ったときのことです。responseについてsubjectを変量効果，time（1,2,4,6,8,10）を固定効果とする切片ランダムモデルで次のように行いました。

> library(nlme)
> Dat<- read.csv("test3.csv")
> Slope<- lme(response ~ time, random= ~1 |subject, method="REML", data=Dat)
> summary(Slope)
(省略）
Random effects:
 Formula: ~1 | subject
        (Intercept)  Residual
StdDev:   0.7424373 0.7297328

Fixed effects: response ~ time 
                Value  Std.Error DF  t-value p-value
(Intercept) 17.823875 0.24980521 69 71.35109       0
time        -0.115773 0.02500512 69 -4.62997       0

この結果，変数timeのt値は出ているのに，p値がバグのような結果となってしまいます。
そこで，次のように切片・傾きランダムモデルで解析すると，最もらしいp値が得られました。

> Slope<- lme(response ~ time, random= ~time|subject, method="REML", data=Dat)
> summary(Slope)
(省略）
Random effects:
 Formula: ~time | subject
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 1.1003544 (Intr)
time        0.1288360 -0.731
Residual    0.5850345       
Fixed effects: response ~ time 
                Value Std.Error DF  t-value p-value
(Intercept) 17.823875 0.3182558 69 56.00487  0.0000
time        -0.115773 0.0398434 69 -2.90570  0.0049

今回の切片ランダムモデルではどうしてp値が出ないのか，理由が分かりません。モデルに問題があるのでしょうか。
ご教示いただければ幸いです。
データを添付ファイルで送りましたが失敗しましたので，添付させてください。

subject	time	response	subject	time	response	subject	time	response
0001	1	17.76	0001	4	18.48	0001	8	17.76
0002	1	17.4	0002	4	18.36	0002	8	16.8
0003	1	19.56	0003	4	18.96	0003	8	14.16
0004	1	17.28	0004	4	16.44	0004	8	17.4
0005	1	15.84	0005	4	15.6	0005	8	15.24
0006	1	16.2	0006	4	16.2	0006	8	17.16
0007	1	17.76	0007	4	16.32	0007	8	16.92
0008	1	16.56	0008	4	16.2	0008	8	15.24
0009	1	18.84	0009	4	18.72	0009	8	17.76
0010	1	17.52	0010	4	18.6	0010	8	17.52
0011	1	18.24	0011	4	18.72	0011	8	16.2
0012	1	17.88	0012	4	17.76	0012	8	17.28
0013	1	18.24	0013	4	16.44	0013	8	15.72
0014	1	17.4	0014	4	16.56	0014	8	16.68
0001	2	18.48	0001	6	18.12	0001	10	17.52
0002	2	17.52	0002	6	16.68	0002	10	16.68
0003	2	19.44	0003	6	17.52	0003	10	16.08
0004	2	16.44	0004	6	16.92	0004	10	16.08
0005	2	16.08	0005	6	15.72	0005	10	15.24
0006	2	16.8	0006	6	16.92	0006	10	16.44
0007	2	17.52	0007	6	17.88	0007	10	17.28
0008	2	16.2	0008	6	16.56	0008	10	16.2
0009	2	19.08	0009	6	18.72	0009	10	18.48
0010	2	17.76	0010	6	18.48	0010	10	17.52
0011	2	18.72	0011	6	18.12	0011	10	17.52
0012	2	18.12	0012	6	17.88	0012	10	16.2
0013	2	17.52	0013	6	16.44	0013	10	15.6
0014	2	17.52	0014	6	17.16	0014	10	18.12

No.20766　Re: 混合効果モデル　　【青木繁伸】　2014/01/22(Wed) 19:07

データが違うか何かで，あなたのと同じ結果にならないけど，説明はつく。

> p値がバグのような結果
> どうしてp値が出ないのか

P 値はちゃんと出ているでしょう

上の場合は両方とも 0.0001 より小さいので 0 と表示されているだけ
> pt(71.35109, 69, lower.tail=FALSE)*2
[1] 2.176332e-66
> pt(4.62997, 69, lower.tail=FALSE)*2
[1] 1.667496e-05
下の場合は片方が 0.0001 より大きいので小数点以下4桁まで表示される。
そのとき，もう一方は小数点以下4桁は 0 なので 0.0000 が表示されるそれだけのこと。
data.frame の表示規則でこんなことが生じる。
> pt(56.00487, 69, lower.tail=FALSE)*2
[1] 2.966261e-59
> pt(2.90570, 69, lower.tail=FALSE)*2
[1] 0.004921145
summary の結果の tTable 要素を見れば，詳しい数値がわかりますよ
> a <- summary(Slope)
> a$tTable
                  Value Std.Error DF    t-value      p-value
(Intercept) 17.55428571 0.3028033 13 57.9725802 4.417642e-17
time         0.05142857 0.1230691 13  0.4178838 6.828520e-01

No.20767　Re: 混合効果モデル　　【赤岳】　2014/01/22(Wed) 19:15

青木先生，
お世話になっております。
そういうことですか。。。
ほんとうにありがとうございます。
また，tTableのコマンドもご教示ありがとうございます。
助かりました。