考虑一下 lme4 中包含的 sleepstudy 数据。Bates 在他关于 lme4的在线书籍中讨论了这一点。在第 3 章中，他考虑了两种数据模型。
M0:Reaction∼1+Days+(1|Subject)+(0+Days|Subject)
和
MA:Reaction∼1+Days+(Days|Subject)
该研究涉及 18 名受试者，在 10 天睡眠不足的情况下进行了研究。在基线和随后几天计算反应时间。反应时间和睡眠剥夺持续时间之间有明显的影响。学科之间也存在显着差异。模型 A 允许随机截距和斜率效应之间存在相互作用的可能性：想象一下，例如，反应时间差的人更容易受到睡眠剥夺的影响。这意味着随机效应呈正相关。
在 Bates 的示例中，格子图没有明显的相关性，模型之间也没有显着差异。然而，为了研究上面提出的问题，我决定采用 sleepstudy 的拟合值，提高相关性并查看两个模型的性能。
从图中可以看出，较长的反应时间与较大的性能损失有关。用于模拟的相关性为 0.58
我根据我的人工数据的拟合值，使用 lme4 中的模拟方法模拟了 1000 个样本。我将 M0 和 Ma 分别拟合并查看了结果。原始数据集有 180 个观测值（18 个受试者各 10 个），模拟数据具有相同的结构。
底线是差别很小。
- 固定参数在两种模型下具有完全相同的值。
- 随机效应略有不同。每个模拟样本有 18 个截距和 18 个斜率随机效应。对于每个样本，这些影响被强制加到 0，这意味着两个模型之间的平均差异为（人为）0。但方差和协方差不同。MA 下的中位协方差为 104，而 M0 下的协方差为 84（实际值，112）。M0 下斜率和截距的方差大于 MA，大概是为了在没有自由协方差参数的情况下获得额外的摆动空间。
- lmer 的 ANOVA 方法给出了一个 F 统计量，用于将 Slope 模型与仅具有随机截距的模型进行比较（由于睡眠剥夺而没有影响）。显然，这个值在两种模型下都非常大，但在 MA 下它通常（但不总是）更大（平均 62 对 55 的平均值）。
- 固定效应的协方差和方差是不同的。
- 大约有一半的时间，它知道 MA 是正确的。用于比较 M0 和 MA 的中值 p 值为 0.0442。尽管存在有意义的相关性和 180 个平衡的观察结果，但只有大约一半的时间会选择正确的模型。
- 两种模型下的预测值不同，但非常细微。预测之间的平均差异为 0，sd 为 2.7。预测值本身的 sd 为 60.9
那么为什么会这样呢？@gung 合理地猜测，未能包括相关的可能性会迫使随机效应不相关。也许应该；但在这个实现中，允许将随机效应关联起来，这意味着数据能够将参数拉向正确的方向，而不管模型如何。错误模型的错误出现在可能性中，这就是为什么您可以（有时）在该级别区分两个模型的原因。混合效应模型基本上是为每个主题拟合线性回归，受模型认为它们应该是什么的影响。与正确模型相比，错误的模型会强制拟合不太合理的值。但归根结底，这些参数是由与实际数据的拟合决定的。
这是我有点笨拙的代码。想法是拟合睡眠研究数据，然后建立一个具有相同参数的模拟数据集，但随机效应的相关性更大。该数据集被输入到 simulation.lmer() 以模拟 1000 个样本，每个样本都适合两种方式。一旦我配对了拟合对象，我就可以提取拟合的不同特征并使用 t 检验或其他方法进行比较。
# Fit a model to the sleep study data, allowing non-zero correlation
fm01 <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=sleepstudy, REML=FALSE)
# Now use this to build a similar data set with a correlation = 0.9
# Here is the covariance function for the random effects
# The variances come from the sleep study. The covariance is chosen to give a larger correlation
sigma.Subjects <- matrix(c(565.5,122,122,32.68),2,2)
# Simulate 18 pairs of random effects
ranef.sim <- mvrnorm(18,mu=c(0,0),Sigma=sigma.Subjects)
# Pull out the pattern of days and subjects.
XXM <- model.frame(fm01)
n <- nrow(XXM) # Sample size
# Add an intercept to the model matrix.
XX.f <- cbind(rep(1,n),XXM[,2])
# Calculate the fixed effects, using the parameters from the sleep study.
yhat <- XX.f %*% fixef(fm01 )
# Simulate a random intercept for each subject
intercept.r <- rep(ranef.sim[,1], each=10)
# Now build the random slopes
slope.r <- XXM[,2]*rep(ranef.sim[,2],each=10)
# Add the slopes to the random intercepts and fixed effects
yhat2 <- yhat+intercept.r+slope.r
# And finally, add some noise, using the variance from the sleep study
y <- yhat2 + rnorm(n,mean=0,sd=sigma(fm01))
# Here is new "sleep study" data, with a stronger correlation.
new.data <- data.frame(Reaction=y,Days=XXM$Days,Subject=XXM$Subject)
# Fit the new data with its correct model
fm.sim <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=new.data, REML=FALSE)
# Have a look at it
xyplot(Reaction ~ Days | Subject, data=new.data, layout=c(6,3), type=c("p","r"))
# Now simulate 1000 new data sets like new.data and fit each one
# using the right model and zero correlation model.
# For each simulation, output a list containing the fit from each and
# the ANOVA comparing them.
n.sim <- 1000
sim.data <- vector(mode="list",)
tempReaction <- simulate(fm.sim, nsim=n.sim)
tempdata <- model.frame(fm.sim)
for (i in 1:n.sim){
tempdata$Reaction <- tempReaction[,i]
output0 <- lmer(Reaction ~ 1 + Days +(1|Subject)+(0+Days|Subject), data = tempdata, REML=FALSE)
output1 <- lmer(Reaction ~ 1 + Days +(Days|Subject), data=tempdata, REML=FALSE)
temp <- anova(output0,output1)
pval <- temp$`Pr(>Chisq)`[2]
sim.data[[i]] <- list(model0=output0,modelA=output1, pvalue=pval)
}