Learn R Programming

simPH (version 1.3.13)

coxsimtvc: Simulate time-interactive quantities of interest from Cox Proportional Hazards models

Description

coxsimtvc simulates time-interactive relative hazards, first differences, and hazard ratios from models estimated with coxph using the multivariate normal distribution. These can be plotted with simGG.

Usage

coxsimtvc(
  obj,
  b,
  btvc,
  qi = "Relative Hazard",
  Xj = NULL,
  Xl = NULL,
  tfun = "linear",
  pow = NULL,
  nsim = 1000,
  from,
  to,
  by = 1,
  ci = 0.95,
  spin = FALSE,
  extremesDrop = TRUE
)

Arguments

obj

a coxph fitted model object with a time interaction.

b

the non-time interacted variable's name.

btvc

the time interacted variable's name.

qi

character string indicating what quantity of interest you would like to calculate. Can be 'Relative Hazard', 'First Difference', 'Hazard Ratio', 'Hazard Rate'. Default is qi = 'Relative Hazard'. If qi = 'First Difference' or qi = 'Hazard Ratio' then you can set Xj and Xl.

Xj

numeric vector of fitted values for b. Must be the same length as Xl or Xl must be NULL.

Xl

numeric vector of fitted values for Xl. Must be the same length as Xj. Only applies if qi = 'First Difference' or qi = 'Hazard Ratio'.

tfun

function of time that btvc was multiplied by. Default is "linear". It can also be "log" (natural log) and "power". If tfun = "power" then the pow argument needs to be specified also.

pow

if tfun = "power", then use pow to specify what power the time interaction was raised to.

nsim

the number of simulations to run per point in time. Default is nsim = 1000.

from

point in time from when to begin simulating coefficient values

to

point in time to stop simulating coefficient values.

by

time intervals by which to simulate coefficient values.

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.

Value

a simtvc object

Details

Simulates time-varying relative hazards, first differences, and hazard ratios using parameter estimates from coxph models. Can also simulate hazard rates for multiple strata.

Relative hazards are found using: $$RH = e^{\beta_{1} + \beta_{2}f(t)}$$ where \(f(t)\) is the function of time.

First differences are found using: $$FD = (e^{(X_{j} - X_{l}) (\beta_{1} + \beta_{2}f(t))} - 1) * 100$$ where \(X_{j}\) and \(X_{l}\) are some values of \(X\) to contrast.

Hazard ratios are calculated using: $$FD = e^{(X_{j} - X_{l}) (\beta_{1} + \beta_{2}f(t))}$$ When simulating non-stratifed time-varying harzards coxsimtvc uses the point estimates for a given coefficient \(\hat{\beta}_{x}\) and its time interaction \(\hat{\beta}_{xt}\) along with the variance matrix (\(\hat{V}(\hat{\beta})\)) estimated from a coxph model. These are used to draw values of \(\beta_{1}\) and \(\beta_{2}\) from the multivariate normal distribution \(N(\hat{\beta},\: \hat{V}(\hat{\beta}))\).

When simulating stratified time-varying hazard rates \(H\) for a given strata \(k\), coxsimtvc uses: $$H_{kxt} = \hat{\beta}_{k0t}e^{\hat{\beta_{1}} + \beta_{2}f(t)}$$ The resulting simulation values can be plotted using simGG.

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.

Golub, Jonathan, and Bernard Steunenberg. 2007. ''How Time Affects EU Decision-Making.'' European Union Politics 8(4): 555-66.

Licht, Amanda A. 2011. ''Change Comes with Time: Substantive Interpretation of Nonproportional Hazards in Event History Analysis.'' Political Analysis 19: 227-43.

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

Run this code
# NOT RUN {
# Load Golub & Steunenberg (2007) Data
data("GolubEUPData")

# Load survival package
library(survival)

# Expand data (not run to speed processing time, but should be run)
GolubEUPData <- SurvExpand(GolubEUPData, GroupVar = 'caseno',
                     Time = 'begin', Time2 = 'end', event = 'event')

# Create time interactions
BaseVars <- c('qmv', 'backlog', 'coop', 'codec', 'qmvpostsea', 'thatcher')
GolubEUPData <- tvc(GolubEUPData, b = BaseVars, tvar = 'end', tfun = 'log')

# Run Cox PH Model
M1 <- coxph(Surv(begin, end, event) ~ qmv + qmvpostsea + qmvpostteu +
                coop + codec + eu9 + eu10 + eu12 + eu15 + thatcher +
                agenda + backlog + qmv_log + qmvpostsea_log + coop_log +
                codec_log + thatcher_log + backlog_log,
            data = GolubEUPData, ties = "efron")

# Create simtvc object for Relative Hazard
Sim1 <- coxsimtvc(obj = M1, b = "qmv", btvc = "qmv_log",
                   tfun = "log", from = 80, to = 2000,
                   Xj = 1, by = 15, ci = 0.99, nsim = 100)

# Create simtvc object for First Difference
Sim2 <- coxsimtvc(obj = M1, b = "qmv", btvc = "qmv_log",
                 qi = "First Difference", Xj = 1,
                 tfun = "log", from = 80, to = 2000,
                 by = 15, ci = 0.95)

# Create simtvc object for Hazard Ratio
Sim3 <- coxsimtvc(obj = M1, b = "backlog", btvc = "backlog_log",
                  qi = "Hazard Ratio", Xj = c(191, 229),
                  Xl = c(0, 0),
                  tfun = "log", from = 80, to = 2000,
                  by = 15, ci = 0.5)
# }
# NOT RUN {
# }

Run the code above in your browser using DataLab