coxsimSpline: Simulate quantities of interest for penalized splines from...
In simPH: Simulate and Plot Estimates from Cox Proportional Hazards Models

Description Usage Arguments Details Value References See Also Examples

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
)

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

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|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) = β[k][1](x)[1+] + β[k][2](x)[2+] + β[k][3](x)[3+] + … + β[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 β[k][1], … β[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.

a simspline object

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.

simGG, survival, strata, and coxph

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