回归参数的置信区间：贝叶斯与经典

机器算法验证 r 回归 贝叶斯 置信区间 常客
2022-01-30 02:40:34

给定两个长度为 n 的数组 x 和 y，我拟合模型 y = a + b*x 并希望计算斜率的 95% 置信区间。这是 (b - delta, b + delta) 以通常的方式找到 b 并且

delta = qt(0.975,df=n-2)*se.slope

se.slope 是斜率的标准误差。从 R 获得斜率的标准误差的一种方法是summary(lm(y~x))$coef[2,2]。

现在假设我写出给定 x 和 y 的斜率的可能性，将其乘以“平坦”先验，并使用 MCMC 技术从后验分布中抽取样本m 。定义

lims = quantile(m,c(0.025,0.975))

我的问题：(lims[[2]]-lims[[1]])/2大约等于上面定义的 delta 吗？

附录下面是一个简单的 JAGS 模型，这两个模型似乎不同。

model {
 for (i in 1:N) {
  y[i] ~ dnorm(mu[i], tau)
  mu[i] <- a + b * x[i]
 }
 a ~ dnorm(0, .00001)
 b ~ dnorm(0, .00001)
 tau <- pow(sigma, -2)
 sigma ~ dunif(0, 100)
}

我在 R 中运行以下命令：

N <- 10
x <- 1:10
y <- c(30.5,40.6,20.5,59.1,52.5,
       96.0,121.4,78.9,112.1,128.4)
lin <- lm(y~x)

#Calculate delta for a 95% confidence interval on the slope
delta.lm <- qt(0.975,df=N-2)*summary(lin)$coef[2,2]

library('rjags')
jags <- jags.model('example.bug', data = list('x' = x,'y' = y,'N' = N),
                   n.chains = 4,n.adapt = 100)
update(jags, 1000)
params <- jags.samples(jags,c('a', 'b', 'sigma'),7500)
lims <- quantile(params$b,c(0.025,0.975))
delta.bayes <- (lims[[2]]-lims[[1]])/2

cat("Classical confidence region: +/-",round(delta.lm, digits=4),"\n")
cat("Bayesian confidence region:  +/-",round(delta.bayes,digits=4),"\n")

并得到：

经典置信区域：+/- 4.6939

贝叶斯置信区域：+/- 5.1605

多次重新运行，贝叶斯置信区域始终比经典置信区域更宽。那么这是由于我选择的先验吗？

4个回答

“问题”在 sigma 上。尝试信息较少的设置

tau ~ dgamma(1.0E-3,1.0E-3)
sigma <- pow(tau, -1/2)

在你的 jags 文件中。然后更新一堆

update(10000)

获取参数，并总结您感兴趣的数量。它应该与经典版本相当吻合。

澄清：更新只是为了确保无论您决定选择哪个先验，您都可以到达您要去的地方，尽管像这种具有扩散先验和随机起始值的模型的链确实需要更长的时间才能收敛。在实际问题中，您会在总结任何内容之前检查收敛性，但收敛性不是我认为的示例中的主要问题。

如果你从 b | 的后部采样 y 并计算 lims（如您定义的那样）它应该与 (b - delta, b + delta) 相同。具体来说，如果你计算 b | y 在一个平坦的先验下，它与 b 的经典采样分布相同。

有关详细信息，请参阅：Gelman 等人。（2003 年）。贝叶斯数据分析。CRC出版社。第 3.6 节

编辑：

Ringold，你观察到的行为与贝叶斯思想是一致的。贝叶斯可信区间 (CI) 通常比经典可信区间宽。原因是，正如您正确猜到的那样，超先验考虑了由于未知参数引起的可变性。

对于像这样的简单场景（不是一般情况）：

贝叶斯 CI > 经验贝叶斯 CI > 经典 CI ；> == 更宽

对于线性高斯模型，最好使用 bayesm 包。它实现了先验的半共轭族，而杰弗里斯先验是该族的一个极限情况。请参阅下面的示例。这些是经典模拟，无需使用 MCMC。

我不记得关于回归参数的可信区间是否与通常的最小二乘置信区间完全相同，但无论如何它们非常接近。

> # required package
> library(bayesm)
> # data
> age <- c(35,45,55,65,75)
> tension <- c(114,124,143,158,166)
> y <- tension
> # model matrix
> X <- model.matrix(tension~age)
> # prior parameters
> Theta0 <- c(0,0)
> A0 <- 0.0001*diag(2)
> nu0 <- 0
> sigam0sq <- 0
> # number of simulations
> n.sims <- 5000
> # run posterior simulations
> Data <- list(y=y,X=X)
> Prior <- list(betabar=Theta0, A=A0, nu=nu0, ssq=sigam0sq)
> Mcmc <- list(R=n.sims)
> bayesian.reg <- runireg(Data, Prior, Mcmc)
> beta.sims <- t(bayesian.reg$betadraw) # transpose of bayesian.reg$betadraw
> sigmasq.sims <- bayesian.reg$sigmasqdraw
> apply(beta.sims, 1, quantile, probs = c(0.025, 0.975))
[,1] [,2]
2.5% 53.33948 1.170794
97.5% 77.23371 1.585798
> # to be compared with: 
> frequentist.reg <- lm(tension~age)

鉴于简单的线性回归在经典和贝叶斯分析与杰弗里的先验之间在分析上是相同的，两者都是解析的，诉诸诸如 MCMC 之类的数值方法来进行贝叶斯分析似乎有点奇怪。MCMC 只是一个数值积分工具，它允许贝叶斯方法用于更复杂的分析难以解决的问题，就像 Newton-Rhapson 或 Fisher Scoring 是用于解决难以解决的经典问题的数值方法一样。

使用 Jeffrey 的先验 p(a,b,s) 与 1/s 成比例的后验分布 p(b|y)（其中 s 是误差的标准差）是一个学生 t 分布，位置为 b_ols，尺度 se_b_ols (" ols”用于“普通最小二乘”估计）和 n-2 个自由度。但是b_ols的抽样分布也是一个student t，位置为b，尺度为se_b_ols，自由度为n-2。因此它们是相同的，除了 b 和 b_ols 已交换，因此在创建区间时，置信区间的“est +- 界限”会反转为可信区间中的“est -+ 界限”。

因此，置信区间和可信区间在分析上是相同的，使用哪种方法无关紧要（前提是没有额外的先验信息） - 所以采用计算成本更低的方法（例如，矩阵求逆较少的方法）。您对 MCMC 的结果表明，与 MCMC 一起使用的特定近似值给出的可信区间与精确的分析可信区间相比太宽了。这可能是一件好事（尽管我们希望近似值更好），因为近似贝叶斯解决方案看起来比精确贝叶斯解决方案更保守。

其它你可能感兴趣的问题
上一篇“混合效应建模”和“潜在增长建模”有什么区别？ 下一篇为什么倾向得分匹配比匹配更好？