如何使用 lme4 将没有随机效应的模型与具有随机效应的模型进行比较？

Question

我可以使用 nlme 包中的 gls() 来构建没有随机效应的 mod1。 然后我可以将使用 AIC 的 mod1 与使用 lme() 构建的 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)

Answer 1

虽然我同意 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)

Answer 2

使用现代 (>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