如何为 MCMCglmm 模型正确编写比例逆 Wishart 先验代码？

Question

我正在运行使用 MCMCglmm() 创建一个具有两个随机效应的多元模型（4 个响应变量）。我目前正在使用逆 Wishart 先验。

###current set up for prior
prior.miw<-list(R=list(V=diag(4), nu=0.002), 
                G=list(G1=list(V=diag(4), 
                nu=0.002,   
                alpha.mu=c(0,0,0,0), 
                alpha.V=diag(4)*1000),
                G2=list(V=diag(4), #need to repeat to deal with second random effect
                nu=0.002,   
                alpha.mu=c(0,0,0,0), 
                alpha.V=diag(4)*1000)))

根据我的模型，看来我真的应该运行缩放逆 Wishart 先验（参见 page 14 of Lemoine 2019）。

我的问题：如何将其变成适合我的 MCMCglmm() 模型的先验？

此外，如果这不是我应该运行宁的那种先验，请随时权衡。

Answer 1

这是一个 two-part 问题：

对于可能性集中在小值的多元随机效应，我应该使用什么先验？（我假设这就是您正在寻找默认逆 Wishart 先验的替代方法的原因）[更适合 CrossValidated]
其中哪些在 MCMCglmm 中可用，我如何在那里实施它们？ [适合 Stack Overflow]

一般技巧是将协方差矩阵分解为多元分量（相关矩阵或逆相关矩阵或其他）和标准差（或逆标准差）的缩放参数向量； Lemoine 建议将 U(0,100) 用于缩放先验，我认为这很糟糕（为什么平坦？我无法找到 Gelman 和 Hill 2007 的精确页面，他们在那里讨论 which 分布用于缩放先验......但如果他们真的建议在方差尺度上进行均匀分布，我会感到有点惊讶......）

更新实际查看了您的代码 (!)：我认为您做的是正确的，只是 nu=0.002 看起来确实很极端；见讨论结束。

这基本上就是 MCMCglmm 所做的，但它使用不同的（IMO 更好）选择来缩放先验。听起来很可怕:

These priors are all from the non-central scaled F-distribution, which implies the prior for the standard deviation is a non-central folded scaled t-distribution (Gelman, 2006).

但这归结为选择四个参数，您真正需要考虑的只有两个。

V：先验均值方差（或先验均值协方差矩阵，如果您有多元随机效应项）。根据课程笔记，“不失一般性，V 可以设置为一个”（或者在多变量模型的情况下，设置为单位矩阵）
alpha.mu：我们几乎总是希望它为零（或者在你的例子中，一个零向量）；这样标准偏差的先验将是学生 t 分布。（alpha.mu != 0 可能有一个用例，但我从来没有运行穿过它。）
alpha.V：V 设置为 1（或单位矩阵），这是协方差矩阵的先验均值。对于您的问题，具有合理尺度的对角矩阵是一个不错的选择
nu：形状参数；作为 nu → ∞ 我们得到标准差的 half-Normal 先验，nu=1 我们得到柯西分布。值越小，尾部越厚（conservative/allowing 更宽的样本更少，但尾部出现奇怪的采样行为的危险也更大）。

对于单变量情况，Hadfield 说 V=1 的 t 先验是

2 * dt(sqrt(v)/sqrt(alpha.V), df = nu, 
       ncp = alpha.mu/sqrt(alpha.V))

v 限制为零。

如果我们设置 alpha.mu=0 像明智的人一样 ncp （non-centrality 参数）参数为零；
前面的 2 使密度加倍，因为我们只查看分布的正半部分（“折叠”）
除以 sqrt(alpha.V) 意味着我们将 t-distribution 缩放 alpha.V
换句话说，我们将通过 sqrt(alpha.V)*abs(rt(1, nu=nu))

我希望 Hadfield 给出了一个 多变量 inverse-Wishart 案例中的参数扩展示例。你在哪里找到你的？

我认为nu=0.002是一个危险的选择。为什么？

让你的先验“尽可能平坦”很诱人，但这有两个潜在的不良后果。

对于某些情况，例如 inverse-Gamma 或 Gamma 分布，较小的形状参数也会在零附近产生较大的尖峰；这就是你上面 R=list(V=diag(4), nu=0.002) 之前的残差情况。如果您的数据提供了信息（对残差的合理猜测），这是可以的，但如果不是，则采样会很糟糕。
对于 parameter-expanded 的情况，我们没有 spike-at-zero 问题（half-Normal、half-t 和 half-Cauchy 为零），但是通过使尾巴变得如此肥大，您可能会允许估算器对真正荒谬（并且在计算上有问题）的值进行采样。（同样，如果您的数据信息量很大，这不是问题。）考虑 nu=1 的 t 分布对应于 Cauchy 分布，其尾巴非常肥大以至于 其均值是无限的 ，然后想象将形状参数缩小 500 倍 (nu=0.002) ...

如何为 MCMCglmm 模型正确编写比例逆 Wishart 先验代码？

How to properly code a scaled inverse Wishart prior for a MCMCglmm model?

r

bayesian

mcmc

multivariate-testing