Estimating Effects After Matching

Identifying the estimand

Before an effect is estimated, the estimand must be specified and clarified. Although some aspects of the estimand depend not only on how the effect is estimated after matching but also on the matching method itself, other aspects must be considered at the time of effect estimation and interpretation. Here, we consider three aspects of the estimand: the population the effect is meant to generalize to (the target population), the effect measure, and whether the effect is marginal or conditional.

The target population. Different matching methods allow you to estimate effects that can generalize to different target populations. The most common estimand in matching is the average treatment effect in the treated (ATT), which is the average effect of treatment for those who receive treatment. This estimand is estimable for matching methods that do not change the treated units (i.e., by weighting or discarding units) and is requested in matchit() by setting estimand = "ATT" (which is the default). The average treatment effect in the population (ATE) is the average effect of treatment for the population from which the sample is a random sample. This estimand is estimable only for methods that allow the ATE and either do not discard units from the sample or explicit target full sample balance, which in MatchIt is limited to full matching, subclassification, and profile matching when setting estimand = "ATE". When treated units are discarded (e.g., through the use of common support restrictions, calipers, cardinality matching, or [coarsened] exact matching), the estimand corresponds to neither the population ATT nor the population ATE, but rather to an average treatment effect in the remaining matched sample (ATM), which may not correspond to any specific target population. See Greifer and Stuart (2021) for a discussion on the substantive considerations involved when choosing the target population of the estimand.

Marginal and conditional effects. A marginal effect is a comparison between the expected potential outcome under treatment and the expected potential outcome under control. This is the same quantity estimated in randomized trials without blocking or covariate adjustment and is particularly useful for quantifying the overall effect of a policy or population-wide intervention. A conditional effect is the comparison between the expected potential outcomes in the treatment groups within strata. This is useful for identifying the effect of a treatment for an individual patient or a subset of the population.

Effect measures. The outcome types we consider here are continuous, with the effect measured by the mean difference; binary, with the effect measured by the risk difference (RD), risk ratio (RR), or odds ratio (OR); and time-to-event (i.e., survival), with the effect measured by the hazard ratio (HR). The RR, OR, and HR are noncollapsible effect measures, which means the marginal effect on that scale is not a (possibly) weighted average of the conditional effects within strata, even if the stratum-specific effects are of the same magnitude. For these effect measures, it is critical to distinguish between marginal and conditional effects because different statistical methods target different types of effects. The mean difference and RD are collapsible effect measures, so the same methods can be used to estimate marginal and conditional effects.

Our primary focus will be on marginal effects, which are appropriate for all effect measures, easily interpretable, and require few modeling assumptions. The “Common Mistakes” section includes examples of commonly used methods that estimate conditional rather than marginal effects and should not be used when marginal effects are desired.

G-computation

To estimate marginal effects, we use a method known as g-computation (Snowden, Rose, and Mortimer 2011) or regression estimation (Schafer and Kang 2008). This involves first specifying a model for the outcome as a function of the treatment and covariates. Then, for each unit, we compute their predicted values of the outcome setting their treatment status to treated, and then again for control, leaving us with two predicted outcome values for each unit, which are estimates of the potential outcomes under each treatment level. We compute the mean of each of the estimated potential outcomes across the entire sample, which leaves us with two average estimated potential outcomes. Finally, the contrast of these average estimated potential outcomes (e.g., their difference or ratio, depending on the effect measure desired) is the estimate of the treatment effect.

When doing g-computation after matching, a few additional considerations are required. First, when we take the average of the estimated potential outcomes under each treatment level, this must be a weighted average that incorporates the matching weights. Second, if we want to target the ATT or ATC, we only estimate potential outcomes for the treated or control group, respectively (though we still generate predicted values under both treatment and control).

G-computation as a framework for estimating effects after matching has a number of advantages over other approaches. It works the same regardless of the form of the outcome model or type of outcome (e.g., whether a linear model is used for a continuous outcome or a logistic model is used for a binary outcome); the only difference might be how the average expected potential outcomes are contrasted in the final step. In simple cases, the estimated effect is numerically identical to effects estimated using other methods; for example, if no covariates are included in the outcome model, the g-computation estimate is equal to the difference in means from a t-test or coefficient of the treatment in a linear model for the outcome. There are analytic approximations to the SEs of the g-computation estimate, and these SEs can incorporate pair/subclass membership (described in more detail below).

