coxsimSpline() R function from [simPH]

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

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


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.

Returns

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)}h(t|X[i])=h[0](t)exp(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+}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(hi) = (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.

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)

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, tools:::Rd_expr_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.

coxsimSpline function