为什么线性混合模型适用于 SAS 和 nlme 而不是 lme4？答案

【问题标题】：Why does a linear mixed model work in SAS and nlme but not lme4?为什么线性混合模型适用于 SAS 和 nlme 而不是 lme4？
【发布时间】：2016-06-30 02:26:09
【问题描述】：

我的数据由对照组中的 20 名受试者和实验组中的 20 名受试者组成。感兴趣的 DV 是在每个参与者上测量的峰值功率的变化分数。还有一个虚拟变量xVarExp，其中仅包含实验组中的受试者的 1。我对个人反应感兴趣，这些数字的差异是总结这一点的统计数据。我也对每个组的手段感兴趣；扩展和控制。

我的数据结构如下：

structure(list(Subject = structure(1:40, .Label = c("Alex", "Ariel", 
"Ashley", "Bernie", "Casey", "Chris", "Corey", "Courtney", "Devon", 
"Drew", "Dylan", "Frances", "Gene", "Jaimie", "Jean", "Jesse", 
"Jo", "Jody ", "Jordan", "Kelly", "Kerry", "Kim", "Kylie", "Lauren", 
"Lee", "Leslie", "Lindsay", "Morgan", "Pat", "Reilly", "Robin", 
"Sage", "Sam", "Sidney", "Terry", "Tristan", "Vic", "Wil", "Wynn", 
"Zane"), class = "factor"), Group = structure(c(1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L), .Label = c("Control", "Exptal"), class = "factor"), 
    xVarExp = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1), DV = c(3.3, -0.8, 2.7, 2.8, 0.6, 5.2, 
    1, 3.4, 1.3, -2.4, 8.5, 3.5, -1.9, 4.3, 1.2, -1.9, -0.6, 
    1.3, -2.6, -1, -3.7, 1.9, 4.6, 2.9, 7.2, -1.7, 4.2, 3.9, 
    -3.2, 9.9, 2.7, -1.7, 7.9, 8.1, 3.8, 2.8, 4.6, 0.8, 2.5, 
    4.1)), .Names = c("Subject", "Group", "xVarExp", "DV"), row.names = c(NA, 
-40L), class = "data.frame")

统计学家是 SAS 用户，并使用以下代码获得了合理的答案：

title "Analyzing change scores";
proc mixed data=import plots(only)=StudentPanel(conditional) alpha=0.1 nobound;
class Subject Group;
model DV=Group/residual outp=pred ;
random xVarExp/subject=Subject;
lsmeans Group/diff=control("Control") cl alpha=0.1;
run;

我开始使用 R 和 lme4，我相信代码是：

Model1 <- lmer(DV ~ Group + (1|Subject/xVarExp), 
             data = RawData)

但是，我收到以下信息：Error: number of levels of each grouping factor must be < number of observations

我设法在与 SAS 输出匹配的 nlme 中使用以下语法使建模工作：

Model2 <- lme(DV ~ Group, 
            random = ~ 1|xVarExp/Subject, data = RawData)

我的问题是：1）为什么该模型在 nlme 中有效，但在 lme4 中无效？和 2) 如何匹配 SAS 语法以使模型进入 lme4？

谢谢！

【问题讨论】：

这实际上是您的完整数据吗？
是的，我只是在 R 中做了 dput(RawData)。有什么遗漏吗？
所以您没有对您的主题进行重复测量，xVarExp 实际上只是 Group 的数字版本？
每个受试者都有一个前后测量，它们之间的区别是变化分数/感兴趣的DV。如果需要，我可以上传 pre 和 post 数据。之前只有一次测量，后期只有一次测量。
您可以通过一些lmerControl 参数绕过lmer 中的这种行为。例如，您可以将观察次数与级别数的检查更改为"ignore" 而不是停止。如control = lmerControl(check.nobs.vs.nlev = "ignore", check.nobs.vs.nRE="ignore")。如果您尝试允许每组不同的残差方差，那么您可能希望您的随机效应类似于(xVarExp-1|Subject)。

标签： r sas lme4 nlme

【解决方案1】：

lme4 包有一些导致错误的内置模型检查。如果您需要使用lmer 拟合不寻常的线性混合模型，您可以更改忽略模型，默认情况下通过lmerControl 中的参数检查该错误。

要允许与您正在拟合的模型中的残差项具有相同数量级别的随机效应，您需要将 check.nobs.vs.nlev 和 check.nobs.vs.nRE 从默认的 "ignore" 更改为 "ignore" @。因此，您希望每组具有不同残差的模型可能看起来像

Model1 <- lmer(DV ~ Group + (xVarExp-1|Subject), 
            data = RawData, control = lmerControl(check.nobs.vs.nlev = "ignore",
                                        check.nobs.vs.nRE="ignore"))

但是，如果您想要的模型允许每组具有不同的残差方差，那么您可以考虑使用 nlme 中的gls。在gls 中，您可以轻松扩展线性模型以允许非常量方差。该模型看起来像

Model2 <- gls(DV ~ Group, data = RawData, weights = varIdent(form = ~1|Group))

这两个模型对固定效应给出了相同的估计值和标准误：

summary(Model1)$coefficients
            Estimate Std. Error  t value
(Intercept)    1.395  0.6419737 2.172986
GroupExptal    1.685  1.0449396 1.612533

summary(Model2)$tTable
            Value Std.Error  t-value    p-value
(Intercept) 1.395 0.6419737 2.172986 0.03608378
GroupExptal 1.685 1.0449396 1.612533 0.11512145

【讨论】：