coxsimPoly {simPH}

Simulate quantities of interest for a range of values for a polynomial nonlinear effect from Cox Proportional Hazards models

Description

coxsimPoly simulates quantities of interest for polynomial covariate effects estimated from Cox Proportional Hazards models. These can be plotted with simGG.

Usage

coxsimPoly(
  obj,
  b = NULL,
  qi = "Relative Hazard",
  pow = 2,
  Xj = NULL,
  Xl = NULL,
  nsim = 1000,
  ci = 0.95,
  spin = FALSE,
  extremesDrop = TRUE
)

Arguments

obj

a coxph class fitted model object with a polynomial coefficient. These can be plotted with simGG.

b

character string name of the coefficient you would like to simulate. To find the quantity of interest using only the polynomial and not the polynomial + the linear terms enter the polynomial created using I, e.g. I(natreg^2) as a string.

qi

quantity of interest to simulate. Values can be "Relative Hazard", "First Difference", "Hazard Ratio", and "Hazard Rate". The default is qi = "Relative Hazard". If qi = "Hazard Rate" and the coxph model has strata, then hazard rates for each strata will also be calculated.

pow

numeric polynomial used in coxph.

Xj

numeric vector of fitted values for b to simulate for.

Xl

numeric vector of values to compare Xj to. If NULL, then it is automatically set to 0.

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 quantities of interest for polynomial covariate effects. For example if a nonlinear effect is modeled with a second order polynomial–i.e. \beta_{1}x_{i} + \beta_{2}x_{i}^{2}–we can draw n simulations from the multivariate normal distribution for both \beta_{1} and \beta_{2}. Then we simply calculate quantities of interest for a range of values and plot the results as before. For example, we find the first difference for a second order polynomial with:

\%\triangle h_{i}(t) = (\mathrm{e}^{\beta_{1}x_{j-1} + \beta_{2}x_{j-l}^{2}} - 1) * 100

where x_{j-l} = x_{j} - x_{l}.

Note, you must use I to create the polynomials.

Value

a simpoly, coxsim 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.

Keele, Luke. 2010. ”Proportionally Difficult: Testing for Nonproportional Hazards in Cox Models.” Political Analysis 18(2): 189-205.

Carpenter, Daniel P. 2002. ”Groups, the Media, Agency Waiting Costs, and FDA Drug Approval.” American Journal of Political Science 46(3): 490-505.

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.simpoly, survival, strata, and coxph

Examples

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

# Load survival package
library(survival)

# Run basic model
M1 <- coxph(Surv(acttime, censor) ~ prevgenx + lethal + deathrt1 +
            acutediz + hosp01  + hhosleng + mandiz01 + femdiz01 +
            peddiz01 + orphdum + natreg + I(natreg^2) +
            I(natreg^3) + vandavg3 + wpnoavg3 +
            condavg3 + orderent + stafcder, data = CarpenterFdaData)

# Simulate simpoly First Difference
Sim1 <- coxsimPoly(M1, b = "natreg", qi = "First Difference",
                pow = 3, Xj = seq(1, 150, by = 5), nsim = 100)

## Not run: 
# Simulate simpoly Hazard Ratio with spin probibility interval
Sim2 <- coxsimPoly(M1, b = "natreg", qi = "Hazard Ratio",
              pow = 3, Xj = seq(1, 150, by = 5), spin = TRUE)

## End(Not run)