Weibull PH模型生存数据模拟报错：integrate返回长度异常

报错核心是多处维度不匹配，导致函数返回向量而非integrate/uniroot要求的标量，具体问题及修复如下：

• 基线协变量索引错误：原代码baseline_cov[,i]%*%(gamma)取的是第i列协变量，而非第i个个体的协变量行，应改为baseline_cov[i,] %*% gamma，否则会得到向量而非标量。
• 生存概率使用错误：invS中调用log(surv_prob)用到了整个向量，而循环中传入的是单个个体的u = surv_prob[i]，应改为log(u)，避免返回向量。
• 不必要的Vectorize：integrate本身会逐点调用目标函数，不需要额外用Vectorize(h)，反而可能引发维度问题，直接传入h即可。

修正后的完整代码 #

### VALUES TO TEST ####
library(MASS) # 加载mvrnorm所在包
N <- 100 # 补充原代码缺失的N定义
vcov <- matrix(c(1,0.5,0.5,1), 2,2) # 补充原代码缺失的vcov定义

survival_probability= runif(N,0,1) #Uniform (0,1)
raneff <- mvrnorm(N, mu=c(0,0), Sigma = vcov)

vcov_covariates <-matrix(c(1,0,0.5, 0,1,0.5, 0.5,0.5,1),3,3)
covariate_vector <- mvrnorm(N, mu = c(0,0,0), vcov_covariates )

baseline_covariates_w1 <- covariate_vector[,1]
baseline_covariates_w2 <- covariate_vector[,2]
X2 <- covariate_vector[,3]
alpha = 0.2#correlation between the processes

lambda = 0.2
beta = c(0.5, 3, 1,1)
gamma = c(1,1)

getSurvTimeWeibullModel2 <-function(sigma.t, alpha, surv_prob, beta, random_effects, gamma, X2, baseline_cov, N, max.FUtime, method="Weibull-PH"){
  #' @description Get's survival time from Weibull distribution
  #' @param sigma.t scale parameter for the Weibull baseline risk function;
  #' @param alpha association parameter
  #' @param surv_prob double between 0 and 1 N by 1 vector
  #' @param beta vector of fixed effects Should contain beta_0, beta_1, beta_2 and beta_3
  #' @param random_effects # a matrix of N by 2
  #' @param gamma #fixed effects for baseline cov in survival model p by 1
  #' @param X2 # vector of 1 x N or N by 1
  #' @param baseline_cov N x p
  #' @param N number of individuals
  #' @return double, survival time a vector of N by 1
  
  invS <- function (t, u, i) {
    timedependentMeasurement <- function(v){
      # 计算与时变标记的关联项
      alpha*((beta[1] + random_effects[i, 1]) + (beta[2] + random_effects[i,2])*v + beta[3]*X2[i] + beta[4]*X2[i]*v)
    }
    
    h <- function (s) {
      TD.i <- timedependentMeasurement(s)
      if(tolower(method)=="weibull-ph"){ # 统一小写避免大小写匹配问题
        exp(log(sigma.t) + (sigma.t - 1) * log(s) + 
              baseline_cov[i,] %*% gamma + TD.i) # 修正协变量索引
      } else {
        stop("Unsupported method") # 增加错误处理
      }
    }
    
    # 去掉Vectorize，直接传入h；用log(u)代替log(surv_prob)
    integrate(h, lower=0, upper = t)$value + log(u)
  }
  
  trueTimes <- numeric(N)
  for (i in 1:N) {
    Root <- uniroot(invS, interval = c(1e-05, max.FUtime), u = surv_prob[i],
                    i = i)$root
    trueTimes[i] <- Root
  }
  return(trueTimes)
}

# 调用函数
getSurvTimeWeibullModel2(sigma.t = 1, alpha, surv_prob=survival_probability, beta, 
                         random_effects=raneff, gamma, X2, 
                         baseline_cov=cbind(baseline_covariates_w1,baseline_covariates_w2), 
                         N=N, max.FUtime=100)

额外补充说明 #

1. 原代码缺失N和vcov的定义，修正时补充了这两个变量，确保代码可直接运行。
2. 增加了method参数的大小写不敏感处理，避免因大小写输入错误引发问题。
3. 给h函数增加了默认分支的错误提示，提升代码鲁棒性。

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