For all these reasons, we use g-computation when possible for all effect estimates, even if there are simpler methods that would yield the same estimates. Using a single workflow (with some slight modifications depending on the context; see below) facilitates implementing best practices regardless of what choices a user makes.

Modeling the Outcome

The goal of the outcome model is to generate good predictions for use in the g-computation procedure described above. The type and form of the outcome model should depend on the outcome type. For continuous outcomes, one can use a linear model regressing the outcome on the treatment; for binary outcomes, one can use a generalized linear model with, e.g., a logistic link; for time-to-event outcomes, one can use a Cox proportional hazards model.

An additional decision to make is whether (and how) to include covariates in the outcome model. One may ask, why use matching at all if you are going to model the outcome with covariates anyway? Matching reduces the dependence of the effect estimate on correct specification of the outcome model; this is the central thesis of Ho et al. (2007). Including covariates in the outcome model after matching has several functions: it can increase precision in the effect estimate, reduce the bias due to residual imbalance, and make the effect estimate “doubly robust”, which means it is consistent if either the matching reduces sufficient imbalance in the covariates or if the outcome model is correct. For these reasons, we recommend covariate adjustment after matching when possible. There is some evidence that covariate adjustment is most helpful for covariates with standardized mean differences greater than .1 (Nguyen et al. 2017), so these covariates and covariates thought to be highly predictive of the outcome should be prioritized in treatment effect models if not all can be included due to sample size constraints.

Although there are many possible ways to include covariates (e.g., not just main effects but interactions, smoothing terms like splines, or other nonlinear transformations), it is important not to engage in specification search (i.e., trying many outcomes models in search of the “best” one). Doing so can invalidate results and yield a conclusion that fails to replicate. For this reason, we recommend only including the same terms included in the propensity score model unless there is a strong a priori and justifiable reason to model the outcome differently.

It is important not to interpret the coefficients and tests of covariates in the outcome model. These are not causal effects and their estimates may be severely confounded. Only the treatment effect estimate can be interpreted as causal assuming the relevant assumptions about unconfoundedness are met. Inappropriately interpreting the coefficients of covariates in the outcome model is known as the Table 2 fallacy (Westreich and Greenland 2013). To avoid this, we only display the results of the g-computation procedure and do not examine or interpret the outcome models themselves.

Estimating Standard Errors and Confidence Intervals

Uncertainty estimation (i.e., of SEs, confidence intervals, and p-values) may consider the variety of sources of uncertainty present in the analysis, including (but not limited to!) estimation of the propensity score (if used), matching (i.e., because treated units might be matched to different control units if others had been sampled), and estimation of the treatment effect (i.e., because of sampling error). In general, there are no analytic solutions to all these issues, so much of the research done on uncertainty estimation after matching has relied on simulation studies. The two primary methods that have been shown to perform well in matched samples are using cluster-robust SEs and the bootstrap, described below.

To compute SEs after g-computation, a method known as the delta method is used; this is a way to compute the SEs of the derived quantities (the expected potential outcomes and their contrast) from the variance of the coefficients of the outcome models. For nonlinear models (e.g., logistic regression), the delta method is only an approximation subject to error (though in many cases this error is small and shrinks in large samples). Because the delta method relies on the variance of the coefficients from the outcome model, it is important to correctly estimate these variances, using either robust or cluster-robust methods as described below.

The Standard Case

For almost all matching methods, whether a caliper, common support restriction, exact matching specification, or \(k\):1 matching specification is used, estimating the effect in the matched dataset is straightforward and involves fitting a model for the outcome that incorporates the matching weights³, then estimating the treatment effect using g-computation (i.e., using marginaleffects::avg_comparisons()) with a cluster-robust SE to account for pair membership. This procedure is the same for continuous and binary outcomes with and without covariates.

There are a few adjustments that need to be made for certain scenarios, which we describe in the section “Adjustments to the Standard Case”. These adjustments include for the following cases: when matching for the ATE rather than the ATT, for matching with replacement, for matching with a method that doesn’t involve creating pairs (e.g., cardinality and profile matching and coarsened exact matching), for subclassification, for estimating effects with binary outcomes, and for estimating effects with survival outcomes. You must read the Standard Case to understand the basic procedure before reading about these special scenarios.

