使用avg_predictions()时flexsurv模型参考水平出现NA值问题

问题描述 #

拟合多种参数生存模型（Weibull、Gamma、广义Gamma等）后，使用marginaleffects包的avg_predictions()计算预测平均生存时间时，因子变量的参考水平（如示例中sex=1）对应的标准误、z值、p值及置信区间均为NA，其他水平的不确定性估计正常显示。

复现代码 #

library(survival)
library(flexsurv)
library(marginaleffects)

data(lung)

# 拟合多种参数生存模型
fmodel1 <- list(
  "Weibull" = flexsurvreg(Surv(time, status) ~ factor(sex), data = lung, dist = "weibull"),
  "Gamma" = flexsurvreg(Surv(time, status) ~ factor(sex), data = lung, dist = "gamma"),
  "Gengamma" = flexsurvreg(Surv(time, status) ~ factor(sex), data = lung, dist = "gengamma"),
  "llogis" = flexsurvreg(Surv(time, status) ~ factor(sex), data = lung, dist = "llogis"),
  "llnorm" = flexsurvreg(Surv(time, status) ~ factor(sex), data = lung, dist = "lognormal"),
  "gompertz" = flexsurvreg(Surv(time, status) ~ factor(sex), data = lung, dist = "gompertz")
)

# 计算平均预测
avg_predictions(fmodel1$Weibull, variables = "sex")

输出结果 #

sex Estimate Std. Error   z Pr(>|z|)    S 2.5 % 97.5 %
   1      331         NA  NA       NA   NA    NA     NA
   2      490       64.5 7.6   <0.001 44.9   363    616

Columns: sex, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
Type:  response 

原因分析 #

这是marginaleffects处理因子参考水平的逻辑导致：默认使用虚拟变量编码时，参考水平的系数被设为0，avg_predictions()在计算不确定性时依赖参数协方差矩阵与预测函数的梯度，但参考水平的梯度计算会因编码方式被跳过，导致相关不确定性指标缺失。

解决方案 #

显式使用datagrid()生成包含所有因子水平的新数据框，强制marginaleffects为每个水平计算完整的预测及不确定性：

# 显式指定所有sex水平的新数据
avg_predictions(fmodel1$Weibull, newdata = datagrid(sex = unique(lung$sex)))

验证方法 #

可先通过flexsurv原生的predict()函数确认模型本身能返回参考水平的标准误：

predict(fmodel1$Weibull, newdata = data.frame(sex = factor(c(1,2))), type = "response", se.fit = TRUE)

若该输出包含两个水平的标准误，说明问题确实出在marginaleffects的默认处理逻辑，使用datagrid()即可解决。

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