如何使用lme4比较没有随机效应的模型和有随机效应的模型？

Question

我可以使用nlme包中的gls()函数来建立没有随机效应的mod1模型。 然后，我可以使用lme()函数来建立包含随机效应的mod2模型，并通过AIC对比mod1和mod2。

mod1 = gls(response ~ fixed1 + fixed2, method="REML", data)
mod2 = lme(response ~ fixed1 + fixed2, random = ~1 | random1, method="REML",data)
AIC(mod1,mod2)

是否有类似于lme4包中gls()的东西，它允许我构建没有随机效应的mod3，并将其与使用包括随机效应的lmer()构建的mod4进行比较？

mod3 = ???(response ~ fixed1 + fixed2, REML=T, data)
mod4 = lmer(response ~ fixed1 + fixed2 + (1|random1), REML=T, data)
AIC(mod3,mod4)

- It Figures

Answer 1

使用现代版本（>1.0）的lme4，你可以直接比较lmer拟合和相应的lm模型，但是你必须使用ML --- 对于没有随机效应的模型，很难想出“REML标准”的合理类比（因为它将涉及数据的线性转换，将所有固定效应设置为零...）

你应该知道，在具有和不具有方差分量的模型之间进行信息论比较存在理论问题：请参见GLMM FAQ获取更多信息。

library(lme4)
fm1 <- lmer(Reaction~Days+(1|Subject),sleepstudy, REML=FALSE)
fm0 <- lm(Reaction~Days,sleepstudy)
AIC(fm1,fm0)
##     df      AIC
## fm1  4 1802.079
## fm0  3 1906.293

我更喜欢以这种格式输出（使用delta-AIC而不是原始AIC值）：

bbmle::AICtab(fm1,fm0)
##     dAIC  df
## fm1   0.0 4 
## fm0 104.2 3 

为了测试，我们模拟没有随机效应的数据（我不得不尝试几个随机数种子才能得到一个实际上估计出组间标准差为零的示例）：

要进行测试，让我们模拟没有随机效应的数据（我不得不尝试几个随机数种子，才能得到一个实际上会将受试者之间的标准偏差估计为零的示例）：

rr <- simulate(~Days+(1|Subject),
               newparams=list(theta=0,beta=fixef(fm1),
                         sigma=sigma(fm1)),
               newdata=sleepstudy,
               family="gaussian",
               seed=103)[[1]]
ss <- transform(sleepstudy,Reaction=rr)
fm1Z <- update(fm1,data=ss)
VarCorr(fm1Z)
##  Groups   Name        Std.Dev.
##  Subject  (Intercept)  0.000  
##  Residual             29.241
fm0Z <- update(fm0,data=ss)
all.equal(c(logLik(fm0Z)),c(logLik(fm1Z)))  ## TRUE

- Ben Bolker

---

我可以使用当前（devel）版本的lme4成功地运行示例。如果是使用你自己的数据，则需要更多信息；请在StackOverflow上提出新问题，或发送电子邮件到[e-mail][首先订阅该列表; 你可以通过搜索找到信息/订阅页面]，在任何情况下都需要提供可重现的示例（以及packageVersion("lme4")的结果）。 - Ben Bolker

当随机模型使用ML结果重新拟合后变为奇异时，应采取什么方法？但是，当使用REML拟合时，它并不是奇异的...有什么建议吗？ - Asier

这并不令人意外，特别是对于高度参数化/不稳定的模型而言。你能否在CrossValidated或[e-mail]上发布一个可重现的示例？ - Ben Bolker

非常感谢，Ben。我会尝试去做......我同意你的看法，我认为这个模型参数太多了。然而奇怪的是，它仍然比其他将因素建模为线性和二次协变量的模型具有更低的AICc。这意味着我需要在可重复的示例中拟合整个数据集。 - Asier

完成了，我在Stack Overflow提供了一个可重现的示例：stackoverflow.com/q/60892398/13099627?sem=2 - Asier - Asier

@BenBolker虽然我同意你提出的是最简单的解决方案，但没有任何随机效应的模型的受限似然函数计算起来相当简单，详见我的回答。 - Jarle Tufto

Answer 2

尽管我同意Ben的观点，即最简单的解决方案是将REML设置为FALSE，但在没有随机效应的模型中，最大REML似然是明确定义的，并且可以通过众所周知的关系相对简单地计算。

在普通轮廓似然函数和受限制似然函数之间。

以下代码模拟数据，其中线性混合模型（LMM）的随机截距的估计方差为0，因此LMM的最大受限对数似然应该等于不包括任何随机效应的模型的受限似然。

通过上述公式计算LM的受限似然，并评估与LMM相同的值。

一个更简单的替代方法是使用glmmTMB。

library(lme4)
#> Loading required package: Matrix
# simulate some toy data for which the LMM ends up at the boundary
set.seed(5)
n <- 100 # the sample size
x <- rnorm(n) 
y <- rnorm(n)
group <- factor(rep(1:10,10))

# fit the LMM via REML
mod1 <- lmer(y ~ x + (1|group), REML=TRUE, control=lmerControl(boundary.tol=1e-8))
#> boundary (singular) fit: see ?isSingular
logLik(mod1)
#> 'log Lik.' -147.8086 (df=4)

# fit a model without random effects and compute its maximum REML log likelihood
mod0 <- lm(y ~ x)
p <- length(coef(mod0)) # number of fixed effect parameters
X <- model.matrix(mod0) # the fixed effect design matrix
sigma.REML <- summary(mod0)$sigma # REMLE of sigma
# the maximum ordinary log likelihood evaluated at the REML estimates
logLik.lm.at.REML <- sum(dnorm(residuals(mod0), 0, sigma.REML, log=TRUE))
# the restricted log likelihood of the model without random effects (via above formula)
logLik.lm.at.REML + p/2*log(2*pi) - 1/2*(- p*log(sigma.REML^2) + determinant(crossprod(X))$modulus)
#> [1] -147.8086
#> attr(,"logarithm")
#> [1] TRUE

library(glmmTMB)
data <- data.frame(y,x,group)
logLik(glmmTMB(y~x, family = gaussian(), data=data, REML=TRUE))
#> 'log Lik.' -147.8086 (df=3)
logLik(glmmTMB(y~x+(1|group), family = gaussian(), data=data, REML=TRUE))
#> 'log Lik.' -147.8086 (df=4)

- Jarle Tufto