Dirichlet回归更换参考类别时优化不收敛问题求助

数据概况 #

我有一个包含行为比例的data.table（df2），数据结构和头部如下：

Model    Cat_id    Day     Active Inactive Maintenance Total propActive propInactive propMaintenance Model2
 1: C.RF1    Cho       1 days  1936   75672    8792        86400 0.02240741    0.8758333      0.10175926  C.RF1
 2: C.RF1    Cho       2 days  1307   78236    6857        86400 0.01512731    0.9055093      0.07936343  C.RF1
 3: C.RF1    Cho       3 days  1360   73784    11256       86400 0.01574074    0.8539815      0.13027778  C.RF1
 4: C.RF1    Cho       4 days  2828   70666    12906       86400 0.03273148    0.8178935      0.14937500  C.RF1
 5: C.RF1    Cho       5 days  2988   74130    9282        86400 0.03458333    0.8579861      0.10743056  C.RF1
 6: C.RF1    Cho       6 days  1809   74477    10114       86400 0.02093750    0.8620023      0.11706019  C.RF1

Classes ‘data.table’ and 'data.frame':  1152 obs. of  11 variables:
$ Model          : Factor w/ 16 levels "C.RF1","C.RF2",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Cat_id         : Factor w/ 12 levels "Cho","George",..: 1 1 1 1 1 1 2 2 2 2 ...
$ Day            : 'difftime' num  1 2 3 4 ...
..- attr(*, "units")= chr "days"
$ Active         : int  1936 1307 1360 2828 2988 1809 2616 2697 2540 3796 ...
$ Inactive       : int  75672 78236 73784 70666 74130 74477 74862 75588 74547 73742 ...
$ Maintenance    : int  8792 6857 11256 12906 9282 10114 8922 8115 9313 8862 ...
$ Total          : int  86400 86400 86400 86400 86400 86400 86400 86400 86400 86400 ...
$ propActive     : num  0.0224 0.0151 0.0157 0.0327 0.0346 ...
$ propInactive   : num  0.876 0.906 0.854 0.818 0.858 ...
$ propMaintenance: num  0.1018 0.0794 0.1303 0.1494 0.1074 ...
$ Model2         : Factor w/ 16 levels "H.RF1","C.RF1",..: 2 2 2 2 2 2 2 2 2 2 ...
- attr(*, ".internal.selfref")=<externalptr> 
- attr(*, "sorted")= chr [1:3] "Model" "Cat_id" "Day"

初始建模 #

我用Dirichlet回归测试行为比例差异，代码运行正常，但默认以C.RF1作为参考类别，所有模型都和它对比：

dirig <- DR_data(df2[, c("propActive","propInactive","propMaintenance")], base=1)
m1 <- DirichReg(dirig ~ Model + Day, data = df2, model = "common")
summary(m1)

参考类别更换后的问题 #

用relevel更换参考类别为C.RF3时，出现优化不收敛错误，其他参考类别无此问题：

df2$Model2 <- relevel(df2$Model, ref="C.RF3")
dirig <- DR_data(df2[, c("propActive","propInactive","propMaintenance")], base=1)
m1 <- DirichReg(dirig ~ Model2 + Day, data = df2, model = "common")
summary(m1)

错误信息：

Error in summary.DirichletRegModel(m1) : 
Optimization did not converge in 264 + 2 iterations and exited with code 8

补充尝试 #

曾将model参数改为"alternative"，但该模型同时将模型第一个类别和行为第一个类别设为参考，只能看到非Active行为的模型对比；修改base参数更换行为参考类别会改变显著性结果，不可行。且更换模型后，原common模型正常的场景也会出现相同收敛错误。

解决方案 #

1. 调整优化迭代参数 #

DirichletReg默认迭代次数有限，可通过control参数增加迭代次数、放宽收敛阈值：

# 增加迭代次数，调整收敛阈值
ctrl <- list(maxit = 1000, rel.tol = 1e-6)
m1 <- DirichReg(dirig ~ Model2 + Day, data = df2, model = "common", control = ctrl)
summary(m1)

2. 检查并调整参考类别的数据 #

出现收敛问题的C.RF3可能存在极端数据（比如某行为比例接近0/1），先排查数据特征：

# 查看C.RF3的行为比例分布
subset(df2, Model == "C.RF3") %>% 
  select(propActive, propInactive, propMaintenance) %>% 
  summary()

如果存在极端值，可对原始计数做轻微收缩处理（加0.5）后重新计算比例，避免完全0/1的情况：

# 调整原始计数，重新计算比例
df2 <- df2 %>% 
  mutate(
    Active_adj = Active + 0.5,
    Inactive_adj = Inactive + 0.5,
    Maintenance_adj = Maintenance + 0.5,
    Total_adj = Active_adj + Inactive_adj + Maintenance_adj,
    propActive_adj = Active_adj / Total_adj,
    propInactive_adj = Inactive_adj / Total_adj,
    propMaintenance_adj = Maintenance_adj / Total_adj
  )
# 用调整后的比例建模
dirig_adj <- DR_data(df2[, c("propActive_adj","propInactive_adj","propMaintenance_adj")], base=1)
m1_adj <- DirichReg(dirig_adj ~ Model2 + Day, data = df2, model = "common")
summary(m1_adj)

3. 不修改参考类别，直接指定对比 #

跳过relevel，用emmeans直接指定参考类别做对比，避免因子顺序影响收敛：

# 拟合原模型
dirig <- DR_data(df2[, c("propActive","propInactive","propMaintenance")], base=1)
m1 <- DirichReg(dirig ~ Model + Day, data = df2, model = "common")
# 用emmeans指定以C.RF3为参考做 pairwise 对比
library(emmeans)
emm <- emmeans(m1, ~ Model, type = "response")
contrast(emm, ref = "C.RF3", method = "pairwise")

4. 更换优化器 #

DirichletReg默认用BFGS优化器，可尝试更换为Nelder-Mead等其他优化器：

ctrl <- list(optimizer = "Nelder-Mead", maxit = 1000)
m1 <- DirichReg(dirig ~ Model2 + Day, data = df2, model = "common", control = ctrl)
summary(m1)

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