JAGS技术问询：参数向量与自变量矩阵相乘及Dirichlet多变量模型拟合

你的数据是100个样本的3物种比例丰度（y矩阵），以及对应的预测变量矩阵x，下面是完整的拟合流程，涵盖数据预处理、JAGS模型编写和R端运行代码：

第一步：数据预处理（修正截距项） #

你提到x的第一列为截距项，但当前生成代码的x是2列随机数，所以需要先给x添加截距列：

# 给x添加截距列（第一列全为1）
x <- cbind(1, x)

现在x变成100×3的矩阵（截距+2个预测变量）。

第二步：编写JAGS模型代码 #

Dirichlet分布的参数alpha必须为正，因此我们用对数线性模型关联预测变量与alpha（通过exp()确保参数为正）。将以下代码保存为dirichlet_model.jag：

model {
  # 1. 每个物种对应一组回归系数，设置弱信息先验
  for(species in 1:3) {          
    for(pred in 1:p) {           
      beta[pred, species] ~ dnorm(0, 0.001)
    }
  }
  
  # 2. 计算每个样本的Dirichlet参数alpha
  for(i in 1:n) {
    for(species in 1:3) {
      # 内积计算单个样本的对数alpha
      log_alpha[i, species] <- inprod(x[i, ], beta[, species])
    }
    # 指数转换得到正的alpha参数
    for(species in 1:3) {
      alpha[i, species] <- exp(log_alpha[i, species])
    }
    
    # 3. 观测模型：每个样本的物种比例服从Dirichlet分布
    y[i, ] ~ dchdirch(alpha[i, ])
  }
}

第三步：R端运行JAGS模型 #

使用rjags包完成模型编译、烧录和抽样：

library(rjags)

# 整理输入数据列表
data_list <- list(
  y = y,
  x = x,
  n = nrow(y),
  p = ncol(x)
)

# 参数初始化函数
init_fun <- function() {
  list(beta = matrix(rnorm(p*3), ncol = 3))
}

# 编译模型（启动3条独立链）
model <- jags.model(
  file = "dirichlet_model.jag",
  data = data_list,
  inits = init_fun,
  n.chains = 3
)

# 烧录迭代（丢弃初始不稳定样本）
update(model, n.iter = 1000)

# 抽样获取后验样本
post_samples <- coda.samples(
  model,
  variable.names = c("beta"),
  n.iter = 5000
)

# 查看后验结果
summary(post_samples)
plot(post_samples)

注意事项： #

• 如果实际数据中存在比例为0的情况，建议给y的每个元素加极小的伪计数（比如y <- y + 1e-6），避免Dirichlet分布出现数值不稳定。
• 先验可以根据需求调整，比如需要更强正则化时，可缩小正态先验的标准差（即增大精度参数）。

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