竞争事件下如何用模拟路径结合Schoenfeld残差检验PH假设

Hey there! I've worked through similar problems before, so let's break down how to replicate that SAS-style simulation-based proportional hazards (PH) test for your Fine-Gray model in R—including generating those simulated paths and the comparison plot you mentioned.

1. Recap: Get Your Baseline Model & Observed Schoenfeld Residuals #

First, let's start with the code you likely already have to fit the Fine-Gray model and extract residuals. I'll use a generic example dataset structure (adjust to match your actual data):

library(cmprsk)

# Example data (replace with your actual dataset)
set.seed(123)
n <- 500
dat <- data.frame(
  time = rexp(n, rate = 0.1),
  status = sample(c(1,2), n, replace = TRUE, prob = c(0.7, 0.3)),
  x = rnorm(n)
)

# Fit Fine-Gray model
fit <- crr(ftime = dat$time, fstatus = dat$status, cov1 = dat$x)

# Extract Schoenfeld residuals
schoen_resid <- residuals(fit, type = "schoenfeld")

# Calculate observed residual process (mean residual at each event time for the target event)
event_times <- sort(unique(dat$time[dat$status == 1]))
obs_process <- sapply(event_times, function(t) {
  mean(schoen_resid[dat$time == t & dat$status == 1, ])
})

2. Generate Simulated Paths Under the PH Null Hypothesis #

SAS's simulation approach uses parametric bootstrapping: we generate 1000 datasets that follow the PH assumption (using the fitted model's parameters), fit a Fine-Gray model to each, and extract the residual process for each simulation.

Here's a function to handle this simulation loop (note: this can take a minute for large datasets/1000 simulations—we can add parallelization later if needed):

simulate_finegray_ph <- function(fit, dat, n_sim = 1000) {
  beta_hat <- fit$coef
  base_cum_haz <- fit$cumhaz
  event_times <- base_cum_haz$time
  cum_haz_target <- base_cum_haz$hazard
  
  # Estimate cumulative hazard for competing event (simplified approach)
  fit_competing <- crr(ftime = dat$time, fstatus = dat$status, cov1 = dat$x, failcode = 2)
  cum_haz_competing <- fit_competing$cumhaz$hazard
  
  sim_processes <- list()
  
  for(i in 1:n_sim) {
    # Generate linear predictor from fitted model
    lin_pred <- as.matrix(dat[, names(beta_hat)]) %*% beta_hat
    
    # Simulate time to target event
    u_target <- runif(nrow(dat))
    t_target <- sapply(u_target, function(u) {
      min(event_times[cum_haz_target * exp(lin_pred) >= -log(u)])
    })
    
    # Simulate time to competing event
    u_competing <- runif(nrow(dat))
    t_competing <- sapply(u_competing, function(u) {
      min(event_times[cum_haz_competing >= -log(u)])
    })
    
    # Determine final time and status
    sim_time <- pmin(t_target, t_competing)
    sim_status <- ifelse(t_target < t_competing, 1, 2)
    
    # Fit model to simulated data
    sim_fit <- crr(ftime = sim_time, fstatus = sim_status, cov1 = dat[, names(beta_hat)])
    
    # Extract residual process (aligned to observed event times)
    sim_resid <- residuals(sim_fit, type = "schoenfeld")
    sim_process <- sapply(event_times, function(t) {
      mean(sim_resid[sim_time == t & sim_status == 1, ], na.rm = TRUE)
    })
    
    sim_processes[[i]] <- sim_process
  }
  
  # Convert list to matrix for easier analysis
  sim_matrix <- do.call(rbind, sim_processes)
  return(list(obs_process = obs_process, sim_matrix = sim_matrix, event_times = event_times))
}

# Run 1000 simulations (adjust n_sim as needed)
set.seed(456)
sim_results <- simulate_finegray_ph(fit, dat, n_sim = 1000)

3. Calculate the Simulation-Based p-Value #

