qgcomp.multinomial.boot: Quantile g-computation for multinomial outcomes

Description

This function estimates a dose-response parameter representing a one quantile increase in a set of exposures of interest. This model estimates the parameters of a marginal structural model (MSM) based on g-computation with quantized exposures. Note: this function allows linear and non-additive effects of individual components of the exposure, as well as non-linear joint effects of the mixture via polynomial basis functions, which increase the computational computational burden due to the need for non-parametric bootstrapping.

Usage

qgcomp.multinomial.boot(
  f,
  data,
  expnms = NULL,
  q = 4,
  breaks = NULL,
  id = NULL,
  weights,
  alpha = 0.05,
  B = 200,
  rr = TRUE,
  degree = 1,
  seed = NULL,
  bayes = FALSE,
  MCsize = nrow(data),
  parallel = FALSE,
  parplan = FALSE,
  ...
)

Value

a qgcompfit object, which contains information about the effect measure of interest (psi) and associated variance (var.psi), as well as information on the model fit (fit) and information on the marginal structural model (msmfit) used to estimate the final effect estimates.

Arguments

f: R style formula
data: data frame
expnms: character vector of exposures of interest
q: NULL or number of quantiles used to create quantile indicator variables representing the exposure variables. If NULL, then gcomp proceeds with un-transformed version of exposures in the input datasets (useful if data are already transformed, or for performing standard g-computation)
breaks: (optional) NULL, or a list of (equal length) numeric vectors that characterize the minimum value of each category for which to break up the variables named in expnms. This is an alternative to using 'q' to define cutpoints.
id: (optional) NULL, or variable name indexing individual units of observation (only needed if analyzing data with multiple observations per id/cluster). Note that qgcomp.multinomial.noboot will not produce cluster-appropriate standard errors. qgcomp.multinomial.boot can be used for this, which will use bootstrap sampling of clusters/individuals to estimate cluster-appropriate standard errors via bootstrapping.
weights: "case weights" - passed to the "weight" argument of multinom
alpha: alpha level for confidence limit calculation
B: integer: number of bootstrap iterations (this should typically be >=200, though it is set lower in examples to improve run-time).
rr: (not used)
degree: polynomial bases for marginal model (e.g. degree = 2 allows that the relationship between the whole exposure mixture and the outcome is quadratic (default = 1).
seed: integer or NULL: random number seed for replicable bootstrap results
bayes: use underlying Bayesian model (arm package defaults). Results in penalized parameter estimation that can help with very highly correlated exposures. Note: this does not lead to fully Bayesian inference in general, so results should be interpreted as frequentist.
MCsize: integer: sample size for simulation to approximate marginal zero inflated model parameters. This can be left small for testing, but should be as large as needed to reduce simulation error to an acceptable magnitude (can compare psi coefficients for linear fits with qgcomp.multinomial.noboot to gain some intuition for the level of expected simulation error at a given value of MCsize). This likely won't matter much in linear models, but may be important with binary or count outcomes.
parallel: use (safe) parallel processing from the future and future.apply packages
parplan: (logical, default=FALSE) automatically set future::plan to plan(multisession) (and set to existing plan, if any, after bootstrapping)
...: arguments to nnet::multinom

Details

Estimates correspond to the average expected change in the probability of an outcome type per quantile increase in the joint exposure to all exposures in `expnms'. Test statistics and confidence intervals are based on a non-parametric bootstrap, using the standard deviation of the bootstrap estimates to estimate the standard error. The bootstrap standard error is then used to estimate Wald-type confidence intervals. Note that no bootstrapping is done on estimated quantiles of exposure, so these are treated as fixed quantities

Examples

Run this code

data("metals") # from qgcomp package
# create categorical outcome from the existing continuous outcome (usually, one will already exist)
metals$ycat = factor(quantize(metals, "y",q=4)$data$y, levels=c("0", "1", "2", "3"), 
                     labels=c("cct", "ccg", "aat", "aag")) 
# restrict to smaller dataset for simplicity
smallmetals = metals[,c("ycat", "arsenic", "lead", "cadmium", "mage35")]

### 1: Define mixture and underlying model ####
mixture = c("arsenic", "lead", "cadmium")
f0 = ycat ~ arsenic + lead + cadmium # the multinomial model 
# (be sure that factor variables are properly coded ahead of time in the dataset)
rr = qgcomp.multinomial.boot(
 f0, 
 expnms = mixture,
 q=4, 
 data = smallmetals, 
 B = 5, # set to higher values in real examples
 MCsize = 100,  # set to higher values in small samples
 )

rr2 = qgcomp.multinomial.noboot(
 f0, 
 expnms = mixture,
 q=4, 
 data = smallmetals
 )
 
 ### 5: Create summary qgcomp object for nice printing ####
 
 summary(rr, tests=c("H")) # include homogeneity test
 
 # 95% confidence intervals
 #confint(rr, level=0.95)
 #rr$breaks # quantile cutpoints for exposures
 # homogeneity_test(rr)
 #joint_test(rr)

qdat = simdata_quantized(
  outcometype="multinomial", 
  n=10000, corr=c(-.9), coef=cbind(c(.2,-.2,0,0), c(.1,.1,.1,.1)), 
  q = 4
)

 rr_sim = qgcomp.multinomial.noboot(
  y~x1+x2+x3+x4, 
  expnms = c("x1", "x2", "x3", "x4"),
  q=4, 
  data = qdat
 )
 
 rr_sim2 = qgcomp.multinomial.boot(
  y~x1+x2+x3+x4, 
  expnms = c("x1", "x2", "x3", "x4"),
  q=4, 
  data = qdat,
  B=1
 )