The Hazard-Response model

The hazard-response model is defined through the system of ODEs (Christen and Rubio 2023):

\[ \begin{cases} h'(t) = \lambda h(t) \left(1 - \dfrac{h(t)}{\kappa}\right) - \alpha q(t) h(t), & h(0) = h_0 \\ q'(t) = \beta q(t) \left( 1- \frac{q(t)}{\kappa} \right) -\alpha q(t) h(t) , & q(0) = q_0 \\ H'(t) = h(t), & H(0) = 0, \end{cases} \]

with \(\lambda>0\), \(\alpha \geq 0\), \(\beta>0\), \(\kappa>0\), \(h_0>0\), and \(q_0>0\). This model represents a particular version of the competitive Lotka-Volterra model (Christen and Rubio 2023). This model formulation assumes that the hazard function \(h(t)\) is in competition with a response \(q(t)\) (associated with the immune system and any treatment or intervention at the population level).

Routines

source("routines.R")

Examples

Hazard and Survival function

# Parameter values
lambda <- 1.8
kappa <- 0.1
alpha <- 6
beta <- 4.8
h0 <- 1e-2
q0 <- 1e-6

# True parameter values for plots
trueparf <- c(lambda, kappa, alpha, beta, kappa, alpha, h0, q0)
# True parameter values for simulations
truepar <- c(lambda, kappa, alpha, beta)

# Hazard function: Solver
paramsHR  <-  c(lambda, kappa, alpha, beta)
initHT0      <- c(h = h0, q = q0, H = 0 )
times = seq(0,36,by = 0.1)
out <- ode(initHT0, times, hazmodHR, paramsHR, method = "lsode", jacfunc = jacODE, jactype = "fullusr")


# Hazard function
plot(as.matrix(out[,c(1,2)]), xlim = c(0,max(times)), ylim = c(0,0.125), 
     type = "l", lwd = 2,
     xlab = "Time", ylab = "Hazard", main = "Hazard-Response ODE", 
     cex.axis = 1.5, cex.lab = 1.5)
points(as.matrix(out[,c(1,3)]), lty=2, type="l", lwd = 2)
legend("bottomright", legend=c("Hazard","Response"), lwd = c(2,2), lty = c(1,2))

# Survival function
plot(as.vector(out[,1]),exp(-as.vector(out[,4])), xlim = c(0,max(times)), ylim = c(0,1), type = "l", lwd = 2,
     xlab = "Time", ylab = "Survival", main = "Hazard-Response ODE", cex.axis = 1.5, cex.lab = 1.5)

Approximated simulation from the Hazard-Response model

Christen and Rubio (2023) proposed a method for approximately simulating from the Hazard-Response model. The following code shows how to simulate from this model using the command simHR. We compare the true hazard model against the Kaplan-Meier estimator of the approximate simulations.

# Approximate simulation
set.seed(1234)
sim <- simHR(n = 5000, par = truepar, tmin = 0, tmax = 150, step = 1e-3)$sim

# Censoring time
cens <- 20

# Vital status
status <- ifelse(sim < cens, 1, 0)

mean(status)
## [1] 0.7482
# Observed times
time <- ifelse(sim < cens, sim, cens)

# Kaplan-Meier estimator for the simulated times
kmsim <- survfit(Surv(time, status) ~ 1)

# Comparison
plot(as.vector(out[,1]),exp(-as.vector(out[,4])), xlim = c(0,max(times)), ylim = c(0,1), type = "l", lwd = 2,
     xlab = "Time", ylab = "Survival", main = "Hazard-Response ODE", cex.axis = 1.5, cex.lab = 1.5)
points(kmsim$time, kmsim$surv, type = "l", col = "gray", lwd = 2, lty = 1, ylim = c(0,1),
     xlab = "Time", ylab = "Survival")
points(kmsim$time, kmsim$surv, type = "l", col = "gray", lwd = 2, lty = 2)
legend("topright", legend = c("True model","Kaplan-Meier"), col = c("black","gray"), lwd = c(2,2), lty = c(1,2))

Christen, J. A., and F. J. Rubio. 2023. “Dynamic Survival Analysis: Modelling the Hazard Function via Ordinary Differential Equations.” Preprint NA: NA–.