为什么线性混合模型适用于 SAS 和 nlme 但不适用于 lme4?
Why does a linear mixed model work in SAS and nlme but not lme4?
我的数据包括 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
我设法在 nlme 中使用下面的语法进行建模,它与 SAS 的输出相匹配:
Model2 <- lme(DV ~ Group,
random = ~ 1|xVarExp/Subject, data = RawData)
我的问题是:1) 为什么该模型在 nlme 中有效但在 lme4 中无效?和 2) 如何匹配 SAS 语法以使模型在 lme4 中运行?
谢谢!
lme4 包有一些导致错误的内置模型检查。如果您需要使用 lmer
拟合不寻常的线性混合模型,您可以通过 lmerControl
中的参数更改忽略默认情况下模型检查该错误。
要允许随机效应与您正在拟合的模型中的残差项具有相同数量的水平,您需要将 check.nobs.vs.nlev
和 check.nobs.vs.nRE
更改为 "ignore"
从默认值 "stop"
。因此,您希望每组具有不同残差的模型可能类似于
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
我的数据包括 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
我设法在 nlme 中使用下面的语法进行建模,它与 SAS 的输出相匹配:
Model2 <- lme(DV ~ Group,
random = ~ 1|xVarExp/Subject, data = RawData)
我的问题是:1) 为什么该模型在 nlme 中有效但在 lme4 中无效?和 2) 如何匹配 SAS 语法以使模型在 lme4 中运行?
谢谢!
lme4 包有一些导致错误的内置模型检查。如果您需要使用 lmer
拟合不寻常的线性混合模型,您可以通过 lmerControl
中的参数更改忽略默认情况下模型检查该错误。
要允许随机效应与您正在拟合的模型中的残差项具有相同数量的水平,您需要将 check.nobs.vs.nlev
和 check.nobs.vs.nRE
更改为 "ignore"
从默认值 "stop"
。因此,您希望每组具有不同残差的模型可能类似于
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