Introduction

qgcomp is a package to implement g-computation for analyzing the effects of exposure mixtures. Quantile g-computation yields estimates of the effect of increasing all exposures by one quantile, simultaneously. This, it estimates a “mixture effect” useful in the study of exposure mixtures such as air pollution, diet, and water contamination.

The implementation in qgcomp is based on a generalization of weighted quantile sums (WQS) regression, which estimates the expected change in an outcome under a hypothetical intervention to increase all exposures in the mixture by one quantile. In the case in which all exposures have linear, additive effects on the quantile scale (i.e. after transforming all exposures into categorical variables defined by quantiles) qgcomp and WQS regression are asymptotically equivalent when all variables have effects in the same direction (the 'directional homogeneity' assumption of WQS). In moderate samples when underlying assumptions are not met exactly, WQS will tend to yield more biased effect estimates due to the assumption of directional homogeneity, and may also have lower precision due to the use of sample splitting. In cases in which the assumptions underlying WQS are met,qgcomp provides valid estimates of WQS 'weights', which estimate the relative contribution of each exposure to the overall mixture effect. Thus, qgcomp will provide valid effect estimates of the entire exposure mixture in general cases, while also allowing deviations from linearity and additivity assumptions.

Using terminology from methods developed for causal effect estimation, quantile g-computation estimates the parameters of a marginal structural model that characterizes the change in the expected potential outcome given a joint intervention on all exposures, possibly conditional on confounders. Under the assumptions of exchangeability, causal consistency, positivity, no interference, and correct model specification, this model yields a causal effect for an intervention on the mixture as a whole. While these assumptions may not be met exactly, they provide a useful roadmap for how to interpret the results of a qgcomp fit, and where efforts should be spent in terms of ensuring accurate model specification and selection of exposures that are sufficient to control co-pollutant confounding.

How to use the qgcomp package

Here we use a running example from the metals dataset from the from the package qgcomp to demonstrate some features of the package and method.

Namely, the examples below demonstrate use of the package for: 1. Fast estimation of exposure effects under a linear model for quantized exposures for continuous (normal) outcomes 2. Estimating conditional and marginal odds/risk ratios of a mixture effect for binary outcomes 3. Adjusting for non-exposure covariates when estimating effects of the mixture 4. Allowing non-linear and non-homogenous effects of individual exposures and the mixture as a whole by including product terms 5. Using qgcomp to fit a time-to-event model to estimate conditional and marginal hazard ratios for the exposure mixture (coming soon)

For analogous approaches to estimating exposure mixture effects, illustrative examples can be seen in the gQWS package help files, which implements WQS regression, and at https://jenfb.github.io/bkmr/overview.html, which describes Bayesian kernel machine regression.

The qgcomp package comes with some example data, which comprise a set of simulated well water exposures and two health outcomes (one continuous, one binary). The exposures are transformed to have mean = 0.0, standard deviation = 1.0. The data are used throughout to demonstrate usage and features of the qgcomp package.

