Description Usage Arguments Details Value See Also Examples
This function estimates a linear doseresponse 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 gcomputation with quantized exposures. Note: this function allows linear and nonadditive effects of individual components of the exposure, as well as nonlinear joint effects of the mixture via polynomial basis functions, which increase the computational computational burden due to the need for nonparametric bootstrapping.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 
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 untransformed version of exposures in the input datasets (useful if data are already transformed, or for performing standard gcomputation) 
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.noboot will not produce clusterappropriate standard errors. Qgcomp.boot can be used for this, which will use bootstrap sampling of clusters/individuals to estimate clusterappropriate standard errors via bootstrapping. 
weights 
"case weights"  passed to the "weight" argument of

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 runtime). 
rr 
logical: if using binary outcome and rr=TRUE, qgcomp.boot will estimate risk ratio rather than odds ratio 
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 ( 
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.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 glm (e.g. family) 
Estimates correspond to the average expected change in the (log) outcome per quantile increase in the joint exposure to all exposures in ‘expnms’. Test statistics and confidence intervals are based on a nonparametric bootstrap, using the standard deviation of the bootstrap estimates to estimate the standard error. The bootstrap standard error is then used to estimate Waldtype confidence intervals. Note that no bootstrapping is done on estimated quantiles of exposure, so these are treated as fixed quantities
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.
qgcomp.noboot
, and qgcomp
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108  set.seed(30)
# continuous outcome
dat < data.frame(y=rnorm(100), x1=runif(100), x2=runif(100), z=runif(100))
# Conditional linear slope
qgcomp.noboot(y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=4, family=gaussian())
# Marginal linear slope (population average slope, for a purely linear,
# additive model this will equal the conditional)
## Not run:
qgcomp.boot(f=y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=4,
family=gaussian(), B=200) # B should be at least 200 in actual examples
# no intercept model
qgcomp.boot(f=y ~ 1+z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=4,
family=gaussian(), B=200) # B should be at least 200 in actual examples
# Note that these give different answers! In the first, the estimate is conditional on Z,
# but in the second, Z is marginalized over via standardization. The estimates
# can be made approximately the same by centering Z (for linear models), but
# the conditional estimate will typically have lower standard errors.
dat$z = dat$z  mean(dat$z)
# Conditional linear slope
qgcomp.noboot(y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=4, family=gaussian())
# Marginal linear slope (population average slope, for a purely linear,
# additive model this will equal the conditional)
qgcomp.boot(f=y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=4,
family=gaussian(), B=200) # B should be at least 200 in actual examples
# Population average mixture slope which accounts for nonlinearity and interactions
qgcomp.boot(y ~ z + x1 + x2 + I(x1^2) + I(x2*x1), family="gaussian",
expnms = c('x1', 'x2'), data=dat, q=4, B=200)
# generally nonlinear/nonaddiive underlying models lead to nonlinear mixture slopes
qgcomp.boot(y ~ z + x1 + x2 + I(x1^2) + I(x2*x1), family="gaussian",
expnms = c('x1', 'x2'), data=dat, q=4, B=200, deg=2)
# binary outcome
dat < data.frame(y=rbinom(50,1,0.5), x1=runif(50), x2=runif(50), z=runif(50))
# Conditional mixture OR
qgcomp.noboot(y ~ z + x1 + x2, family="binomial", expnms = c('x1', 'x2'),
data=dat, q=2)
#Marginal mixture OR (population average OR  in general, this will not equal the
# conditional mixture OR due to noncollapsibility of the OR)
qgcomp.boot(y ~ z + x1 + x2, family="binomial", expnms = c('x1', 'x2'),
data=dat, q=2, B=3, rr=FALSE)
# Population average mixture RR
qgcomp.boot(y ~ z + x1 + x2, family="binomial", expnms = c('x1', 'x2'),
data=dat, q=2, rr=TRUE, B=3)
# Population average mixture RR, indicator variable representation of x2
# note that I(x==...) operates on the quantilebased category of x,
# rather than the raw value
res = qgcomp.boot(y ~ z + x1 + I(x2==1) + I(x2==2) + I(x2==3),
family="binomial", expnms = c('x1', 'x2'), data=dat, q=4, rr=TRUE, B=200)
res$fit
plot(res)
# now add in a nonlinear MSM
res2 = qgcomp.boot(y ~ z + x1 + I(x2==1) + I(x2==2) + I(x2==3),
family="binomial", expnms = c('x1', 'x2'), data=dat, q=4, rr=TRUE, B=200,
degree=2)
res2$fit
res2$msmfit # correct point estimates, incorrect standard errors
res2 # correct point estimates, correct standard errors
plot(res2)
# Log risk ratio per one IQR change in all exposures (not on quantile basis)
dat$x1iqr < dat$x1/with(dat, diff(quantile(x1, c(.25, .75))))
dat$x2iqr < dat$x2/with(dat, diff(quantile(x2, c(.25, .75))))
# note that I(x>...) now operates on the untransformed value of x,
# rather than the quantized value
res2 = qgcomp.boot(y ~ z + x1iqr + I(x2iqr>0.1) + I(x2>0.4) + I(x2>0.9),
family="binomial", expnms = c('x1iqr', 'x2iqr'), data=dat, q=NULL, rr=TRUE, B=200,
degree=2)
res2
# using parallel processing
qgcomp.boot(y ~ z + x1iqr + I(x2iqr>0.1) + I(x2>0.4) + I(x2>0.9),
family="binomial", expnms = c('x1iqr', 'x2iqr'), data=dat, q=NULL, rr=TRUE, B=200,
degree=2, parallel=TRUE, parplan=TRUE)
# weighted model
N=5000
dat4 < data.frame(id=seq_len(N), x1=runif(N), x2=runif(N), z=runif(N))
dat4$y < with(dat4, rnorm(N, x1*z + z, 1))
dat4$w=runif(N) + dat4$z*5
qdata = quantize(dat4, expnms = c("x1", "x2"), q=4)$data
# first equivalent models with no covariates
qgcomp.noboot(f=y ~ x1 + x2, expnms = c('x1', 'x2'), data=dat4, q=4, family=gaussian())
qgcomp.noboot(f=y ~ x1 + x2, expnms = c('x1', 'x2'), data=dat4, q=4, family=gaussian(),
weights=w)
set.seed(13)
qgcomp.boot(f=y ~ x1 + x2, expnms = c('x1', 'x2'), data=dat4, q=4, family=gaussian(),
weights=w)
# using the correct model
set.seed(13)
qgcomp.boot(f=y ~ x1*z + x2, expnms = c('x1', 'x2'), data=dat4, q=4, family=gaussian(),
weights=w, id="id")
(qgcfit < qgcomp.boot(f=y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat4, q=4,
family=gaussian(), weights=w))
qgcfit$fit
summary(glm(y ~ z + x1 + x2, data = qdata, weights=w))
## End(Not run)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.