To get a p-value, we compare the "extremeness" of our observed residual process to the simulated ones. A common test statistic is the maximum absolute value of the residual process:

# Calculate observed test statistic
obs_stat <- max(abs(sim_results$obs_process))

# Calculate test statistic for each simulated path
sim_stats <- apply(sim_results$sim_matrix, 1, function(x) max(abs(x), na.rm = TRUE))

# Compute p-value (proportion of simulated stats >= observed stat)
p_value <- mean(sim_stats >= obs_stat)
cat("Simulation-based PH test p-value:", round(p_value, 2), "\n")

This should give you a p-value similar to the 0.55 you got from SAS.

4. Plot Observed vs. Top 20 Simulated Paths #

Let's use ggplot2 to create the comparison plot you mentioned:

library(ggplot2)
library(tidyr)

# Format observed data for plotting
obs_df <- data.frame(
  time = sim_results$event_times,
  value = sim_results$obs_process,
  type = "Observed",
  group = "Observed"
)

# Format top 20 simulated paths
sim_top20_df <- as.data.frame(sim_results$sim_matrix[1:20, ]) %>%
  mutate(sim_id = factor(1:20)) %>%
  pivot_longer(cols = -sim_id, names_to = "time", values_to = "value") %>%
  mutate(
    time = as.numeric(time),
    type = "Simulated",
    group = sim_id
  )

# Combine datasets
plot_df <- rbind(obs_df, sim_top20_df)

# Create the plot
ggplot(plot_df, aes(x = time, y = value, color = type, group = group)) +
  geom_line(aes(alpha = type)) +
  scale_alpha_manual(values = c(Observed = 1, Simulated = 0.4)) +
  scale_color_manual(values = c(Observed = "#e74c3c", Simulated = "#95a5a6")) +
  labs(
    x = "Event Time",
    y = "Mean Schoenfeld Residual",
    title = "Observed vs. Simulated Residual Processes",
    subtitle = "Top 20 Simulated Paths Under PH Null Hypothesis"
  ) +
  theme_minimal() +
  theme(legend.position = "top")

Pro Tip: Speed Up Simulations #

If 1000 simulations are too slow, use parallel computing with foreach and doParallel:

library(foreach)
library(doParallel)

cl <- makeCluster(4) # Use 4 cores (adjust to your machine)
registerDoParallel(cl)

# Rewrite the simulation loop to use foreach
sim_processes <- foreach(i = 1:1000, .packages = "cmprsk") %dopar% {
  beta_hat <- fit$coef
  base_cum_haz <- fit$cumhaz
  event_times <- base_cum_haz$time
  cum_haz_target <- base_cum_haz$hazard
  
  fit_competing <- crr(ftime = dat$time, fstatus = dat$status, cov1 = dat$x, failcode = 2)
  cum_haz_competing <- fit_competing$cumhaz$hazard
  
  lin_pred <- as.matrix(dat[, names(beta_hat)]) %*% beta_hat
  
  u_target <- runif(nrow(dat))
  t_target <- sapply(u_target, function(u) {
    min(event_times[cum_haz_target * exp(lin_pred) >= -log(u)])
  })
  
  u_competing <- runif(nrow(dat))
  t_competing <- sapply(u_competing, function(u) {
    min(event_times[cum_haz_competing >= -log(u)])
  })
  
  sim_time <- pmin(t_target, t_competing)
  sim_status <- ifelse(t_target < t_competing, 1, 2)
  
  sim_fit <- crr(ftime = sim_time, fstatus = sim_status, cov1 = dat[, names(beta_hat)])
  sim_resid <- residuals(sim_fit, type = "schoenfeld")
  sim_process <- sapply(event_times, function(t) {
    mean(sim_resid[sim_time == t & sim_status == 1, ], na.rm = TRUE)
  })
  
  sim_process
}

stopCluster(cl)
sim_matrix <- do.call(rbind, sim_processes)

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