coxsimSpline: Simulate quantities of interest for penalized splines from Cox Proportional Hazards models

Description

coxsimSpline simulates quantities of interest from penalized splines using multivariate normal distributions.

Usage

coxsimSpline(obj, bspline, bdata, qi = "Relative Hazard", Xj = 1, Xl = 0,
  nsim = 1000, ci = 0.95, spin = FALSE, extremesDrop = TRUE)

Arguments

obj

a coxph class fitted model object with a penalized spline. These can be plotted with simGG.

bspline

a character string of the full pspline call used in obj. It should be exactly the same as how you entered it in coxph. You also need to enter a wh

bdata

a numeric vector of splined variable's values.

quantity of interest to simulate. Values can be "Relative Hazard", "First Difference", "Hazard Ratio", and "Hazard Rate". The default is qi = "Relative Hazard". Think carefully before using

numeric vector of fitted values for b to simulate for.

numeric vector of values to compare Xj to. Note if qi = "Relative Hazard" or "Hazard Rate" only Xj is relevant.

nsim

the number of simulations to run per value of Xj. Default is nsim = 1000.

the proportion of simulations to keep. The default is ci = 0.95, i.e. keep the middle 95 percent. If spin = TRUE then ci is the confidence level of the shortest probability interval. Any value from 0 through 1 may be

spin

logical, whether or not to keep only the shortest probability interval rather than the middle simulations. Currently not supported for hazard rates.

extremesDrop

logical whether or not to drop simulated quantity of interest values that are Inf, NA, NaN and $> 1000000$ for spin = FALSE or $> 800$ for spin = TRUE. These values are difficult to plot

Value

a simspline object

Details

Simulates relative hazards, first differences, hazard ratios, and hazard rates for penalized splines from Cox Proportional Hazards models. These can be plotted with simGG. A Cox PH model with one penalized spline is given by: $$h(t|\mathbf{X}_{i})=h_{0}(t)\mathrm{e}^{g(x)}$$

where $g(x)$ is the penalized spline function. For our post-estimation purposes $g(x)$ is basically a series of linearly combined coefficients such that:

$$g(x) = \beta_{k_{1}}(x)_{1+} + \beta_{k_{2}}(x)_{2+} + \beta_{k_{3}}(x)_{3+} + \ldots + \beta_{k_{n}}(x)_{n+}$$

where $k$ are the equally spaced spline knots with values inside of the range of observed $x$ and $n$ is the number of knots.

We can again draw values of each $\beta_{k_{1}},\: \ldots \beta_{k_{n}}$ from the multivariate normal distribution described above. We then use these simulated coefficients to estimates quantities of interest for a range covariate values. For example, the first difference between two values $x_{j}$ and $x_{l}$ is:

$$%\triangle h_{i}(t) = (\mathrm{e}^{g(x_{j}) - g(x_{l})} - 1) * 100$$ {FD(h[i](t)) = (exp(g(x[j]) - g(x[l])) - 1) * 100}

Relative hazards and hazard ratios can be calculated by extension.

Currently coxsimSpline does not support simulating hazard rates form multiple stratified models.

References

Luke Keele, "Replication data for: Proportionally Difficult: Testing for Nonproportional Hazards In Cox Models", 2010, http://hdl.handle.net/1902.1/17068 V1 [Version].

King, Gary, Michael Tomz, and Jason Wittenberg. 2000. ''Making the Most of Statistical Analyses: Improving Interpretation and Presentation.'' American Journal of Political Science 44(2): 347-61.

Liu, Ying, Andrew Gelman, and Tian Zheng. 2013. ''Simulation-Efficient Shortest Probability Intervals.'' Arvix. http://arxiv.org/pdf/1302.2142v1.pdf.

Examples

Run this code

# Load Carpenter (2002) data
data("CarpenterFdaData")

# Load survival package
library(survival)

# Run basic model
# From Keele (2010) replication data
M1 <- coxph(Surv(acttime, censor) ~  prevgenx + lethal + deathrt1 +
           acutediz + hosp01  + pspline(hospdisc, df = 4) +
           pspline(hhosleng, df = 4) + mandiz01 + femdiz01 + peddiz01 +
           orphdum + natreg + vandavg3 + wpnoavg3 +
           pspline(condavg3, df = 4) + pspline(orderent, df = 4) +
           pspline(stafcder, df = 4), data = CarpenterFdaData)

# Simulate Relative Hazards for orderent
Sim1 <- coxsimSpline(M1, bspline = "pspline(stafcder, df = 4)",
                    bdata = CarpenterFdaData$stafcder,
                    qi = "Hazard Ratio",
                    Xj = seq(1100, 1700, by = 10),
                    Xl = seq(1099, 1699, by = 10), spin = TRUE)

# Simulate Hazard Rates for orderent
Sim2 <- coxsimSpline(M1, bspline = "pspline(orderent, df = 4)",
                       bdata = CarpenterFdaData$orderent,
                       qi = "Hazard Rate",
                       Xj = seq(2, 53, by = 3), nsim = 100)

Run the code above in your browser using DataCamp Workspace