library(qgcomp)
library(knitr)
data("metals", package="qgcomp")
head(metals)
##       arsenic     barium     cadmium    calcium   chloride   chromium
## 1 -0.14247764  0.2382489  0.06918218  0.8723956 -0.2254397  0.2001649
## 2 -0.22391247 -0.9986707 -0.20965619 -0.2722012 -0.1266698  0.2278990
## 3 -0.02313404 -0.5244978 -0.24330700  1.8023805 -0.1022276  0.4675177
## 4 -0.62555286 -0.3314166 -0.02263188 -0.6656564  0.4876452 -0.1042095
## 5 -0.15937578 -0.1211859 -0.18205546 -0.8087310 -0.2581438  0.9556171
## 6 -0.56329940  0.1975043  0.05797721 -1.0910304 -0.1674069  0.6789692
##        copper       iron         lead  magnesium   manganese    mercury
## 1  0.10449837 -0.1497640  0.058250135 -0.7107938 -0.05500511 -0.7321676
## 2 -0.64702980 -0.2051216 -0.114339703  0.3620378 -0.33681125  0.8230223
## 3  0.48890800  0.7763546 -0.227492111 -0.9253601  0.31202800 -1.5662318
## 4 -0.06923684 -0.2121876 -0.135158599  0.1474715 -0.39678244  0.1348865
## 5 -0.81346265 -0.2205297  0.333810130  0.1474715 -0.32391298  1.6935092
## 6  0.42913190  1.1449590 -0.001203963 -1.0151688 -0.44985876  0.9643995
##       nitrate    nitrite         ph   selenium     silver     sodium
## 1 -0.27148155 -0.4705337  0.3457076 -0.1571526  0.3421772 -0.3976339
## 2 -0.44566685 -0.8811492  0.6266144 -0.3304766  1.5884644  0.2721376
## 3 -0.30969203 -1.4922476 -0.4970128  0.6337753  1.6978706 -0.3560701
## 4  0.04515444  1.4350805  1.0479746 -0.9619823 -0.9106139  1.2102928
## 5 -0.36964859  0.3343747  0.2052542 -0.8430120  0.2845793  0.6521498
## 6 -0.02550827  1.0223160 -2.6038139 -0.9206695 -0.5999680 -0.4011965
##       sulfate total_alkalinity total_hardness       zinc mage35
## 1 -0.11699667        0.1834647      0.6830483 -0.1811071      0
## 2 -0.14129067        0.7665308     -0.1483585  0.5000657      0
## 3 -0.17027646        0.8248374      1.4748644 -0.1441533      0
## 4  0.05568995        2.0142922     -0.5838574 -0.1604838      0
## 5 -0.11751844        1.1863384     -0.7422206 -0.2281649      0
## 6 -0.17422680       -1.5890562     -1.2205620  0.8078975      0
##             y disease_state
## 1  0.57233831             0
## 2 -1.12962925             0
## 3  0.62308138             0
## 4  0.11517818             0
## 5  0.04546277             0
## 6  0.43851586             0

Example 1: linear model

# we save the names of the mixture variables in the variable "Xnm"
Xnm <- c(
    'arsenic','barium','cadmium','calcium','chromium','copper',
    'iron','lead','magnesium','manganese','mercury','selenium','silver',
    'sodium','zinc'
)
covars = c('nitrate','nitrite','sulfate','ph', 'total_alkalinity','total_hardness')



# Example 1: linear model
# Run the model and save the results "qc.fit"
system.time(qc.fit <- qgcomp.noboot(y~.,dat=metals[,c(Xnm, 'y')], family=gaussian()))
## Including all model terms as exposures of interest
##    user  system elapsed 
##   0.020   0.002   0.028
#   user  system elapsed 
#  0.011   0.002   0.018 

# contrasting other methods with computational speed
# WQS regression
#system.time(wqs.fit <- gwqs(y~NULL,mix_name=Xnm, data=metals[,c(Xnm, 'y')], family="gaussian"))
#   user  system elapsed 
# 35.775   0.124  36.114 

# Bayesian kernel machine regression (note that the number of iterations here would 
#  need to be >5,000, at minimum, so this underestimates the run time by a factor
#  of 50+
#system.time(bkmr.fit <- kmbayes(y=metals$y, Z=metals[,Xnm], family="gaussian", iter=100))
#   user  system elapsed 
# 81.644   4.194  86.520 


#first note that qgcomp is very fast

# View results: scaled coefficients/weights and statistical inference about
# mixture effect
qc.fit
## Scaled effect size (positive direction, sum of positive coefficients = 0.425)
##   calcium    barium   mercury manganese   arsenic  chromium      iron 
##  0.727104  0.076455  0.067479  0.046667  0.030523  0.018701  0.016592 
##   cadmium      lead    copper 
##  0.013656  0.002462  0.000362 
## 
## Scaled effect size (negative direction, sum of negative coefficients = -0.167)
## magnesium      zinc    silver    sodium  selenium 
##    0.3906    0.2473    0.1801    0.1489    0.0331 
## 
## Mixture slope parameters (Delta method CI):
## 
##      Estimate Std. Error Lower CI Upper CI t value  Pr(>|t|)
## psi1 0.257774   0.073702  0.11332  0.40223  3.4975 0.0005178