Here, we demonstrate the faster analytic approach to estimating confidence intervals; for the bootstrap approach, see the section “Using Bootstrapping to Estimate Confidence Intervals” below.

First, we will perform full matching on the propensity score for the ATT. Remember, all matching methods use this exact procedure or a slight variation, so this section is critical even if you are using a different matching method.

#Optimal full matching on the PS for the ATT
mF <- matchit(A ~ X1 + X2 + X3 + X4 + X5 + 
                 X6 + X7 + X8 + X9, data = d,
               method = "full", estimand = "ATT")
mF

## A `matchit` object
##  - method: Optimal full matching
##  - distance: Propensity score             - estimated with logistic regression
##  - number of obs.: 2000 (original), 2000 (matched)
##  - target estimand: ATT
##  - covariates: X1, X2, X3, X4, X5, X6, X7, X8, X9

#Extract matched data
md <- match.data(mF)

head(md)

##   A      X1      X2      X3       X4 X5      X6      X7      X8       X9      Y_C Y_B     Y_S distance weights subclass
## 1 0  0.1725 -1.4283 -0.4103 -2.36059  1 -1.1199  0.6398 -0.4840 -0.59385  0.07104   0  278.46  0.08461 0.20795      150
## 2 0 -1.0959  0.8463  0.2456 -0.12333  1 -2.2687 -1.4491 -0.5514 -0.31439  0.15619   0  330.63  0.01855 0.04475      119
## 3 0  0.1768  0.7905 -0.8436  0.82366  1 -0.2221  0.2971 -0.6966 -0.69516 -0.85180   1  369.94  0.22210 0.88379      288
## 4 0 -0.4595  0.1726  1.9542 -0.62661  1 -0.4019 -0.8294 -0.5384  0.20729 -2.35184   0   91.06  0.04180 0.17676      303
## 5 1  0.3563 -1.8121  0.8135 -0.67189  1 -0.8297  1.7297 -0.6439 -0.02648  0.68058   0  182.73  0.43291 1.00000       86
## 6 0 -2.4313 -1.7984 -1.2940  0.04609  1 -1.2419 -1.1252 -1.8659 -0.56513 -5.62260   0 2563.73  0.04998 0.12626      158

Typically one would assess balance and ensure that this matching specification works, but we will skip that step here to focus on effect estimation. See vignette("MatchIt") and vignette("assessing-balance") for more information on this necessary step. Because we did not use a caliper, the target estimand is the ATT.

We perform all analyses using the matched dataset, md, which, for matching methods that involve dropping units, contains only the units retained in the sample.

First, we fit a model for the outcome given the treatment and (optionally) the covariates. It’s usually a good idea to include treatment-covariate interactions, which we do below, but this is not always necessary, especially when excellent balance has been achieved. You can also include the propensity score (usually labeled distance in the match.data() output), which can add some robustness, especially when modeled flexibly (e.g., with polynomial terms or splines) (Austin 2017); see here for an example.

#Linear model with covariates
fit1 <- lm(Y_C ~ A * (X1 + X2 + X3 + X4 + X5 + 
                        X6 + X7 + X8 + X9),
           data = md, weights = weights)

Next, we use marginaleffects::avg_comparisons() to estimate the ATT.

avg_comparisons(fit1, variables = "A",
                vcov = ~subclass,
                newdata = subset(A == 1))

## 
##  Estimate Std. Error    z Pr(>|z|)    S 2.5 % 97.5 %
##      1.75      0.479 3.65   <0.001 11.9  0.81   2.69
## 
## Term: A
## Type:  response 
## Comparison: mean(1) - mean(0)
## Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted

Let’s break down the call to avg_comparisons(): to the first argument, we supply the model fit, fit1; to the variables argument, the name of the treatment ("A"); to the vcov argument, a formula with subclass membership (~subclass) to request cluster-robust SEs; and to the newdata argument, a version of the matched dataset containing only the treated units (subset(A == 1)) to request the ATT. Some of these arguments differ depending on the specifics of the matching method and outcome type; see the sections below for information.

