Estimate propensity scores for multivariate continuous exposure by assuming joint normal conditional densities. Simultaneously estimate stabilized inverse probability of treatment weights (IPTW) using joint normal density for marginal distribution of exposure.
mvGPS(D, C, common = FALSE, trim_w = FALSE, trim_quantile = 0.99)numeric matrix of dimension \(n\) by \(m\) designating values of the exposures
either a list of numeric matrices of length \(m\) of dimension
\(n\) by \(p_j\) designating values of the confounders for each exposure
value or if common is TRUE a single matrix of of dimension \(n\) by
\(p\) that represents common confounders for all exposures.
logical indicator for whether C is a single matrix of common confounders for all exposures. default is FALSE meaning C must be specified as list of confounders of length \(m\).
logical indicator for whether to trim weights. default is FALSE
numeric scalar used to specify the upper quantile to trim weights if applicable. default is 0.99
list of score and wts, where score is the mvGPS score values and wts are the corresponding stabilized inverse probability of treatment weights
Generalized propensity scores (GPS) were proposed by hirano_continuous;textualmvGPS and imai_causalGPS;textualmvGPS to extend propensity scores to handle continuous exposures. The GPS is constructed using the conditional density of the exposure given a set of confounders. These original methods and the subsequent literature have focused on the case of a single continuous exposure where the GPS could be estimated using normal densities, kernel smoothing kennedy2017mvGPS, gradient boosting zhu_boostingmvGPS, and empirical likelihoods fong2018mvGPS. In this package we provide an extension to this literature to allow for multivariate continuous exposures to be estimated.
Assume that we have a set of continuous exposures, \(D\), of length
m, i.e., \(\mathbf{D}=D_{1}, \dots, D_{m}\) collected on \(n\)
units. Further, we assume that there exists a set of confounders
\(\mathbf{C}_{1},\dots,\mathbf{C}_{m}\) for each
exposure of length \(p_{j}\) for \(j=1,\dots,m\). The confounders are
related to both the exposures and the outcome of interest such. Note that
the confounders need not be identical for all exposures.
In our function we therefore expect that the argument D is a numeric
matrix of dimension \(n\times m\), and that C is a list of length
\(m\) where each element is a matrix of dimension \(n\times p_{j}\). For
the case where we assume that all exposures have common confounders we
set common=TRUE and C must be a matrix of dimension
\(n\times p\).
We define the multivariate generalized propensity score, mvGPS, as $$mvGPS=f_{\mathbf{D}\mid \mathbf{C}_{1},\dots,\mathbf{C}_{m}}$$ where \(f\) represents a joint multivariate conditional density function. For our current development we specify \(f\) as multivariate normal, i.e., $$\mathbf{D}\mid \mathbf{C}_{1},\dots,\mathbf{C}_{m}\sim N_{m}(\boldsymbol{\mu}, \Sigma).$$
Factorizing this joint density we can rewrite the mvGPS as the product of univariate conditional densities, i.e., $$mvGPS=f_{D_{m}\mid \mathbf{C}_{m}, D_{m-1},\dots,D_{1}}\times\cdots\times f_{D_{1}\mid\mathbf{C}_{1}}.$$ We use this factorized version in our implementation, with parameters for each distribution estimated through least squares.
Following robins2000;textualmvGPS, we use the mvGPS to construct stabilized inverse probability of treatment (IPTW) weights. These have been shown to balance confounders and return unbiased estimated of the dose-response. Weights are constructed as $$w=\frac{f_{\mathbf{D}}}{f_{\mathbf{D}\mid \mathbf{C}_{1},\dots,\mathbf{C}_{m}}},$$ where the marginal density \(f_{\mathbf{D}}\) of the exposures is assumed to be multivariate normal as well.
Writing the weights using completely factorized densities we have $$w=\frac{f_{D_{m}\mid D_{m-1},\dots, D_{1}}\times\cdots\times f_{D_{1}}}{f_{D_{m}\mid \mathbf{C}_{m}, D_{m-1},\dots,D_{1}}\times\cdots\times f_{D_{1}\mid\mathbf{C}_{1}}}.$$
Often when using weights based on the propensity score, practitioners are
concerned about the effect of extreme weights. It has been shown that an
effective way to protect extreme weights is to trim them at a particular
percentile crump2009dealing,lee2011weightmvGPS. We allow
users to specify whether trimmed weights should be returned and at which
threshold. To trim weights set trim_w=TRUE and specify the desired
percentile as trim_quantile=q. Note that trimming is applied at
both the upper and lower percentile thresholds, i.e.,
$$w^{*}=w 1_{\{Q(w, 1-q)\le w \le Q(w, q)\}}+Q(w, 1-q) 1_{\{w < Q(w, 1-q)\}} + Q(w, q) 1_{\{w > Q(w, q)\}}$$
# NOT RUN {
#generating confounded bivariate exposure
sim_dt <- gen_D(method="u", n=200, rho_cond=0.2, s_d1_cond=2, s_d2_cond=2, k=3,
C_mu=rep(0, 3), C_cov=0.1, C_var=1, d1_beta=c(0.5, 1, 0), d2_beta=c(0, 0.3, 0.75), seed=06112020)
D <- sim_dt$D
C <- sim_dt$C
#generating weights and mvGPS
out_mvGPS <- mvGPS(D=D, C=list(C[, 1:2], C[, 2:3]))
# can apply trimming with default 99th percentile
out_mvGPS_trim <- mvGPS(D=D, C=list(C[, 1:2], C[, 2:3]), trim_w=TRUE)
#or assume all exposures have equivalent confounders
out_mvGPS_common <- mvGPS(D=D, C=C, common=TRUE)
# }
Run the code above in your browser using DataLab