Survival Models

Data Structure

data {
  int<lower=0> N;
  vector<lower=0>[N] time;      // Observed/censored time
  array[N] int<lower=0,upper=1> event;  // 1=event, 0=censored
  matrix[N, K] X;               // Covariates
}

Exponential Model

Stan

parameters {
  real alpha;           // Log baseline hazard
  vector[K] beta;
}
model {
  alpha ~ normal(0, 2);
  beta ~ normal(0, 1);

  for (n in 1:N) {
    real lambda = exp(alpha + X[n] * beta);
    if (event[n] == 1)
      target += exponential_lpdf(time[n] | lambda);
    else
      target += exponential_lccdf(time[n] | lambda);  // Survival
  }
}

JAGS (with censoring)

model {
  for (i in 1:N) {
    is.censored[i] ~ dinterval(t[i], t.cen[i])
    t[i] ~ dexp(lambda[i])
    log(lambda[i]) <- alpha + inprod(X[i,], beta[])
  }
  alpha ~ dnorm(0, 0.25)
  for (k in 1:K) { beta[k] ~ dnorm(0, 1) }
}

Weibull Model

Stan (AFT Parameterization)

parameters {
  real alpha;                    // Intercept (log scale)
  vector[K] beta;
  real<lower=0> shape;           // Weibull shape
}
model {
  alpha ~ normal(0, 5);
  beta ~ normal(0, 2);
  shape ~ exponential(1);

  for (n in 1:N) {
    real mu = alpha + X[n] * beta;
    if (event[n] == 1)
      target += weibull_lpdf(time[n] | shape, exp(mu));
    else
      target += weibull_lccdf(time[n] | shape, exp(mu));
  }
}

JAGS

model {
  for (i in 1:N) {
    is.censored[i] ~ dinterval(t[i], t.cen[i])
    t[i] ~ dweib(shape, lambda[i])
    log(lambda[i]) <- alpha + inprod(X[i,], beta[])
  }
  shape ~ dgamma(1, 0.001)
  alpha ~ dnorm(0, 0.01)
  for (k in 1:K) { beta[k] ~ dnorm(0, 0.01) }
}

Log-Normal Model

Stan

parameters {
  real alpha;
  vector[K] beta;
  real<lower=0> sigma;
}
model {
  for (n in 1:N) {
    real mu = alpha + X[n] * beta;
    if (event[n] == 1)
      target += lognormal_lpdf(time[n] | mu, sigma);
    else
      target += lognormal_lccdf(time[n] | mu, sigma);
  }
}

Piecewise Exponential (Cox-like)

Stan

data {
  int<lower=0> N;
  int<lower=0> J;               // Number of intervals
  vector[J] cuts;               // Cut points
  matrix[N, J] d;               // Time in each interval
  array[N] int<lower=0,upper=1> event;
  array[N] int<lower=1,upper=J> interval;  // Event interval
  matrix[N, K] X;
}
parameters {
  vector[J] log_baseline;       // Log baseline hazard per interval
  vector[K] beta;
}
model {
  log_baseline ~ normal(0, 2);
  beta ~ normal(0, 1);

  for (n in 1:N) {
    real log_hazard = log_baseline[interval[n]] + X[n] * beta;

    // Contribution from all intervals
    for (j in 1:J)
      target += -d[n,j] * exp(log_baseline[j] + X[n] * beta);

    // Event contribution
    if (event[n] == 1)
      target += log_hazard;
  }
}

Frailty Model (Random Effects)

Stan

data {
  int<lower=0> N;
  int<lower=0> G;               // Number of groups
  array[N] int<lower=1,upper=G> group;
  // ... rest of survival data
}
parameters {
  real alpha;
  vector[K] beta;
  real<lower=0> shape;
  vector[G] frailty_raw;        // Non-centered
  real<lower=0> sigma_frailty;
}
transformed parameters {
  vector[G] frailty = sigma_frailty * frailty_raw;
}
model {
  sigma_frailty ~ exponential(1);
  frailty_raw ~ std_normal();

  for (n in 1:N) {
    real mu = alpha + X[n] * beta + frailty[group[n]];
    // ... Weibull likelihood with censoring
  }
}

Generated Quantities

generated quantities {
  // Hazard ratio for 1-unit increase in X[,1]
  real HR = exp(beta[1]);

  // Median survival at X=0
  real median_survival = exp(alpha) * pow(log(2), 1/shape);

  // Survival function at time t=1
  array[N] real S_1;
  for (n in 1:N)
    S_1[n] = exp(-pow(1 / exp(alpha + X[n] * beta), shape));
}

Interpretation

• Weibull shape: <1 decreasing hazard, =1 constant (exponential), >1 increasing
• HR: Hazard ratio (multiplicative effect on hazard)
• AFT: Accelerated failure time (multiplicative effect on survival time)