If, in addition to the effect estimate, we want the average estimated potential outcomes, we can use marginaleffects::avg_predictions(), which we demonstrate below. Note the interpretation of the resulting estimates as the expected potential outcomes is only valid if all covariates present in the outcome model (if any) are interacted with the treatment.

avg_predictions(fit1, variables = "A",
                vcov = ~subclass,
                newdata = subset(A == 1))

## 
##  A Estimate Std. Error     z Pr(>|z|)     S 2.5 % 97.5 %
##  0     2.14      0.382  5.61   <0.001  25.6  1.39   2.89
##  1     3.89      0.229 16.99   <0.001 212.7  3.44   4.34
## 
## Type:  response 
## Columns: A, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high

We can see that the difference in potential outcome means is equal to the average treatment effect computed previously⁴. All of the arguments to avg_predictions() are the same as those to avg_comparisons().

Adjustments to the Standard Case

This section explains how the procedure might differ if any of the following special circumstances occur.

#Subclassification on the PS for the ATT
mS <- matchit(A ~ X1 + X2 + X3 + X4 + X5 + 
                 X6 + X7 + X8 + X9, data = d,
               method = "subclass", estimand = "ATT")

#Extract matched data
md <- match.data(mS)

fitS <- lm(Y_C ~ subclass * (A * (X1 + X2 + X3 + X4 + X5 + 
                                    X6 + X7 + X8 + X9)),
           data = md)

avg_comparisons(fitS,
                variables = "A",
                vcov = "HC3",
                newdata = subset(A == 1))

## 
##  Estimate Std. Error    z Pr(>|z|)    S 2.5 % 97.5 %
##      1.65      0.364 4.54   <0.001 17.4  0.94   2.37
## 
## Term: A
## Type:  response 
## Comparison: mean(1) - mean(0)
## Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted

#Logistic regression model with covariates
fit2 <- glm(Y_B ~ A * (X1 + X2 + X3 + X4 + X5 + 
                        X6 + X7 + X8 + X9),
            data = md, weights = weights,
            family = quasibinomial())

#Compute effects; RR and confidence interval
avg_comparisons(fit2,
                variables = "A",
                vcov = ~subclass,
                newdata = subset(A == 1),
                comparison = "lnratioavg",
                transform = "exp")

## 
##  Estimate Pr(>|z|)    S 2.5 % 97.5 %
##      1.57   <0.001 52.0  1.41   1.75
## 
## Term: A
## Type:  response 
## Comparison: ln(mean(1) / mean(0))
## Columns: term, contrast, estimate, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted

library("survival")

#Cox Regression for marginal HR
coxph(Surv(Y_S) ~ A, data = md, robust = TRUE, 
      weights = weights, cluster = subclass)

## Call:
## coxph(formula = Surv(Y_S) ~ A, data = md, weights = weights, 
##     robust = TRUE, cluster = subclass)
## 
##   coef exp(coef) se(coef) robust se  z      p
## A 0.45      1.57     0.05      0.04 12 <2e-16
## 
## Likelihood ratio test=63  on 1 df, p=2e-15
## n= 2000, number of events= 2000

#get_matches() after matching with replacement
gm <- get_matches(mR)

#Austin & Cafri's (2020) SE estimator
fs <- coxph(Surv(Y_S) ~ A, data = gm, robust = TRUE, 
      weights = weights, cluster = subclass)
Vs <- fs$var
ks <- nlevels(gm$subclass)

fi <- coxph(Surv(Y_S) ~ A, data = gm, robust = TRUE, 
      weights = weights, cluster = id)
Vi <- fi$var
ki <- length(unique(gm$id))

fc <- coxph(Surv(Y_S) ~ A, data = gm, robust = TRUE, 
      weights = weights)
Vc <- fc$var
kc <- nrow(gm)

#Compute the variance and sneak it back into the fit object
fc$var <- (ks/(ks-1))*Vs + (ki/(ki-1))*Vi - (kc/(kc-1))*Vc

fc

Using Bootstrapping to Estimate Confidence Intervals

