在R的HMSC 3.3-7中计算RRR模型WAIC时遇维度不匹配错误

问题描述 #

在为群落生态学数据实现HMSC模型并采用降秩回归(RRR)方法时，计算RRR模型的WAIC出现如下错误：

Error in X %*% Beta : non-conformable arguments

推测是拟合后的HMSC对象中X与Beta的维度不匹配，复现《Joint Species Distribution Modelling: With Applications in R》中的示例代码后仍出现相同问题，可复现代码如下：

library(Hmsc)

## MCMC test parameters
thin = 1
samples = 100
transient = 50
nChains = 2
verbose = 0


## First Block
ny = 20
nc = 10
ns = 10
XData = data.frame(matrix(rnorm(ny*nc), ncol = nc, nrow = ny))

## Second Block
X.FULL = as.matrix(XData)
IX.FULL = cbind(rep(1, ny), X.FULL)
nc.FULL = nc + 1
beta.FULL = matrix(rnorm(ns*nc.FULL), ncol = ns, nrow = nc.FULL)
eps = matrix(rnorm(ns*ny), ncol = ns, nrow = ny)
L.FULL = (IX.FULL %*% beta.FULL)
Y.FULL = L.FULL + eps

## Third Block
pc = princomp(X.FULL)
X.PC = pc$scores[,1]
IX.PC = cbind(rep(1,ny), X.PC)
nc.PC = 2
beta.PC = matrix(rnorm(ns*nc.PC), ncol = ns,
                 nrow = nc.PC)
L.PC = (IX.PC %*% beta.PC)
Y.PC = L.PC + eps

## Fourth Block
wRRR = matrix(c(1,-1,1,-1,1,-1,1,-1,1,-1))
X.RRR = X.FULL %*% wRRR
IX.RRR = cbind(rep(1, ny), X.RRR)
nc.RRR = 2
beta.RRR = matrix(rnorm(ns*nc.RRR), ncol = ns,
                  nrow = nc.RRR)
L.RRR = (IX.RRR %*% beta.RRR)
Y.RRR = L.RRR + eps

## Fifth Block
Y = list(Y.FULL, Y.PC, Y.RRR)
beta = list(beta.FULL, Y.PC, beta.RRR)

## Sixth Block
models = list()
for(dataset in 1:3){
    tmp = list()
    for (model in 1:3){
      switch(model,{
        m = Hmsc(Y = Y[[dataset]], XData = XData,
                 XFormula = ~., distr = "normal")
      },
      {
        pc = princomp(XData)
        XData.PC = data.frame(pc$scores[,1])
        m = Hmsc(Y = Y[[dataset]], XData = XData.PC,
                 XFormula = ~., distr = "normal")
      },
      {
      m = Hmsc(Y = Y[[dataset]], XData = XData, XFormula = ~1,
               XRRRData = XData, XRRRFormula =~.-1, ncRRR=1,
               distr = "normal")
      }
      )
      m = sampleMcmc(m,
                     thin = thin,
                     samples = samples,
                     transient = transient,
                     nChains = nChains,
                     verbose = verbose)
      tmp[[model]] = m
  }
models[[dataset]] = tmp
}

## WAIC Block
WAIC = matrix(NA, nrow = 3, ncol = 3)
for(dataset in 1:3){
  for (model in 1:3) {
    WAIC[dataset, model] = computeWAIC(models[[dataset]]
                                       [[model]])
    }
}

错误原因分析 #

这个错误的核心是旧版本HMSC包的computeWAIC函数在处理RRR模型时，未正确区分固定效应（X）和RRR效应（XRRR）的参数结构：

• 你的RRR模型中，XFormula = ~1 构造的固定效应设计矩阵X是20行×1列（仅包含截距项）
• RRR部分的参数BetaRRR是1行×10列（对应ncRRR=1和10个物种）
• 但旧版computeWAIC错误地尝试将X与合并后的固定效应+RRR参数矩阵（2行×10列）执行乘法，导致矩阵维度不兼容，触发non-conformable arguments错误。

解决方法 #

1. 升级HMSC包（推荐） #

该bug在HMSC的后续版本中已被修复，直接更新包即可解决：

install.packages("Hmsc")

更新后重新运行原代码，WAIC计算即可正常执行。

2. 手动实现WAIC计算（临时方案） #

如果暂时无法升级包，可以自定义适配RRR模型的WAIC计算函数：

computeWAIC_RRR <- function(model) {
  post_len <- length(model$postList)
  # 初始化对数似然矩阵
  log_lik <- matrix(0, nrow = model$ny, ncol = model$ns)
  
  for (i in 1:post_len) {
    # 提取当前迭代的参数
    BetaFixed <- model$postList[[i]]$Beta
    BetaRRR <- model$postList[[i]]$BetaRRR
    XRRR <- model$XRRRLatent
    sigma <- model$postList[[i]]$sigma
    
    # 计算线性预测（固定效应+RRR效应）
    lin_pred <- model$X %*% BetaFixed + XRRR %*% BetaRRR
    # 累加对数似然
    log_lik <- log_lik + dnorm(model$Y, mean = lin_pred, sd = sigma, log = TRUE)
  }
  
  # 计算平均对数似然
  avg_log_lik <- log_lik / post_len
  # 计算lppd（点预测密度对数和）
  lppd <- sum(log(rowMeans(exp(avg_log_lik))))
  # 计算pWAIC（有效参数数）
  pWAIC <- sum(apply(avg_log_lik, 1, var))
  # 计算WAIC
  WAIC <- -2 * (lppd - pWAIC)
  
  return(WAIC)
}

然后在WAIC计算块替换原函数：

WAIC[dataset, model] = computeWAIC_RRR(models[[dataset]][[model]])

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