使用JM包拟合联合模型时遇生存时间无效错误求助

问题背景 #

我尝试使用JM包拟合联合纵向-生存模型：

• 纵向部分：用nlme::lme构建随机斜率模型
• 生存部分：用survival::coxph构建Cox模型

调用jointModel(..., method="spline-PH-GH")时触发如下错误：

Error en survreg(Surv(Time, d) ~ ., data = dat):
Invalid survival times for this distribution

可复现代码 #

library(nlme)
library(survival)
library(dplyr)
library(JM)

# 1) 纵向子模型
lmeFit <- lme(
  bmi        ~ time_days*group + sexo + fumador + diabetes + hipercol + hta + educ + cohort,
  random     = ~ time_days | id2,
  data       = xxx,
  method     = "REML",
  na.action  = na.exclude,
  control    = lmeControl(opt = "optim")
)

# 2) 转换为每个受试者一行的生存模型数据集
cox_xxx <- xxx %>%
  group_by(id2) %>%
  summarise(
    # 基线BMI（time_days == 0时）
    bmi           = first(bmi[time_days == 0]),
    # 最后一次观测的时间和事件状态
    time_days_cox = max(time_days),
    event         = event[which.max(time_days)],
    # 基线协变量
    sexo          = first(sexo),
    group         = first(group),
    fumador       = first(fumador),
    diabetes      = first(diabetes),
    hipercol      = first(hipercol),
    hta           = first(hta),
    educ          = first(educ),
    cohort        = first(cohort),
    .groups       = "drop"
  )

# 3) 生存子模型（每个受试者一行）
coxfit <- coxph(
  Surv(time_days_cox, event) ~ group*sexo + fumador + diabetes + hipercol + hta + educ + cohort,
  data  = cox_xxx,
  x     = TRUE,
  model = TRUE
)

# 4) 联合模型
jointFit2 <- jointModel(
  lmeObject  = lmeFit,
  survObject = coxfit,
  timeVar    = "time_days",
  method     = "spline-PH-GH"
)

排查与疑问 #

查阅相关内容后发现，该错误源于JM包内部调用survreg(..., dist="weibull"或"lognormal")，推测是数据中存在一个或多个time_days_cox为0或负值的情况。

请问JM包中是否有内置方法可以绕过该错误，同时不破坏时间的比例性？我尝试添加参数inits=c(1,1,0,0)但无效。

示例数据 #

xxx_ex <- 
structure(list(id2 = structure(c(2L, 7L, 4L, 5L, 9L, 1L, 3L, 
8L, 6L, 10L), levels = c("01_110186", "01_110382", "01_116661", 
"01_220163", "01_221785", "01_223497", "01_224759", "01_225705", 
"01_226139", "19_3441"), class = "factor"), cohort = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L), levels = c("1", "2"), class = "factor"), 
    group = structure(c(1L, 2L, 1L, 3L, 2L, 1L, 1L, 1L, 1L, 1L
    ), levels = c("1", "2", "3"), class = "factor"), seg = c(49.39, 
    64.1, 45.88, 77.69, 43.52, 40.76, 66.36, 61.04, 63.64, 68.13
    ), fecha = structure(c(17981, 18920, 11296, 19054, 15385, 
    15358, 16282, 19111, 11205, 19571), class = "Date"), visit_no = c(2L, 
    6L, 1L, 9L, 3L, 1L, 4L, 9L, 1L, 12L), baseline = structure(c(15867, 
    9967, 11296, 10465, 11873, 15358, 14139, 10645, 11205, 12941
    ), class = "Date"), time_days = c(2114, 8953, 0, 8589, 3512, 
    0, 2143, 8466, 0, 6630), time_years = c(5.78781656399726, 
    24.5119780971937, 0, 23.5154004106776, 9.61533196440794, 
    0, 5.86721423682409, 23.1786447638604, 0, 18.1519507186858
    ), event = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), bmi = c(23.335, 
    25.209, 28.082, 21.218, 22.656, 22.231, 28.387, 28.441, 23.14, 
    17.259), sexo = structure(c(2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 
    2L, 1L), levels = c("1", "2"), class = "factor"), fumador = structure(c(1L, 
    2L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 1L), levels = c("1", "2", 
    "3"), class = "factor"), diabetes = structure(c(2L, 2L, 2L, 
    2L, 2L, 2L, 1L, 2L, 2L, 2L), levels = c("1", "2"), class = "factor"), 
    hipercol = structure(c(2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 
    2L), levels = c("1", "2"), class = "factor"), hta = structure(c(2L, 
    2L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 2L), levels = c("1", "2"), class = "factor"), 
    educ = structure(c(3L, 3L, 1L, 3L, 2L, 3L, 3L, 1L, 3L, 1L
    ), levels = c("0", "1", "3"), class = "factor")), row.names = c(NA, 
-10L), class = "data.frame")

cox_xxx_ex <- 
structure(list(id2 = structure(c(6L, 7L, 10L, 3L, 5L, 2L, 8L, 
9L, 4L, 1L), levels = c("01_117186", "01_220973", "01_221445", 
"01_221580", "01_225818", "19_2411", "19_257", "19_3603", "19_4183", 
"19_886"), class = "factor"), bmi = c(20.895, 26.743, 26.854, 
26.939, 27.732, 20.673, 23.386, 27.531, 19.141, 20.898), time_days_cox = c(4955, 
2719, 5444, 6669, 2690, 7496, 7273, 2129, 0, 4236), event = c(0L, 
1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L), sexo = structure(c(1L, 2L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L), levels = c("1", "2"), class = "factor"), 
    group = structure(c(3L, 1L, 2L, 3L, 3L, 3L, 3L, 2L, 2L, 4L
    ), levels = c("1", "2", "3", "4"), class = "factor"), fumador = structure(c(1L, 
    2L, 2L, 3L, 2L, 2L, 3L, 2L, 3L, 3L), levels = c("1", "2", 
    "3"), class = "factor"), diabetes = structure(c(1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L), levels = "2", class = "factor"), 
    hipercol = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L), levels = "2", class = "factor"), hta = structure(c(2L, 
    1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L), levels = c("1", "2"), class = "factor"), 
    educ = structure(c(1L, 3L, 1L, 1L, 3L, 1L, 1L, 2L, 3L, 1L
    ), levels = c("0", "1", "3"), class = "factor"), cohort = structure(c(2L, 
    2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L), levels = c("1", "2"), class = "factor")), row.names = c(NA, 
-10L), class = "data.frame")

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