AIC，方差分析错误：模型并非都适合相同数量的观察，模型并非都适合相同大小的数据集

机器算法验证 r 混合模式 aic
2022-03-21 15:15:28

我有这样的模型：

require(nlme)

set.seed(123)
n <- 100
k <- 5
cat <- as.factor(rep(1:k, n))
cat_i <- 1:k # intercept per kategorie
x <- rep(1:n, each = k)
sigma <- 0.2
alpha <- 0.001
y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma)
plot(x, y)

m1 <- lm(y ~ x)
summary(m1)

m2 <- lm(y ~ cat + x)
summary(m2)

m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit)
summary(m3)

现在我正在尝试评估模型中是否应该存在随机效应。所以我使用 AIC 或 anova 比较模型，我得到以下错误：

> AIC(m1, m2, m3)
   df       AIC
m1  3 1771.4696
m2  7 -209.1825
m3  4 -154.0245
Warning message:
In AIC.default(m1, m2, m3) :
  models are not all fitted to the same number of observations  
> anova(m2, m3)
Error in anova.lmlist(object, ...) : 
  models were not all fitted to the same size of dataset

如您所见，在这两种情况下，我都使用相同的数据集。我找到了两种补救措施，但我认为它们并不令人满意：

  1. 添加method = "ML"到 lme() 调用- 不确定更改方法是否是个好主意。
  2. 改为使用lmer()。令人惊讶的是，尽管 lmer() 使用 REML 方法，但它仍然有效。但是我不喜欢这个解决方案，因为lmer()它不显示系数的 p 值 - 我喜欢使用旧的lme()。

你知道这是否是一个错误，我们如何解决这个问题？

2个回答

快速搜索表明它是可能的（尽管我不得不承认我认为它不是）并且它不是错误......只是 R 中的方法被隐藏并导致看起来“意外”的另一种情况'，但 RTFM 人群说，“它在文档中。” 无论如何......您的解决方案是anova将lme作为第一个参数和lm模型作为第二个（如果您愿意，还有第三个）参数。如果这看起来很奇怪，那是因为它有点奇怪。原因是当您调用 时，仅当第一个参数是对象时才调用anova该方法。否则，它会调用（它又会调用;如果你深入研究，你会明白为什么）。有关您希望如何拨打电话的详细信息anova.lmelmeanova.lmanova.lmlistanova.lmanova在这种情况下，请为anova.lme. 您会看到您可以将其他模型与lme模型进行比较，但它们必须位于第一个参数以外的位置。显然，也可以anova使用该函数拟合模型，gls而无需过多关心模型参数的顺序。但我不知道足够的细节来确定这是否是一个好主意，或者它究竟意味着什么（看起来可能很好，但你的电话）。从那个链接比较lm，lme似乎有据可查并被引用为一种方法，所以我会在那个方向犯错，是你吗。

祝你好运。

这绝对是奇特的。作为第一个想法：在模型具有不同的固定效应结构（例如）的情况下进行模型比较时，我们最好使用m2ML“改变”。（它会将它与相乘，其中）有趣的是，它使用它可以工作，这让我相信它可能不是一个错误。似乎它几乎强制执行“良好做法”。m3MLREMLykkX=0method="ML"

话虽如此，让我们看看引擎盖下：

 methods(AIC)  
 getAnywhere('AIC.default')

 A single object matching ‘AIC.default’ was found
 It was found in the following places
   registered S3 method for AIC from namespace stats
   namespace:stats with value

 function (object, ..., k = 2) 
 {
     ll <- if ("stats4" %in% loadedNamespaces()) 
         stats4:::logLik
     else logLik
     if (!missing(...)) {
         lls <- lapply(list(object, ...), ll)
         vals <- sapply(lls, function(el) {
             no <- attr(el, "nobs") #THIS IS THE ISSUE!
             c(as.numeric(el), attr(el, "df"), if (is.null(no)) NA_integer_ else no)
         })
         val <- data.frame(df = vals[2L, ], ll = vals[1L, ])
         nos <- na.omit(vals[3L, ])
         if (length(nos) && any(nos != nos[1L])) 
             warning("models are not all fitted to the same number of observations")
         val <- data.frame(df = val$df, AIC = -2 * val$ll + k * val$df)
             Call <- match.call()
             Call$k <- NULL
         row.names(val) <- as.character(Call[-1L])
         val
     }
     else {
         lls <- ll(object)
         -2 * as.numeric(lls) + k * attr(lls, "df")
     }     
 }

在您的情况下，您可以看到：

  lls <- lapply(list(m2,m3), stats4::logLik)
  attr(lls[[1]], "nobs")
  #[1] 500
  attr(lls[[2]], "nobs")
  #[1] 498

显然logLik正在做某事（也许？）出乎意料...？不，不是真的，如果您查看 , 的文档logLik，?logLik您会看到它明确说明：

 There may be other attributes depending on the method used: see
 the appropriate documentation.  One that is used by several
 methods is ‘"nobs"’, the number of observations used in estimation
 (after the restrictions if ‘REML = TRUE’)

这使我们回到了最初的观点，您应该使用ML.

用 CS 中的一句俗话说：“这不是错误；这是（真正的）功能！”

编辑：（只是为了解决您的评论：）假设您使用lmer这个时间适合另一个模型：

m3lmer <- lmer(y ~ x + 1|cat)

并且您执行以下操作：

lls <- lapply(list(m2,m3, m3lmer), stats4::logLik)
attr(lls[[3]], "nobs")
#[1] 500
 attr(lls[[2]], "nobs")
#[1] 498

两者之间似乎存在明显差异，但实际上并非如 Gavin 所解释的那样。只是为了说明显而易见的：

 attr( logLik(lme(y ~ x, random = ~ 1|cat, na.action = na.omit, method="ML")),
 "nobs")
#[1] 500

我认为在方法论方面发生这种情况是有充分理由的。lme确实会尝试为您理解 LME 回归，而lmer在进行模型比较时，它会立即退回到它的 ML 结果。我认为在这个问题上没有错误，lme每个lmer包背后都有不同的理由。

另请参阅 Gavin Simposon 的评论，对发生的事情进行更深入的解释anova()（实际上发生了同样的事情AIC）

其它你可能感兴趣的问题
上一篇机器学习在数论中的应用 下一篇统一的先验如何从最大似然和后验模式得出相同的估计？