coxsimSpline {simPH}

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.

bdata

a numeric vector of the splined variable's values.

qi

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 qi = "Hazard Rate". You may be creating very many simulated values which can be very computationally intensive to do. Adjust the number of simulations per fitted value with nsim.

Xj

numeric vector of fitted values for b to simulate for.

Xl

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.

ci

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 used.

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 simGG and may prevent spin from finding the central interval.

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.

Value

a simspline object

References

Gandrud, Christopher. 2015. simPH: An R Package for Illustrating Estimates from Cox Proportional Hazard Models Including for Interactive and Nonlinear Effects. Journal of Statistical Software. 65(3)1-20.

Luke Keele, "Replication data for: Proportionally Difficult: Testing for Nonproportional Hazards In Cox Models", 2010, doi: 10.7910/DVN/VJAHRG 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. https://arxiv.org/pdf/1302.2142v1.pdf.

See Also

simGG, survival, strata, and coxph

Examples

# 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)

## Not run: 
# 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)


## End(Not run)
# 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)