使用mvord运行双变量固定效应有序Probit模型时遇正定错误

针对你遇到的Error in chol.default(derivs_for_se$H) : the leading minor of order 373 is not positive definite错误，本质是模型的Hessian矩阵非正定，通常由参数过多、多重共线性或数值优化不稳定导致，以下是适配mvord包的具体解决方案：

1. 重构固定效应（优先处理） #

你将ResponseId作为固定效应纳入模型，若ResponseId的类别数较多（比如373个水平），会直接导致参数数量暴增，引发Hessian矩阵数值不稳定。

替代方案：改用随机效应 #

将ResponseId设为随机效应，大幅减少待估参数：

fit <- mvord(
  formula = MMO2(Safety,Comfort) ~ Separation_from_Vehicle_Lane + Bike_Yielding_Priority + 
    Art_Ads_in_Center + Uni_Directional + Bike_Volumes + Facility_Style + X4_Legged 
    + Q23_2 + Q3 + Q4 + Q6 + Q32 + Q33 + Q5 + Q23_1 + Q23_3 + Q23_4 + Q23_5 + 
    Q23_6 + Q23_7 + Q23_8 + Violations + Errors + Positive_behavior + Roundabout_Number + (1 | ResponseId),  
  data = safety_comfort_full,
  error.structure = cor_general(~ 1),
  link = mvprobit(),
  threshold.constraints = NULL,
  control = mvord.control()
)

若必须保留固定效应： #

• 检查ResponseId的水平数：length(unique(safety_comfort_full$ResponseId))，合并样本量极小的类别；
• 确认是否存在重复的ResponseId观测，考虑聚合个体层面数据后再建模。

2. 排查并消除多重共线性 #

自变量间高度相关会导致Hessian矩阵奇异，可通过方差膨胀因子（VIF）检测：

# 近似检测（将有序因变量转为数值做线性回归）
lm_temp <- lm(as.numeric(Safety) ~ Separation_from_Vehicle_Lane + Bike_Yielding_Priority + 
    Art_Ads_in_Center + Uni_Directional + Bike_Volumes + Facility_Style + X4_Legged 
    + Q23_2 + Q3 + Q4 + Q6 + Q32 + Q33 + Q5 + Q23_1 + Q23_3 + Q23_4 + Q23_5 + 
    Q23_6 + Q23_7 + Q23_8 + Violations + Errors + Positive_behavior + Roundabout_Number, 
  data = safety_comfort_full)
car::vif(lm_temp)

• 移除VIF>10的变量；
• 对Q23系列这类同源变量，可先做因子分析降维，用提取的公因子替代原始变量。

3. 调整数值优化参数 #

修改mvord.control()中的优化设置，提升模型收敛稳定性：

fit <- mvord(
  # 其他参数不变
  control = mvord.control(maxit = 1000, rel.tol = 1e-6, optimizer = "nlminb")
)

• maxit：增加迭代次数，确保模型充分收敛；
• rel.tol：降低收敛阈值；
• optimizer：切换优化器（nlminb比默认的optim更适合复杂模型）。

4. 检查数据质量 #

• 检查自变量方差：移除方差为0的变量（所有观测值相同）：

  var_check <- sapply(safety_comfort_full[, c("Separation_from_Vehicle_Lane", "Bike_Yielding_Priority", ...)], var)
  drop_vars <- names(var_check[var_check == 0])
  safety_comfort_full <- safety_comfort_full[, !names(safety_comfort_full) %in% drop_vars]
  

• 检查因变量水平：确保Safety和Comfort的有序因子水平数≥2；
• 处理缺失值：用sum(is.na(safety_comfort_full))排查，删除或插补缺失观测。

5. 逐步构建模型 #

从最简模型开始逐步添加变量，定位引发错误的具体变量：

# 最简模型
fit_simple <- mvord(
  formula = MMO2(Safety,Comfort) ~ Separation_from_Vehicle_Lane + Bike_Yielding_Priority,
  data = safety_comfort_full,
  error.structure = cor_general(~ 1),
  link = mvprobit()
)

若最简模型运行正常，逐个添加剩余变量，每次运行后检查是否报错，找到问题变量后针对性处理。

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