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 )
obj |
a |
bspline |
a character string of the full |
bdata |
a numeric vector of the splined variable's values. |
qi |
quantity of interest to simulate. Values can be
|
Xj |
numeric vector of fitted values for |
Xl |
numeric vector of values to compare |
nsim |
the number of simulations to run per value of |
ci |
the proportion of simulations to keep. The default is
|
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 |
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.
# 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)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.