求助：在R的rjags包中为Weibull AFT贝叶斯生存模型添加时变协变量

Got it, let's tackle adding time-dependent covariates to your Weibull AFT Bayesian survival model in rjags. The key here is restructuring your data and adjusting the model to account for how covariates change over time—here's a step-by-step breakdown with code:

核心思路：时变协变量的区间拆分 #

Time-dependent covariates like age can't be plugged in like fixed variables because their values change as time passes. The standard fix is to split each individual's follow-up into time intervals where the covariate stays constant. For example:

• If someone starts at age 30, is observed for 8 years, and disposes at year 8, you might split their data into yearly intervals, each with their age that year (30, 31, ..., 37).
• Only the final interval will have an event marker (if they weren't censored).

Step 1: Restructure Your Data into Counting Process Format #

We'll use the survival package's survSplit function to split your data into start-stop intervals. Assuming age increases linearly over time (a safe default if you don't have repeated age measurements), we'll calculate the age for each interval's midpoint:

# Load required packages
library(dplyr)
library(survival)
library(rjags)

# Convert raw data to counting process format
# Your data.s should have: id, duration, event, ageatstart, marital, education, income
data_count <- survSplit(
  Surv(duration, event) ~ .,
  data = data.s,
  cut = seq(0, max(data.s$duration), by = 1), # Split into 1-year intervals (adjust as needed)
  start = "start_time",
  end = "end_time",
  event = "event_interval"
)

# Calculate time-dependent age for each interval (midpoint age)
data_count <- data_count %>%
  mutate(age_timevar = ageatstart + (start_time + end_time)/2)

# Check the restructured data
head(data_count)

Step 2: Update the JAGS Model String #

The model now needs to iterate over intervals (not individuals) and calculate the probability of an event occurring in each interval. Here's the adjusted model, including your fixed covariates plus the time-dependent age:

mod_string_wei_timevar <- "model{
  # Loop over each time interval (not each individual)
  for(k in 1 : n_intervals) {
    # Length of the current interval
    interval_length[k] <- end_time[k] - start_time[k]
    
    # Probability of event in this interval (derived from Weibull survival function)
    event_interval[k] ~ dbern(p_event[k])
    p_event[k] <- 1 - exp( - ( (end_time[k]/lambda[k])^alpha - (start_time[k]/lambda[k])^alpha ) )
    
    # Lambda includes time-dependent age + fixed covariates
    lambda[k] <- exp( beta[1] + beta[2]*age_timevar[k] + beta[3]*marital[k] + beta[4]*education[k] + beta[5]*income[k] )
  }
  
  # Priors for coefficients and shape parameter
  for (i in 1:5) { # 1=intercept, 2=age_timevar, 3=marital, 4=education,5=income
    beta[i] ~ dnorm(0.0, 1.0/1.0e5) # Weakly informative prior
  }
  alpha ~ dgamma(1, 0.0001) # Prior for Weibull shape parameter
  
  # Optional: Calculate median survival for a reference case
  # Example: age=40, marital=1, education=2, income=50000
  lambda_ref <- exp( beta[1] + beta[2]*40 + beta[3]*1 + beta[4]*2 + beta[5]*50000 )
  median_ref <- lambda_ref * (log(2))^(1/alpha)
}"

Key Model Changes: #

• We now loop over intervals (k) instead of individuals (i), since one person can have multiple rows.
• p_event[k] calculates the probability of an event happening in the interval, using the Weibull survival function's properties.
• lambda[k] includes your time-dependent age variable (age_timevar[k]) alongside fixed covariates.

Step 3: Run the Updated Model #

Adjust your data setup and model execution to match the new structure:

# Prepare data list for JAGS
n_intervals <- nrow(data_count)
data_list <- list(
  n_intervals = n_intervals,
  start_time = data_count$start_time,
  end_time = data_count$end_time,
  event_interval = data_count$event_interval,
  age_timevar = data_count$age_timevar,
  marital = data_count$marital,
  education = data_count$education,
  income = data_count$income
)

# Initial values
set.seed(72)
params = c("beta", "alpha", "median_ref") # Include any derived quantities you want
inits1 = function() {
  list(
    beta = rnorm(5, -8.0, 5.0), # 5 beta parameters (adjust prior mean/sd as needed)
    alpha = rgamma(1, 3.0, 4.0)
  )
}

# Run model (increase adapt/burn-in for better convergence)
chain.num = 3
adapt.num = 1000 # Critical for complex models—let JAGS tune samplers
mod.syd_timevar = jags.model(
  textConnection(mod_string_wei_timevar),
  data = data_list,
  inits = inits1,
  n.chains = chain.num,
  n.adapt = adapt.num
)

# Burn-in period
burn.count = 1000
update(mod.syd_timevar, burn.count)

# Sample from posterior
thin.num = 2 # Reduce autocorrelation
iteration.num = 2000
mod.syd_timevar.sim = coda.samples(
  model = mod.syd_timevar,
  variable.names = params,
  n.iter = iteration.num,
  thin = thin.num
)

# Process results
mod.syd_timevar.csim = do.call(rbind, mod.syd_timevar.sim)

# Check convergence and summarize
summary(mod.syd_timevar.sim)
plot(mod.syd_timevar.sim) # Look for stable chains
gelman.diag(mod.syd_timevar.sim) # Ensure potential scale reduction factor is <1.1

Quick Tips for Success #

• Interval granularity: Adjust the by argument in survSplit to balance accuracy and computation time (e.g., 6-month intervals instead of yearly if needed).
• Covariate centering: Center numerical variables like age and income (subtract their mean) to make beta coefficients easier to interpret and improve convergence.
• Convergence checks: Time-dependent models are more complex—always verify chains are stable with plots and Gelman diagnostics.

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