如何使用 BIC 比较混合效应和固定效应广义线性模型？

机器算法验证 混合模式 广义线性模型 aic
2022-04-03 06:36:11

假设我有两个广义线性模型，一个只有固定效应，另一个有固定和随机效应，我如何比较使用 AIC/BIC 哪个模型最简约？我只能用 glmer 拟合混合效应模型，因为它在尝试拟合固定效应模型时会返回错误。我只能用 glm 拟合固定效应模型，因为这不允许随机效应。这给我留下了两组模型，一组装有 glmer，另一组装有 glm。

AIC 能否比较不同类型的模型？声明（我自己强调）：

这取决于。AIC 是对数似然的函数。如果两种类型的模型以相同的方式计算对数似然（即包括相同的常数），那么如果模型是嵌套的，那么可以。

我有理由确定 glm() 和 lmer() 不使用可比较的对数似然性。

关于嵌套模型的观点也有待讨论。有人说 AIC 仅对嵌套模型有效，因为这就是理论的呈现/工作方式。其他人将其用于各种比较。

DRAFT r-sig-mixed-models FAQ指出：

我可以将 AIC 用于混合模型吗？如何计算随机效应的自由度数？

是的，谨慎行事。

本页随后列出了一些注意事项和建议，例如使用修改后的 AIC 计算或REML=FALSE.

作为该领域的新手，我发现阅读有关使用 AIC/BIC 比较模型和使用混合效应模型进行群体效应的文章非常令人大开眼界，并且也是直观的“常识”方法。如何将这两种方法放在一起让我觉得我有些落入了裂缝，似乎没有一个明确的答案。当然，如果我遵循 AIC/BIC 简约方法，那么我应该将固定效应模型与混合效应模型与 AIC 进行比较，看看随机效应是否值得存在？

我错过了一些明显的东西吗？通常如何比较混合效应和固定效应广义线性模型？

1个回答

据我所知，您可以比较glmer()和glm()模型的可能性，至少对于family=binomial（尚未针对其他家庭进行测试）。如果方差分量估计为零，那么可能性应该是相同的，这显然是这种情况。下面是一个例子来说明这一点：

dat <- structure(list(id = c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 6L, 
6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 
9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L, 
12L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L, 14L, 14L, 14L, 
14L, 14L, 15L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 16L, 17L, 
17L, 17L, 17L, 17L, 18L, 18L, 18L, 18L, 18L, 19L, 19L, 19L, 19L, 
19L, 20L, 20L, 20L, 20L, 20L), xi = c(0, 0, 0, 0, 0, -1, -1, 
-1, -1, -1, -1, -1, -1, -1, -1, 0.8, 0.8, 0.8, 0.8, 0.8, -0.9, 
-0.9, -0.9, -0.9, -0.9, 0.7, 0.7, 0.7, 0.7, 0.7, 0.1, 0.1, 0.1, 
0.1, 0.1, -1.7, -1.7, -1.7, -1.7, -1.7, 0.3, 0.3, 0.3, 0.3, 0.3, 
-2.8, -2.8, -2.8, -2.8, -2.8, 2.7, 2.7, 2.7, 2.7, 2.7, -0.1, 
-0.1, -0.1, -0.1, -0.1, -0.2, -0.2, -0.2, -0.2, -0.2, 2, 2, 2, 
2, 2, -0.6, -0.6, -0.6, -0.6, -0.6, 1.1, 1.1, 1.1, 1.1, 1.1, 
0.2, 0.2, 0.2, 0.2, 0.2, -0.4, -0.4, -0.4, -0.4, -0.4, 2, 2, 
2, 2, 2, -1.1, -1.1, -1.1, -1.1, -1.1), xij = c(1.1, 1.1, 0.2, 
0.9, 0.4, -2.1, -0.4, -0.7, 0, 0.8, -0.4, 0.2, -1, 0, -1.2, 1.1, 
1.9, 0.9, -1.4, -0.8, -0.3, -0.7, 0.7, -1.2, 1.1, -1.5, 0.3, 
-1.7, -2, 0.2, 2, -0.5, -1.2, -0.2, -2.3, -0.6, -0.6, -1.6, -0.4, 
-1.5, -0.5, 0.8, 0.1, -0.3, -0.7, 0.7, 0.3, -0.4, 0.4, 0.5, -0.8, 
0.6, 0.3, 0.6, 0.2, -0.8, 0, -2.3, 0.5, 0, 0.9, 0.6, 2.2, 0.6, 
-0.3, 0.3, 0.5, -2.2, 2, -0.6, -0.7, -0.3, -0.7, 1.7, -0.7, -0.3, 
0.6, -0.9, -1.9, -0.5, 1.6, -0.5, 0.4, 1.1, 0.5, -1.8, 1.2, 1.7, 
-1.1, 0.2, -0.6, -1.1, 2.1, 0.4, 0.9, 0.5, -2, 1.6, 0.1, 0.7), 
    yi = c(0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 
    0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L)), .Names = c("id", 
"xi", "xij", "yi"), row.names = c(NA, -100L), class = "data.frame")

library(lme4)

res0 <- glm(yi ~ xi + xij, data=dat, family=binomial)
summary(res0)

res1 <- glmer(yi ~ xi + xij + (1 | id), data=dat, family=binomial)
summary(res1)

logLik(res0)
logLik(res1)
anova(res1, res0)

最后三行产生：

> logLik(res0)
'log Lik.' -29.96427 (df=3)
> logLik(res1)
'log Lik.' -29.96427 (df=4)
> 
> anova(res1, res0)
Data: dat
Models:
res0: yi ~ xi + xij
res1: yi ~ xi + xij + (1 | id)
     Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
res0  3 65.929 73.744 -29.964   59.929                        
res1  4 67.929 78.349 -29.964   59.929     0      1          1

因此，（对数）似然是相同的，因为id水平方差分量估计为零。因此，混合效应模型的 AIC 值比预期的大 2 个点（因为模型多了一个参数）。

不过需要注意的一件事：默认glmer()值为nAGQ=1，这意味着使用拉普拉斯近似。让我们使用“适当的”自适应正交：

res1 <- glmer(yi ~ xi + xij + (1 | id), data=dat, family=binomial, nAGQ=7)
logLik(res0)
logLik(res1)
anova(res1, res0)

这产生：

>     logLik(res0)
'log Lik.' -29.96427 (df=3)
>     logLik(res1)
'log Lik.' -29.96427 (df=4)
>     anova(res1, res0)
Error in anova.merMod(res1, res0) : 
  GLMMs with nAGQ>1 have log-likelihoods incommensurate with glm() objects

方差分量仍然估计为零，并且（对数）似然是相同的。但是anova()会吐出一个错误，表明这些模型不应该相互比较。

其它你可能感兴趣的问题
上一篇最著名的集合相似性度量？ 下一篇在 PCA 中，是否有系统的方法来删除变量以最大化两个群体的隔离？