One advantage of quantile g-computation over other methods that estimate “mixture effects” (the effect of changing all exposures at once), is that it is very computationally efficient. Contrasting methods such as WQS (gWQS package) and Bayesian Kernel Machine regression (bkmr package), quantile g-computation can provide results many orders of magnitude faster. For example, the example above ran 3000X faster for quantile g-computation versus WQS regression, and we estimate the speedup would be several hundred thousand times versus Bayesian kernel machine regression.

Quantile g-computation yields fixed weights in the estimation procedure, similar to WQS regression. However, note that the weights from qgcomp.noboot can be negative or positive. When all effects are linear and in the same direction (“directional homogeneity”), quantile g-computation is equivalent to weighted quantile sum regression in large samples.

The overall mixture effect from quantile g-computation (psi1) is interpreted as the effect on the outcome of increasing every exposure by one quantile, possibly conditional on covariates. Given the overall exposure effect, the weights are considered fixed and so do not have confidence intervals or p-values.

Example 2: conditional odds ratio, marginal odds ratio in a logistic model

This example introduces the use of a binary outcome in qgcomp via the qgcomp.noboot function, which yields a conditional odds ratio or the qgcomp.boot, which yields a marginal odds ratio or risk/prevalence ratio. These will not equal each other when there are non-exposure covariates (e.g. confounders) included in the model because the odds ratio is not collapsible (both are still valid). Marginal parameters will yield estimates of the population average exposure effect, which is often of more interest due to better interpretability over conditional odds ratios. Further, odds ratios are not generally of interest when risk ratios can be validly estimated, so qgcomp.boot will estimate the risk ratio by default for binary data (set rr=FALSE to allow estimation of ORs when using qgcomp.boot).

qc.fit2 <- qgcomp.noboot(disease_state~., expnms=Xnm, 
          data = metals[,c(Xnm, 'disease_state')], family=binomial(), 
          q=4)
qcboot.fit2 <- qgcomp.boot(disease_state~., expnms=Xnm, 
          data = metals[,c(Xnm, 'disease_state')], family=binomial(), 
          q=4, B=10,# B should be 200-500+ in practice
          seed=125, rr=FALSE)
qcboot.fit2b <- qgcomp.boot(disease_state~., expnms=Xnm, 
          data = metals[,c(Xnm, 'disease_state')], family=binomial(), 
          q=4, B=10,# B should be 200-500+ in practice
          seed=125, rr=TRUE)


# Compare a qgcomp.noboot fit:
qc.fit2
## Scaled effect size (positive direction, sum of positive coefficients = 1.69)
##  mercury  arsenic  calcium     zinc   silver   copper  cadmium selenium 
##   0.3049   0.1845   0.1551   0.1153   0.0828   0.0730   0.0645   0.0199 
## 
## Scaled effect size (negative direction, sum of negative coefficients = -0.733)
##    barium      lead  chromium      iron manganese magnesium    sodium 
##    0.2679    0.2571    0.2304    0.1267    0.0542    0.0481    0.0156 
## 
## Mixture log(OR) (Delta method CI):
## 
##      Estimate Std. Error Lower CI Upper CI Z value Pr(>|z|)
## psi1  0.95579    0.46656 0.041347   1.8702  2.0486   0.0405
# and a qgcomp.boot fit:
qcboot.fit2
## Mixture log(OR) (bootstrap CI):
## 
##      Estimate Std. Error Lower CI Upper CI Z value Pr(>|z|)
## psi1  0.95579    0.44328 0.086982   1.8246  2.1562  0.03107
# and a qgcomp.boot fit, where the risk/prevalence ratio is estimated, 
#  rather than the odds ratio:
qcboot.fit2b
## Mixture log(RR) (bootstrap CI):
## 
##      Estimate Std. Error Lower CI Upper CI Z value Pr(>|z|)
## psi1  0.77156    0.32901  0.12673   1.4164  2.3451  0.01902

