使用R语言simsurv包混合Weibull模型遇端点符号错误求助

核心错误：协变量与系数的混淆 #

你犯了一个关键的参数传递错误：把作为协变量系数的frailty_beta错误地放进了协变量数据框x中，同时未通过betas参数明确指定协变量的效应大小。这导致simsurv默认所有协变量的系数为0，完全偏离了Stata中的模型设定，进而触发uniroot求解失败的错误。

在你的代码中，z是个体水平的frailty变量，frailty_beta=0.5是该变量的系数——正确的做法是将z放入x，并通过betas参数指定其系数，而非把系数本身作为协变量传入。

修正后的代码 #

# Load relevant packages
library(survival)
library(simsurv)

# Set sample size for the truth (i.e., a large number)
ss_truth <- 100000

# Sample age and frailty
age <- rnorm(ss_truth, 60, 6)
z <- rnorm(ss_truth, 0, 1)  # 个体的frailty值
t_lim_truth <- 80

# Other parameters
dist <- "weib"
lambda1 <- 0.6
lambda2 <- 0.7
gamma1 <- 1.8
gamma2 <- 0.5
pmix <- 0.3
frailty_beta <- 0.5  # frailty变量的系数

# Set seed for reproducibility
set.seed(1)

# Adjust for mixture versus non-mixture models
pmixopt <- pmix
mixopt <- TRUE

# Adjust length of lambda for mixture versus non-mixture models
lambda_in <- if(is.na(lambda2)) lambda1 else c(lambda1, lambda2)

# Adjust length of gamma for mixture versus non-mixture models
gamma_in <- if(is.na(gamma2)) gamma1 else c(gamma1, gamma2)

# Generate "truth" - 修正协变量与参数传递
td <- simsurv(
  dist = dist,
  lambda = lambda_in,
  gamma = gamma_in,
  pmix = pmixopt,
  maxt = t_lim_truth,
  x = data.frame(frailty = z),  # 仅传入个体协变量值
  betas = c(frailty = frailty_beta),  # 明确指定协变量系数
  mixture = mixopt
)

为什么原代码会触发uniroot错误？ #

当gamma2=0.5（对应递减的风险率）时，Weibull分布的生存函数下降极慢。加上你未正确应用frailty效应，部分个体的生存概率在maxt=80时仍远高于随机生成的生存概率阈值u，导致uniroot函数无法找到符号相反的端点（两端均为正），从而报错。修正协变量效应后，frailty会改变个体的风险水平，让生存时间的分布更符合预期，解决uniroot的求解问题。

备选方案：若仍报错，调整uniroot搜索区间 #

如果修正后仍遇到相同错误，可修改simsurv的内部函数，扩大uniroot的搜索上限：

# 重写simweib函数，调整uniroot的upper参数
simweib <- function (x, betas, lambda, gamma, pmix, maxt, mixture) {
  if (!mixture) {
    linear.predictor <- as.matrix(x) %*% betas
    lambda_i <- lambda * exp(linear.predictor)
    u <- stats::runif(n = nrow(x))
    t <- (-log(u)/lambda_i)^(1/gamma)
    event <- as.integer(t <= maxt)
    t <- pmin(t, maxt)
    data.frame(event = event, time = t)
  } else {
    linear.predictor <- as.matrix(x) %*% betas
    lambda_i <- lambda * exp(linear.predictor)
    gamma_i <- gamma
    pmix_i <- pmix
    u <- stats::runif(n = nrow(x))
    v <- stats::runif(n = nrow(x))
    t <- rep(NA, nrow(x))
    event <- rep(NA, nrow(x))
    for (i in seq_len(nrow(x))) {
      if (v[i] <= pmix_i) {
        lambda_i1 <- lambda_i[1, i]
        gamma_i1 <- gamma_i[1]
        rootfn_surv <- function(t, survival, lambda, gamma) {
          exp(-lambda * t^gamma) - survival
        }
        surv <- u[i]
        # 扩大uniroot的upper区间到1e5
        t_i <- stats::uniroot(rootfn_surv, survival = surv, lambda = lambda_i1, 
                              gamma = gamma_i1, lower = 0, upper = 1e5)$root
      } else {
        lambda_i2 <- lambda_i[2, i]
        gamma_i2 <- gamma_i[2]
        rootfn_surv <- function(t, survival, lambda, gamma) {
          exp(-lambda * t^gamma) - survival
        }
        surv <- u[i]
        # 扩大uniroot的upper区间到1e5
        t_i <- stats::uniroot(rootfn_surv, survival = surv, lambda = lambda_i2, 
                              gamma = gamma_i2, lower = 0, upper = 1e5)$root
      }
      event[i] <- as.integer(t_i <= maxt)
      t[i] <- pmin(t_i, maxt)
    }
    data.frame(event = event, time = t)
  }
}

# 替换simsurv包中的simweib函数
assignInNamespace("simweib", simweib, ns = "simsurv")

# 重新生成数据
td <- simsurv(
  dist = dist,
  lambda = lambda_in,
  gamma = gamma_in,
  pmix = pmixopt,
  maxt = t_lim_truth,
  x = data.frame(frailty = z),
  betas = c(frailty = frailty_beta),
  mixture = mixopt
)

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