Evaluating Different Covariate Balancing Methods: A Monte Carlo Simulation

2.1 Randomized Controlled Trials, Exact Matching, and Propensity Score Matching

In many scientific fields, randomized controlled trials (RCTs) are regarded as the most reliable method of estimating causal effects and are frequently incorporated into real-world policy-making processes (Banerjee and Duflo 2011; Gneezy and List 2014; Leigh 2018).

However, conducting RCTs on every policy is not feasible due to various factors, such as budgetary constraints and resistance arising from individual or organizational interests. According to Holland (1986), some “treatments” of interest cannot be randomly assigned because it would be unethical. Moreover, experiments only offer a limited range of outcome distributions, making it challenging to identify parameters (e.g. the distribution of treatment effects) that depend on joint distribution. Therefore, we cannot eliminate using observational data for policy analysis and evaluation.

Nonetheless, in observational data, the treatment and comparison groups are likely to be affected differently by factors (covariates) other than the treatment, unlike those in experimental data. Consequently, these factors (covariates) may impact the outcome difference between the treatment and comparison groups when one uses unadjusted observational data.

Matching is one technique for dealing with sample selection biases caused by non-randomized data that are inherent in observational data analysis. We create pairs by selecting observation units from the treatment and comparison groups with similar values of covariates other than the treatment. As a result of this matching, the effects of covariates other than the treatment are, on average, similar in the two groups. For example, in exact matching (EM), we create pairs by taking observation units with the same values for all covariates from the treatment and comparison groups, which can reduce the bias caused by the covariate imbalance between the treatment and comparison groups to zero. This allows us to estimate the effect of the treatment on the outcome more accurately. The important point is that matching adjusts only for bias due to observed covariates and does not correct for bias derived from unobserved covariates.

However, simple matching, such as EM, becomes infeasible as the covariates increase. Indeed, EM is infeasible with even one continuous covariate. PSM can effectively deal with this problem, referred to as the curse of dimensionality, which occurs as more covariates are taken into account in the covariate matrix X, causing a sharp decline in the probability of achieving exact matches for all units. In other words, inexact matching, rather than EM, is required when there are many discrete covariates and/or at least one continuous covariate. To perform inexact matching, some distance metric must be used to measure the similarity of the treated and untreated units concerning the values of the observable variables. Propensity scores and Mahalanobis distance are two examples of such distance metrics. In PSM, we estimate and quantify the “propensity” of each observation unit to receive the treatment by parametric estimators of binary response models, such as the probit and logit^[2] regression, by using observed covariates. The predicted value from the regression is the propensity score, i.e. the conditional probability that observation unit i(i = 1,⋯N) is assigned a treatment (W _i  = 1; W _i  = 0 for non-treatment) under given observed covariates x _i (x _i  ∈ X _i ). Put simply, the propensity score is a single value for each unit of observation that summarizes the effect of multiple covariates on treatment assignment.

There are several PSM methods. However, because King and Nielsen (2019) criticize PSM with 1-to-1 greedy matching, we limit our discussion to this type of PSM. In nearest-neighbor caliper matching, we create pairs for a given observation unit i(i ∈ I₁) in the treatment group by taking the observation unit j(j ∈ I₀) in the comparison group, given that the difference between the propensity scores of i and j is less than or equal to a predetermined acceptable value (caliper) ε. Suppose C is the set of observation units of the comparison group; P _i and P _j are the propensity scores of the observation units in the treatment and comparison groups, respectively; and C(i) is the set of observation units j in the comparison group that is matched to the observation unit i of the treatment group whose propensity score is P _i . For observation unit i in the treatment group, observation unit j in the comparison group is selected if the absolute difference in propensity scores between i and j satisfies the following conditions:

As a result of such matching, the overall effect of covariates is approximately homogenized between the treatment and matched comparison groups, allowing for randomization although individual covariates may differ considerably (for details and variations of PSM methods, see Rosenbaum and Rubin (1983), Dehejia and Wahba (1999, 2002, Caliendo and Kopeinig (2008), Rosenbaum (2010), Guo and Fraser (2014), and Imbens and Rubin (2015), among others). Although various matching algorithms have been developed, PSM is one of the most extensively used. As of September 4, 2021, 28,330 academic papers referenced PSM in the Scopus literature search. In contrast, only 426 and 428 academic papers mentioned MDM and CEM, respectively.

2.2 Weaknesses of PSM With 1-to-1 Greedy Matching

King and Nielsen (2019) contend that PSM with 1-to-1 greedy matching has serious problems. According to them, the weakness of PSM is that it is modeled after a fully randomized controlled trial, in which treatment assignment is determined solely by the probability that each observation unit receives the treatment, independent of the covariates. In other words, when treatment is assigned via randomization, there is no relationship between covariates and treatment status in the population. As a result, the observation unit pairs are not necessarily similar in terms of their covariates even though the difference in propensity scores of the paired observation units is less than or equal to a specific predetermined acceptance value and is minimal. PSM effectively addresses the curse of dimensionality by allowing researchers to match observation units that are different in terms of individual covariates. However, because of this feature, PSM with 1-to-1 greedy matching worsens the covariate balance at a certain point and increases bias in estimating treatment effects.

2.3 Advantages of MDM and CEM

In a fully blocked RCT, treatments are randomly assigned after observation units with similar attributes are grouped or stratified based on the value of covariates that may affect the outcome. Thus, if ideally implemented, a fully blocked RCT will balance the covariates in the treatment and comparison groups. In other words, a fully blocked randomized experimental design will block the treated and comparison groups on the observed covariates exactly at the start. Therefore, the imbalance in these experiments is always zero by design without having to remove any observation units from the sample. In contrast, in a completely randomized experimental design, treatment assignment is random in regard to the observed covariates because it only depends on the scalar probability of treatment for all units. Randomness does not mean there is zero imbalance in any one sample (except by rare coincidence or in asymptotic samples or when modified to create randomized block designs). PSM approximates a fully randomized experiment whereas CEM and MDM approximate a totally blocked experiment. King and Nielsen (2019) argue that this means that CEM and MDM can achieve lower degrees of imbalance, model dependence, and bias than PSM (see also Iacus, King, and Porro 2011).

MDM calculates distance d(i, j), i.e. Mahalanobis distance, between any observation unit in the treatment group and all observation units in the comparison group, which is defined as follows:

where u and v are the values of the covariates used for matching observation unit i in the treatment group and observation unit j in the comparison group in the case of single-nearest-neighbor matching. ∑ is the sample variance-covariance matrix of the covariates used for matching. The observation unit j of the comparison group with the smallest distance d(i, j) is selected and matched to the observation unit i of the treatment group. The Mahalanobis distance increases with the degree of covariate difference. Therefore, each pair will have similar covariate values if the comparison units near the treated units are located within the Mahalanobis distance. In contrast, PSM matches units with similar propensity scores. Because the covariate distribution is reduced into a single dimension by propensity scores, it follows that in PSM, the unit pairs are not necessarily similar regarding their covariates.

As its name implies, CEM is a matching method in the EM family. EM can effectively remove errors introduced by covariate imbalance between treatment and comparison groups with a small number of discrete covariates. However, EM suffers from the curse of dimensionality and tends to reduce the number of observation units to be matched because real-world data usually contain multiple covariates and continuous variables. CEM attempts to overcome this curse of dimensionality differently than PSM. In CEM’s algorithm, each covariate is grouped by a coarse index, which matches the observation units in the treatment and comparison groups. The original covariates used to create the coarse indices are retained in the matched data and used in subsequent analyses (e.g. a second-step analysis that estimates a parametric model using weights provided by CEM; Blackwell et al. 2009; Iacus, King, and Porro 2011, 2012).

CEM’s algorithm allows for matching on coarse indices such as school systems. For example, the covariate of education level is represented by a continuous variable of years of education. However, in the US institutional context, years of education can also be grouped into five coarse indices: 0–6 years for grade school students, 7–8 years for middle school students, 9–12 years for high school students, 13–16 years for college students, and 17 years or more for graduate school students.

3.1 Responses From other Researchers

The results of the simulation analyses by King and Nielsen (2019) suggest that MDM and CEM, which approximate a fully blocked randomized experiment, achieve better covariate balance than PSM, which approximates a fully randomized controlled experiment. King and Nielsen (2019) also suggest that the causal effect estimates from MDM and CEM are less biased than the estimates from PSM.

Inspired by the preprint version of their study, Ripollone et al. (2018) investigated whether King and Nielsen’s (2019) findings can be replicated when PSM is applied to real observational data rather than simulation data. They used US insurance claims data from the Pharmaceutical Assistance Contract for the Elderly and the Medicaid Analytic eXtract databases to test whether further pruning of matched sets would lead to a rise in covariate imbalance. Their analysis results suggest that covariate imbalance worsened when matched sets were gradually pruned in descending order of propensity score distance calipers (0.01–0.05). However, they also found that the imbalance in the matched data set improved after it was applied with commonly used propensity score calipers to define an acceptable match. Indeed, when they used standard propensity score calipers (e.g. 0.025), the data pruning stopped in the lowest region of the imbalance trend. Ripollone et al. (2018) concluded that when standard propensity score calipers are used on the aforementioned observational data, PSM does not increase covariate imbalance.

Guo, Fraser, and Chen (2020) agree with King and Nielsen (2019) that the PSM aims to approximate a fully randomized controlled experiment and that, in many cases, the matching methods that approximate fully blocked randomized experiments are more desirable for studies using observational data. However, in empirical research, priority is typically given to ensure that the matched sample remains representative of the research population of interest. In this context, losing a substantive proportion of observations can be problematic. Guo, Fraser, and Chen (2020) argue that PSM is not necessarily inferior to MDM and CEM in maintaining observation units in the sample.

Furthermore, Guo, Fraser, and Chen (2020) contend that blocking can lead to incomplete matching. Indeed, appropriate blocking measures are necessary to implement a matching design approximating a fully blocked randomized experiment. When appropriate blocking measures are not evident, blocked matching can result in imprecise estimates. After conducting simulation analysis with two conditions of data generation, i.e. selection on the observables and selection on the unobservables, Guo, Fraser, and Chen (2020) state that they cannot conclude that the alternative methods King and Nielsen (2019) recommend more effectively reduce finite-sample bias than PSM.

3.2 Critical Review

There are several concerns regarding the analyses by King and Nielsen (2019) and their critics. King and Nielsen (2019) draw 100 comparison units and 100 treated units from uniform distributions with limits of 0 and 5 (Uniform(0, 5)) and 1 and 6 (Uniform(1, 6)), respectively, for each of the two covariates, X_i1 and X_i2, of the equations below. The observations within the overlapping [1, 5] × [1, 5] squares are considered the data from a completely random experiment. In contrast, the observations outside these squares introduce an imbalance in the data. The following equation generates the outcome variable Y _i :

where ε∼N(0, 1), i.e. ε follows a standard normal distribution N(0, 1). T _i is the treatment variable.

King and Nielsen (2019) generated 100 simulated data sets. They estimated propensity scores using a logistic regression model (T _i  = X_i1 + X_i2 + ϵ _i ) and performed 1-to-1 greedy matchings without replacement. They also performed 1-to-1 matching using the Mahalanobis distance. They then estimated 512 estimation models, which include X₁,X₂, and all possible combinations of their three quadratic terms and four cubic terms for each data set and each of the 85 calipers. Thus, the simulation is performed in a setting with no unobserved confounders (the assumption of “unconfoundedness” is achieved), but the correct functional form (specification of estimating equation) is not known.

King and Nielsen (2019) aim to illustrate how the combination of model dependence, analyst discretion, and PSM can cause biases in this somewhat extreme setup. Finding out how data are generated is the goal of the research process in observational data analysis. However, using one model as if it were known makes little sense when our understanding of the data-generating process is limited. The analyst therefore faces model dependence, whereby two or more models that fit the data essentially equally generate empirically different causal estimates. The analyst chooses one or more results to publish in response to this ambiguity, but this qualitative decision made without any constraints from a range of possible estimates is probably biased. King and Nielsen (2019) intended to demonstrate how the bias can be worsened in this situation by deleting the worst-matched observation pairs according to the absolute propensity score distance.

King and Nielsen (2019) calculated the variances of the coefficients of T (true value is 2) obtained from this simulation for each number of observations pruned due to the single-nearest-neighbor matching with a caliper. They were then plotted in a figure. The larger the variance, the greater the variation in the estimated coefficient of T, meaning that the change in the estimated coefficient of T depends on the specification of the estimating equation for the outcome variable Y _i . King and Nielsen’s (2019) analysis shows that in MDM matching, the variance of the estimated coefficient of T decreases as the number of removed observations increases whereas in PSM matching, the variance of the estimated coefficient of T continues to increase as the number of removed observations increases beyond a certain point.

To illustrate the model dependence and bias in the causal effect estimates that can result from the analyst selecting the specification of the estimating equation arbitrarily, King and Nielsen (2019) also present a figure plotting the estimated coefficient of T for each number of removed observations. Their assumption in this figure is that researchers choose the largest value among the estimated coefficients from the various models, hoping the causal effect is large. Again, the simulation results show that in MDM matching, the bias of the estimated coefficient of T decreases as the number of removed observations increases whereas in PSM matching, the bias of the estimated coefficient of T continues to increase as the number of removed observations increases beyond a certain point.

King and Nielsen (2019) finally attempted to replicate their results on the real observational data used in published papers. The results again show that the covariate balance worsens as the number of removed observations increases in PSM matching. In contrast, the covariate balance improves as the number of removed observations increases in MDM and CEM.

King and Nielsen’s (2019) findings nevertheless raise three concerns. The first concern relates to the assumption of unconfoundedness. Indeed, in many studies, researchers using observational data make the unconfoundedness assumption, that is, the assumption that potential outcomes are conditionally independent of treatment assignment. However, in reality, it is difficult to assume that unobserved confounding is zero. Heckman et al. (1996, 1998 decomposed the conventional measure of selection bias in observational studies into three components: those corresponding to (1) failure of common support condition (bias arising from regions of nonoverlapping support), (2) failure to weight the data appropriately (bias due to different distributions of observed characteristics), and (3) selection bias precisely defined (residual selection on unobserved characteristics). Matching and weighting estimators can help eliminate the first and second components, but matching and weighting do not necessarily eliminate the third component. In light of the above, it would be more practical and interesting to investigate how and which matching techniques could or could not help reduce the estimated bias of causal effects in a circumstance in which the assumption of unconfoundedness is violated (i.e. in a setting where unobserved confounding factors exist).

The second concern relates to the presentation of the bias in the estimated coefficients of T. In practice, it is highly uncommon, although it may be possible, for estimation to be performed using all possible combinations of quadratic and cubic terms of the covariates. Instead, the estimation will often be conducted after initial selection of a relatively limited number of specifications of the estimating equations following theories and previous studies. As a result, the resulting estimated coefficients are not necessarily the largest among all possible specifications of the estimating equations. Given this point, it would probably be more realistic and practical to plot the average estimated coefficients rather than to plot the maximum estimated coefficients obtained from various models for each number of removed observations.

The third concern relates to their methodological approach. King and Nielsen (2019) investigate whether the covariate balance and bias in causal effect estimates improve or worsen as more observations are removed through matching. However, in empirical research, it is usually a top priority to make sure that the matched sample continues to reflect the research population of interest. The loss of a sizeable number of observations can be problematic in this situation. Therefore, researchers who apply matching to observational data often attempt to improve the covariate balance while maintaining the number of observations as much as feasible. They then seek to reduce the bias in the causal effect estimates due to the improved covariate balance. From the perspective of applied research, examining how the bias in the causal effect estimates changes as the covariate balance varies would be more beneficial.

Regarding the first concern (i.e. the assumption of unconfoundedness), Guo, Fraser, and Chen (2020) conduct simulations in a more realistic setting. The settings Guo, Fraser, and Chen (2020) use are particularly beneficial because they perform simulations with and without observed/unobserved confoundings. (We will describe the detailed simulation setup of Guo, Fraser, and Chen (2020) in the next section.).

For the second concern (i.e. the presentation of the bias in the estimated coefficients of the treatment variable [T]), we also have concerns with the method Guo, Fraser, and Chen (2020) use. Indeed, the mean bias in Guo, Fraser, and Chen (2020) is the difference in real values between the true value of the treatment variable’s coefficient and the estimated coefficients. The upward and downward biases cancel each other out in this calculation; therefore, the true magnitude of the bias is unknown. The mean bias should be calculated using the absolute difference between the true value of the coefficient of the treatment variable and the estimated coefficients to determine the actual size of the bias more accurately. Ripollone et al. (2018) use real observational data for their simulation. Therefore, they could not examine whether PSM would increase bias in estimating treatment effects when PSM successively pruned observations because the true value of the treatment variable’s coefficient was unknown for them.

Finally, for the third concern, Guo, Fraser, and Chen (2020) also have limitations similar to those of King and Nielsen (2019). Guo, Fraser, and Chen (2020) only present the average number of remaining observations after the simulations, the average estimated coefficients of the treatment variable, and their biases. In addition, they perform matchings with only one width of the caliper. Guo, Fraser, and Chen (2020) did not examine the relationship between changes in covariate balance and changes in biases in the estimates of causal effects.

4.1 Relationship Between the Covariate Balance and the Number of Removed Observations

In the case of continuous covariates, the standardized difference is the difference in means between two groups divided by the estimated standard deviation of the same covariate (e.g. Austin 2009). Therefore, standardized differences are comparable for covariates measured on different scales. The following equation defines the absolute standardized difference, d, in the case of continuous covariates:

where x‾T and x‾C are the sample means of the covariate X for the treatment (T) and comparison (C) groups before adjustment (weighting or matching), x‾TA and x‾CA are the sample means of the covariate X for the treatment (T) and comparison (C) groups after adjustment, and sT2, sC2, sTA2, and sCA2 are the corresponding sample variances of the covariate X for both groups.

We derive the absolute standardized difference using the following formula when the covariate is a binary variable:

where P _T , P _TA and P _C , P _CA are the proportions of the binary variables in the treatment and comparison groups, respectively.

First, let us examine whether the simulations in this paper replicate the consequence King and Nielsen (2019) mention, according to which PSM with 1-to-1 greedy matching without replacement worsens the covariate balance beyond a certain point.

In Figure 4, we plot the overall MASD of all covariates after 1-to-1 matching without replacement for each number of observations removed by matching, using the sample generated in Setting 1 with no unobserved confounders and no selection bias. Regarding PSM, we generated propensity scores by logit regression using all covariates except for Z.

Figure 4:

Relationship between MASD and the number of removed observations (1).

Figure 4 suggests that MDM (1-to-1) and CEM (k-to-k) effectively improve covariate balance: both achieved less than 10 % MASD in the matched units. MDM achieved less than 10 % MASD with a mean of 2105 matched units. CEM improves covariate balance to the same extent as MDM but leads to a far greater reduction in units. Indeed, CEM requires pruning of more than 2000 units to achieve an MASD of less than 10 %. In contrast, PSM (1-to-1) only achieves an MASD in the 30 % range with a mean of 1649 matched units. This result presumably occurs because the instrumental variable Z is not used in the calculation of the propensity scores whereas Z is included in the calculation of the overall MASD of all covariates after matching. Consistent with the conjecture, Figure 5 shows that the overall MASD of all covariates after PSM also becomes less than 10 % with a mean of 1067 matched units when the instrumental variable Z is included in the propensity score model.

Figure 5:

Relationship between MASD and the number of removed observations (2).

The phenomenon King and Nielsen (2019) mention is present not just in PSM but also in MDM and CEM. Indeed, as Figures 4 and 5 shows, the covariate balance worsens slightly in all matching methods as the removal of observations by matching advances. Figures 6–9 shows that the same phenomenon occurs even when we use 1-to-5 matching for PSM with and without Z in the propensity score model and j-to-k matching for CEM with and without auto cut points for CEM. These results suggest that matching and weighting under extreme caliper settings do not improve covariate balance, regardless of the matching and weighting methods.

Figure 6:

Relationship between MASD and the number of removed observations (3).

Figure 7:

Relationship between MASD and the number of removed observations (4).

Figure 8:

Relationship between MASD and the number of removed observations (5).

Figure 9:

Relationship between MASD and the number of removed observations (6).

Iacus, King, and Porro (2012), who developed CEM, recommend checking covariate balance with the multivariate imbalance measure (MIM), which they also developed. Therefore, we calculated MIMs after matching and plotted them for each number of removed observations.

Figures 10–15 shows that when they are evaluated using a MIM, there is no substantive difference in the covariate balancing performance of PSM and MDM. In contrast, CEM’s performance is better than those of PSM and MDM. However, a phenomenon King and Nielsen (2019) observe still seems present in the data sets generated by Setting 1, albeit slightly, for all matching methods, even when tested on MIM. The phenomenon can also be observed in the data sets generated by Settings 2 (i.e. selection on observables) and 3 (i.e. selection on unobservables); see Appendix A (Figures A1-1–A1-12; A2-1–A2-12).

Figure 10:

Relationship between MIM and the number of removed observations (1).

Figure 11:

Relationship between MIM and the number of removed observations (2).

Figure 12:

Relationship between MIM and the number of removed observations (3).

Figure 13:

Relationship between MIM and the number of removed observations (4).

Figure 14:

Relationship between MIM and the number of removed observations (5).

Figure 15:

Relationship between MIM and the number of removed observations (6).

The relationship between the covariate balance measures (the MASD and MIM) and the number of removed units after matching seems simply to suggest that implementing matching to the extreme does not result in the most improved covariate balance, regardless of matching methods.

4.2 Relationship Between the Bias in the Treatment Effect Estimate and the Covariate Balance

4.2.1 Setting 1

Next, we examine the relationship between the MAE and RMSE of the estimated coefficient of the treatment variable W and the overall MASDs of the covariates. We examined 34 estimators in our simulations. Table 1 provides a summary of the examined estimators. Tables 2–4 show the simulation results from Settings 1–3, which are presented in ascending order of RMSE. Each error’s effect on RMSE is in proportion to the size of the squared error. Therefore, RMSE can detect larger errors more easily than MAE, in which each error influences in proportion to the absolute value of the error. (Tables B1–B3 in Appendix B present the mean, min, max, and the 5th, 25th, 50th, 75th, and 95th percentile values of MAE of treatment effects.)

Table 1:

Summary of the examined estimators.

ID	Estimation	Estimates	Propensity score (PS) calculation	Matching	Calipers/cut-points
1	OLS	Simple difference in means
2	OLS	Two-step (doubly robust estimation)
3	EB	Simple difference in means
4	EB	Two-step (doubly robust estimation)
5	IPW	IPW estimates	Without Z (instrumental variable)
6	IPW	Two-step (doubly robust estimation)	Without Z (instrumental variable)
7	IPW	IPW estimates	With Z (instrumental variable)
8	IPW	Two-step (doubly robust estimation)	With Z (instrumental variable)
9	PSM	PSM estimates	Without Z (instrumental variable)	1-to-1	No caliper
10	PSM	PSM estimates	With Z (instrumental variable)	1-to-1	No caliper
11	PSM	PSM estimates	Without Z (instrumental variable)	1-to-5	No caliper
12	PSM	PSM estimates	With Z (instrumental variable)	1-to-5	No caliper
13	PSM	Two-step (doubly robust estimation)	Without Z (instrumental variable)	1-to-1	0.05, 0.1, 0.15, 0.2, 0.25
14	PSM	Two-step (doubly robust estimation)	With Z (instrumental variable)	1-to-1	0.05, 0.1, 0.15, 0.2, 0.25
15	PSM	Two-step (doubly robust estimation)	Without Z (instrumental variable)	1-to-5	0.05, 0.1, 0.15, 0.2, 0.25
16	PSM	Two-step (doubly robust estimation)	With Z (instrumental variable)	1-to-5	0.05, 0.1, 0.15, 0.2, 0.25
17	MDM	Simple difference in means		1-to-1	0.3, 0.975, 1.65, 2.325, 3
18	MDM	Two-step (doubly robust estimation)		1-to-1	0.3, 0.975, 1.65, 2.325, 3
19	CEM	Simple difference in means		k-to-k	21–25
20	CEM	Two-step (doubly robust estimation)		k-to-k	21–25
21	CEM	Simple difference in means		j-to-k	21–25
22	CEM	Two-step (doubly robust estimation)		j-to-k	21–25
23	CEM	Simple difference in means		k-to-k	Auto (13)
24	CEM	Two-step (doubly robust estimation)		k-to-k	Auto (13)
25	CEM	Simple difference in means		j-to-k	Auto (13)
26	CEM	Two-step (doubly robust estimation)		j-to-k	Auto (13)
27	CEM (weighting)	Simple difference in means		k-to-k	21–25
28	CEM (weighting)	Two-step (doubly robust estimation)		k-to-k	21–25
29	CEM (weighting)	Simple difference in means		j-to-k	21–25
30	CEM (weighting)	Two-step (doubly robust estimation)		j-to-k	21–25
31	CEM (weighting)	Simple difference in means		k-to-k	Auto (13)
32	CEM (weighting)	Two-step (doubly robust estimation)		k-to-k	Auto (13)
33	CEM (weighting)	Simple difference in means		j-to-k	Auto (13)
34	CEM (weighting)	Two-step (doubly robust estimation)		j-to-k	Auto (13)

Table 2:

Simulation results from setting 1.

ID	Estimation	Estimates	PS calculation	Matching	Calipers/cut-points	Mean number of matched units (treated/untreated) [*]	Mean abs. std. diff. (MASD) after adj.	Mean multivariate imbalance measure (MIM) after adj.	Mean abs. error (MAE) of treatment effects	Root mean squared error (RMSE) of treatment effects	Rank
											MASD	MIM	MAE	RMSE
5	IPW	IPW estimates	Without Z				0.4625		0.0323	0.0405	30		1	1
11	PSM	PSM estimates	Without Z	1-to-5	No caliper		0.3563		0.0422	0.0527	26		2	2
16	PSM	Two-step	With Z	1-to-5	0.05, 0.1, 0.15, 0.2, 0.25	1995 (1240/755)	0.4086	0.9758	0.0437	0.0548	29	12	3	3
2	OLS	Parametric					0.5883		0.0446	0.0557	33		4	4
15	PSM	Two-step	Without Z	1-to-5	0.05, 0.1, 0.15, 0.2, 0.25	2293 (1241/1052)	0.5398	0.9753	0.0450	0.0563	32	11	5	5
6	IPW	Two-step	Without Z				0.4625		0.0474	0.0593	30		6	6
14	PSM	Two-step	With Z	1-to-1	0.05, 0.1, 0.15, 0.2, 0.25	1067 (533/533)	0.0548	0.9388	0.0490	0.0614	9	9	7	7
13	PSM	Two-step	Without Z	1-to-1	0.05, 0.1, 0.15, 0.2, 0.25	1649 (825/825)	0.3776	0.9671	0.0497	0.0623	28	10	8	8
9	PSM	PSM estimates	Without Z	1-to-1	No caliper		0.3619		0.0566	0.0713	27		9	9
26	CEM	Two-step		j-to-k	Auto (13)	763 (382/381)	0.1638	0.9243	0.0592	0.0744	14	7	10	10
8	IPW	Two-step	With Z				0.3382		0.0625	0.0790	24		11	11
7	IPW	IPW estimates	With Z				0.3382		0.0644	0.0814	24		12	12
24	CEM	Two-step		k-to-k	Auto (13)	544 (272/272)	0.0373	0.8927	0.0689	0.0859	5	5	14	13
32	CEM (weighting)	Two-step		k-to-k	Auto (13)		0.3045		0.0689	0.0859	22		14	13
34	CEM (weighting)	Two-step		j-to-k	Auto (13)		0.2549		0.0688	0.0863	16		13	15
17	MDM	Diff. in means		1-to-1	0.3, 0.975, 1.65, 2.325, 3	2105 (1053/1053)	0.0612	0.9794	0.0692	0.0868	10	13	16	16
18	MDM	Two-step		1-to-1	0.3, 0.975, 1.65, 2.325, 3	2105 (1053/1053)	0.0612	0.9794	0.0700	0.0879	10	13	17	17
33	CEM (weighting)	Diff. in means		j-to-k	Auto (13)		0.2549		0.0729	0.0912	16		18	18
23	CEM	Diff. in means		k-to-k	Auto (13)	544 (272/272)	0.0373	0.8927	0.0731	0.0913	5	5	19	19
31	CEM (weighting)	Diff. in means		k-to-k	Auto (13)		0.3045		0.0731	0.0913	22		19	19
25	CEM	Diff. in means		j-to-k	Auto (13)	763 (382/381)	0.1638	0.9243	0.1101	0.1277	14	7	24	21
12	PSM	PSM estimates	With Z	1-to-5	No caliper		0.0818		0.1056	0.1327	12		23	22
3	EB	Diff. in means					0.0000		0.1055	0.1332	1		21	23
4	EB	Two-step					0.0000		0.1055	0.1332	1		22	24
10	PSM	PSM estimates	With Z	1-to-1	No caliper		0.1133		0.1503	0.1900	13		25	25
21	CEM	Diff. in means		j-to-k	21–25	115 (58/57)	0.0511	0.5384	0.1578	0.1988	7	3	27	26
22	CEM	Two-step		j-to-k	21–25	115 (58/57)	0.0511	0.5384	0.1577	0.1989	7	3	26	27
19	CEM	Diff. in means		k-to-k	21–25	108 (54/54)	0.0248	0.5113	0.1609	0.2027	3	1	30	28
27	CEM (weighting)	Diff. in means		k-to-k	21–25		0.3000		0.1609	0.2027	20		30	28
20	CEM	Two-step		k-to-k	21–25	108 (54/54)	0.0248	0.5113	0.1608	0.2027	3	1	28	30
28	CEM (weighting)	Two-step		k-to-k	21–25		0.3000		0.1608	0.2027	20		28	30
29	CEM (weighting)	Diff. in means		j-to-k	21–25		0.2936		0.1626	0.2053	18		32	32
30	CEM (weighting)	Two-step		j-to-k	21–25		0.2936		0.1626	0.2054	18		33	33
1	OLS	Diff. in means					0.5883		0.3704	0.3718	33		34	34

1. [*] Due to rounding, the sum of treated and untreated units does not necessarily correspond to the number of matched units.

Table 3:

Simulation results from setting 2.

ID	Estimation	Estimates	PS calculation	Matching	Calipers/cut-points	Mean number of matched units (treated/untreated) [*]	Mean abs. std. diff. (MASD) after adj.	Mean multivariate imbalance measure (MIM) after adj.	Mean abs. error (MAE) of treatment effects	Root mean squared error (RMSE) of treatment effects	Rank
											MASD	MIM	MAE	RMSE
1	OLS	Diff. in means					0.5883		0.0264	0.0329	33		1	1
16	PSM	Two-step	With Z	1-to-5	0.05, 0.1, 0.15, 0.2, 0.25	1995 (1240/755)	0.4086	0.9758	0.0401	0.0502	29	12	2	2
2	OLS	Parametric					0.5883		0.0409	0.0511	33		3	3
15	PSM	Two-step	Without Z	1-to-5	0.05, 0.1, 0.15, 0.2, 0.25	2293 (1241/1052)	0.5398	0.9753	0.0413	0.0516	32	11	4	4
6	IPW	Two-step	Without Z				0.4625		0.0435	0.0543	30		5	5
14	PSM	Two-step	With Z	1-to-1	0.05, 0.1, 0.15, 0.2, 0.25	1067 (533/533)	0.0548	0.9388	0.0449	0.0563	9	9	6	6
13	PSM	Two-step	Without Z	1-to-1	0.05, 0.1, 0.15, 0.2, 0.25	1649 (825/825)	0.3776	0.9671	0.0455	0.0571	28	10	7	7
26	CEM	Two-step		j-to-k	Auto (13)	763 (382/381)	0.1638	0.9243	0.0543	0.0682	14	7	8	8
8	IPW	Two-step	With Z				0.3382		0.0573	0.0724	24		9	9
25	CEM	Diff. in means		j-to-k	Auto (13)	763 (382/381)	0.1638	0.9243	0.0591	0.0737	14	7	10	10
7	IPW	IPW estimates	With Z				0.3382		0.0596	0.0752	24		11	11
17	MDM	Diff. in means		1-to-1	0.3, 0.975, 1.65, 2.325, 3	2105 (1053/1053)	0.0612	0.9794	0.0604	0.0758	10	13	12	12
24	CEM	Two-step		k-to-k	Auto (13)	544 (272/272)	0.0373	0.8927	0.0631	0.0787	5	5	14	13
32	CEM (weighting)	Two-step		k-to-k	Auto (13)		0.3045		0.0631	0.0787	22		14	13
34	CEM (weighting)	Two-step		j-to-k	Auto (13)		0.2549		0.0631	0.0791	16		13	15
33	CEM (weighting)	Diff. in means		j-to-k	Auto (13)		0.2549		0.0640	0.0802	16		16	16
23	CEM	Diff. in means		k-to-k	Auto (13)	544 (272/272)	0.0373	0.8927	0.0642	0.0802	5	5	18	17
31	CEM (weighting)	Diff. in means		k-to-k	Auto (13)		0.3045		0.0642	0.0802	22		18	17
18	MDM	Two-step		1-to-1	0.3, 0.975, 1.65, 2.325, 3	2105 (1053/1053)	0.0612	0.9794	0.0641	0.0806	10	13	17	19
3	EB	Diff. in means					0.0000		0.0967	0.1221	1		20	20
4	EB	Two-step					0.0000		0.0967	0.1221	1		21	21
12	PSM	PSM estimates	With Z	1-to-5	No caliper		0.0818		0.1061	0.1328	12		22	22
22	CEM	Two-step		j-to-k	21–25	115 (58/57)	0.0511	0.5384	0.1445	0.1823	7	3	23	23
21	CEM	Diff. in means		j-to-k	21–25	115 (58/57)	0.0511	0.5384	0.1452	0.1831	7	3	24	24
20	CEM	Two-step		k-to-k	21–25	108 (54/54)	0.0248	0.5113	0.1474	0.1858	3	1	25	25
28	CEM (weighting)	Two-step		k-to-k	21–25		0.3000		0.1474	0.1858	20		25	25
19	CEM	Diff. in means		k-to-k	21–25	108 (54/54)	0.0248	0.5113	0.1476	0.1860	3	1	27	27
27	CEM (weighting)	Diff. in means		k-to-k	21–25		0.3000		0.1476	0.1860	20		27	27
30	CEM (weighting)	Two-step		j-to-k	21–25		0.2936		0.1490	0.1882	18		29	29
29	CEM (weighting)	Diff. in means		j-to-k	21–25		0.2936		0.1492	0.1885	18		30	30
10	PSM	PSM estimates	With Z	1-to-1	No caliper		0.1133		0.1571	0.1984	13		31	31
5	IPW	IPW estimates	Without Z				0.4625		0.4692	0.4707	30		32	32
9	PSM	PSM estimates	Without Z	1-to-1	No caliper		0.3619		0.4772	0.4824	27		33	33
11	PSM	PSM estimates	Without Z	1-to-5	No caliper		0.3563		0.4807	0.4834	26		34	34

1. [*] Due to rounding, the sum of treated and untreated units does not necessarily correspond to the number of matched units.

Table 4:

Simulation results from setting 3.

ID	Estimation	Estimates	PS calculation	Matching	Calipers/cut-points	Mean number of matched units (treated/untreated) [*]	Mean abs. std. diff. (MASD) after adj.	Mean multivariate imbalance measure (MIM) after adj.	Mean abs. error (MAE) of treatment effects	Root mean squared error (RMSE) of treatment effects	Rank
											MASD	MIM	MAE	RMSE
25	CEM	Diff. in means		j-to-k	Auto (13)	768 (382/386)	0.1660	0.9242	0.0810	0.0988	14	7	1	1
5	IPW	IPW estimates	Without Z				0.4599		0.1202	0.1265	30		2	2
11	PSM	PSM estimates	Without Z	1-to-5	No caliper		0.3539		0.1253	0.1350	26		3	3
9	PSM	PSM estimates	Without Z	1-to-1	No caliper		0.3600		0.1253	0.1408	27		4	4
2	OLS	Parametric					0.5883		0.1469	0.1551	33		5	5
15	PSM	Two-step	Without Z	1-to-5	0.05, 0.1, 0.15, 0.2, 0.25	2293 (1241/1052)	0.5360	0.9756	0.1516	0.1602	32	11	6	6
6	IPW	Two-step	Without Z				0.4599		0.1610	0.1701	30		12	7
13	PSM	Two-step	Without Z	1-to-1	0.05, 0.1, 0.15, 0.2, 0.25	1652 (826/826)	0.3740	0.9680	0.1601	0.1705	28	10	11	8
16	PSM	Two-step	With Z	1-to-5	0.05, 0.1, 0.15, 0.2, 0.25	2004 (1240/763)	0.4082	0.9762	0.1622	0.1705	29	12	13	9
23	CEM	Diff. in means		k-to-k	Auto (13)	549 (275/275)	0.0361	0.8918	0.1541	0.1741	5	5	7	10
31	CEM (weighting)	Diff. in means		k-to-k	Auto (13)		0.3098		0.1541	0.1741	22		7	10
33	CEM (weighting)	Diff. in means		j-to-k	Auto (13)		0.2592		0.1554	0.1753	16		9	12
17	MDM	Diff. in means		1-to-1	0.3, 0.975, 1.65, 2.325, 3	2109 (1054/1054)	0.0608	0.9790	0.1590	0.1772	10	13	10	13
26	CEM	Two-step		j-to-k	Auto (13)	768 (382/386)	0.1660	0.9242	0.1768	0.1911	14	7	14	14
14	PSM	Two-step	With Z	1-to-1	0.05, 0.1, 0.15, 0.2, 0.25	1079 (540/540)	0.0543	0.9399	0.1890	0.1984	9	9	18	15
34	CEM (weighting)	Two-step		j-to-k	Auto (13)		0.2592		0.1811	0.1995	16		15	16
24	CEM	Two-step		k-to-k	Auto (13)	549 (275/275)	0.0361	0.8918	0.1816	0.1998	5	5	16	17
32	CEM (weighting)	Two-step		k-to-k	Auto (13)		0.3098		0.1816	0.1998	22		16	17
18	MDM	Two-step		1-to-1	0.3, 0.975, 1.65, 2.325, 3	2109 (1054/1054)	0.0608	0.9790	0.1999	0.2170	10	13	21	19
8	IPW	Two-step	With Z				0.3373		0.2050	0.2185	24		25	20
7	IPW	IPW estimates	With Z				0.3373		0.2047	0.2191	24		24	21
12	PSM	PSM estimates	With Z	1-to-5	No caliper		0.0813		0.1939	0.2251	12		19	22
3	EB	Diff. in means					0.0000		0.2040	0.2339	1		22	23
4	EB	Two-step					0.0000		0.2040	0.2339	1		23	24
21	CEM	Diff. in means		j-to-k	21–25	113 (57/56)	0.0519	0.5259	0.1955	0.2419	7	3	20	25
19	CEM	Diff. in means		k-to-k	21–25	106 (53/53)	0.0246	0.5030	0.2064	0.2542	3	1	26	26
27	CEM (weighting)	Diff. in means		k-to-k	21–25		0.3092		0.2064	0.2542	20		26	26
22	CEM	Two-step		j-to-k	21–25	113 (57/56)	0.0519	0.5259	0.2115	0.2584	7	3	29	28
29	CEM (weighting)	Diff. in means		j-to-k	21–25		0.3006		0.2105	0.2593	18		28	29
20	CEM	Two-step		k-to-k	21–25	106 (53/53)	0.0246	0.5030	0.2126	0.2607	3	1	30	30
28	CEM (weighting)	Two-step		k-to-k	21–25		0.3092		0.2126	0.2607	20		30	30
30	CEM (weighting)	Two-step		j-to-k	21–25		0.3006		0.2166	0.2656	18		32	32
1	OLS	Diff. in means					0.5883		0.2705	0.2724	33		34	33
10	PSM	PSM estimates	With Z	1-to-1	No caliper		0.1101		0.2327	0.2742	13		33	34

1. [*] Due to rounding, the sum of treated and untreated units does not necessarily correspond to the number of matched units.

Table 2 presents the simulation results obtained from Setting 1, where the assumption of unconfoundedness is satisfied. Overall, there is no substantive difference in the ranking results based on MAE and RMSE. IPW and PSM are generally lower in MAE and RMSE than CEM, MDM, and EB. IPW estimates (without Z) (ID 5) achieved the lowest bias, with an average bias of about 6 % in MAE. PSM estimates (without Z) (ID 11) attained the next lowest bias, with an average bias of about 8 % in MAE. However, the two-step procedure, i.e. doubly robust estimation, is not necessarily less effective in reducing bias than IPW/PSM estimates. In fact, Table 2 shows that half of the top 10 in the RMSE ranking are two-step estimates (IDs 16, 15, 6, 14, and 13), except IPW estimates in the first place (ID 5), OLS estimates in the fourth place (ID 2), and PSM estimates in the second and ninth places (IDs 11 and 9). In Setting 1, where the assumption of unconfoundedness is satisfied, the OLS estimation result (ID 2) from unmatched/unweighted samples is equivalent to those obtained from samples applied with PSM, and its estimated bias in terms of RMSE is lower than the estimated bias from samples applied with MDM and CEM.

Regarding the instrumental variable Z, Table 2 suggests that not including Z in the propensity score calculation (e.g. IDs 5, 11, 15, and 6) is, on average, more likely to result in lower bias than including Z in the propensity score calculation (e.g. IDs 16 and 14) in terms of MAE and RMSE. The results are consistent with those of Bhattacharya and Vogt (2007) and Wooldridge (2016). For 1-to-1/k-to-k matching and 1-to-many/j-to-k matching, the latter (e.g. IDs 11, 16, and 15) appears to achieve lower bias in terms of MAE and RMSE than the former (e.g. IDs 14, 13, and 9).

EB (IDs 3 and 4) performs very well in improving covariate balance. After covariate distribution adjustment using EB, the MASD between the treatment and comparison groups is nearly zero. The problem is that despite its covariate-balancing ability, EB does not very effectively reduce bias in terms of MAE and RMSE: estimation with an EB-weighted sample shows a mean bias of about 20 % for MAE in the simple difference in means and the two-step procedures.

MDM performs as well as other matching/weighting methods in improving covariate balance. However, MDM is less effective than IPW, PSM, and CEM regarding bias reduction. Indeed, estimates from samples obtained by MDM (IDs 17 and 18), which have less than 10 % MASD between the treatment and comparison groups, show an average bias of about 14 % when measured by MAE.

CEM (matching), which uses j-to-k match and auto cut-points (13), is next to PSM regarding bias reduction. It shows a mean bias of about 12 % for MAE (ID 26). The simulation results from Setting 1 also suggest that in the case of CEM, weighting is slightly inferior to matching in terms of bias reduction: it results in a mean bias of about 14 % for MAE (e.g. IDs 32 and 34).

Overall, comparing RMSE and MASD rankings in Table 2, it appears that the improvement in covariate balance measured by MASD does not necessarily lead to an improvement in bias as measured by RMSE. A similar trend emerges even when we compare RMSE and MIM rankings, though only with post-matching data because the CEM command does not accept weights calculated with methods other than cem.

Then, what is the relationship between RMSE and MASD in each sample obtained by weighting or matching? Does the improvement in covariate balance measured by MASD again not necessarily lead to an improvement in bias measured by RMSE in each sample? We examine the relationship between changes in covariate balance as measured by MASD and changes in bias as measured by RMSE by plotting them on figures. For simplicity, the figures compare the bias of the simple difference in means and parametric OLS estimators using unweighted and unmatched samples with the bias of estimations using samples applied with EB, IPW, PSM, MDM, and CEM. We created figures for all 32 estimators (except for the simple difference in means and parametric OLS estimators using raw data). The number of figures for the three settings amounts to 96. Therefore, we present the figures with the lowest mean bias for each of EB, IPW, PSM, MDM, and CEM in the main text, and we present the remaining figures in Appendix C (Figures C1-1–C1-32; C2-1–C2-32; C3-1–C3-32).

In Setting 1, the biases of a simple difference in means and OLS estimation on raw data are generally stable regardless of the magnitude of covariate imbalance (see Figures 16 and 20). In contrast, the covariate imbalance is stable at almost zero in the sample obtained with EB weighting, but the bias in terms of RMSE is highly variable and large. The same trend emerges for the two-step procedure, with almost the same size of bias (Appendix C, Figure C1-2).

Figure 16:

The bias in the treatment effect estimate and the covariate balance: ID 1 versus ID 3.

In contrast, Figure 17 shows that in the sample obtained from IPW, the bias in terms of RMSE is generally constant regardless of covariate balance, as in the case of a simple difference in means and OLS estimation on raw data. For IPW, we observe a trend in which the bias increases slightly when the covariate balance is achieved to nearly the fullest extent. The same tendency is observed for the two-step procedure, and the size of the bias is slightly larger (Appendix C, Figure C1-4).

Figure 17:

The bias in the treatment effect estimate and the covariate balance: ID 1 versus ID 5.

Figure 18 suggests that a similar pattern emerges in the PSM. Again, the bias in terms of RMSE is generally constant regardless of the covariate balance in the sample obtained from PSM, but the bias increases slightly when the maximum covariate balance is achieved. The same tendency is observed in other PSM estimates (Appendix C, Figures C1-7–C1-14), but we observe no substantive increase in bias proportional to the improvement in covariate balance for ID 10 and ID 12 (Appendix C, Figures C1-8 and C1-10). Although biases in average PSM estimates for ID 10 and ID 12, obtained when Z (instrumental variable) was included in the propensity score calculation, decrease in proportion to the improvement in covariate balance, the biases are almost twice as large as those obtained from other PSM estimates.

Figure 18:

The bias in the treatment effect estimate and the covariate balance: ID 1 versus ID 11.

In contrast, Figure 19 shows that when MDM is applied, biases in the estimates of the coefficients of W improve in conjunction with the improvement in covariate balance. However, the biases in the estimated coefficients of W obtained from samples applied with MDM are generally higher than those of the estimates obtained from IPW and PSM. The same result is observed for the two-step procedure, with almost the same size of bias (Appendix C, Figure C1-16).

Figure 19:

The bias in the treatment effect estimate and the covariate balance: ID 1 versus ID 17.

Finally, Figure 20 shows that a similar pattern is observed in CEM, too. Here again, the bias in terms of RMSE is generally constant regardless of the covariate balance of the samples obtained from the CEM. Moreover, the bias shows a slight tendency to increase when the covariate balance is maximized. It should be noted, however, that such a trend is not always observed in the case of CEM, depending on the estimation setup. In other words, in k-to-k matching using five sets of cut points (Figures C1-17 and C1-18; IDs 19, 20) and weighting after k-to-k matching using auto cut points (Figures C1-29 and C1-30; IDs 31, 32), there is almost no tendency for the mean bias to increase even when the covariate balance is maximized. Unfortunately, they did not achieve the lowest bias obtained from the CEM.

Figure 20:

The bias in the treatment effect estimate and the covariate balance: ID 2 versus ID 26.

Our examination of the relationship between changes in covariate balance as measured by MASD and changes in bias as measured by RMSE suggests that beyond a certain point, the improvement in covariate balance can no longer help reduce bias in the estimated effects of treatment regardless of the weighting/matching method.

4.2.2 Setting 2

Even when unconfoundedness could be reasonably assumed, there is always the risk of confounding bias in estimates of causal effects when using observational data. In Setting 2, selection is on observables, in which stochastic dependence between W and u occurs because of the correlation between Z and u, leading to selection bias. Because Z is an observable confounder, conditioning on Z controls the dependence between W and u and mitigates the selection bias.

Table 3 presents the simulation results obtained from Setting 2. Again, there is no substantive difference in the ranking results based on MAE and RMSE. Overall, IPW and PSM are lower in MAE and RMSE than CEM, MDM, and EB. The lowest bias was achieved by a simple difference in means from raw data (ID 1), with an average bias of about 5 % in MAE. The OLS estimation results (ID 2) from unmatched/unweighted samples are again equivalent to or better than those obtained from samples applied with EB, IPW, PSM, MDM, and CEM. Biases in the OLS coefficient estimates from unmatched/unweighted samples are lower than the estimated bias from samples applied with EB, IPW, MDM, and CEM.

Here again, the two-step procedure, i.e. doubly robust estimation, is not necessarily less effective in reducing bias than IPW/PSM estimates: we attained the next-lowest bias using PSM estimates (with Z) from the two-step procedure (ID 16), with an average bias of about 8 % in MAE. In addition, Table 3 shows that 7 of the top 10 in the RMSE ranking are two-step estimates (IDs 16, 15, 6, 14, 13, 26, 8). As mentioned above, in Setting 2, selection is on observables: Z is an observable confounder, and conditioning on Z controls the dependence between W and u and mitigates selection bias. Put differently, the unobservable confounding factor is indirectly controlled by controlling for the observable confounders because of the correlation between observable and unobservable confounders (Z and u; Caliendo, Mahlstedt, and Mitnik 2017; Groenwold and Klungel 2015). Thanks to the indirect control of the unobservable confounding factor through conditioning on Z, the two-step procedure in Setting 2 seems no less effective than the IPW/PSM estimates in reducing bias.

Table 3 suggests, in contrast to Table 2, that including Z in the propensity score calculation (e.g. IDs 16, 14, and 8) does not necessarily result in higher bias than not including Z in the propensity score calculation (e.g. IDs 15, 6, and 13) in terms of MAE and RMSE. Apparently, the results are inconsistent with those of Bhattacharya and Vogt (2007) and Wooldridge (2016). However, Bhattacharya and Vogt (2007) suggest that the propensity score method becomes more inconsistent as the instrumental variable’s strength increases, i.e. as the instrument becomes a better predictor of treatment assignment. In Setting 2, the coefficient of Z is 0.5 and the correlation between Z and u is 0.4, suggesting that Z is neither a very strong nor a very weak instrument. In their simulation, Bhattacharya and Vogt (2007) set the coefficients of the instrument as 0.2 and 0.8 for the weak and strong instruments, respectively. Our Z is a moderately weak instrument compared to Bhattacharya and Vogt’s (2007) setting. Due to the moderately weak correlation between Z and u, it is presumed that including Z in the propensity score calculation does not necessarily result in a greater bias in terms of MAE and RMSE than not including Z in the propensity score calculation. For 1-to-1/k-to-k matching and 1-to-many/j-to-k matching, again, the latter (e.g. IDs 16 and 15) appears to achieve lower bias in terms of MAE and RMSE compared to the former (e.g. IDs 14 and 13).

When we compare RMSE and MASD/MIM rankings in Table 3, we again observe that the improvement in covariate balance measured by MASD/MIM does not necessarily lead to an improvement in bias as measured by RMSE. Also, by plotting the relationship between changes in MASD and changes in RMSE for each of EB, IPW, PSM, MDM, and CEM on figures, we find the following results, which are quite similar to those from Setting 1 and are consistent with the findings from the comparison of RMSE and MASD/MIM rankings in Table 3 (see Figures 21–25 for figures with the lowest mean bias for each estimator and Figures C2-1–C2-32 in Appendix C for all figures):

1. the biases of OLS estimators on raw data are generally stable regardless of the magnitude of covariate imbalance in the data (Figures 21 and 22);

2. the covariate imbalance is stable at almost zero in the sample obtained with EB weighting, but the bias in terms of RMSE is highly variable and large (Figure 21);

3. the biases in terms of RMSE are generally constant for IPW, PSM, and CEM regardless of covariate balance in the sample (Figures 22, 23, and 25);

4. all estimators show that the bias tends to increase slightly when the covariate balance is achieved to nearly the fullest extent (Figures 21–25) except for several estimates from EB (Figures C2-1 and C2-2), IPW (Figure C2-3), PSM (Figures C2-8 and C2-10), and CEM (Figures C2-17 and C2-18).

5. even when maximizing the covariate balance leads to little or no increase in the bias, unfortunately, it does not achieve the lowest bias that can be attained from each estimator (Figures C2-3, C2-8, C2-10, C2-17, and C2-18).

Figure 21:

The bias in the treatment effect estimate and the covariate balance: ID 1 versus ID 3.

Figure 22:

The bias in the treatment effect estimate and the covariate balance: ID 2 versus ID 6.

Figure 23:

The bias in the treatment effect estimate and the covariate balance: ID 2 versus ID 16.

Figure 24:

The bias in the treatment effect estimate and the covariate balance: ID 1 versus ID 17.

Figure 25:

The bias in the treatment effect estimate and the covariate balance: ID 2 versus ID 26.

4.2.3 Setting 3

Setting 3 has the issue of selection bias because of the stochastic dependence between W and u caused by the correlation between u and v, an unobservable confounder. Therefore, in Setting 3 (i.e. selection on unobservables), it is impossible to condition on v and mitigate the selection bias. As expected, including unobservable confounding factors causes an increase in the mean absolute bias in the coefficient estimates of W.

Table 4 exhibits the simulation results obtained from Setting 3. Here, we achieved the lowest bias with a simple difference in means from samples applied with CEM (ID 25) though with a larger average bias of about 16 % in MAE compared to those from Settings 1 and 2 due to selection on unobservables. Interestingly, the ranking results based on MAE and RMSE show a notable difference: in the MAE ranking, we observed a simple difference in means based on samples obtained from CEM and MDM (IDs 23, 31, 33, and 17) in the 7th to 10th places, but in the RMSE ranking, they dropped to the 10th to 13th places. Instead, doubly robust estimators based on IPW and PSM samples (IDs 6, 13, and 16) moved up to 7th–9th in the RMSE ranking. The results in Table 4 suggest that in Setting 3, estimation with samples obtained by CEM and MDM tends to produce larger outlier errors.

The OLS estimation results (ID 2) from unmatched/unweighted samples are again equivalent to those obtained from samples applied with PSM. Biases in the OLS coefficient estimates from unmatched/unweighted samples are lower than the estimated bias from samples applied with EB and MDM.

Regarding the two-step procedure, doubly robust estimators are not among the top four in the RMSE ranking. In fact, in Setting 3, the two-step procedure does not seem as effective as a difference in means after CEM (ID 25) and IPW/PSM estimates (IDs 5, 11, and 9) in reducing bias. These results seem consistent with the fact that the two-step procedure cannot be doubly robust in Setting 3, i.e. selection on unobservables. Indeed, in Setting 3, the propensity score model is incorrect, and the parametric model is misspecified.

Table 4 suggests that not including Z in the propensity score calculation (e.g. IDs 5, 11, 9, 15, 6, and 13) results in lower bias than including Z in the propensity score calculation (e.g. ID 16) in terms of MAE and RMSE, which is consistent with Bhattacharya and Vogt (2007) and Wooldridge (2016). Furthermore, for 1-to-1/k-to-k matching and 1-to-many/j-to-k matching, again, the latter (e.g. IDs 25 and 11) seems to achieve lower bias in terms of MAE and RMSE than the former (e.g. IDs 9 and 13).

Comparing RMSE and MASD/MIM rankings in Table 4, here again, we find that the improvement in covariate balance measured using MASD/MIM does not necessarily lead to an improvement in bias as measured using RMSE.

The relationship between changes in MASD and changes in RMSE for each of EB, IPW, PSM, MDM, and CEM demonstrates the following results (see Figures 26–30 for figures with the lowest mean bias for each estimator and Figures C3-1–C3-32 in Appendix C for all figures):

1. the biases of OLS estimators on raw data are generally stable regardless of the magnitude of covariate imbalance in the data (Figures 26 and C3-2);

2. the covariate imbalance is stable at almost zero in the sample obtained with EB weighting, but the bias in terms of RMSE is highly variable and large (Figures 26 and C3-2);

3. the biases in terms of RMSE are generally constant for IPW, PSM, and CEM regardless of covariate balance in the sample (Figures 27, 28, and 30);

4. a simple difference in means after MDM (Figure 29) suggests that the bias in terms of RMSE increases as the covariate balance improves whereas the bias in the estimates from the two-step procedure (Figure C3-16) is generally constant;

5. all estimators show that the bias tends to increase slightly when the covariate balance is achieved to nearly the fullest extent (Figures 26–30) except for several estimates from EB (Figures C3-1 and C3-2), PSM (Figure C3-10), MDM (Figure C3-16), and CEM (Figures C3-22, C3-24, C3-29 and C3-30);

6. however, again, even when maximizing the covariate balance leads to little or no increase in the ATET bias, it does not achieve the lowest bias that can be attained from each estimator.

Figure 26:

The bias in the treatment effect estimate and the covariate balance: ID 1 versus ID 3.

Figure 27:

The bias in the treatment effect estimate and the covariate balance: ID 1 versus ID 5.

Figure 28:

The bias in the treatment effect estimate and the covariate balance: ID 1 versus ID 11.

Figure 29:

The bias in the treatment effect estimate and the covariate balance: ID 1 versus ID 17.

Figure 30:

The bias in the treatment effect estimate and the covariate balance: ID 1 versus ID 25.

These findings are consistent with those from Settings 1 and 2 and MASD/MIM rankings in Table 4.