The bootstrap is an alternative to the delta method for estimating confidence intervals for estimated effects. See the section Bootstrapping above for details. Here, we’ll demonstrate two forms of the bootstrap: 1) the standard bootstrap, which involve resampling units and performing matching and effect estimation within each bootstrap sample, and 2) the cluster bootstrap, which involves resampling pairs after matching and estimating the effect in each bootstrap sample. For both, we will use functionality in the boot package. It is critical to set a seed using set.seed() prior to performing the bootstrap in order for results to be replicable.

boot_fun <- function(data, i) {
  boot_data <- data[i,]
  
  #Do 1:1 PS matching with replacement
  m <- matchit(A ~ X1 + X2 + X3 + X4 + X5 + 
                 X6 + X7 + X8 + X9,
               data = boot_data,
               replace = TRUE)
  
  #Extract matched dataset
  md <- match.data(m, data = boot_data)
  
  #Fit outcome model
  fit <- glm(Y_B ~ A * (X1 + X2 + X3 + X4 + X5 + 
                 X6 + X7 + X8 + X9),
             data = md, weights = weights,
             family = quasibinomial())
  
  ## G-computation ##
  #Subset to treated units for ATT; skip for ATE
  md1 <- subset(md, A == 1)
  
  #Estimated potential outcomes under treatment
  p1 <- predict(fit, type = "response",
                newdata = transform(md1, A = 1))
  Ep1 <- mean(p1)
  
  #Estimated potential outcomes under control
  p0 <- predict(fit, type = "response",
                newdata = transform(md1, A = 0))
  Ep0 <- mean(p0)
  
  #Risk ratio
  return(Ep1 / Ep0)
}

library("boot")
set.seed(54321)
boot_out <- boot(d, boot_fun, R = 199)

boot_out

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = d, statistic = boot_fun, R = 199)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    1.347  0.1417      0.1937

boot.ci(boot_out, type = "perc")

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 199 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot_out, type = "perc")
## 
## Intervals : 
## Level     Percentile     
## 95%   ( 1.144,  1.891 )  
## Calculations and Intervals on Original Scale
## Some percentile intervals may be unstable

mNN <- matchit(A ~ X1 + X2 + X3 + X4 + X5 + 
                 X6 + X7 + X8 + X9, data = d)
mNN

## A `matchit` object
##  - method: 1:1 nearest neighbor matching without replacement
##  - distance: Propensity score             - estimated with logistic regression
##  - number of obs.: 2000 (original), 882 (matched)
##  - target estimand: ATT
##  - covariates: X1, X2, X3, X4, X5, X6, X7, X8, X9

md <- match.data(mNN)

#Unique pair IDs
pair_ids <- levels(md$subclass)

#Unit IDs, split by pair membership
split_inds <- split(seq_len(nrow(md)), md$subclass)

cluster_boot_fun <- function(pairs, i) {
  
  #Extract units corresponding to selected pairs
  ids <- unlist(split_inds[pairs[i]])
  
  #Subset md with block bootstrapped indices
  boot_md <- md[ids,]
  
  #Fit outcome model
  fit <- glm(Y_B ~ A * (X1 + X2 + X3 + X4 + X5 + 
                 X6 + X7 + X8 + X9),
             data = boot_md, weights = weights,
             family = quasibinomial())
  
  ## G-computation ##
  #Subset to treated units for ATT; skip for ATE
  md1 <- subset(boot_md, A == 1)
  
  #Estimated potential outcomes under treatment
  p1 <- predict(fit, type = "response",
                newdata = transform(md1, A = 1))
  Ep1 <- mean(p1)
  
  #Estimated potential outcomes under control
  p0 <- predict(fit, type = "response",
                newdata = transform(md1, A = 0))
  Ep0 <- mean(p0)
  
  #Risk ratio
  return(Ep1 / Ep0)
}

library("boot")
set.seed(54321)
cluster_boot_out <- boot(pair_ids, cluster_boot_fun,
                         R = 199)

cluster_boot_out

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = pair_ids, statistic = cluster_boot_fun, R = 199)
## 
## 
## Bootstrap Statistics :
##     original   bias    std. error
## t1*    1.588 0.001319      0.1265

boot.ci(cluster_boot_out, type = "perc")

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 199 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = cluster_boot_out, type = "perc")
## 
## Intervals : 
## Level     Percentile     
## 95%   ( 1.356,  1.857 )  
## Calculations and Intervals on Original Scale
## Some percentile intervals may be unstable

