coxsimInteract() R function from [simPH]

Simulate quantities of interest for linear multiplicative interactions from Cox Proportional Hazards models

coxsimInteract simulates quantities of interest for linear multiplicative interactions using multivariate normal distributions. These can be plotted with simGG.


coxsimInteract(
  obj,
  b1,
  b2,
  qi = "Marginal Effect",
  X1 = NULL,
  X2 = NULL,
  means = FALSE,
  expMarg = TRUE,
  nsim = 1000,
  ci = 0.95,
  spin = FALSE,
  extremesDrop = TRUE
)

Arguments

obj: a coxph class fitted model object with a linear multiplicative interaction.
b1: character string of the first constitutive variable's name. Note b1 and b2 must be entered in the order in which they are entered into the coxph model.
b2: character string of the second constitutive variable's name.
qi: quantity of interest to simulate. Values can be "Marginal Effect", "First Difference", "Hazard Ratio", and "Hazard Rate". The default is qi = "Hazard Ratio". If qi = "Hazard Rate" and the coxph model has strata, then hazard rates for each strata will also be calculated.
X1: numeric vector of fitted values of b1 to simulate for. If qi = "Marginal Effect" then only X2 can be set. If you want to plot the results, X1 should have more than one value.
X2: numeric vector of fitted values of b2 to simulate for.
means: logical, whether or not to use the mean values to fit the hazard rate for covaraiates other than b1 b2 and b1*b2. Note: it does not currently support models that include polynomials created by I. Note: EXPERIMENTAL. lines are not currently supported in simGG if means = TRUE.
expMarg: logical. Whether or not to exponentiate the marginal effect.
nsim: the number of simulations to run per value of X. Default is nsim = 1000.
ci: the proportion of middle 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 siminteract class object

Details

Simulates marginal effects, first differences, hazard ratios, and hazard rates for linear multiplicative interactions. Marginal effects are calculated as in Brambor et al. (2006) with the addition that we take the exponent, so that it resembles a hazard ratio. You can choose not to take the exponent by setting the argument expMarg = FALSE. For an interaction between variables $X$ and $Z$ the marginal effect for $X$ is:

ME_{X} = e^(\beta_{X} + \beta_{XZ}Z)ME[X] = exp(\beta[X] +\beta[XZ]Z)

Note that for First Differences the comparison is not between two values of the same variable but two values of the constitute variable and 0 for the two variables.

Examples


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

# Load survival package
library(survival)

# Run basic model
M1 <- coxph(Surv(acttime, censor) ~ lethal*prevgenx,
            data = CarpenterFdaData)

# Simulate Marginal Effect of lethal for multiple
# values of prevgenx
Sim1 <- coxsimInteract(M1, b1 = "lethal", b2 = "prevgenx",
                       X2 = seq(2, 115, by = 5), spin = TRUE)

## Not run:

# Change the order of the covariates to make a more easily
# interpretable relative hazard graph.
M2 <- coxph(Surv(acttime, censor) ~ prevgenx*lethal +
             orphdum, data = CarpenterFdaData)

# Simulate Hazard Ratio of lethal for multiple
# values of prevgenx
Sim2 <- coxsimInteract(M2, b1 = "prevgenx", b2 = "lethal",
                    X1 = seq(2, 115, by = 2),
                    X2 = c(0, 1),
                    qi = "Hazard Ratio", ci = 0.9)

# Simulate First Difference
Sim3 <- coxsimInteract(M2, b1 = "prevgenx", b2 = "lethal",
                       X1 = seq(2, 115, by = 2),
                       X2 = c(0, 1),
                       qi = "First Difference", spin = TRUE)

# Simulate Hazard Rate
Sim4 <- coxsimInteract(M2, b1 = "prevgenx", b2 = "lethal",
                       X1 = 90, X2 = 1, qi = "Hazard Rate",
                       means = TRUE)
## End(Not run)

# Example with a categorical variable
# Download data
data(hmohiv)

# Create category lables
hmohiv$drug <- factor(hmohiv$drug, labels = c('not treated', 'treated'))

M3 <- coxph(Surv(time,censor) ~ drug*age, data = hmohiv)

# Note: Use relevant coefficient name as shown in model summary, e.g.
# 'drugtreated'.
Sim5 <- coxsimInteract(M3, b1 = "drugtreated", b2 = 'age', X2 = 20:54)

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.

Brambor, Thomas, William Roberts Clark, and Matt Golder. 2006. ''Understanding Interaction Models: Improving Empirical Analyses.'' Political Analysis 14(1): 63-82.

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.

coxsimInteract function