WAIC 计算期间的警告:如何进行?

机器算法验证 r 贝叶斯 模型选择 aic 哀伤
2022-04-19 14:53:17

我正在使用Rwaic()中的包中的函数计算 WAIC(广泛适用或 Watanabe-Akaike 信息标准) 'loo'。这样做时,我看到如下警告消息:

Warning messages: 1: 118 (5.1%) p_waic estimates greater than 0.4. We recommend trying loo() instead.

我该怎么办?我注意到有几个在线示例/教程忽略了这个警告(例如这里 - 这是一个可重现的示例)。5.1%(在我的例子中)或 3.3%(在这个例子中)是一个危险的高数字吗?在比较两个模型时,这可能会如何影响我的推理?如果在我的比较中 WAIC 值真的不同(比如超过 40 个 WAIC 单位不同)怎么办?我可以相信一个模型是首选的推论吗?

改为使用loo()(用于留一法交叉验证),我得到一个不同的警告:

Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

我的模型使用伯努利误差,因此不能错误指定误差项。在我的特定模型类别中,线性预测器也不太可能被严重错误指定。

警告,血腥细节如下

我的模型是一个潜变量模型(见这里) 共同模拟 15 种鸟类在约 100 个地点的分布。该模型不包括协变量,因此可以被认为是基于模型的排序。地点 j 的物种 i 的发生概率在概率尺度上建模为特定于物种的截距加上几个 (1-3) 潜在因素的影响,这些潜在因素的行为类似于 GLM 中的协变量。每个潜在因子在每个站点都有一个特定的值(根据数据估计),并且每个物种都有一个特定于物种的因子负载,该负载控制其跨站点的丰度。因子载荷矩阵乘以其转置矩阵是物种间完整 (15x15) 方差-协方差矩阵的近似值。我想用两个潜在因素来拟合模型(因为这会产生易于在二维中可视化的排序)。然而,将潜在因子的数量更改为一或三个不会提高 p_waic 估计值大于 0.4 的百分比;这让我相信我无法通过为线性预测器选择更好的规范来改善这个问题。

我进行模型选择的原因不是要挖掘大量可能的协变量,而是要对两个代表有意义不同的生物学假设/场景的替代模型进行推断。我想知道证据是否强烈支持一种假设而不是另一种假设。

对于那些感兴趣的人,JAGS 中的完整模型规范,包括对数似然的计算,如下所示:

##### Model specification #####
LV_C <- function() {
  ## Data Level ##
  for(i in 1:n) { # sites
    for(j in 1:p){ # species
      eta[i,j] <- inprod(lv.coefs[j,2:(num.lv+1)],lvs[i,]) # LV part of the linear predictor
      Z[i,j] ~ dnorm(lv.coefs[j,1] + eta[i,j], 1) # This line and the next one implement a probit link
      y[i,j] ~ dbern(step(Z[i,j]))
      loglikelihood[(i-1)*(p)+j] <- y[i,j]*log(phi(lv.coefs[j,1] + eta[i,j])) + (1 - y[i,j])*log(1 - phi(lv.coefs[j,1] + eta[i,j]))
    }
  }


  ## Latent variables ## 
  for(i in 1:n) { for(k in 1:num.lv) { lvs[i,k] ~ dnorm(0,1) } } # Says what Latent Variable values are

  ## Process level and priors ##
  for(j in 1:p) { lv.coefs[j,1] ~ dnorm(0,0.01) } ## Separate species intercepts

  for(i in 1:(num.lv-1)) { for(j in (i+2):(num.lv+1)) { lv.coefs[i,j] <- 0 } } ## Constraints to 0 on upper diagonal
  for(i in 1:num.lv) { lv.coefs[i,i+1] ~ dunif(0,20) } ## Sign constraints on diagonal elements
  for(i in 2:num.lv) { for(j in 2:i) { lv.coefs[i,j] ~ dunif(-20,20) } } ## Free lower diagonals
  for(i in (num.lv+1):p) { for(j in 2:(num.lv+1)) { lv.coefs[i,j] ~ dunif(-20,20) } } ## All other elements
}


#  Set up for Run

jags.data = list(y=ssm2, n=dim(ssm2)[1], p=dim(ssm2)[2],
                 num.lv=2)

params = c('lvs', 'lv.coefs', 'loglikelihood')

p <- dim(ssm2)[2]
n <- dim(ssm2)[1]
inits <- function(jjj) {
  Tau <- rWishart(1,p+1,diag(p))[,,1]
  Sigma <- solve(Tau)
  Z <- abs(t(mvtnorm::rmvnorm(n,rep(0,p),Sigma)))
  Z <- ifelse(as.matrix(ssm2), Z, -1 * Z)

  list(Z=Z)
}
1个回答

我很高兴您关心这些诊断信息。我负责那个特定的警告信息。WAIC 没有很好的方法来诊断可靠性,0.4 阈值是根据经验选择的。PSIS-LOO 中的 Pareto k 诊断要好得多。您没有提到具体的 k 值,但如果它们大于 1,那么理论和实验表明误差可以是任意大的(参见 Aki Vehtari、Andrew Gelman 和 Jonah Gabry(2017)。Pareto 平滑重要性抽样。https : //arxiv.org/abs/1507.02646)。WAIC 和 PSIS-LOO 是相连的,因此如果 PSIS-LOO 失败,那么 WAIC 会失败得更多(Aki Vehtari、Andrew Gelman 和 Jonah Gabry (2017)。使用留一法交叉验证和 WAIC 进行实用贝叶斯模型评估。在统计中和计算,27(5):1413–1432,https://arxiv.org/abs/1507.04544)。因此,即使一个 p_waic 可以有任意大的误差,那么总的 WAIC 也可以有任意大的误差,不应该相信最终结果。

由于您在 WAIC 上有很大差异,因此差异可能与更准确/稳健的方法(例如 k-fold-CV)相似,但如果不计算更准确/稳健的结果,您无法确定。在您的情况下,WAIC 和 PSIS-LOO 很可能会失败,因为您有一个灵活的模型,并且一些观察结果非常有影响力,也就是说,完整的数据后验和留一法后验是如此不同,以至于使用完整的后验由于重要性采样 LOO 中的提案分布失败。WAIC 改用泰勒级数近似,这在有高度影响的观察的情况下也会失败(这解释了为什么 p_waic 被用作临时诊断)。