Example 3: adjusting for covariates, plotting estimates

In the following code we run a maternal age-adjusted linear model with qgcomp (family = "gaussian"). Further, we plot

qc.fit3 <- qgcomp.noboot(y ~ mage35 + arsenic + barium + cadmium + calcium + chloride + 
                           chromium + copper + iron + lead + magnesium + manganese + 
                           mercury + selenium + silver + sodium + zinc,
                         expnms=Xnm,
                         metals, family=gaussian(), q=4)
qc.fit3
## Scaled effect size (positive direction, sum of positive coefficients = 0.434)
##   calcium    barium   mercury manganese   arsenic  chromium   cadmium 
##   0.71389   0.08325   0.06767   0.05064   0.03706   0.02049   0.01150 
##      iron      lead 
##   0.01052   0.00497 
## 
## Scaled effect size (negative direction, sum of negative coefficients = -0.178)
## magnesium      zinc    silver    sodium  selenium    copper 
##    0.3726    0.2325    0.1820    0.1636    0.0312    0.0182 
## 
## Mixture slope parameters (Delta method CI):
## 
##      Estimate Std. Error Lower CI Upper CI t value Pr(>|t|)
## psi1 0.256057   0.073911   0.1112  0.40092  3.4644 0.000584
plot(qc.fit3)
qcboot.fit3 <- qgcomp.boot(y ~ mage35 + arsenic + barium + cadmium + calcium + chloride + 
                           chromium + copper + iron + lead + magnesium + manganese + 
                           mercury + selenium + silver + sodium + zinc,
                         expnms=Xnm,
                         metals, family=gaussian(), q=4, B=10,# B should be 200-500+ in practice
                         seed=125)
qc.fit3
## Scaled effect size (positive direction, sum of positive coefficients = 0.434)
##   calcium    barium   mercury manganese   arsenic  chromium   cadmium 
##   0.71389   0.08325   0.06767   0.05064   0.03706   0.02049   0.01150 
##      iron      lead 
##   0.01052   0.00497 
## 
## Scaled effect size (negative direction, sum of negative coefficients = -0.178)
## magnesium      zinc    silver    sodium  selenium    copper 
##    0.3726    0.2325    0.1820    0.1636    0.0312    0.0182 
## 
## Mixture slope parameters (Delta method CI):
## 
##      Estimate Std. Error Lower CI Upper CI t value Pr(>|t|)
## psi1 0.256057   0.073911   0.1112  0.40092  3.4644 0.000584
plot(qcboot.fit3)

plot of chunk unnamed-chunk-4plot of chunk unnamed-chunk-4

From the first plot we see weights from qgcomp.noboot function, which include both positive and negative effect directions. When the weights are all on a single side of the null, these plots are easy to in interpret since the weight corresponds to the proportion of the overall effect from each exposure. WQS uses a constraint in the model to force all of the weights to be in the same direction - unfortunately such constraints lead to biased effect estimates. The qgcomp package takes a different approach and allows that “weights” might go in either direction, indicating that some exposures may beneficial, and some harmful, or there may be sampling variation due to using small or moderate sample sizes (or, more often, systematic bias such as unmeasured confounding). The “weights” in qgcomp correspond to the proportion of the overall effect when all of the exposures have effects in the same direction, but otherwise they correspond to the proportion of the effect in a particular direction, which may be small (or large) compared to the overall “mixture” effect.NOTE: that the left and right sides of the plot should not be compared with each other because the length of the bars corresponds to the effect size only relative to other effects in the same direction. The darkness of the bars corresponds to the overall effect size - in this case the bars on the right (positive) side of the plot are darker because the overall “mixture” effect is positive. Thus, the shading allows one to make informal comparisons across the left and right sides: a large, darkly shaded bar indicates a larger independent effect than a large, lightly shaded bar.