Moderation Analysis

Moderation analysis involves determining whether a treatment effect differs across levels of another variable. The use of matching with moderation analysis is described in Green and Stuart (2014). The goal is to achieve balance within each subgroup of the potential moderating variable, and there are several ways of doing so. Broadly, one can either perform matching in the full dataset, requiring exact matching on the moderator, or one can perform completely separate analyses in each subgroup. We’ll demonstrate the first approach below; see the blog post “Subgroup Analysis After Propensity Score Matching Using R” by Noah Greifer for an example of the other approach.

There are benefits to using either approach, and Green and Stuart (2014) find that either can be successful at balancing the subgroups. The first approach may be most effective with small samples, where separate propensity score models would be fit with greater uncertainty and an increased possibility of perfect prediction or failure to converge (Wang et al. 2018). The second approach may be more effective with larger samples or with matching methods that target balance in the matched sample, such as genetic matching (Kreif et al. 2012). With genetic matching, separate subgroup analyses ensure balance is optimized within each subgroup rather than just overall. The chosen approach should be that which achieves the best balance, though we don’t demonstrate assessing balance here to maintain focus on effect estimation.

The full dataset approach involves pooling information across subgroups. This could involve estimating propensity scores using a single model for both groups but exact matching on the potential moderator. The propensity score model could include moderator-by-covariate interactions to allow the propensity score model to vary across subgroups on some covariates. It is critical that exact matching is done on the moderator so that matched pairs are not split across subgroups.

We’ll consider the binary variable X5 to be the potential moderator of the effect of A on Y_C. Below, we’ll estimate a propensity score using a single propensity score model with a few moderator-by-covariate interactions. We’ll perform nearest neighbor matching on the propensity score and exact matching on the moderator, X5.

mP <- matchit(A ~ X1 + X2 + X5*X3 + X4 + 
                X5*X6 + X7 + X5*X8 + X9, data = d,
              exact = ~X5, method = "nearest")
mP

## A `matchit` object
##  - method: 1:1 nearest neighbor matching without replacement
##  - distance: Propensity score             - estimated with logistic regression
##  - number of obs.: 2000 (original), 882 (matched)
##  - target estimand: ATT
##  - covariates: X1, X2, X5, X3, X4, X6, X7, X8, X9

Although it is straightforward to assess balance overall using summary(), it is more challenging to assess balance within subgroups. The easiest way to check subgroup balance would be to use cobalt::bal.tab(), which has a cluster argument that can be used to assess balance within subgroups, e.g., by cobalt::bal.tab(mP, cluster = "X5"). See the vignette “Appendix 2: Using cobalt with Clustered, Multiply Imputed, and Other Segmented Data” on the cobalt website for details.

If we are satisfied with balance, we can then model the outcome with an interaction between the treatment and the moderator.

mdP <- match.data(mP)

fitP <- lm(Y_C ~ A * X5, data = mdP, weights = weights)

To estimate the subgroup ATTs, we can use avg_comparisons(), this time specifying the by argument to signify that we want treatment effects stratified by the moderator.

avg_comparisons(fitP, variables = "A",
                vcov = ~subclass,
                newdata = subset(A == 1),
                by = "X5")

## 
##  X5 Estimate Std. Error    z Pr(>|z|)    S 2.5 % 97.5 %
##   0     2.21      0.670 3.29   <0.001 10.0 0.894   3.52
##   1     2.18      0.569 3.83   <0.001 12.9 1.065   3.29
## 
## Term: A
## Type:  response 
## Comparison: mean(1) - mean(0)
## Columns: term, contrast, X5, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted

We can see that the subgroup mean differences are quite similar to each other. Finally, we can test for moderation using another call to avg_comparisons(), this time using the hypothesis argument to signify that we want to compare effects between subgroups:

avg_comparisons(fitP, variables = "A",
                vcov = ~subclass,
                newdata = subset(A == 1),
                by = "X5",
                hypothesis = "pairwise")

## 
##  Estimate Std. Error      z Pr(>|z|)   S 2.5 % 97.5 %
##    0.0275      0.879 0.0313    0.975 0.0  -1.7   1.75
## 
## Term: 0 - 1
## Type:  response 
## Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high

