R语言mvord包多变量有序Probit模型报错及公式适配问询

我用R的mvord包拟合多变量有序Probit模型，有两个有序因变量transit.commute.freq.pre.agg和transit.commute.freq.agg，各自对应不同的预测变量集。想联合估计模型，同时为每个因变量单独指定公式（类似GJRM包的分方程设定），但传入公式列表时出现错误：

Error in x$terms %||% attr(x, "terms") %||% stop("no terms component nor attribute") :
no terms component nor attribute

简化公式后仍报错：

commute_21_23 <- mvord(
  formula = list(
    transit.commute.freq.pre.agg ~ as.factor(h_cma),
    transit.commute.freq.agg ~ as.factor(age.gen)
  ), 
  data = emp_21_23_500
)

已确认因变量是有序因子且无缺失值，现提出两个问题：

1. mvord包是否支持为多变量有序模型的每个结果变量单独指定公式？
2. 若不支持，有没有输入重构方法或更适配的替代包（无法用GJRM，因为它不支持有序响应）？

示例数据：

structure(list(transit.commute.freq.agg=structure(3:1,levels=c("Non-transit users","Occasional transit users", "Regular transit users", "Frequent transit users"),class=c("ordered", "factor")),transit.commute.freq.pre.agg=structure(c(2L,4L,1L),levels =c("Non-transit users","Occasional transit users","Regular transit users","Frequent transit users"), class = c("ordered", "factor")),h_cma = structure(c(6L,1L,5L), levels =c("Toronto", "Montreal", "Vancouver","Ottawa-Gatineau","Calgary","Hamilton"),class="factor")),row.names = c("1", "2", "3"), class = "data.frame")

1. mvord包的公式支持情况 #

mvord包不支持直接传入公式列表为每个因变量单独指定预测变量，但可以通过主公式+mvformula参数的方式实现类似需求：

• 先定义包含所有因变量和所有可能预测变量的主公式
• 再用mvformula参数的列表，为每个因变量明确指定对应的预测变量子集

正确用法示例：

# 主公式：包含所有因变量和预测变量
main_formula <- cbind(transit.commute.freq.pre.agg, transit.commute.freq.agg) ~ h_cma + age.gen

# 用mvformula指定每个因变量对应的预测变量
commute_21_23 <- mvord(
  formula = main_formula,
  mvformula = list(
    transit.commute.freq.pre.agg ~ h_cma,  # 第一个因变量仅使用h_cma
    transit.commute.freq.agg ~ age.gen     # 第二个因变量仅使用age.gen
  ),
  data = emp_21_23_500,
  family = cumulative(probit)  # 指定有序Probit模型类型
)

注意：mvformula列表中每个元素的左侧必须与主公式中的因变量名完全匹配。

2. 替代方案 #

如果mvord的写法不符合需求，可考虑以下工具：

• brms包：支持贝叶斯框架下的多变量有序模型，能直接为每个因变量独立指定公式，同时联合估计变量间的相关性，示例代码：

  library(brms)
  
  model <- brm(
    bf(transit.commute.freq.pre.agg ~ h_cma) +
    bf(transit.commute.freq.agg ~ age.gen) +
    set_rescor(TRUE),  # 开启因变量间的相关性估计
    data = emp_21_23_500,
    family = cumulative(probit)
  )
  

• MASS包的polr函数：如果不需要联合估计因变量间的相关性，可分别拟合单变量有序Probit模型，但无法得到联合模型的相关性参数。

内容的提问来源于stack exchange，提问作者It's lighthouse