从 mcmc 对象中提取系数答案

【问题标题】：Extracting coefficients from mcmc object从 mcmc 对象中提取系数
【发布时间】：2019-04-09 00:12:50
【问题描述】：

从对象中提取东西一直是 R 对我来说最令人困惑的方面之一。我已经使用 rjags 拟合了贝叶斯线性回归模型，并具有以下 mcmc 对象：

summary(m_csim)
Iterations = 1:150000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 150000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

            Mean        SD  Naive SE Time-series SE
BR2     0.995805 0.0007474 1.930e-06      3.527e-06
BR2adj  0.995680 0.0007697 1.987e-06      3.633e-06
b[1]   -5.890842 0.1654755 4.273e-04      1.289e-02
b[2]    1.941420 0.0390239 1.008e-04      1.991e-03
b[3]    1.056599 0.0555885 1.435e-04      5.599e-03
sig2    0.004678 0.0008333 2.152e-06      3.933e-06

2. Quantiles for each variable:

            2.5%       25%       50%       75%    97.5%
BR2     0.994108  0.995365  0.995888  0.996339  0.99702
BR2adj  0.993932  0.995227  0.995765  0.996229  0.99693
b[1]   -6.210425 -6.000299 -5.894810 -5.784082 -5.55138
b[2]    1.867453  1.914485  1.940372  1.967466  2.02041
b[3]    0.942107  1.020846  1.057720  1.094442  1.16385
sig2    0.003321  0.004082  0.004585  0.005168  0.00657

为了提取系数的意思我做了b = colMeans(mod_csim)[3:5]。我想计算可信区间，所以我也需要提取 0.025 和 0.975 分位数。我如何以编程方式做到这一点？

【问题讨论】：

通常有这样的方法。不确定rjags，但是coef和fixef一般都是用的，要么在模型对象上，要么在模型对象的摘要上。也有不错的包，比如tidybayes。
如果您包含一个简单的reproducible example，其中包含可用于测试和验证可能解决方案的示例输入和所需输出，则更容易为您提供帮助。
我会从str(m_csim) 开始，看看它到底是什么。然后，请记住 summary functions 在打印之前多次执行自己的计算并运行 str(summary(m_csim))。

标签： r mcmc jags

【解决方案1】：

你可能正在寻找

model_summary_object <- summary(m_csim)
model_summary_object$quantiles[,c('2.5%','97.5%')]

【讨论】：

【解决方案2】：

您可能可以直接提取分位数。正如其他人指出的那样，您可以调用str(m_csim)，并且可以使用str(m_csim, max.level=1) 限制str 调用的输出，并继续在max.level= 参数中添加一个，直到您看到看起来像分位数的东西。

我喜欢做的是将 MCMC 输出转换为data.frame，这样更容易使用。我使用 jagsUI 而不是 rjags，但我经常做类似的事情

mcmc_df <- as.data.frame(as.matrix(MY_MCMC_OBJECT$samples))

注意：它可能与 rjags 有点不同，但我相信你可以稍微挖掘一下。

好处：然后我可以为每个 mcmc_df$PARAMETER 访问一个向量，并创建一个分位数矩阵

mcmc_quants <- apply(mcmc_df, 2, quantile, p=c(0.025, 0.25, 0.5, 0.75, 0.975))

或任何你想要的分位数。

【讨论】：

【解决方案3】：

我希望我没有超出我的知识范围，但我想从“一般”而不是专门针对 rjags 的角度来回答。 m_csim 是一个对象，可以在其上使用许多方法。您已使用摘要方法查看某些内容。正如人们评论的那样，可能有一种 coef 方法。然而，正如其他人评论的那样（当我在回复时！），使用 str() 查看对象包含的内容是查看对象中包含哪些信息以及如何处理它的最佳方法。如果使用 str() 不仅没有显示如何找到系数，而且还没有显示有关置信区间的足够信息以让您找到所需的 CI，我会感到非常惊讶。

【讨论】：