5

我一直在尝试将重复测量模型从 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

自由度是相似的,但不是固定效应的测试,我不知道这是从哪里来的。有人知道我在这里做错了什么吗?

4

2 回答 2

1

您的 R 代码在多个方面与 SAS 代码不同。其中一些是可修复的,但我无法修复所有方面以重现 SAS 分析。

  1. R 代码拟合具有随机mice效应的混合效应模型,而 SAS 代码拟合允许残差之间存在相关性的广义线性模型,但没有随机效应(因为没有RANDOM陈述)。在 R 中,您必须使用同一个包中的gls函数。nlme

  2. 在 R 代码中,同一组内的所有观察值都具有相同的方差,而在 SAS 代码中,您有一个非结构化协方差矩阵,即每组内的每个时间点都有自己的方差。您可以使用weights=varIdent(form=~1|group*time).

  3. 在 R 代码中,无论组如何,每只鼠标的相关矩阵都是相同的。在 SAS 代码中,每个组都有自己的相关矩阵。这是我不知道如何在 R 中重现的部分。

我必须指出,R 模型似乎更有意义 - SAS 估计的方差和相关性太多(顺便说一句,您可以使用语句的RandRCORR选项看到有意义的排列)。repeated

于 2014-09-12T22:05:45.973 回答
0

“在 R 代码中,无论组如何,每只鼠标的相关矩阵都是相同的。在 SAS 代码中,每个组都有自己的相关矩阵。这是我不知道如何在 R 中重现的部分。” - 尝试:correlation=corSymm(~1|group*time)

于 2019-02-18T05:11:50.073 回答