来自 mgcv 的后验覆盖率答案

【问题标题】：Coverage of posterior from mgcv来自 mgcv 的后验覆盖率
【发布时间】：2019-12-05 21:53:17
【问题描述】：

我正在使用 R 中的 mgcv 包拟合空间二项式模型，并希望模拟预测点的后验分布（代码如下）。我一直在使用模拟数据来测试后验的覆盖属性。我发现，当总体流行率约为 0.5 (50%) 时，覆盖率很差（约 35% 的真实值位于 95% 的后验区间内），但随着您远离 0.5，这种情况会有所改善。例如，当平均患病率为 1% 时，约 97% 位于 95% 的后验范围内。我想我的问题是：

这是在这种方法中使用 GAM/mgcv 的固有限制吗？
我对后验的贝叶斯解释不正确吗？
我是否在代码中指定错误？
有更好的方法吗？（我尝试使用spaMM 包来拟合具有空间相关随机效应（使用拉普拉斯近似）的模型，该模型效果更好。毫无疑问，MCMC 方法会更好，但地统计方法在缩放数量时有局限性模型/预测点，所以我很想使用mgcv。

非常欢迎任何想法/cmets！

干杯，休

library(mgcv)
library(RandomFields)
library(raster)

# Simluate some data
set.seed(1981)
mean <- 0
model <- RMexp(var=0.5, scale=50)
simu <- RandomFields::RFsimulate(model, x=1:256, 
                                 y=1:256, RFoptions(spConform=FALSE))

# Convert to raster
simu_raster <- raster(nrows = 256, ncol = 256, xmn=0, xmx=1, ymn=0, ymx=1)
simu_raster[] <- as.vector(simu)

# Add mean and onvert to probability
log_odds_raster <- mean + simu_raster 
prev_raster <- exp(log_odds_raster) / (1 + exp(log_odds_raster))

# simulate 1000 candidate sampling points
candidate_points <- coordinates(prev_raster)[sample(1:nrow(coordinates(prev_raster)), 1000),]

# Sample 100 of those and take binomial sample of 100 individuals per location 
sampled_points_idx <- sample(1:nrow(candidate_points), 100)
sampled_points <- as.data.frame(candidate_points[sampled_points_idx,])
sampled_points$n_pos <- rbinom(100, 100, extract(prev_raster, sampled_points))
sampled_points$n_neg <- 100 - sampled_points$n_pos

# Fit spatial GAM
spatial_mod <- gam(cbind(n_pos, n_neg) ~ s(x, y), 
                   data = sampled_points,
                   family="binomial")

# check k and plot observed v predicted
gam.check(spatial_mod)

# Simulate 1000 draws from the posterior at every non-sampled location
prediction_data <- as.data.frame(candidate_points[-sampled_points_idx,])
prediction_data$prev <- extract(prev_raster, prediction_data)
Cg <- predict(spatial_mod, prediction_data, type = "lpmatrix")
sims <- rmvn(1000, mu = coef(spatial_mod), V = vcov(spatial_mod, unconditional = TRUE))
fits <- Cg %*% t(sims)
fits_prev <- exp(fits) / (1 + exp(fits))

# For every prediction point, see whether the true/simulated prevalence
# lies within the posterior with correct accuracy. i.e. 95% of the time, 
# the true value should lie within the 95% BCI. 
BCI_95 <- apply(fits_prev, 1, FUN=function(x){quantile(x, prob = c(0.025, 0.975))})
within_BCI <- c()
for(i in 1:nrow(prediction_data)){
  within_BCI <- c(within_BCI, (prediction_data$prev[i] >= BCI_95[1,i] &
                                  prediction_data$prev[i] <= BCI_95[2,i]))
}
mean(within_BCI)

【问题讨论】：

标签： r spatial gam mgcv

【解决方案1】：

更新：

多次运行此实验后，均值和覆盖率之间的关系不像上面的示例那样极端。

另外，在上面，我使用了默认的k（即s(x,y)）并使用gam.check，这表明k足够高。但是如果你使用更高的 k（即s(x, y, k=100)），让样条曲线更加摆动，那么预测区间自然会更宽（即更不确定）并且覆盖范围会提高。覆盖范围仍然变化很大，但要好得多。

很想听听别人的想法。

【讨论】：