In boennecd/dynamichazard: Dynamic Hazard Models using State Space Models

knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>", fig.path = "README-", dpi = 124, 
  error = FALSE)

dynamichazard

The goal of dynamichazard is to estimate time-varying effects in survival analysis. The time-varying effects are estimated with state space models where the coefficients follow a given order random walk. The advantageous of using state space models is that you can extrapolate beyond the last observed time period. For more details, see the ddhazard vignette at https://cran.r-project.org/web/packages/dynamichazard/vignettes/ddhazard.pdf.

The particle filter and smoother methods can estimate more general models then the random walk model. See the /examples directory for some examples.

Installation

You can install dynamichazard from github with:

# install.packages("remotes")
remotes::install_github("boennecd/dynamichazard")

You can also download the package from CRAN by calling:

install.packages("dynamichazard")

Example - ddhazard

I will use the aids data set from the JMbayes package. The data set is from a a randomized clinical trial between two drugs for HIV or aids patients. The event is when the patient die. Patient are right censored at the end of the study. The data set is longitudinal/panel format with rows for patient. Though, the only longitudinal variable is the CD4 count (T-cells count) which is presumably affected by the drug. Thus, I will not use it in the model. The other the other columns of interest are:

AZT is one of the two enrollment criteria. It indicates whether the patient was enrolled due to intolerance to the drug zidovudine or whether the drug failed prior to the study start.
prevOI is the other enrollment criteria. Patients are enrolled either due AIDS diagnosis or two CD4 counts of 300 or fewer. The variable indicates which is the case.
drug is either ddC or ddI depending on which of the two drugs the patient is randomly assigned to.
gender.

The analysis is given below with comments:

library(dynamichazard)
library(JMbayes) # Contain the aids data set

# We remove the data we dont neeed
aids <- aids[aids$Time == aids$stop, ]
aids <- aids[, !colnames(aids) %in% c("Time", "death", "obstime", "CD4")]

# A look at head of data
head(aids)
max(aids$stop)                  # Last observation time
max(aids$stop[aids$event == 1]) # Last person with event

# Fit model with extended Kalman filter
fit <- ddhazard(
  Surv(stop, event) ~ ddFixed_intercept() + ddFixed(AZT) + gender + 
    ddFixed(drug) + ddFixed(prevOI),
  aids,
  model = "exponential",          # piecewise constant exponentially distributed 
                                  # arrivals times
  by = .5,                        # Length of time intervals in state space 
                                  # model
  max_T = 19,                     # Last period we observe when modeling
  Q = .1^2,                       # Covariance matrix for state equation in 
                                  # first iteration
  Q_0 = 1,                        # Covariance matrix for the prior
  control = ddhazard_control(
    eps = .001,                   # tolerance for EM-algorithm
    LR = .9,                      # Learning rate
    n_max = 100                   # Max number iterations in EM
    ))

# Plot the estimates. Dashed lines are 95% confidence bounds
plot(fit)

# Bootstrap the estimates
set.seed(87754771)
boot_out <- ddhazard_boot(fit, R = 1000) # R is number of bootstrap samples

# Plot bootstrap estimates. Dashed lines are 2.5% and 97.5% quantiles of the 
# bootstrap estimates. Transparent lines are bootstrap estimates
plot(fit, ddhazard_boot = boot_out)

Bootstrapping only slightly changes the confidence bounds. An example of a paper analyzing the CD4 count can be found in [@guo2004]. They also fit a static model (time-invariant coefficients) of the survival times with an exponential model. The estimates are comparable with those above as expected.

Example - particle filter and smoother

A particle filter and smoother is also included in the package. The computational complexity of these methods match those of the extended Kalman filter but with a much larger constant. Below, I fit a model for the aids data. We only use a time-varying effect for gender this time.

set.seed(20170907)
pf_fit <- PF_EM(
  Surv(stop, event) ~ ddFixed_intercept() + ddFixed(AZT) + gender + 
    ddFixed(drug) + ddFixed(prevOI),
  aids, model = "exponential",
  by = .5,  max_T = 19, Q = .5^2, Q_0 = 1,
  control = PF_control(
    n_threads = 4, # I use a quad-core machine
    # set number of particles
    N_fw_n_bw = 1000, N_first = 1000,
    N_smooth = 1, # Does not matter with Brier_O_N_square
    smoother = "Brier_O_N_square", # Select smoother
    eps = 1e-5, n_max = 1000, est_a_0 = FALSE, nu = 6L, 
    covar_fac = 1.2, averaging_start = 150L)
  #, trace = 1 # comment back to get feedback during estimation
  )

# Compare estimates of Q
pf_fit$Q / .5
fit$Q
plot(pf_fit)
plot(pf_fit$log_likes) # log-likelihoods
logLik(pf_fit)

# better estimate of final log-likelihood
logLik(PF_forward_filter(pf_fit, N_fw = 10000, N_first = 10000))

For more details, see the "Particle filters in the dynamichazard package" vignette at https://cran.r-project.org/web/packages/dynamichazard/vignettes/Particle_filtering.pdf