MatchThem:: Matching and Weighting after Multiple Imputation
Farhad Pishgar
Johns Hopkins University School of Medicine, BaltimoreNoah Greifer
Johns Hopkins Bloomberg School of Public Health, BaltimoreClémence Leyrat
London School of Hygiene & Tropical Medicine, LondonElizabeth Stuart
Johns Hopkins Bloomberg School of Public Health, Baltimore2021-08-17
Abstract
Balancing the distributions of the confounders across the exposure levels in an observational study through matching or weighting is an accepted method to control for confounding due to these variables when estimating the association between an exposure and outcome and reducing the degree of dependence on certain modeling assumptions. Despite the increasing popularity in practice, these procedures cannot be immediately applied to datasets with missing values. Multiple imputation of the missing data is a popular approach to account for missing values while preserving the number of units in the dataset and accounting for the uncertainty in the missing values. However, to the best of our knowledge, there is no comprehensive matching and weighting software that can be easily implemented with multiply imputed datasets. In this paper, we review this problem and suggest a framework to map out the matching and weighting of multiply imputed datasets to 5 actions as well as the best practices to assess balance in these datasets after matching and weighting. We also illustrate these approaches using a companion package for R, MatchThem.Introduction
Researchers often seek to estimate the effect of a treatment, exposure, or policy on an outcome but may be unable to randomly assign participants into the groups to be compared. The inability to randomize can lead to differences between the distributions of participant characteristics between exposure groups (known as covariate imbalance), which is the source of confounding bias in a naïve estimate of the exposure effect. When enough confounders—causes of exposure status and the outcome of interest—have been observed, one strategy for reducing this bias is to equate the confounder distributions across the exposure groups by matching or weighting units prior to estimating the exposure effect. Ideally, after matching or weighting, the exposure groups will be balanced, and a simple or covariate-adjusted estimate of the difference in average outcomes between the exposure groups will be unbiased for the true exposure effect (Stuart 2010). Matching and weighting can also enhance the robustness to misspecification of any outcome models used to estimate the exposure effect (Daniel E. Ho et al. 2007).
Despite increasing popularity in practice, matching and weighting methods cannot be immediately applied to datasets with missing values. There are several solutions to address the problem of missing data in causal effect estimation, but a standard and relatively easy-to-use strategy is multiple imputation of the missing data, which preserves the number of units in the dataset while accounting for some of the uncertainty in the missing values (Cham and West 2016). Multiple imputation involves filling in the missing data points using estimates of their values, repeating the process multiple times with randomness incorporated into each prediction to arrive at a set of multiple complete datasets containing the imputed values. The difficulty of analyzing multiply imputed data is that any analysis must be carried out within each imputed dataset, and the results pooled together using specific combining rules to arrive at a single set of estimates. Because matching and weighting are iterative, multi-step procedures, it is not straightforward to implement an analysis using these methods without extensive programming.
The MatchThem R package, which we introduce here, offers an analysis pipeline for estimating exposure effects using matching and weighting with multiply imputed data. The functions MatchThem offers blend seamlessly with functions used in other R packages for matching, weighting, and the generation and analysis of multiply imputed data. The aims of the present paper are to briefly review the issues around matching and weighting with multiply imputed data (section 2), to describe the structure and functionality of MatchThem (section 3), to describe the steps involved in implementing the best practices for these procedures (section 4), and to demonstrate the typical use of the MatchThem R package (section 5).
1.1. Notation
Let \(i = {1, 2, 3, ..., n}\) index the \(n\) units in a dataset, in which the causal effects of a binary exposure indicator (\(z\)) on an outcome (\(y\)) in the presence of a set of potential confounders (\(X = \{x_{1}, x_{2}, x_{3}, ... \}\)) are to be estimated (such that \(z_{i} = 0\) indicates that unit \(i\) is assigned to the control group and \(z_{i} = 1\) indicates that the unit \(i\) is assigned to the treated group) (Figure 1).
The typical context in which matching and weighting are used is one where data have been collected from an observational study in which the exposure is not randomized, yielding systematic differences between exposed and unexposed units on a set of measured potential confounders (often referred to as covariates). The situation we consider here is one in which the values of some of the covariates or the outcome are missing for some units in the observed dataset (we do not consider situations in which the exposure status is missing as the methods described herein may not apply to such scenarios). In order to account for the missingness in the covariates and outcome, the missing values are multiply imputed, creating \(m\) complete datasets. Though we briefly explain the procedure of multiple imputation in section 2, here, we focus on the procedures following imputation; see (White, Royston, and Wood 2011) and (Azur et al. 2011) for accessible introductions to multiple imputation for medical researchers.
1.2. Software requirements
The MatchThem package works with the R statistical software and programming language and can be installed in R (with version \(\geqslant 3.6.0\)) running on different platforms. MatchThem can be installed from the Comprehensive R Archive Network by executing the following commands in the R software console (the MatchThem package depends on the MatchIt (Daniel E. Ho et al. 2011) and WeightIt (Greifer 2020b) packages; these lines will install those packages too):
install.packages("MatchThem")
library(MatchThem)Matching, weighting, and missing data
2.1. Matching
Matching and weighting are methods to equate the distributions of the covariates between exposure groups (Stuart 2010). Matching does so by duplicating, selecting, or dropping units from the dataset in such a way that the resulting exposure groups have similar covariate distributions. Typically, matching relies on a distance measure constructed from the covariates to pair similar units between exposure groups, which then form the resulting matched sample; a popular distance measure is the propensity score difference, the propensity score being the predicted probability of being in the exposed group given the covariates. Propensity scores can be used in a variety of matching algorithms (Daniel E. Ho et al. 2011; Williamson et al. 2012), though other distance measures can be used as well (as there have been some recent concerns about the use of propensity scores for matching (King and Nielsen 2019)). Matching produces a set of matching weights (often 1 for those retained and 0 for those dropped) and matched pair membership, which can be incorporated into a regression of the outcome on the exposure to estimate the exposure effect in the matched sample. If the balance is achieved across the exposure groups in the matched sample, then bias in the exposure effect estimate will be reduced.
Matching is implemented in a number of R packages, but the
MatchIt package provides access to a variety of matching
methods for complete (i.e., non-missing) data. The MatchIt
function matchit() performs the requested form of matching
on the supplied dataset, returning an object from which matching weights
and matched pair membership can be extracted for use in effect
estimation. MatchThem interfaces with MatchIt to
extend matchit() for use with multiply imputed data.
2.2. Weighting
Weighting is another way to achieve balance and reduce bias in the estimate of an exposure effect. Weights for each unit can be estimated so that the distributions of the covariates are the same across the exposure groups in the weighted samples. The weights then function like survey weights and can be used in a weighted regression of the outcome on the exposure to estimate the exposure effect. A common way of estimating weights is to use a function of the estimated propensity score, a procedure known as inverse probability weighting (IPW), though there have been some developments that bypass estimating the propensity score to estimate the weights directly (Hainmueller 2012; Zubizarreta 2015).
The R package WeightIt implements a variety of weighting
methods and functions similarly to MatchIt. The
WeightIt’s function, weightit(), performs the
requested form of weighting and returns an object containing the
estimated weights. MatchThem interfaces with WeightIt
to extend weightit() for use with multiply imputed
data.
2.3. Assessing covariate balance
After matching or weighting, one must assess the degree to which the balancing method was successful at achieving covariate balance in the exposure groups. This involves using numerical and graphical criteria to compare the distributions of covariates across the groups (Austin 2009). If the balance is not achieved, the matching or weighting specification should be changed, and the procedure performed again until a satisfactory balance is found. The cobalt package provides tools for assessing balance after matching and weighting and has tools for summarizing balance in multiply imputed data. MatchThem interfaces with cobalt to facilitate balance checking as part of a complete analysis pipeline. We refer readers to the cobalt documentation for further explanation of cobalt’s capabilities in order to maintain focus on the structure and functionality of MatchThem, but we will include the use of cobalt in the demonstration of the analysis pipeline in section 5.
2.4. Missing data
A major obstacle for most matching and weighting procedures is that they cannot be performed in a straightforward way for units with missing values in the covariates because these procedures either search control and treat groups for units with similar covariate values or rely on the predictions from a model with the covariates as the predictors (i.e., the propensity scores), which cannot be computed in the presence of missing data.
Complete-case analysis, i.e., excluding units with missing values in the potential confounders or outcome, is a simple and naïve approach for handling missing data. However, the complete-case analysis may not be a valid option in all instances; the assumption of missingness completely-at-random, described below, is required to justify complete-case analysis and is often violated in observational data. In addition, it is possible that dropping units with any missing values may yield a dataset with few remaining units (Pigott 2001). Another approach is to replace the missing values with an arbitrary constant and include indicators for missingness as additional covariates in \(X\), though this can also allow bias to remain (Knol et al. 2010). The preferred method to address the problem of missing data that preserves the number of units in the dataset and often yields unbiased effect estimates is to impute the missing values using multiple imputation (Leyrat et al. 2019).
Multiple imputation refers to the procedure of substituting the missing values with a set of plausible values that reflect the uncertainty in predicting the true unobserved values, which results in \(m\) imputed (filled-in) datasets (Sterne et al. 2009). Multiple imputation is justified when the mechanism behind the missingness is ignorable, i.e., given the observed data, units with missing data represent a random subset of the dataset (‘missing-completely-at-random’ in Rubin’s language (Rubin 1987)) or when the probability that a value is missing relies on values of other observed variables, but not on the missing value itself or unobserved factors (‘missing-at-random’ in Rubin’s language (Rubin 1987)).
Several multiple imputation methods are described in the literature,
and multiple statistical packages can be used to generate multiply
imputed data. The general framework of these methods is the same: impute
the missing values to produce \(m\)
datasets, analyze the imputed datasets separately, and pool the results
obtained in each imputed dataset using standard combining rules to
arrive at a single estimate for the sample (Sterne et al. 2009; Rubin 1987). A popular
method of multiple imputation is multiple imputation with chained
equations (MICE), which involves iteratively fitting models to predict
the missing values and is implemented in the mice R
package. The mice package contains the functions
mice() to impute the missing values, with() to
run a supplied analysis model on each imputed dataset, and
pool() to pool the results of the models to arrive at a
single set of coefficient estimates and standard errors, facilitating
the creation and analysis of multiply imputed data in a single analysis
pipeline requiring minimal programming. We refer the reader to the
mice documentation for further details (van Buuren and Groothuis-Oudshoorn 2011).
2.5. Matching and weighting multiply imputed datasets
Given the limitations of conducting a complete-case analysis, multiply imputing missing data before matching or weighting has become a popular alternative for use with missing data. There has been some research examining the performance and optimal use of matching and weighting with multiply imputed data, with a focus on the correct sequence of actions involved. There are two main approaches that have been identified:
The within approach: In this approach, matching or weighting is performed within each imputed dataset, using the observed and imputed covariate values, and the exposure effects estimated in each of the \(m\) matched or weighted datasets are pooled together (Leyrat et al. 2019).
The across approach: In this approach, propensity scores are averaged across the imputed datasets, and, using this averaged measure, matching or weighting is performed in the imputed datasets. Finally, the estimated exposure effects obtained from analyzing the matched or weighted datasets are pooled together (Mitra and Reiter 2016).
The across approach has been demonstrated to have inferior statistical performance as compared to the within approach in many common scenarios (Leyrat et al. 2019; de Vries and Groenwold 2016), though early research favored its use (Mitra and Reiter 2016). In particular, the across approach seems most effective when the outcomes are not used to impute the missing covariate values (de Vries and Groenwold 2016). Although some recommend avoiding the inclusion of the outcome variable during or prior to matching and weighting with propensity scores (Rubin 2001), statistical evidence favors the use of the outcome variable in multiple imputation of covariates (Leyrat et al. 2019). In addition, the across approach is not compatible with matching and weighting methods that do not involve propensity scores, such as coarsened exact matching (de Vries and Groenwold 2016), Mahalanobis distance matching, and entropy balancing (Hainmueller 2012). Both approaches are implemented in MatchThem to facilitate comparison between them, though the within approach is the default in MatchThem functions and is the approach we recommend.
It should be noted that the across approach described by Mitra and Reiter (2016) differs slightly from that described here; in their procedure, the averaged propensity scores are used to estimate the causal effect in a single dataset consisting of just the observed exposure and outcome values, which are assumed to be non-missing. The procedure described here is in the spirit of the original method but allows for the presence of imputed outcomes and the use of imputed covariates in the effect estimation. When there is no missingness in the outcome and covariates are not used in the effect estimation, the two versions of this approach coincide.
Package contents and structure
MatchThem provides functions and S3 classes to
facilitate the use of matching and weighting with multiply imputed data
and the estimation of exposure effects and their uncertainty (i.e.,
standard errors), which requires special care when done with matched or
weighted multiply imputed data. In particular, MatchThem
extends the functionality of MatchIt and WeightIt for
matching and weighting to multiply imputed data and the functionality of
mice for the analysis of multiply imputed data to matched and
weighted data. MatchThem provides wrappers for functions in
these packages to create a smooth workflow requiring minimal
programming. Table 1 contains a summary of the
functions and classes contained in MatchThem.
| Function | Input | Output | Extends | Description |
matchthem() |
mids object |
mimids object |
MatchIt::matchit() |
Performs the requested form of matching on the imputed datasets. |
weightthem() |
mids object |
wimids object |
WeightIt::weightit() |
Performs the requested form of weighting on the imputed datasets. |
complete() |
mimids or wimids object |
data.frame\(^{1}\) |
mice::complete() |
Extracts one or more imputed datasets from the supplied input along with the outputs of the matching or weighting procedure. |
with() |
mimids or wimids object |
mimira object |
mice::with() |
Runs the supplied analysis model on each imputed dataset, incorporating the outputs of the matching or weighting procedure. |
pool() |
mimira object |
mimipo object\(^{2}\) |
mice::pool() |
Pools the coefficients and standard errors estimated across the imputed datasets to a single set of coefficient and standard error estimates. |
The MatchThem functions matchthem() and
weightthem() are wrappers for
MatchIt::matchit() and WeightIt::weightit()
that apply them to each imputed dataset, supplied in the form of a
"mids" object, the output of a call to
mice::mice(), which contains the multiply imputed datasets
(matchthem() and weightthem() are also
compatible with "amelia" objects from the Amelia
package, but they are first transformed into "mids" objects
before matching or weighting is performed on them).
matchthem() and weightthem() apply the
corresponding functions to the imputed datasets using the requested
approach, storing the outputs along with the original imputed data in a
"mimids" or "wimids" object, which extend
mice’s "mids" class to additionally contain
the matching and weighting output. The "mimids" and
"wimids" classes have a number of methods that extend
mice’s functions for analyzing "mids" objects;
in particular, MatchThem offers complete(),
with(), and pool(), which function similarly
to their equivalents in mice. MatchThem also contains
methods for cbind(), print(),
summary(), and plot() with
"mimids" and "wimids" objects. We describe the
syntax of these functions below.
3.1. matchthem()
The syntax for matchthem() is as follows:
matchthem(formula, datasets,
approach = "within",
method = "nearest",
distance = "logit",
distance.options = list(),
discard = "none",
reestimate = FALSE, ...)The formula argument corresponds to the model formula,
which relates the exposure (on the left-hand side) to the covariates.
The datasets argument corresponds to the
"mids" or "amelia" object containing the
multiply imputed datasets. The approach argument, with
options "within" and "across", corresponds to
the approach used. The method argument corresponds to the
method of matching used, which, as of now, can be one of the nearest
neighbor matching ("nearest"), optimal full matching
("full"), propensity score subclassification
("subclass"), optimal pair matching
("optimal"), exact matching ("exact"),
coarsened exact matching ("cem"), and genetic matching
("genetic"). The distance argument corresponds
to the method used to define the distances between units; it can be
"mahalanobis" for Mahalanobis distance matching or a method
of estimating propensity scores (the default, "glm",
estimates propensity scores using logistic regression). Note that only
the methods that involve propensity scores are allowed with the across
approach; as of now, these include "nearest",
"full", "subclass", "optimal",
and "genetic", and only when propensity scores are
requested. The other arguments, including those supplied to
..., control other aspects of the matching procedure and
are, along with formula, method, and distance, passed directly to
matchit().
3.2. weightthem()
The syntax for weightthem() is as follows:
weightthem(formula, datasets,
approach = "within",
method = "ps", ...)The formula, datasets, and
approach arguments have the same meaning as those for
matchthem(). The method argument controls the
weighting method used, which, as of now, can be one of logistic
regression propensity score weighting ("ps"), covariate
balancing propensity score weighting ("cbps") and its
nonparametric variety ("npcbps"), generalized boosted
modeling propensity score weighting ("gbm"), entropy
balancing ("ebal"), SuperLearner propensity score weighting
("super"), optimization-based weighting
("optweight"), energy balancing ("energy"),
and Bayesian additive regression tree propensity score weighting
("bart"). Note that only methods that involve propensity
scores can be used with the across approach; as of now, these include
"ps", "cbps", "gbm",
"super", and "bart". Arguments supplied to
... are passed to weightit() to control
details of the weight estimation.
3.3. complete()
The syntax for complete.mimids() is as follows (the
syntax for complete.wimids() is identical):
complete.mimids(data, action = 1,
include = FALSE, mild = FALSE, all = TRUE, ...)
complete.wimids(data, action = 1,
include = FALSE, mild = FALSE, all = TRUE, ...)complete() extracts the multiply imputed and matched or
weighted datasets from a "mimids" or "wimids"
object, yielding a "data.frame" containing the imputed data
and the outputs of the matching or weighting function. These additional
outputs include the estimated weights (which are produced with both
matching and weighting), matched pair membership, and estimated
propensity scores (if used). This function extends
complete.mids() from the mice package. The
data argument corresponds to a "mimids" or
"wimids" object, the output of a call to
matchthem() or weightthem(). The
action argument controls the format of the output; it can
be supplied as a number to extract a single dataset corresponding to
that imputation number or a string, such as "long", to
produce an object (of one of several types) containing all of the
imputed datasets; the mild argument provides another way to
control output. The all argument controls whether all units
should be included in the output or just units with a weight greater
than zero (i.e., units that have not been discarded or that were left
unmatched). Though the datasets produced by complete() can
be analyzed separately and the estimates manually combined using the
appropriate combining rules, the with() and
pool() methods facilitate proper analysis of the matched or
weighted imputed data.
3.4. with()
The syntax for with.mimids() and
with.wimids() are as follows:
with.mimids(data, expr, cluster, ...)
with.wimids(data, expr, ...)The with() method for "mimids" and
"wimids" objects extends the with() method for
"mids" objects provided in mice. It works by
extracting the imputed datasets one-by-one along with the matching or
weighting outputs and applies the modeling function supplied to
expr (e.g., a call to glm(),
survey::svyglm(), or survival::coxph() with
the outcome model formula included) to each of the datasets. No data
argument needs to be supplied to the modeling function because the
imputed datasets are automatically supplied by with().
with() automatically incorporates the estimated weights in
the estimation function when possible. The output of a call to with is a
"mimira" object, which contains the outputs of the models
run on each imputed dataset and extends the corresponding
"mira" class from mice.
In general, for generalized linear outcome models, the
svyglm() function in the survey
package should be used as it correctly accounts for the weights and
produces approximately correct standard errors, whereas the standard
errors resulting from a normal call to glm() will be
inaccurate. with() automatically constructs and supplies
the svydesign object containing the data and weights, so
neither need to be supplied. When matching is used, and a modeling
function from the survey package is supplied, information about
pair membership is also supplied to appropriately account for the
clustering in the standard error estimation as recommended by (Austin 2011) (this functionality can be
controlled using the cluster argument).
3.5. pool()
The syntax for pool() is as follows:
pool(object, dfcom = NULL)pool() takes in a "mimira" object (supplied
to the object argument) to pool the models and provide a
single set of coefficient estimates and information required to compute
their standard errors. The dfcom argument controls the
degrees of freedom used for the tests of the coefficients and confidence
intervals, which typically is close to the number of units in the
original dataset. Because matching and weighting can yield datasets with
different numbers of units remaining, the default is to use the smallest
degrees of freedom from the supplied models if possible; otherwise, a
large value is used to approximate a \(z\)-test. MatchThem re-exports
pool() as a generic with methods for "mira"
objects (the output of mice::with.mids()) and
"mimira" objects as mice::pool() is not a
generic. The output of MatchThem::pool.mimira() is a
"mimipo" object, which can be used with the methods
available in mice for "mipo" objects (e.g.,
summary() and print()).
3.6. Methods for "mimids" and "wimids"
objects
MatchThem also contains methods for cbind(),
print(), summary(), and plot()
with "mimids" and "wimids" objects.
cbind() allows one to add variables from an external
dataset, not included in the original "mids" object, that
one might wish to be involved in the effect estimation model, such as an
outcome not involved in the imputation or a variable collected after the
multiple imputation occurred. The print(),
summary(), and plot() methods simply apply the
corresponding function to the matchit or
weightit objects contained within the "mimids"
or "wimids" object. The MatchIt and
WeightIt documentation details their functionality. An
additional argument, n, determines to which imputed dataset
the function should be applied (e.g.,
print.mimids(x, n = 1) prints the output of the call to
matchit() used on the first imputed dataset).
Suggested workflow
The suggested workflow for pre-processing imputed datasets with matching or weighting and then analyzing them to estimate exposure effects using MatchThem is as follows (Figure 2):
Imputing the Missing Data in the Dataset: Data are imputed using functions in mice or Amelia. Data imputed using another package can be coerced to a
"mids"object by the mice functionas.mids()for use with MatchThem functions (van Buuren and Groothuis-Oudshoorn 2011; Honaker, King, and Blackwell 2011).Matching or Weighting the Imputed Datasets: Matching or weighting are performed using
matchthem()orweightthem(), respectively, on the imputed datasets. The functions perform the matching or weighting within each imputed dataset using the specified approach.Assessing Balance on the Matched or Weighted Datasets: Functions in the cobalt R package can be used to assess the balance to ensure that the resulting bias is small across imputed datasets. The
bal.tab()andlove.plot()functions in the cobalt package can be used directly on the output ofmatchthem()andweightthem()(Greifer 2020a). If the balance is not achieved, step 2 should be repeated with different approaches or methods until it is.Analyzing the Matched or Weighted Datasets: Using
with()function from MatchThem package, causal effects and their standard errors are estimated in each of the matched or weighted imputed datasets. Robust standard errors should be used with weighting and most matching methods and are available through integration with the survey package (Lumley 2004).Pooling the Causal Effect Estimates: The
pool()function from the package is used to pool the obtained causal effect estimates and standard errors from each dataset using Rubin’s rules.
Example
In this section, we review the suggested workflow for matching and
weighting multiply imputed datasets, using an example. The research
question in this context is whether osteoporosis at baseline is
associated with increased odds of developing knee osteoarthritis in the
follow-up or not (Figure 1). We will use the
osteoarthritis dataset (included in the MatchThem
package):
library(MatchThem)
data('osteoarthritis')The osteoarthritis dataset contains data on 7
characteristics (age: AGE, gender: SEX, body
mass index: BMI, racial background: RAC,
smoking status: SMK, osteoporosis at baseline:
OSP, and knee osteoarthritis in the follow-up:
KOA) of 2,585 individuals. The dataset contains missing
data in BMI, RAC, SMK, and
KOA variables. We assume the missing values are missing at
random, justifying the use of multiple imputation.
summary(osteoarthritis)5.1. Imputing the missing data in the dataset
We use the mice package to impute the missing data in the
osteoarthritis dataset. See the mice package
reference manual for more details about this step (van Buuren and Groothuis-Oudshoorn 2011). We
use 5 imputations for illustration, though more imputations are always
better.
library(mice)
imputed.datasets <- mice(osteoarthritis, m = 5)The code above produces the 5 imputed datasets and saves them in the
imputed.datasets object (of class "mids").
This object will be supplied to MatchThem functions to perform
matching and weighting in the imputed datasets.
5.2. Matching and weighting the imputed datasets
5.2.1. Matching the imputed datasets
In this example, we use matchthem() to match the
multiply imputed datasets, imputed.datasets, using the
"within" matching approach with nearest neighbor matching
on the propensity score, a caliper of 5% of the standard deviation of
the propensity score, and 2:1 unexposed-to-exposed matching ratio for
matching. The syntax is the same as it is for
MatchIt::matchit(), except that the
imputed.datasets is supplied to the datasets
argument (whereas matchit() takes a
"data.frame" for its data argument) and an
argument to approach is supplied to select the approach to
be used.
matched.datasets <- matchthem(OSP ~ AGE + SEX + BMI + RAC + SMK,
datasets = imputed.datasets,
approach = 'within',
method = 'nearest',
caliper = 0.05,
ratio = 2)
# Matching Observations | dataset: #1 #2 #3 #4 #5After the matching is performed in the 5 imputed datasets, the output
will be saved in the matched.datasets object (of
"mimids" class). The "mimids" object contains
the original imputed data and the output of the calls to
matchit() applied to each imputed dataset.
5.2.2. Weighting the imputed datasets
Here, we use weightthem() to weight the imputed
datasets, imputed.datasets, using the "across"
weighting approach with logistic regression propensity score weighting
targeting the average treatment effect in the matched sample (ATM)
estimand (which mimics the target population resulting from matching
with a caliper (Li and Greene 2013)). The
syntax is the same as it is for WeightIt::weightit() except
that the imputed.datasets is supplied to the
datasets argument (whereas weightit() takes a
"data.frame" for its data argument) and an
argument to approach is supplied to select the approach to
be used.
weighted.datasets <- weightthem(OSP ~ AGE + SEX + BMI + RAC + SMK,
datasets = imputed.datasets,
approach = 'across',
method = 'ps',
estimand = 'ATM')
# Estimating distances | dataset: #1 #2 #3 #4 #5
# Estimating weights | dataset: #1 #2 #3 #4 #5The weighted.datasets object (of "wimids"
class) contains the original imputed data and the output of the calls to
weightit() applied to each imputed dataset.
5.3. Assessing balance on the matched and weighted datasets
Functions in the cobalt package are compatible with
"mimids" and "wimids" objects, and the degree
to which balance was achieved in the matched and weighted datasets of
these objects can be assessed using the cobalt functions
bal.tab(), bal.plot(), and
love.plot(). Here, we illustrate the use of
bal.tab() to compute absolute standardized mean differences
(ASMDs) and Kolmogorov-Smirnov (KS) statistics for each covariate. The
code below produces the largest ASMD and KS statistic after matching for
each covariate across all the imputed datasets, indicating the worst
balance across the datasets (Greifer
2020a).
library(cobalt)
bal.tab(matched.datasets, stats = c('m', 'ks'),
imp.fun = 'max')
# Balance summary across all imputations
# Type Max.Diff.Adj Max.KS.Adj
# distance Distance 0.0107 0.0300
# AGE Contin. 0.0296 0.0450
# SEX_2 Binary 0.0021 0.0021
# BMI Contin. 0.0333 0.0557
# RAC_0 Binary 0.0021 0.0021
# RAC_1 Binary 0.0118 0.0118
# RAC_2 Binary 0.0139 0.0139
# RAC_3 Binary 0.0021 0.0021
# SMK Binary 0.0128 0.0128
# Average sample sizes across imputations
# Control Treated
# All 2106. 479.
# Matched (ESS) 738.17 467.2
# Matched (Unweighted) 810.2 467.2
# Unmatched 1295.8 11.8This information shows that the covariates are well-balanced in the osteoporosis negative and positive groups after matching as the largest ASMD and KS statistics for all covariates across the imputed datasets are close to zero. The sample size information below indicates that some units were left unmatched and dropped from the sample. The displayed values are averaged across the imputed datasets.
We can produce the same balance table for the weighted datasets:
bal.tab(weighted.datasets, stats = c('m', 'ks'),
imp.fun = 'max')
# Balance summary across all imputations
# Type Max.Diff.Adj Max.KS.Adj
# prop.score Distance 0.0040 0.0433
# AGE Contin. 0.0209 0.0376
# SEX_2 Binary 0.0002 0.0002
# BMI Contin. 0.0070 0.0466
# RAC_0 Binary 0.0002 0.0002
# RAC_1 Binary 0.0019 0.0019
# RAC_2 Binary 0.0027 0.0027
# RAC_3 Binary 0.0006 0.0006
# SMK Binary 0.0138 0.0138
# Average effective sample sizes across imputations
# Control Treated
# Unadjusted 2106. 479.
# Adjusted 954.44 478.3As with the matched datasets, the covariates are well balanced in the weighted datasets as demonstrated by the low values of the largest ASMD and KS statistics across the datasets. For more information on the available options for assessing balance with multiply imputed data, we refer the reader to the cobalt documentation (Greifer 2020a).
5.4. Analyzing the matched and weighted datasets
The exposure effect within each imputed dataset can be estimated
using with(), which applies the supplied outcome model to
each of the matched or weighted datasets and stores the output in a
"mimira" object. We illustrate the use of
with() below to estimate the difference in the log odds of
knee osteoarthritis between osteoporosis groups in the matched imputed
datasets:
library(survey)
matched.models <- with(matched.datasets,
svyglm(KOA ~ OSP, family = quasibinomial()),
cluster = TRUE)The models fit in each matched dataset are saved in the
matched.models object (of "mimira" class). We
can run the same code with the weighted imputed datasets:
weighted.models <- with(weighted.datasets,
svyglm(KOA ~ OSP, family = quasibinomial()))Results are saved in the weighted.models object (of
"mimira" class, note that in the calls to
svyglm(), no svydesign() or
weights arguments need to be specified as these are
automatically supplied by with())
5.5. Pooling the causal effect estimates
In order to arrive at a single set of coefficient and standard error
estimates from the imputed datasets, we must pool the estimated models
using pool(). We demonstrate this below on the
matched.models object containing the models we fit
above.
matched.results <- pool(matched.models)The output of the pool() is saved in the
matched.results object, which is of "mimipo"
class. We can run summary() to arrive at a final set of
estimates for the matched data:
summary(matched.results, conf.int = TRUE)
# term estimate std.error statistic df p.value 2.5 % 97.5 %
# 1 (Intercept) -0.2045680 0.08781317 -2.3295816 47.88960 0.024092 -0.38114 -0.027997
# 2 OSP1 -0.1255041 0.14138543 -0.8876733 53.09776 0.378720 -0.40908 0.158067The displayed results show that there is not sufficient evidence of an association between osteoporosis and knee osteoarthritis development in the follow-up in this sample (beta = -0.13 \[-0.41 -- 0.16\], odds ratio = 0.88 \[0.66 -- 1.17\]).
We can run pool() and then summary() on the
model fits in the weighted datasets to arrive at a similar table of
results:
weighted.results <- pool(weighted.models)
summary(weighted.results, conf.int = TRUE)
# term estimate std.error statistic df p.value 2.5 % 97.5 %
# 1 (Intercept) -0.1602278 0.06932782 -2.311162 225.930 0.02172525 -0.29684 -0.023616
# 2 OSP1 -0.1817342 0.13059573 -1.391579 63.768 0.16888557 -0.44265 0.079179Again, there is no evidence for an association between osteoporosis and knee osteoarthritis development in this sample.
Summary
Matching and weighting are popular methods used to balance the distributions of potential confounders across the exposure levels in an observational study. However, these procedures cannot be immediately applied to datasets with missing values. Multiple imputation of the missing data can be an effective approach to account for missing values while preserving the number of units in the dataset and accounting for the uncertainty in the imputation of the missing values. In this paper, we described the functionality of MatchThem, which interfaces with MatchIt, WeightIt, and mice to provide a smooth, straightforward workflow for estimating exposure effects in multiply imputed data using matching or weighting. MatchThem contains functions and classes that encourage the use of best practices in analyzing matched or weighted multiply imputed data without requiring extensive programming by the user. Given the ubiquity of missing data in observational studies, we hope MatchThem will facilitate smart choices and reliable analyses by researchers attempting to estimate causal effects from observational data.