Using qgcomp.boot also allows us to assess linearity of the total exposure effect (the second plot). Similar output is available for WQS (gWQS package), though WQS results will generally be less interpretable when exposure effects are non-linear (see below how to do this with qgcomp.boot.

The plot for the qcboot.fit3 object (using g-computation with bootstrap variance) gives predictions at the joint intervention levels of exposure. It also displays a smoothed (graphical) fit. Generally, we cannot overlay the data over this plot since the regression line corresponds to a change in potentially many exposures at once. Hence, it is useful to explore non-linearity by fitting models that allow for non-linear effects.

Example 4: non-linearity (and non-homogeneity)

Let's close with one more feature of qgcomp (and qgcomp.boot): handling non-linearity. Here is an example where we use a feature of the R language for fitting models with interaction terms. We use y~. + .^2 as the model formula, which fits a model that allows for quadratic term for every predictor in the model.

Similar approaches could be used to include interaction terms between exposures, as well as between exposures and covariates.

qcboot.fit4 <- qgcomp(y~. + .^2,
                         expnms=Xnm,
                         metals[,c(Xnm, 'y')], family=gaussian(), q=4, B=10, seed=125)
plot(qcboot.fit4)

plot of chunk unnamed-chunk-5

Note that allowing for a non-linear effect of all exposures induces an apparent non-linear trend in the overall exposure effect. The smoothed regression line is still well within the confidence bands of the marginal linear model (by default, the overall effect of joint exposure is assumed linear, though this assumption can be relaxed via the 'degree' parameter in qgcomp.boot, as follows:

qcboot.fit5 <- qgcomp(y~. + .^2,
                         expnms=Xnm,
                         metals[,c(Xnm, 'y')], family=gaussian(), q=4, degree=2, B=10, seed=125)
plot(qcboot.fit5)

plot of chunk unnamed-chunk-6

Ideally, the smooth fit will look very similar to the model prediction regression line.

Example 5: Comparing model fits and further exploring non-linearity

Exploring a non-linear fit in settings with multiple exposures is challenging. One way to explore non-linearity, as demonstrated above, is to to include all 2-way interaction terms (including quadratic terms, or “self-interactions”). Sometimes this approach is not desired, either because the number of terms in the model can become very large, or because some sort of model selection procedure is required, which risks inducing over-fit (biased estimates and standard errors that are too small). Short of having a set of a priori non-linear terms to include, we find it best to take a default approach (e.g. taking all second order terms) that doesn't rely on statistical significance, or to simply be honest that the search for a non-linear model is exploratory and shouldn't be relied upon for robust inference. Methods such as kernel machine regression may be good alternatives, or supplementary approaches to exploring non-linearity.

NOTE: qgcomp necessarily fits a regression model with exposures that have a small number of possible values, based on the quantile chosen. By package default, this is q=4, but it is difficult to fully examine non-linear fits using only four points, so we recommend exploring larger values of q, which will change effect estimates (i.e. the model coefficient implies a smaller change in exposures, so the expected change in the outcome will also decrease).

Here, we examine a couple one strategy for default and exploratory approaches to mixtures that can be implemented in qgcomp using a smaller subset of exposures (iron, lead, cadmium), which we choose via the correlation matrix. High correlations between exposures may result from a common source, so small subsets of the mixture may be useful for examining hypotheses that relate to interventions on a common environmental source or set of behaviors. Note that we can still adjust for the measured exposures, even though only 3 our exposures of interest are considered as the mixture of interest. Note that we will require a new R package to help in exploring non-linearity: splines. Note that qgcomp.boot must be used in order to produce the graphics below, as qgcomp.boot does not calculate the necessary quantities.

Graphical approach to explore non-linearity in a correlated subset of exposures using splines

library(splines)
# find all correlations > 0.6 (this is an arbitrary choice)
cormat = cor(metals[,Xnm])
idx = which(cormat>0.6 & cormat <1.0, arr.ind = TRUE)
newXnm = unique(rownames(idx)) # iron, lead, and cadmium


qc.fit6lin <- qgcomp.boot(y ~ iron + lead + cadmium + 
                         mage35 + arsenic + magnesium + manganese + mercury + 
                         selenium + silver + sodium + zinc,
                         expnms=newXnm,
                         metals, family=gaussian(), q=8, B=10)

qc.fit6nonlin <- qgcomp.boot(y ~ bs(iron) + bs(cadmium) + bs(lead) +
                         mage35 + arsenic + magnesium + manganese + mercury + 
                         selenium + silver + sodium + zinc,
                         expnms=newXnm,
                         metals, family=gaussian(), q=8, B=10, degree=2)

qc.fit6nonhom <- qgcomp.boot(y ~ bs(iron)*bs(lead) + bs(iron)*bs(cadmium) + bs(lead)*bs(cadmium) +
                         mage35 + arsenic + magnesium + manganese + mercury + 
                         selenium + silver + sodium + zinc,
                         expnms=newXnm,
                         metals, family=gaussian(), q=8, B=10, degree=3)
# it helps to place the plots on a common y-axis, which is easy due to dependence of the qgcomp plotting functions on ggplot
pl.fit6lin <- plot(qc.fit6lin, suppressprint = TRUE)
pl.fit6nonlin <- plot(qc.fit6nonlin, suppressprint = TRUE)
pl.fit6nonhom <- plot(qc.fit6nonhom, suppressprint = TRUE)

pl.fit6lin + coord_cartesian(ylim=c(-0.75, .75)) + ggtitle("Linear fit: mixture of iron, lead, and cadmium")
pl.fit6nonlin + coord_cartesian(ylim=c(-0.75, .75)) + ggtitle("Linear fit: mixture of iron, lead, and cadmium")
pl.fit6nonhom + coord_cartesian(ylim=c(-0.75, .75)) + ggtitle("Non-linear, non-homogeneous fit: mixture of iron, lead, and cadmium")

plot of chunk unnamed-chunk-7plot of chunk unnamed-chunk-7plot of chunk unnamed-chunk-7

Note one restriction on exploring non-linearity: while we can use flexible function such as splines for individual exposures, the overall fit is limited via the degree parameter to polynomial functions (here a quadratic polynomial fits the non-linear model well, and a cubic polynomial fits the non-linear/non-homogenous model well - though this is an informal argument and does not account for the wide confidence intervals). We note here that only 10 bootstrap iterations are used to calculate confidence intervals, which is far too low.

Statistical approach explore non-linearity in a correlated subset of exposures using splines

The graphical approaches don't give a clear picture of which model might be preferred, but we can compare the model fits using AIC, AICC, or BIC (all various forms of information criterion that weigh model fit with over-parameterization). Both of these criterion suggest the linear model fits best (lowest AIC and BIC), which suggests that the apparently non-linear fits observed in the graphical approaches don't improve prediction of the health outcome, relative to the linear fit, due to the increase in variance associated with including more parameters.

AIC(qc.fit6lin$fit)
## [1] 755.6437
AIC(qc.fit6nonlin$fit)
## [1] 760.5661
AIC(qc.fit6nonhom$fit)
## [1] 772.6342
BIC(qc.fit6lin$fit)
## [1] 813.2352
BIC(qc.fit6nonlin$fit)
## [1] 842.8397
BIC(qc.fit6nonhom$fit)
## [1] 965.9773

Example 6: time-to-event analysis

Example coming soon. See help file for qgcomp.cox.noboot for usage.

References

Alexander P. Keil, Jessie P. Buckley, Katie M. O’Brien, Kelly K. Ferguson, Shanshan Zhao, Alexandra J. White. A quantile-based g-computation approach to addressing the effects of exposure mixtures. https://arxiv.org/abs/1902.04200

Acknowledgements

The development of this package was supported by NIH Grant RO1ES02953101