零膨胀泊松MCMCglmm先验分布维度错误排查求助

背景 #

刚接触贝叶斯统计，在运行包含大量零值的节肢动物类群丰度多响应MCMCglmm模型时，遇到先验分布维度不匹配的问题。原模型包含52个响应变量、7个固定效应和1个随机效应，为调试简化为3个响应变量、3个固定效应的最小可复现示例。

简化测试代码 #

datum <- read.csv("2024_data_raw_MCMCglmm.csv")

datum$Site <- factor(datum$Site) # 随机效应变量

abund_vars <- c("Lycosidae","Gnaphosidae","Pholcidae") # 响应变量
abund_vars_sc <- paste(abund_vars,"s",sep="_") # 标准化变量
datum2 <- datum

fmla.MMLM1 <- cbind(Lycosidae, Gnaphosidae, Pholcidae) ~
  trait:(ShrubVolume + X.Leafy + P_cover_shrubs) - 1
# trait对应固定效应

fam.test <- rep("zipoisson", 3) # 每个变量使用零膨胀泊松分布

prior.model.1 <- list(R = list(V=diag(3)/3, nu=0.003),
                      G = list(G1=list(V=diag(3)/3, nu = 0.003)))
## 先验模型存在问题，运行时报错！##

t2 <- system.time(
  MMLM1.fit <- MCMCglmm(fmla.MMLM1,
                        random=~ us(trait):Site,
                        rcov=~ us(trait):units,
                        prior = prior.model.1,
                        data = datum2,
                        family = fam.test,
                        nitt= 10000, burnin= 2000, thin=10)
)

报错信息 #

Error in priorformat(if (NOpriorG) { : 
  V is the wrong dimension for some prior$G/prior$R elements
Timing stopped at: 0.21 0 0.22

疑问 #

原本认为V的维度应该是diag(3)对应3个响应变量，但不清楚错误原因。更困惑的是，当family设为"poisson"时，相同先验可以正常运行，零膨胀模型是否改变了先验的要求？

解决方案 #

零膨胀模型（如zipoisson）的先验维度要求和普通泊松模型不同——零膨胀模型对每个响应变量包含两个参数：计数部分的泊松参数，以及零膨胀部分的逻辑斯蒂参数。因此，对于3个响应变量，先验矩阵的维度应该是6x6，而非3x3。

具体修改 #

调整先验中的V矩阵维度为响应变量数量的2倍：

prior.model.1 <- list(R = list(V=diag(6)/6, nu=0.003),
                      G = list(G1=list(V=diag(6)/6, nu = 0.003)))

原因解释 #

当使用多响应零膨胀模型时，MCMCglmm会将每个响应变量的两个组件（计数+零膨胀）视为联合分布的一部分，因此协方差矩阵的维度是响应变量数量的2倍。而普通泊松模型每个响应变量只有一个参数，3个变量对应3x3矩阵即可满足要求。

内容的提问来源于stack exchange，提问作者Joseph Coolidge