R（lme）和SAS（proc mix）中线性混合模型的不同结果

Question

我尝试使用nlme包中的lme函数来拟合随机系数模型，示例代码如下：

model_example <- lme(CFB~ BASE + Dose + ADY + Dose:ADY + BASE:ADY,
                   random = ~ ADY| SubjectID, 
                   correlation = corSymm(form = ~ 1 | SubjectID),
                   data = example_data,
                   method = "REML", na.action = na.exclude,
                   control = lmeControl(msMaxIter = 2000, opt = "optim"))

基于此代码，我们尝试拟合具有非结构化相关性的随机截距和斜率模型。

同时，我们也用同样的数据在SAS中拟合模型，代码如下：

proc mixed data example_data MAXITER= 2000;
    class Dose SubjectID;
    model CFB = Dose BASE ADY BASE *ADY Dose*ADY/SOLUTION cl DDFM= BW;
                random Int ADY/type=un subject=SubjectID G GCORR VCORR;
run;

SAS代码也适合具有非结构化相关性的随机截距和斜率模型，但是，这两个模型给我们带来了不同的结果

而棘手的部分是，当我们删除R代码中的相关部分时（根据我对帮助文档的理解，不存在组内相关性），那么我们得到与具有非结构化相关性的SAS代码相同的结果.

model_example <- lme(CFB~ BASE + Dose + ADY + Dose:ADY + BASE:ADY,
                   random = ~ ADY| SubjectID, 
                   #correlation = corSymm(form = ~ 1 | SubjectID),
                   data = example_data,
                   method = "REML", na.action = na.exclude,
                   control = lmeControl(msMaxIter = 2000, opt = "optim"))

大家有什么想法吗？

Answer 1

我认为您可能误解了

corSymm

的作用；它适合非结构化相关矩阵组内，这可能不是您想要的（并且未包含在 SAS 规范中......）也许您将其与组内复合对称性混淆了（即，每个一组内的观察对是同等相关的），这是由对组的截距随机效应引起的（这已经在你的模型中）？

来自

?corSymm

：

form

：形式为“~ t”或“~ t |”的单边公式克’， 指定时间协变量“t”和可选的分组 因子 ‘g’ ... 默认为 ‘~ 1’，对应于 using 数据中观测值的顺序作为协变量，以及 没有团体。 [已添加强调]

这里的文档并不完全明确，但这意味着

corSymm = ~ 1 | g

将与协变量

1 .. n

拟合随机效应，协变量是观测值的组内索引。我们很快就会看到它的作用。

使用和不使用

corSymm

来拟合模型（模拟数据的代码位于问题末尾）：

library(nlme)
m2 <- lme(CFB ~ BASE + Dose + ADY + Dose:ADY + BASE:ADY,
          random = ~ 1 + ADY | SubjectID,
          data = dd, method = "REML")
m2B <- update(m2, correlation = corSymm(form = ~ 1 | SubjectID))

查看对数似然：

logLik(m2)
## 'log Lik.' -1614.338 (df=10)
 logLik(m2B)
## 'log Lik.' -1593.981 (df=55)

m2B

（具有相关性）还有 45 个参数 - 我的示例每组有 10 个观测值，因此相关矩阵为 10×10 并且需要 (10*9)/2 = 45 个参数 - 以及更大的对数似然 (因为模型要复杂得多，所以可以更精确地拟合数据）。

如果我们查看

print(m2B)

或

summary(m2B)

，我们可以看到该模型确实包含完整的 10×10 估计相关矩阵：

Correlation Structure: General
 Formula: ~1 | SubjectID 
 Parameter estimate(s):
 Correlation: 
   1      2      3      4      5      6      7      8      9     
2  -0.143                                                        
3   0.182  0.102                                                 
4   0.012  0.061  0.181                                          
5  -0.138 -0.055  0.060 -0.176                                   
6   0.112 -0.027  0.050 -0.101 -0.224                            
7  -0.177 -0.158 -0.127 -0.134 -0.141  0.108                     
8   0.026  0.256  0.188  0.167 -0.196  0.169  0.191              
9   0.017  0.148  0.124 -0.069  0.068 -0.136 -0.103  0.054       
10  0.009  0.096  0.117  0.112 -0.065  0.189  0.086  0.216  0.110

模拟数据

set.seed(101)
dd <- data.frame(Dose = rnorm(1000), SubjectID = gl(100, 10),
                 ADY = rnorm(1000), BASE = factor(rep(1:2, 500)),
                 step = gl(10,100))

## use ::: to avoid having to load lme4 before running ...
dd$CFB <- lme4:::simulate.formula(~ BASE + Dose + ADY + Dose:ADY + BASE:ADY + (1 + ADY | SubjectID),
                   newdata = dd,
                   newparams = list(beta = rep(0, 6), theta = rep(1, 3),
                                    sigma = 1))[[1]]