我一直在尝试将重复测量模型从 SAS 转换为 R,因为合作者将进行分析但没有 SAS。我们处理 4 组,每组 8 到 10 只动物,然后每只动物有 5 个时间点。模拟数据文件可在此处https://drive.google.com/file/d/0B-WfycVUQyhaVGU2MUpuQkg4Mk0/edit?usp=sharing作为 Rdata 文件和此处https://drive.google.com/file/d/ 0B-WfycVUQyhaR0JtZ0V4VjRkTk0/edit?usp=共享为excel文件:
原始 SAS 代码 (1) 是:
proc mixed data=essai.data_test method=reml;
class group time mice;
model param = group time group*time / ddfm=kr;
repeated time / type=un subject=mice group=group;
run;
这使 :
Type 3 Tests des effets fixes
DDL DDL Valeur
Effet Num. Res. F Pr > F
group 3 15.8 1.58 0.2344
time 4 25.2 10.11 <.0001
group*time 12 13.6 1.66 0.1852
我知道 R 不像 SAS 那样处理自由度,所以我首先尝试获得类似于 (2) 的结果:
proc mixed data=essai.data_test method=reml;
class group time mice;
model param = group time group*time;
repeated time / type=un subject=mice group=group;
run;
我在这里找到了一些提示,将重复测量混合模型公式从 SAS 转换为 R,并且在指定复合对称相关矩阵时,这非常有效。但是,对于一般相关矩阵,我无法获得相同的结果。
在 SAS 中使用(2),我得到以下结果:
Type 3 Tests des effets fixes
DDL DDL Valeur
Effet Num. Res. F Pr > F
group 3 32 1.71 0.1852
time 4 128 11.21 <.0001
group*time 12 128 2.73 0.0026
使用以下 R 代码:
options(contrasts=c('contr.sum','contr.poly'))
mod <- lme(param~group*time, random=list(mice=pdDiag(form=~group-1)),
correlation = corSymm(form=~1|mice),
weights = varIdent(form=~1|group),
na.action = na.exclude, data = data, method = "REML")
anova(mod,type="marginal")
我得到:
numDF denDF F-value p-value
(Intercept) 1 128 1373.8471 <.0001
group 3 32 1.5571 0.2189
time 4 128 10.0628 <.0001
group:time 12 128 1.6416 0.0880
自由度是相似的,但不是固定效应的测试,我不知道这是从哪里来的。有人知道我在这里做错了什么吗?