绘制贝叶斯 beta 回归模型预测的置信区间

Question

我有下面的样本数据和代码，如果你能帮助我绘制贝叶斯 beta 回归模型的可靠预测区间，我将不胜感激。

library(ggplot2)
library(plotly)
library(zoib)

data("GasolineYield", package = "zoib")

re.md <- zoib(yield ~ temp | 1 | 1, data=GasolineYield, 
              joint = FALSE, random=1, EUID=GasolineYield$batch, 
              zero.inflation = FALSE, one.inflation = FALSE, 
              n.iter=3200, n.thin=15, n.burn=200)

pred <- pred.zoib(re.md, data.frame(temp = seq(100, 600, 0.01)))

df <- data.frame(temp = seq(100, 600, 0.01), 
                 yield = (pred$pred[[1]][, 201] + pred$pred[[2]][, 201])/2)

ggplotly( 
ggplot() + 
geom_point(data = GasolineYield, 
aes(x = temp, y = yield, fill = batch), 
size = 4, shape = 21) + 
xlim(100, 600) + 
geom_line(data = df, aes(y = yield, x = temp), col="red") + 
theme_classic())

Answer 1

我对贝叶斯统计没有什么经验（虽然我很想进入它），但我相信这就是你所追求的：

df1 <- data.frame(temp = seq(100, 600, 0.01), 
                  pred$summary)
ggplotly( 
  ggplot() + 
    geom_point(data = GasolineYield, 
               aes(x = temp, y = yield, fill = batch), 
               size = 4, shape = 21) + 
    xlim(100, 600) + 
    geom_line(data = df1, aes(y = mean, x = temp), col="red") + 
    geom_ribbon(data = df1, aes(ymin= X2.5., ymax = X97.5., x = temp), alpha = 0.3) +
    theme_classic())

来自?pred.zoib的帮助：

summary if TRUE (the default), a basic summary on each posterior predictive value, including mean, SD, min, max, med, 2.5% and 97.5% quantiles, are provided.

这与您正在绘制的有点不同，因为摘要中的平均值实际上是：

rowSums(pred$pred[[1]])/ncol(pred$pred[[1]]

可视化差异：

df <- data.frame(temp = seq(100, 600, 0.01), 
                 yield = (pred$pred[[1]][, 201] + pred$pred[[2]][, 201])/2)

ggplotly( 
  ggplot() + 
    geom_point(data = GasolineYield, 
               aes(x = temp, y = yield, fill = batch), 
               size = 4, shape = 21) + 
    xlim(100, 600) + 
    geom_line(data = df1, aes(y = mean, x = temp), col="red") + 
    geom_ribbon(data = df1, aes(ymin= X2.5., ymax = X97.5., x = temp), alpha = 0.3) +
    geom_line(data = df, aes(y = yield, x = temp), col="blue") + 
    theme_classic())

一些额外的注意事项：

all.equal(rowSums(pred$pred[[1]])/ncol(pred$pred[[1]]), df1$mean)
#output
TRUE

all.equal(apply(pred$pred[[1]], 1, quantile, probs = 0.025), df1$X2.5.)
#output
TRUE

all.equal(apply(pred$pred[[1]], 1, quantile, probs = 0.975), df1$X97.5.)
#output
TRUE

max、min等也是如此

我不确定 pred$pred[[2]] 代表什么，但您可以使用上述方法为其生成摘要并绘制如下：

df2 <- data.frame(temp = seq(100, 600, 0.01), 
              mean = apply(pred$pred[[2]], 1, mean),
              X97.5. = apply(pred$pred[[2]], 1, quantile, probs = 0.975),
              X2.5. = apply(pred$pred[[2]], 1, quantile, probs = 0.025))

让我们绘制两者（小心我的 R 在使用 ggplotly 执行此操作时变得无响应）：

  ggplot() + 
    geom_point(data = GasolineYield, 
               aes(x = temp, y = yield, fill = batch), 
               size = 4, shape = 21) + 
    xlim(100, 600) + 
    geom_line(data = df1, aes(y = mean, x = temp), col="red") + 
    geom_ribbon(data = df1, aes(ymin= X2.5., ymax = X97.5., x = temp), alpha = 0.3) +
    geom_line(data = df2, aes(y = mean, x = temp), col="blue") + 
    geom_ribbon(data = df2, aes(ymin= X2.5., ymax = X97.5., x = temp), alpha = 0.3)+
    theme_classic()

绘制贝叶斯 beta 回归模型预测的置信区间

Plotting confidence intervals of predictions from a Bayesian beta regression model

regression

r

bayesian

prediction

ggplot2