使用GAMLSS拟合右删截ZAGA模型时遇Ops.Surv错误求助

我有一个数据集：

• 响应变量y取值范围0-100，其中80%取值为0；
• 解释变量age取值范围18-70。

为适配零膨胀特性，我选择用GAMLSS拟合零膨胀伽马模型（ZAGA）。由于y大于100不符合实际需求，我需要对ZAGA分布进行右删截（将大于100的概率集中在100处，而非截断），但执行拟合代码时触发错误：

Ops.Surv(y, mass.p) : Invalid operation on a survival time

我尝试的代码如下：

# 变量定义
y <- c(0,   0,   0,  35,   0,   0,   0,   0, 100, 100,   0,   0,   0,   0,  40,  50,   0,   0,   0,  40,   0,   0,  45 ,  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0, 100,   0,   0,
0 ,  0,   0,  15,   0,   0,   0,   0,   0,   0,  45,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
0,   0, 100,   0,   0,  35, 100,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0 ,  0,   0,   0,   0,  50,   0,  70 , 50,   0,   0,  75,   0,   0, 100,   0,  83,   0,   0,
0 ,  0, 100)
age <- c(60, 47, 63, 44, 54, 60, 54, 40, 58, 57, 49, 63, 53, 46, 52, 41, 39, 41, 51, 33, 54, 66, 47, 65, 44, 59, 32, 42, 46, 58, 31, 56, 25, 49, 36, 53, 61, 40, 40, 38, 47, 64, 57 ,42, 50, 58, 61, 54, 54, 58, 60, 60, 52, 46, 67 ,52, 50, 57, 56, 34, 64, 61, 45, 62, 60, 53, 50, 36, 51, 48, 28, 66, 57, 60, 54, 58, 52, 61 ,62, 52, 57, 47, 42, 49, 58, 42, 41, 54, 57, 38, 41, 59, 55, 30, 60, 46, 57, 52, 56, 56, 54, 64, 65, 55, 39, 61, 66, 64, 43, 57, 58, 59, 53, 43, 50, 44, 62, 50, 58, 48)

# 加载包并准备删截变量
library('gamlss')
library('survival')
library('gamlss.cens')
ysurv <- Surv(y, y!=100, type="right")
gen.cens(ZAGA, type = "right")

# 拟合模型触发错误
Cens <- gamlss(ysurv ~ pb(age), nu.fo = ~ pb(age), family = ZAGArc)

改用正态删截模型或无删截ZAGA模型均可正常运行。

该错误源于Surv对象（来自survival包）与ZAGA这类零膨胀混合分布的兼容性问题：

• ZAGA分布是点质量（0值）+ 伽马连续分布的混合结构；
• Surv对象的删截逻辑是为单一连续分布设计的，gamlss.cens生成的ZAGArc分布在处理Surv对象时，内部会尝试对生存时间对象和零膨胀点质量概率进行运算，而Surv不支持此类操作，因此报错。

使用GAMLSS内置的cen()函数标记删截状态，替代Surv对象，该函数更适配gamlss的自定义混合分布：

library('gamlss')
library('gamlss.cens')

# 1. 创建删截指示变量：标记y=100为右删截，其余为未删截
cens_status <- ifelse(y == 100, "right", "none")

# 2. 生成右删截版ZAGA分布
gen.cens(ZAGA, type = "right")

# 3. 用cen()包装响应变量，拟合模型
Cens <- gamlss(cen(y, cens_status) ~ pb(age), 
               nu.fo = ~ pb(age), 
               family = ZAGArc)

# 查看模型结果
summary(Cens)

• cen()函数是GAMLSS专门用于处理删截数据的工具，能正确识别零膨胀分布的混合结构，避免Surv对象带来的运算冲突；
• 该方法将y=100的样本视为右删截，模型会自动将伽马分布在100到无穷区间的概率累积到100处，符合需求；
• 若需要验证删截效果，可使用plot()或predict()函数查看拟合结果的分布特征。

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