我正在使用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)
}