21

我可以使用 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)
4

2 回答 2

33

使用现代 (>1.0) 版本，您可以在拟合和相应模型lme4之间进行直接比较，但您必须使用 ML --- 很难为没有随机的模型提出“REML 标准”的合理模拟效果（因为它将涉及将所有固定效果设置为零的数据的线性变换......）lmerlm

您应该知道，有和没有方差分量的模型之间的信息论比较存在理论问题：有关更多信息，请参阅GLMM 常见问题解答。

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
于 2014-06-03T16:48:20.807 回答
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)
于 2021-11-24T12:32:37.423 回答