个人认为,下面的code应该没问题。anova table也好像和你的一样。
但是,我还没想清楚检验simple effect应该怎么计算自由度?
???好像contrast包的结果都用42有点太多了???
是不是应该用residual df=22或者satterthwaite/kenwardroger近似?
<br />
dat=data.frame(<br />
y=c(910.22,809.32,1006.72,884.1,951.2,782,969.67,1066.03,<br />
806.24,701.72,881.67,813.44,1404,1201,1052.19,1044.26,<br />
1132.69,1033.91,1129.53,1023.88,1360.61,1582,1060.12,<br />
1104,1024,964.3,1125,1304.5,1184.97,1004.18,985,1161.12,<br />
996,763.57,1014,1041.34,1540,2877,1541.8,1254.49,1416.76,<br />
1353.62,1405.27,1719.33,1558,1854,1316.22,1472.39),<br />
a=gl(2,12),<br />
b=as.factor(rep(1:2, each = 24)),<br />
rat=as.factor(rep(1:24,2)))<br />
library(nlme)<br />
lmef=lme(y~a*b,dat,~1|rat)<br />
anova(lmef)<br />
contrast(lmef,list(a='1',b='1'),list(a='2',b='1'))<br />
contrast(lmef,list(a='1',b='2'),list(a='2',b='2'))<br />
contrast(lmef,list(a='1',b='1'),list(a='1',b='2'))<br />
contrast(lmef,list(a='2',b='1'),list(a='2',b='2'))<br />
results:
> anova(lmef)
numDF denDF F-value p-value
(Intercept) 1 22 857.4545 <.0001
a 1 22 28.3354 <.0001
b 1 22 22.9274 0.0001
a:b 1 22 4.5574 0.0442
> contrast(lmef,list(a='1',b='1'),list(a='2',b='1'))
lme model parameter contrast
Contrast S.E. Lower Upper t df Pr(>|t|)
-295.4883 101.8443 -495.0994 -95.87723 -2.9 42 0.0059
> contrast(lmef,list(a='1',b='2'),list(a='2',b='2'))
lme model parameter contrast
Contrast S.E. Lower Upper t df Pr(>|t|)
-561.7417 101.8443 -761.3528 -362.1306 -5.52 42 0
> contrast(lmef,list(a='1',b='1'),list(a='1',b='2'))
lme model parameter contrast
Contrast S.E. Lower Upper t df Pr(>|t|)
-165.4708 88.19085 -338.3217 7.380063 -1.88 42 0.0676
> contrast(lmef,list(a='2',b='1'),list(a='2',b='2'))
lme model parameter contrast
Contrast S.E. Lower Upper t df Pr(>|t|)
-431.7242 88.19085 -604.5751 -258.8733 -4.9 42 0
>