survival-models

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);
  }
}