As expected, the difference between the subgroup treatment effects is small and nonsignificant, so there is no evidence of moderation by X5.

Reporting Results

It is important to be as thorough and complete as possible when describing the methods of estimating the treatment effect and the results of the analysis. This improves transparency and replicability of the analysis. Results should at least include the following:

• a description of the outcome model used (e.g., logistic regression, a linear model with treatment-covariate interactions and covariates, a Cox proportional hazards model with the matching weights applied)
• the way the effect was estimated (e.g., using g-computation or as the coefficient in the outcome model)
• the way SEs and confidence intervals were estimated (e.g., using robust SEs, using cluster-robust SEs with pair membership as the cluster, using the BCa bootstrap with 4999 bootstrap replications and the entire process of matching and effect estimation included in each replication)
• R packages and functions used in estimating the effect and its SE (e.g., glm() in base R, avg_comparisons() in marginaleffects, boot() and boot.ci() in boot)
• The effect and its SE and confidence interval

All this is in addition to information about the matching method, propensity score estimation procedure (if used), balance assessment, etc. mentioned in the other vignettes.

1. Failing to include weights

Several methods involve weights that are to be used in estimating the treatment effect. With full matching and stratification matching (when analyzed using MMWS), the weights do the entire work of balancing the covariates across the treatment groups. Omitting weights essentially ignores the entire purpose of matching. Some cases are less obvious. When performing matching with replacement and estimating the treatment effect using the match.data() output, weights must be included to ensure control units matched to multiple treated units are weighted accordingly. Similarly, when performing k:1 matching where not all treated units receive k matches, weights are required to account for the differential weight of the matched control units. The only time weights can be omitted after pair matching is when performing 1:1 matching without replacement. Including weights even in this scenario will not affect the analysis and it can be good practice to always include weights to prevent this error from occurring. There are some scenarios where weights are not useful because the conditioning occurs through some other means, such as when using the direct subclass strategy rather than MMWS for estimating marginal effects after stratification.

2. Failing to use robust or cluster-robust standard errors

Robust SEs are required when using weights to estimate the treatment effect. The model-based SEs resulting from weighted least squares or maximum likelihood are inaccurate when using matching weights because they assume weights are frequency weights rather than probability weights. Cluster-robust SEs account for both the matching weights and pair membership and should be used when appropriate. Sometimes, researchers use functions in the survey package to estimate robust SEs, especially with inverse probability weighting; this is a valid way to compute robust SEs and will give similar results to sandwich::vcovHC().¹⁰

3. Interpreting conditional effects as marginal effects

The distinction between marginal and conditional effects is not always clear both in methodological and applied papers. Some statistical methods are valid only for estimating conditional effects and they should not be used to estimate marginal effects (without further modification). Sometimes conditional effects are desirable, and such methods may be useful for them, but when marginal effects are the target of inference, it is critical not to inappropriately interpret estimates resulting from statistical methods aimed at estimating conditional effects as marginal effects. Although this issue is particularly salient with binary and survival outcomes due to the general noncollapsibility of the OR, RR, and HR, this can also occur with linear models for continuous outcomes or the RD.

The following methods estimate conditional effects for binary or survival outcomes (with noncollapsible effect measures) and should not be used to estimate marginal effects:

• Logistic regression or Cox proportional hazards model with covariates and/or the propensity score included, using the coefficient on treatment as the effect estimate
• Conditional logistic regression after matching
• Stratified Cox regression after matching
• Averaging stratum-specific effect estimates after stratification, including using Mantel-Haenszel OR pooling
• Including pair or stratum fixed or random effects in a logistic regression model, using the coefficient on treatment as the effect estimate

In addition, with continuous outcomes, conditional effects can be mistakenly interpreted as marginal effect estimates when treatment-covariate interactions are present in the outcome model. If the covariates are not centered at their mean in the target population (e.g., the treated group for the ATT, the full sample for the ATE, or the remaining matched sample for an ATM), the coefficient on treatment will not correspond to the marginal effect in the target population; it will correspond to the effect of treatment when the covariate values are equal to zero, which may not be meaningful or plausible. G-computation is always the safest way to estimate effects when including covariates in the outcome model, especially in the presence of treatment-covariate interactions.