mgcv后部覆盖

Question

我正在用R中的mgcv包拟合一个空间二项式模型，并想要模拟预测点的后验分布(代码如下)。我一直在使用模拟数据来测试后验的覆盖特性。我发现，当总流行率在0.5左右(50%)时，覆盖率很低(大约35%的真实值在95%的后间隔内)，但当你离开0.5时，这种情况会有所改善。例如，当平均患病率为1%时，~97%在95%的后部。我想我的问题是：

1. ，这是在这种方法中使用GAMs/mgcv的固有限制吗？
2. 是我对后验错误的贝叶斯解释吗？
3. 是我的代码中的错误内容吗？
4. 有更好的方法吗？(我尝试使用spaMM软件包来拟合一个具有空间相关随机效应的模型(用Laplace近似)，它做得稍微好一些。毫无疑问，MCMC方法会更好，但是地理统计方法在缩放模型/预测点的数量时有局限性，所以我想使用mgcv.

。

任何想法/意见将是非常欢迎！

干杯，休

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)                

Answer 1

更新：

通过多次运行这个实验，平均值和覆盖率之间的关系不像上面的例子那样极端。

另外，我使用了默认的k (即s(x,y))，使用gam.check表示k足够高。但是，如果使用较高的k(即s(x, y, k=100))，允许样条更加摇摆，那么预测间隔自然会更宽(即更不确定)，覆盖范围也会提高。覆盖范围仍然是可变的，但它要好得多。

想听听别人的想法。