plot_surv: Plots Quantiles of the Conditional Hazards

In VAJointSurv: Variational Approximation for Joint Survival and Marker Models

Plots Quantiles of the Conditional Hazards

Description

Plots Quantiles of the Conditional Hazards

Usage

plot_surv(
  time_fixef,
  time_rng,
  x_range,
  fixef_vary,
  vcov_vary,
  frailty_var,
  ps = c(0.025, 0.5, 0.975),
  log_hazard_shift = 0,
  associations,
  xlab = "Time",
  ylab = "Hazard",
  ders = NULL,
  newdata = NULL,
  ...
)

Arguments

time_fixef	the time-varying fixed effects. See .e.g. poly_term. This is for the baseline hazard. Note that many basis expansions have boundary knots. It is important that these are set to cover the full range of survival times including time zero for some expansions.
time_rng	an expansion or a list of expansions for the time-varying random effects of the markers. See marker_term.
x_range	two dimensional numerical vector with the range the hazard should be plotted in.
fixef_vary	fixed effect coefficients for time_fixef.
vcov_vary	covariance matrix for the expansion or expansions in time_rng.
frailty_var	variance of the frailty.
ps	quantiles to plot.
log_hazard_shift	possible shift on the log hazard.
associations	association parameter for each time_rng or possible multiple parameters for each time_rng if ders is supplied.
xlab, ylab, ...	arguments passed to matplot.
ders	a list with integer vectors for how the survival outcome is linked to the markers. 0 implies present values, -1 is integral of, and 1 is the derivative. NULL implies the present value of the random effect for all markers.
newdata	data.frame with data for the weights if any.

Value

A list containing data for plotting.

Examples

# load in the data
library(survival)
data(pbc, package = "survival")

# re-scale by year
pbcseq <- transform(pbcseq, day_use = day / 365.25)
pbc <- transform(pbc, time_use = time / 365.25)

# create the marker terms
m1 <- marker_term(
  log(bili) ~ 1, id = id, data = pbcseq,
  time_fixef = bs_term(day_use, df = 5L),
  time_rng = poly_term(day_use, degree = 1L, raw = TRUE, intercept = TRUE))
m2 <- marker_term(
  albumin ~ 1, id = id, data = pbcseq,
  time_fixef = bs_term(day_use, df = 5L),
  time_rng = poly_term(day_use, degree = 1L, raw = TRUE, intercept = TRUE))

# base knots on observed event times
bs_term_knots <-
  with(pbc, quantile(time_use[status == 2], probs = seq(0, 1, by = .2)))

boundary <- c(bs_term_knots[ c(1, length(bs_term_knots))])
interior <- c(bs_term_knots[-c(1, length(bs_term_knots))])

# create the survival term
s_term <- surv_term(
  Surv(time_use, status == 2) ~ 1, id = id, data = pbc,
  time_fixef = bs_term(time_use, Boundary.knots = boundary, knots = interior))

# expansion of time for the fixed effects in the survival term
time_fixef <- s_term$time_fixef
# expansion of time for the random effects in the marker terms
time_rng <- list(m1$time_rng, m2$time_rng)
# no frailty
frailty_var <- matrix(0L,1)
# var-covar matrix for time-varying random effects
vcov_vary <- c(0.9658, 0.0954, -0.1756, -0.0418, 0.0954, 0.04, -0.0276,
               -0.0128, -0.1756, -0.0276, 0.1189, 0.0077, -0.0418, -0.0128,
               0.0077, 0.0057) |> matrix(4L)
# coefficients for time-varying fixed effects
fixef_vary <- c(1.0495, -0.2004, 1.4167, 1.255, 2.5007, 4.8545, 4.7889)
# association parameters
associations <- c(0.8627, -3.2358, 0.1842)
# constant shift on the log-hazard scale
log_hazard_shift <- -4.498513
# specify how the survival outcome is linked with markers
ders = list(0L, c(0L, -1L))

# plot the hazard with pointwise quantiles
plot_surv(
  time_fixef = time_fixef,
  time_rng = time_rng,
  x_range = c(0, 5), vcov_vary = vcov_vary, frailty_var = frailty_var,
  ps = c(.25, .5, .75), log_hazard_shift = log_hazard_shift,
  fixef_vary = fixef_vary, associations = associations, ders = ders)

VAJointSurv documentation built on Aug. 14, 2022, 9:05 a.m.