Using Propensity Score Methods to Approximate Factorial Experimental Designs to Analyze the Relationship Between Two Variables and an Outcome

Abstract

Researchers have become increasingly interested in programs’ main and interaction effects of two variables (A and B, e.g., two treatment variables or one treatment variable and one moderator) on outcomes. A challenge for estimating main and interaction effects is to eliminate selection bias across A-by-B groups. I introduce Rubin’s causal model to approximate factorial experimental designs for studies with partial randomization and nonrandomization. I apply a Monte Carlo simulation to evaluate several propensity score applications. The findings suggest the following two applications for reducing bias and mean square error of parameter estimates when analyzing the relationship of two variables and an outcome: (a) inverse of propensity score weighting based on one multinomial propensity score model and (b) subclassification based on two binary propensity score models. As a demonstration, I examine whether the effects of the Head Start program, compared to other center-based care, for improving children’s reading achievement vary by child care quality.

Keywords

causal inference factorial experimental designs propensity score

Introduction

Researchers and policy makers have become increasingly interested in both the main and interaction effects of multiple variables in program evaluations. One might want to know not only the main effects of a new math curriculum (A), for example, or a teacher professional training program (B) on students’ math achievement (Y) but also whether the effects of the new math curriculum taught by teachers who also received the teacher professional training are better than the effects of the curriculum taught by teachers who did not receive the training. As another example, it might be valuable to know not only the effects of Head Start (A) on child outcomes in general but also the influence of child care quality (B).

These two examples illustrate scenarios that concern the interaction of two independent treatment variables, (A × B); that is, the effect of one treatment variable also depends on the level of another treatment variable. In the first case, teacher professional development (B) might influence the effect of the math curriculum (A); in the second case, it is quality (B) that might influence the effect of Head Start (A). Or, stated another way, in the first case, the effect of teacher professional development (B) might vary by the curriculum condition (A, new curriculum or not); and, in the second case, the effect of child care quality (B) might differ by the types of care (A, Head Start or not). These are situations in which two independent treatment variables are malleable. One can imagine changing how much professional development teachers receive or trying to improve a program’s quality. There is also a scenario in which these kinds of variables, that is, ones that influence treatment effects, cannot be manipulated, but we might be interested in their influences, nevertheless. For example, in a study that examines how the effects of Head Start (A) vary by race/ethnicity (B), exposure to Head Start can be manipulated, while race/ethnicity cannot. In these situations, we consider these nonmalleable variables, for example, race/ethnicity, to be moderators. A moderator affects the direction and/or magnitude of the relationship between the treatment variable and the outcome variable, and the moderation effects are usually tested through an examination of the interaction term (A × B) of the treatment variable and the moderator (Baron & Kenny, 1986). Note that some researchers (e.g., Kraemer, Kiernan, Essex, & Kupfer, 2008; Kraemer, Stice, Kazdin, Offord, & Kupfer, 2001) suggest that the moderators should precede the treatment variables, however, the Baron and Kenny’s approach (e.g., Baron & Kenny, 1986; Frazier, Tix, & Barron, 2004) does not require this criterion for moderation analysis, which means that the moderators could be posttreatment variables, for example, the program-related subgroup variables (e.g., Hill, Brooks-Gunn, & Waldfogel, 2003; Lochman, Boxmeyer, Powell, Roth, & Windle, 2006; Peck, 2003; Schochet & Burghardt, 2007). This article adopts the Baron and Kenny’s approach for moderator analysis.

The above examples illustrate two commonly seen situations in which two variables, that is, two independent treatment variables, and the main treatment indicator and a moderator of some sort, work together to affect an outcome. Another commonly seen situation is mediation, whereby the treatment variable, which occurs earlier, changes the outcome through changing the mediator, which occurs later (Baron & Kenny, 1986; Kraemer et al., 2001). Detecting a mediation effect is of interest in evaluation research, and researchers have proposed various approaches for causal mediation analysis (Hong & Nomi, 2012; Imai, Keele, & Tingley, 2010). Nevertheless, this article does not concern causal mediation analysis but focuses on the analysis of main and interaction effects of two independent treatment variables and the moderation effect under partial randomization or nonrandomization.

The purpose of this study, based on the counterfactual model (Holland, 1986; Rosenbaum, 2002; Rubin, 1974), is to evaluate several propensity score methods by using a Monte Carlo simulation to compare several propensity score methods and to identify appropriate approaches to reduce bias and mean square error (MSE) of parameter estimates of the main and interaction effects of two variables. In the rest of the first section, I demonstrate the challenge in analyzing the effects of two treatment variables and moderation effects.

The biggest challenge in analyzing two treatment variables and moderation effects is the ability to produce unbiased estimates of both the main and interaction effects of two variables. In the situation in which the two variables/factors (A and B) are both dichotomous and malleable, ideally, we can use a 2 × 2 factorial experimental design to estimate the main and interaction effects of A and B. Figure 1 (Shadish, Cook, & Campbell, 2002, p. 264) demonstrates a 2 × 2 factorial experimental design, where subjects are randomly assigned to four groups: (1) Level 1 of Factor A and Level 1 of Factor B (A1B1), (2) Level 1 of Factor A and Level 2 of Factor B (A1B2), (3) Level 2 of Factor A and Level 1 of Factor B (A2B1), and (4) Level 2 of Factor A and Level 2 of Factor B (A2B2). Due to randomization, the main and interaction effects of two factors can be unbiasedly estimated. However, when random assignment of both variables is not feasible, which includes situations for which only one or none of the two variables concerns random assignment, the students in the four cells (groups) in Figure 1 may systematically differ from each other and produce biased estimates of both the main and interaction effects.

Figure 1.

2 × 2 Factorial design. Adapted from figure 8.1 by Shadish, Cook, and Campbell (2002, p. 264).

To illustrate that partial randomization (only one variable is randomly assigned) may produce a biased estimate of the interaction effect, consider a hypothetical example (Table 1). The purpose of the study in the example is to examine the effects of two variables (treatment status and race/ethnicity) on math achievement. This example involves partial randomization, in which the subjects were randomly assigned to the treatment condition (e.g., new math curriculum), while race/ethnicity (White vs. non-White) cannot be manipulated or randomly assigned. The total sample included 16 participants (8 White and 8 non-White). The White sample had higher socioeconomic status (SES) than did the non-White sample, and the average SES for the total sample is 0.

Table 1.

Data for a Hypothetical Example of a Partial Randomization.

Subject	Treatment status (A)	White (B)	SES	Math (Y)	Matched sample^a
1	0	0	−3	94
2	0	0	−2	96
3	0	0	0	100	*
4	0	0	1	102	*
5	0	1	0	110	*
6	0	1	1	112	*
7	0	1	1	112
8	0	1	2	114
9	1	0	−3	96
10	1	0	−2	99
11	1	0	0	105	*
12	1	0	1	108	*
13	1	1	0	115	*
14	1	1	1	118	*
15	1	1	1	118
16	1	1	2	121

Note. SES = socioeconomic status. ^aMatched sample is indicated by *. Matching on SES across four groups (Treatment × White).

Eight subjects were randomly assigned to the treatment condition. Math achievement was collected after a 1-year intervention. The true productivity function of the math outcome was assumed to be the following: Math = 100 + 10 × White + 2 × SES + 5 × Treatment + 1 × SES × Treatment + 0 × White × Treatment. In this equation, the average treatment effect is assumed to be 5 points and the moderation effect of SES is assumed to be 1 point, but there is no moderation effect of race/ethnicity, that is, the treatment effect did not differ between the White and the non-White participants. One advantage of randomized experiments is that the treatment effect estimates do not rely on statistical models and can be directly estimated by the treatment control mean difference.

Table 2 presents the means and treatment effects using the full sample for Table 1. Note that the overall treatment effect is 5 points, which is an unbiased estimate. For subgroup analysis in regard to race/ethnicity, the treatment effects for non-White and White are 4 and 6 points, respectively. Due to the random assignment of participants to the treatment condition, the treatment effects for subgroups are unbiased estimates, as well. However, the correct interpretation is that the treatment effect for the non-White (lower SES) group is 4 points and, for the White (higher SES) people, 6 points. It would not be accurate, however, if one claimed that the moderation effect of race is 2 points or that the treatment effect is 2 points higher for Whites than for non-Whites. The 2-point treatment effect difference is not due to race/ethnicity itself but, rather, due to SES. It is because the SES of Whites is 2 points higher than that of non-Whites that the treatment effect for Whites appeared 2 points higher than that for non-Whites. If race/ethnicity could be manipulated and participants randomly assigned to White and non-White groups, then there would be no SES difference between the race/ethnicity subgroups, and the treatment effect difference between Whites and non-Whites (i.e., moderation effect) could be estimated unbiasedly. This example illustrates that, in analyzing moderation effects under partial randomization, it is possible to get biased moderation effect estimates when there are critical differences between the subgroups (SES difference between Whites and non-Whites in this example), even if the participants were randomly assigned to treatment conditions.¹

Table 2.

Means and Treatment Effects for the Full Sample in Table 1.

	Control (0)	Treatment (1)	Treatment effect
Non-White (0)	98	102	4
White (1)	112	118	6
Total sample	105	110	5

Note. N = 16.

For the above simple example, we can make the baseline variable (SES) equivalent across four groups (non-White control, non-White treatment, White control, and White treatment) by matching on SES. The participants denoted by the asterisk (*) in Table 1 consist of the matched sample, in which each of the four groups includes two participants with SES equal to 0 and 1. The means and treatment effects for the matched sample (N = 8) are presented in Table 3. The overall treatment effect for this sample is 5.5 points, which is an unbiased estimate, and the treatment effects for two subgroups are both 5.5 points, which also is an unbiased estimate. The moderation effect, that is, the treatment effect difference between the Whites and non-Whites is 0, which is an unbiased estimate. This simple matching analysis illustrates that it is essential to keep the baseline variables equivalent across all four groups to estimate the interaction effect under partial randomization.

Table 3.

Means and Treatment Effects for the Matched Sample in Table 1.

	Control (0)	Treatment (1)	Treatment effect
Non-White (0)	101	106.5	5.5
White (1)	111	116.5	5.5
Total sample	106	111.5	5.5

Note. N = 8.

The above example also indicates that it is possible to achieve unbiased estimates of the interaction effects under partial randomization when analyzing the effects of two variables. Within the past decade or so, statistical methods for analyzing the effects of two variables have been developed. In experimental subgroup analysis, researchers have varyingly extended propensity score methods (Rosenbaum & Rubin, 1983), in an attempt to minimize selection bias, and have identified well-matched subsets of treatment and control groups within which to estimate subgroup treatment effects (Hill et al., 2003; Lochman et al., 2006; Peck, 2003; Schochet & Burghardt, 2007). When the subgroup variable is binary, the basic practice is to match participants between the treatment and control groups within each of the two subgroups. These separate matches may make the baselines equivalent for the treatment and control groups for each subgroup; thus, they may produce unbiased estimates of treatment effects for that subgroup (for further discussion about matching within subgroups, see Peck, 2013). However, without making an effort to make the baselines equivalent among all four groups (Figure 1), we cannot produce an unbiased estimate of the interaction effect (i.e., whether the difference in treatment effects between two subgroups is due solely to the intervention and subgroup factors), as illustrated in Tables 1 and 2.

Other researchers, for example, Imai and van Dyk (2004), have generalized propensity score methods to study the main effects of two treatment variables using subclassification. However, we still lack a systematic study to examine all of the methods noted previously for analyzing the main and interaction effects of two variables on the outcome.

In the second section, I present Rubin’s causal model to approximate factorial design. In the third section, I review and identify applicable propensity score methods. Then, in the fourth section, I present the Monte Carlo simulation used to evaluate these methods. I use a case study, in the fifth section, to demonstrate how to use the suggested propensity score methods to conduct causal analysis. Finally, the sixth section concludes the article by providing a discussion of applications.

Causal Framework for Analyzing the Effects of Two Treatment Variables or Moderator Effects on the Outcome

Using the counterfactual model (Holland, 1986; Rosenbaum, 2002; Rubin, 1974), I present the potential outcomes for individuals under the situation that both treatment A and treatment B (or moderator) are binary variables (Figure 2). This is analogous to a 2 × 2 factorial experimental design (Figure 1). The main effect of treatment A is the marginal mean difference Y(1,.) − Y(0,.); the main effect of treatment B or moderator is the marginal mean difference, Y(.,1) − Y(.,0). The interaction effect of treatment A and treatment B, or moderator effect (i.e., the difference in treatment effect between two levels of moderator), is [Y(1,1) − Y(0,1)] − [Y(1,0) − Y(0,0)].

Figure 2.

Potential outcomes of two variables.

Alternatively, this counterfactual model can be illustrated in Figure 3. Note that, in the figure, this path diagram looks similar to the one that Raudenbush (2011, p. 22) used to illustrate the potential outcome of mediation analysis (Case 2: Treatment → Mediating Treatment → Y), but there are notable differences. In Raudenbush’s diagram, there is a causal link (arrow) from treatment A to treatment B, which, in turn, causes the outcome. In this diagram, however, no causal relationship between treatment A and treatment B (or moderator) is specified, that is, treatment A and treatment B (or moderator) are independent, there is no arrow from treatment A to treatment B (or moderator), and it is treatment A and treatment B (or moderator) together that cause the outcome (Y). Although the mechanism is different, we can use a similar counterfactual model to make a causal inference.

Figure 3.

Potential outcomes and paths of two variables.

Applicable Propensity Score Methods to Approximate Factorial Experimental Designs

In this section, I extend Imai and van Dyk’s (2004) approach to analyzing the effects of two binary variables on the outcome, using a stratification method based on the two propensity score functions. Then, I introduce other potential propensity score methods for analyzing two binary variables.

Imai and van Dyk (2004) used Rubin’s causal model to generalize a propensity score application for the analysis of the effects of two continuous treatment variables on the outcome. To do so, they relied on two standard assumptions (the stable unit treatment value assumption; Rubin, 1980, 1990) and the strong ignorability of the treatment assignment assumption (Rosenbaum & Rubin, 1983) as well as a third assumption of a uniquely parameterized propensity function. The third assumption implies that the propensity score function can be effectively summarized by the parameter that can be calculated from Gaussian regression using the covariates (see Imai & van Dyk, 2004, for more details about this assumption). The two propensity score functions for the two continuous treatment variables are estimated based on two independent Gaussian linear regression models. The two estimated propensity score functions are used to subclassify the data into several groups. Figure 4 illustrates the 3 × 3 subclassification based on two propensity score functions (Imai & van Dyk, 2004).

Figure 4.

3 × 3 Subclassification based on two propensity score functions. Adapted from figure 4 by Imai and van Dyk (2004, p. 861). Data are subclassified into three groups (lower third, middle third, and upper third) based on each of the two propensity score functions, respectively. Each cell of the 3 × 3 table represents a group based on two propensity score functions jointly.

Data are classified into three groups (lower third, middle third, and upper third) based on each of the two propensity score functions. Each cell of the 3 × 3 table represents a group based on two propensity score functions jointly. The overall treatment effect point estimate and its standard error estimate are a weighted average across nine groups (Imbens & Rubin, 2009). In addition to a 3 × 3 subclassification, Imai and van Dyk (2004) presented simulation results for a 2 × 2 and a 4 × 4 subclassification. In general, more groups produce less bias, as in the analysis of a single treatment variable (Lunceford & Davidian, 2004).² Note that the subclassification method based on two propensity score functions could be used for analyses in which both two variables are binary or one variable is binary and another variable is continuous.

Other Propensity Score Methods for Analyzing Multiple Groups

Imai and van Dyk’s (2004) subclassification requires the estimation of two independent propensity score functions based on treatment A and treatment B (or moderator), respectively. There are, however, alternative approaches to using propensity scores for analyzing multiple groups. First, we can collapse two dimensions (2 × 2) of design to one dimension (4 × 1) of design and use the propensity score methods that are appropriate for one-dimensional multiple groups analysis to analyze the effects of two variables. We convert a 2 × 2 design to a 4 × 1 design, that is, a design that has one treatment variable with four levels. For example, Figure 1 presents a 2 × 2 factorial design in which there are two factors (Factors A and B), and each factor has two levels. The conventional outcome analysis involves the use of the two factors and their interaction term as predictors in the statistical model. We can combine these two factors, however, to create one new variable with four levels (A1B1, A1B2, A2B1, and A2B2) and use this new variable as the predictor in the statistical model. By doing so, we simplify our analysis of two predictors to the analysis of one predictor. This is particularly helpful when we need to estimate the propensity scores of the predictors.

Second, the propensity scores can be estimated using a multinomial logistic regression model, for example, Group A1B1 in Figure 1 serves as the reference group, and all the other groups are compared with this reference group. The propensity scores for a participant are the probabilities of being in a particular group, given the covariates.

Third, various propensity score methods can be applied to estimate the group differences, which can be easily converted to the estimates of the main and interaction effects of two variables. These propensity score applications include subclassification (Rosenbaum & Rubin, 1984), optimal nonbipartite matching (Lu, Greevy, Xu, & Beck, 2011; Lu & Rosenbaum, 2004), inverse of propensity score weighting (Imbens, 2000), and propensity score as covariates.

Monte Carlo Simulation

The previous sections provided a discussion of various applications of propensity score methods in analyzing the effects of two variables on the outcome. These applications have multiple dimensions: (a) estimates of propensity scores (using one multinomial propensity score model or two binary propensity score models); (b) application of propensity scores (subclassification, weighting, matching, or as a covariate); and (c) covariate adjustment (using the covariates that are included in the model to estimate propensity scores to estimate the effects of two variables on the outcome). I evaluate these applications of propensity score methods by conducting a Monte Carlo simulation for two scenarios: (a) partial randomization, which concerns random assignment of treatment A, while treatment B (or moderator) is not randomly assigned; and (b) nonrandomization, that is, neither treatment A nor treatment B (or moderator) concerns random assignment, and treatment A and treatment B (or moderator) are artificially correlated through the other variables.

One example of a partial randomization can be seen in a design to study the main and interaction effects of Head Start (treatment A) and child care quality (treatment B). Students are randomly assigned to a treatment group (Head Start) or control group; however, they are not randomly assigned to child care quality groups. In nonrandomization, the students are neither randomly assigned to a treatment group (Head Start) or control group nor randomly assigned to child care quality groups.

In the Monte Carlo simulation, I create the data with known effects of the two variables under different assumptions. I then use different propensity score methods to estimate the main and interaction effects of two variables through an analysis of the simulated data. By comparing the effect estimates with the true effects, we can determine which propensity score applications work better to approximate factorial experiments to produce less-biased effect estimates. To save space, I put the models that produce data for the simulation in Appendix A. The following sections focus on the approaches to analysis and results of the simulation.

Approaches to Analyzing the Main and Interaction Effects of Two Variables on the Outcome

Partial randomization

I first analyze the full sample using an ordinary least square (OLS) model (Model 1), where W ₁ and W ₂ are two covariates, A is an indicator of treatment A, B is an indicator of treatment B (or moderator), and β _A , β _B , and β _A _{× B} are the main and interaction effects of variables A and B. Model 1 is a misspecified model because it omits the cubic term of a covariate (W ₁) as compared to the model used to produce the data (Model A1 in Appendix A), which simulates the reality that researchers might not know the correct model. These results serve as a reference for comparison with propensity score methods.

Y_{i} = β_{0} + β_{1} {(W_{1})}_{i} + β_{2} {(W_{2})}_{i} + β_{A} A_{i} + β_{B} B_{i} + β_{A \times B} A_{i} B_{i} + e_{i} .

I then apply propensity score methods to analyze the data. Table 4 presents various propensity score applications. First, the multinomial logistic regression model is run to estimate the generalized propensity score (Imbens, 2000) using covariates W ₁ and W ₂. The generalized propensity score is the conditional probability of receiving treatment t, given pretreatment covariate W, that is, $r (t, W) = p r (T = t | W)$ , where W includes W ₁ and W ₂. Thus, each participant will have four generalized propensity scores associated with four groups, Y(0,0), Y(0,1), Y(1,0), and Y(1,1) (see Figure 2.

Table 4.

Propensity Score Applications in a Monte Carlo Simulation Under Partial Randomization.

Propensity score model	Data	Application of propensity score	Covariate adjustment^a	Analysis model
One multinomial propensity score model	1. Full sample	Subclassification^b	No	ANOVA (subclassification)
		Subclassification^b	Yes	OLS (subclassification)
		Weighting^c	No	ANOVA (weighting)
		Weighting^c	Yes	OLS (weighting)
		As covariate^b	No	ANCOVA (as covariate)
		As covariate^b	Yes	OLS (as covariate)
	2. Common support sample (overlap among four groups)	Subclassification^b	No	ANOVA (subclassification)
		Subclassification^b	Yes	OLS (subclassification)
		Weighting^c	No	ANOVA (weighting)
		Weighting^c	Yes	OLS (weighting)
	3. Overlap sample of matching one group with all other three groups using Data Item 2 above	Matching^b	No	ANOVA
		Matching^b	Yes	OLS
Two binary propensity score models	4. Common support sample (overlap among four groups)	Subclassification^d	No	ANOVA (subclassification)
Two binary propensity score models	4. Common support sample (overlap among four groups)	Subclassification^d	Yes	OLS (subclassification)

Note. ANOVA = analysis of variance; OLS = ordinary least square; ANCOVA = analysis of covariance. ^aCovariate adjustment is to adjust for covariates W ₁ and W ₂. ^bSubclassification, matching, and covariate are all based on the probabilities of being in the reference group. ^cWeighting is by the inverse of the probabilities of being the group of treatment that participants actually received. ^dSubclassification and matching are based on two propensity score functions.

The generalized propensity score is estimated using the multinomial logistic regression as given subsequently, where the reference group (e.g., Y(0,0) in Figure 2) is coded as 4:

log i t {r (t, W)} = ln (\frac{P (T = t | W_{1}, W_{2})}{P (T = 4 | W_{1}, W_{2})}) = α_{0 t} + α_{1 t} W_{1} + α_{2 t} W_{2} .

The distributions of the generalized propensity scores of being in the reference group (e.g., the probabilities of being in Group Y(0,0) in Figure 2) for the four groups are examined for overlap. A sample with common support can be obtained by excluding participants whose probabilities of being in the reference group are larger than the minimum value of four maximum probabilities or smaller than the maximum value of the four minimum probabilities (called “common support sample,” Item 2 under Data in Table 4). The formal set notation of the common support sample (C) can be expressed as follows:

C = \{\begin{matrix} i | p r_{i}^{k} (T = 4) \in [L, U], \\ L = M a x (M i n (p r_{i}^{1} (T = 4)), M i n (p r_{i}^{2} (T = 4)), M i n (p r_{i}^{3} (T = 4)), M i n (p r_{i}^{4} (T = 4))), \\ U = M i n (M a x (p r_{i}^{1} (T = 4)), M a x (p r_{i}^{2} (T = 4)), M a x (p r_{i}^{3} (T = 4)), M a x (p r_{i}^{4} (T = 4))), \\ k = 1, 2, 3, 4 \end{matrix}\},

where $p r_{i}^{k} (T = 4)$ is the propensity score of being in Group 4 for individual i, whose actual group is k, k = 1, 2, 3, and 4. L and U are the lower and upper bounds of the propensity scores in the common support sample, L is the maximum value of the minimum propensity scores in four actual groups, and U is the minimum value of the maximum propensity scores in the four actual groups. Use of the common support sample is intended to reduce estimate bias and to avoid extreme values (close to 0 or 1) of the propensity score that can cause very large variance.³

This common support sample is further used for the nearest neighbor one-to-one matching by matching the reference group with all the other three groups, respectively, based on the probabilities of being in the reference group. I then have three matched samples, and each matched sample contains one sample of the reference group. I choose the common sample of being in the reference group among the three matched samples. The final matching sample includes the common sample of being in the reference group and its matching sample in all other three groups (Item 3 under Data in Table 4).

For the full sample and common support sample, to make results comparable with a 3 × 3 subclassification based on two propensity score functions, I subclassify data into nine groups based on the probabilities of being in the reference group and obtain the weighted average parameter estimates across nine groups. I also use the inverse of the generalized propensity score as weight (Imbens, 2000) to analyze data using the full sample and common supported sample, respectively. The weight is $\frac{1}{r (T_{i} = t, W)}$ , where T_i denotes the treatment that participant i actually received. Further, using data from the full sample, I use the probabilities of being in the reference group as covariates to estimate the effects of two variables. All of the estimates of the effects of two variables are conducted using OLS (Model 1) and analysis of variance (ANOVA; Model 3), that is, the model without a covariate.

Y_{i} = β_{0} + β_{A} A_{i} + β_{B} B_{i} + β_{A \times B} A_{i} B_{i} + e_{i} .

Finally, two independent binary logistic regression models, based on covariates W ₁ and W ₂, are used to estimate the propensity scores for treatment A and treatment B (or moderator), respectively. The common support sample (Item 4 under Data in Table 4) is obtained by including the overlap sample among the four groups based on propensity scores for variables A and B through the following steps: (a) the first common support sample is the overlap sample between two levels of treatment A based on the propensity scores of treatment A, (b) the second common support sample is the overlap sample between two levels of treatment B based on the propensity scores of treatment B, and (c) the final common support sample is the overlap sample between the first common support sample and second common support sample. As illustrated in Figure 4, I apply a 3 × 3 subclassification based on two propensity score functions and estimate the weighted treatment effect across nine groups (Imai & van Dyk, 2004). Both the OLS (Model 1) and ANOVA models are used for data analysis.

Non-randomization

I conduct analyses similar to those in the partial experiment scenario. I first analyze the full sample using an OLS model (Model 4), where β _A , β _B , and β _A _{× B} are the main and interaction effects of variables A and B. Similar to Model 1, Model 4 is a misspecified model because it omits a cubic term of a covariate and two interaction terms of covariates and variables A and B (see Expression A2 in Appendix A), which simulates the reality that researchers might not know the correct model. These results serve as a reference for a comparison with propensity score methods.

Y_{i} = β_{0} + β_{1} {(W_{1})}_{i} + β_{2} {(W_{2})}_{i} + β_{3} {(W_{3})}_{i} + β_{A} A_{i} + β_{B} B_{i} + β_{A \times B} A_{i} B_{i} + e_{i} .

I then apply propensity score methods to analyze the data. Table 5 presents various propensity score applications. Both the OLS (Model 4) and ANOVA models are used for data analysis.

Table 5.

Propensity Score Applications in a Monte Carlo Simulation Under Nonrandomization.

Propensity score model	Data	Application of propensity score	Covariate adjustment^a	Analysis model
One multinomial propensity score model	1. Common support sample (overlap among four groups)	Subclassification^b	No	ANOVA (subclassification)
		Subclassification^b	Yes	OLS (subclassification)
		Weighting^c	No	ANOVA (weighting)
		Weighting^c	Yes	OLS (weighting)
	2. Overlap sample of matching one group with all other three groups using Data Item 1 above	Matching^d	No	ANOVA
		Matching^d	Yes	OLS
Two binary propensity score models	3. Common support sample (overlap among four groups)	Subclassification^d	No	ANOVA (subclassification)
Two binary propensity score models	3. Common support sample (overlap among four groups)	Subclassification^d	Yes	OLS (subclassification)

Note. ANOVA = analysis of variance; OLS = ordinary least square. ^aCovariate adjustment is to adjust for covariates W ₁, W ₂, and W ₃. ^bSubclassification, matching, and covariates are all based on the probabilities of being in the reference group. ^cWeighting is by the inverse of the probabilities of being the group of treatment that participants actually received. ^dSubclassification and matching are based on two propensity score functions.

I replicate 1,000 times to produce and analyze data. The estimates of bias, percentage bias, and MSE of the parameters (θ) can be calculated as follows:

\overset{\land}{B i a s} (\hat{θ}) = \frac{1}{K} \sum_{k = 1}^{K} {\hat{θ}}_{k} - θ, w h e r e k = 1, 2, \dots, K . K i s t h e n u m b e r o f v a l i d r e p l i c a t i o n s .

Percentage bias is defined as: $(\frac{\overset{\land}{B i a s} (\hat{θ})}{θ}) \times 100$

\overset{\land}{M S E} (\hat{θ}) = \frac{1}{K} \sum_{k = 1}^{K} {({\hat{θ}}_{k} - θ)}^{2} = \overset{\land}{var} (\hat{θ}) + {(\overset{\land}{B i a s} (\hat{θ}))}^{2} .

Bias is a measure of the difference between the effect estimate and the true effect. A positive sign of bias indicates an overestimate of the effect, and a negative sign of bias indicates an underestimate of the effect. The percentage bias indicates the percentage of bias in terms of the true effect. MSE is a measure of the accuracy of the effect estimate. We would hope to have small bias, percentage bias, and MSE for our estimates. I estimate bias, percentage bias, and MSE for β _A , β _B , and β _A _{× B}, respectively. I evaluate the performance of various methods using the percentage bias for the assessment of bias and MSE for the assessment of accuracy (Burton, Altman, Royston, & Holder, 2006).

Simulation Results

Tables 6 and 7 present the average percentages of the final analysis samples in the full sample for partial randomization and nonrandomization simulations, respectively. The common support sample for both propensity score models is quite large. Specifically, it is larger than 85% in the partial randomization simulation and larger than 68% in the nonrandomization simulation. It is not surprising, however, that the matching samples are small because I used one-to-one matching and allowed the sample allocation between two levels of treatment A and treatment B (or moderator) to vary (i.e., the proportion of Level 1 of treatment A ranges from 0.2 to 0.8 to simulate both balanced and nonbalanced experimental designs in the partial randomization simulation). It is smaller than 30% in the partial randomization simulation and smaller than 20% in the nonrandomization simulation.

Table 6.

Average Sample Sizes of Various Analyses in a Partial Randomization Simulation.

Propensity score model	Data	Proportion of Level 1 of treatment A

		0.2	0.3	0.4	0.5	0.6	0.7	0.8
One multinomial propensity score model	1. Full sample	100	100	100	100	100	100	100
	2. Common support sample (overlap among four groups)	87	88	89	90	90	90	90
	3. Overlap sample of matching one group with all the other three groups using Data Item 2 above	8	15	23	28	27	22	15
Two binary propensity score models	4. Common support sample (overlap among four groups based on two propensity scores)	86	87	88	88	88	87	85

Note. Entries are the average percentages of the final analysis samples in the full sample (N = 8,000). A total of 1,000 replications.

Table 7.

Average Sample Sizes of Various Analyses in a Nonrandomization Simulation.

Propensity score model	Data	Proportion of Level 1 of treatment A and treatment B (or moderator)

		P_a = 0.5, P_b = 0.5	P_a = 0.5, P_b = 0.3	P_a = 0.3, P_b = 0.3	P_a = 0.7, P_b = 0.3	P_a = 0.7, P_b = 0.7
NA	1. Full sample	100	100	100	100	100
One multinomial propensity score model	2. Common support sample (overlap among four groups)	82	80	74	81	76
One multinomial propensity score model	3. Overlap sample of matching one group with all other three groups using Data Item 1 above	19	15	14	10	16
Two binary propensity score models	4. Common support sample (overlap among four groups based on two propensity scores)	77	73	68	69	68

Note. Entries are the average percentages of the final analysis samples in the full sample (N = 8,000). A total of 1,000 replications.

In a partial randomization simulation, the results of seven types of sample allocation to the treatment group (Table 6) have very small variation. Thus, I report their average results. Table 8 presents the results for percentage bias and MSE by three parameters separately in a partial randomization simulation. The conventional OLS estimates, by including both variables and their interaction term using the full sample, produce large bias (46.1%, −23.0%, and −49.7%) and MSE (21.86, 5.62, and 25.89) for the coefficients of the interaction term, treatment A and treatment B (or moderator). This implies that, even in an experiment in which one treatment variable is randomly assigned, when the treatment effect is heterogeneous (i.e., the treatment effect varies by the covariate, W ₁), the conventional OLS estimate could produce biased estimates of the effects of two variables and their interaction (the average treatment effect estimate that uses OLS without the interaction term, however, would be unbiased).

Table 8.

Percentage Bias and MSE in a Partial Randomization Simulation.

Propensity score model	Data	Analysis model	%Bias $2 p t ({\hat{β}}_{A \times B})$	%Bias $2 p t ({\hat{β}}_{A})$	%Bias $2 p t ({\hat{β}}_{B})$	MSE $2 p t ({\hat{β}}_{A \times B})$	MSE $2 p t ({\hat{β}}_{A})$	MSE $2 p t ({\hat{β}}_{B})$
NA	1. Full sample	OLS	46.1	−23.0	−49.7	21.86	5.62	25.89
One multinomial propensity score model	1. Full sample	ANOVA (subclassification)	3.4	−1.6	24.8	5.48	2.98	9.50
		OLS (subclassification)	2.6	−0.9	−4.7	3.56	1.90	1.77
		ANOVA (weighting)	6.6	−3.2	40.2	11.20	5.67	20.99
		OLS (weighting)	7.4	−3.6	4.0	3.17	1.50	1.82
		ANCOVA (including propensity score as covariate)	46.5	−23.5	−19.6	22.69	6.03	5.23
		OLS (including propensity score as covariate)	45.9	−22.9	−23.9	21.70	5.54	6.90
	2. Common support sample (overlap among four groups)	ANOVA (subclassification)	2.4	−2.3	8.8	2.30	1.18	3.07
		OLS (subclassification)	1.7	−1.6	−2.5	1.62	0.80	0.83
		ANOVA (weighting)	1.1	−1.6	6.4	6.07	2.83	3.82
		OLS (weighting)	1.3	−1.3	1.0	2.10	1.16	1.28
	3. Overlap sample of matching one group with all other three groups using Data 2 above	ANOVA	0.4	3.8	0.7	1.99	1.18	1.27
		OLS	0.3	4.1	0.2	1.71	1.08	0.87
Two binary propensity score models	4. Common support sample (overlap among four groups based on two propensity scores)	ANOVA (subclassification)	6.0	−3.0	27.5	1.45	0.60	8.84
Two binary propensity score models		OLS (subclassification)	5.9	−2.9	−3.8	1.14	0.46	0.51

Note. ANOVA = analysis of variance; OLS = ordinary least square; ANCOVA = analysis of covariance; MSE = mean square error. A total of 1,000 replications. %Bias is the percentage bias $((\overset{\land}{B i a s} (\hat{θ}) θ) \times 100)$ .

For one multinomial propensity score model, the application of the common support sample has smaller percentage bias and MSE than does the application of the full sample, and the OLS estimate has smaller percentage bias and MSE than does an ANOVA. These results are consistent with previous propensity scoring studies of the use of the common support sample and the inclusion of covariates in the outcome model for double robustness (e.g., D’Agostino, 1998; Lunceford & Davidian, 2004). In addition, including the propensity score as a covariate in the analysis has smaller effects on bias and MSE reduction than do the other propensity score applications (e.g., subclassification and weighting). For two binary propensity score models, subclassification has similar good performance in MSE reduction as do subclassification, weighing, and matching for one multinomial propensity score model when OLS regression is used. However, subclassification for two binary propensity score models produces slightly larger percentage bias than do the other methods.

In a nonrandomization simulation, the results of five combinations of sample allocations have some variations. Thus, I report their individual results and average absolute values. I rely on the average results, however, to evaluate the performance of various methods. Tables 9 through 11 present the results for percentage bias and MSE for the coefficients of the interaction term, treatment A, and treatment B (or moderator), respectively. The two approaches that have smaller average percentage bias and MSE across three parameter estimates than do the conventional OLS analysis are (a) inverse of propensity score weighting based on one multinomial propensity score model and (b) subclassification based on two binary propensity score models.

Table 9.

Percentage Bias and MSE for ${\hat{β}}_{A \times B}$ in a Nonrandomization Simulation.

Propensity score model	Data	Analysis model	%Bias						MSE

			Proportion of Level 1 of two variables					Average^a	Proportion of Level 1 of two variables					Average^b

			P_a = 0.5, P_b = 0.5	P_a = 0.5, P_b = 0.3	P_a = 0.3, P_b = 0.3	P_a = 0.7, P_b = 0.3	P_a = 0.7, P_b = 0.7		P_a = 0.5, P_b = 0.5	P_a = 0.5, P_b = 0.3	P_a = 0.3, P_b = 0.3	P_a = 0.7, P_b = 0.3	P_a = 0.7, P_b = 0.7
NA	1. Full sample	OLS	21.0	−4.0	−29.3	23.9	74.3		30.5	4.86	0.80	9.40	6.56		56.07	15.54
One multinomial propensity score model	2. Common support sample	ANOVA (subclassification)	30.5	29.3	24.1	39.3	34.3	31.5	11.32	10.82	7.80	18.27	14.78	12.60
		OLS (subclassification)	48.3	46.8	45.5	51.0	47.5	47.8	25.07	23.59	22.64	28.39	25.72	25.08
		ANOVA (weighting)	3.4	0.7	5.6	2.5	2.9	3.0	15.98	23.63	19.10	23.10	20.21	20.40
		OLS (weighting)	4.7	3.0	4.5	6.1	6.7	5.0	4.35	4.93	5.05	5.71	4.95	5.00
	3. Overlap sample of matching one group with all other three groups	ANOVA	32.3	35.8	28.6	42.0	28.1	33.4	11.77	14.48	10.21	20.51	9.71	13.34
		OLS	53.1	56.3	58.8	53.3	49.7	54.2	29.43	33.36	36.62	31.00	26.46	31.37
Two binary propensity score models	4. Common support sample (overlap among four groups based on two propensity scores)	ANOVA (subclassification)	7.6	5.7	6.2	3.2	7.5	6.0	1.75	2.34	2.98	2.83	2.87	2.55
Two binary propensity score models		OLS (subclassification)	7.5	6.1	6.2	3.3	8.2	6.3	1.70	2.42	2.85	2.64	3.07	2.53

Note. ANOVA = analysis of variance; OLS = ordinary least square; MSE = mean square error. ^aEntries are the means of the absolute value of percentage bias among five combinations of sample allocations. ^bEntries are the means of the MSE among five combinations of sample allocations. A total of 1,000 replications for each parameter combination. %Bias is the percentage bias $((\overset{\land}{B i a s} (\hat{θ}) / θ) \times 100)$ .

In sum, the simulation results from partial randomization and nonrandomization scenarios suggest two good propensity score applications to reduce bias and MSE of parameter estimates in analyzing two variables: (a) inverse of propensity score weighting based on one multinomial propensity score model and (b) subclassification based on two binary propensity score models. Further, the use of a common support sample and covariate adjustment in outcome analysis will reduce bias and MSE.

Case Study: The Main and Differential Effects of Head Start

The Monte Carlo simulation suggested two appropriate applications of propensity score methods for analyzing the effects of two binary variables on the outcome. In this case study, I focus on these two approaches. Here, the inverse of propensity score weighting method (Imbens, 2000) can be used to analyze two multicategorical variables, while Imai and van Dyk’s (2004) method can be used for an analysis that involves continuous variables.

I demonstrate the applications of Imbens’ (2000) and Imai and van Dyk’s (2004) methods for the examination of the effects of the Head Start program as compared with other center-based care in improving children’s reading achievement, as differing by child care quality (as a binary variable and continuous variable, respectively). Because the Head Start program is a particular center-based care, this case study is also used to determine whether the care in the Head Start program is more effective than the care in the other center-based settings.

Sample and Measures

I use data from a nonrandomization study in which neither of the two variables (child care type and child care quality) concerned random assignment and their effect estimates suffer from selection bias. The data set used for this case study is the Early Childhood Longitudinal Study–Birth cohort (ECLS-B). ECLS-B began as a randomly selected nationally representative sample of 14,000 children born in 2001 in the United States. ECLS-B collects information on a rich array of individual-, household-, teacher-, and child care–level measures. I use the Item Response Theory (IRT) scale score for reading skill measured at kindergarten as the outcome variable, and I focus on observed global classroom/care quality. The quality of Head Start’s and other center-based care was assessed through the Early Childhood Environmental Rating Scale–Revised (ECERS-R; Harms, Clifford, & Cryer, 1998).

The ECLS-B administered 37 of the 43 items included in the original ECERS-R, which covered the following six subscales: (a) furnishing and display, (b) personal care routines, (c) listening and talking, (d) learning activities, (e) interaction, and (f) program structure. Each of the 37 items was rated on a 7-point Likert-type scale. I use the overall ECERS-R score in the analysis. The covariates used in this study include the mental IRT scale score measured at 2 years old, reading and math IRT scales measured at prekindergarten, gender, race/ethnicity (White, Black, Hispanic, and other race/ethnicity), birth weight, age at outcome assessment, English spoken at home, disability status, special education status, health problem, family’s SES, mother’s education, income, poverty level, welfare receipt, home violence, household structure, and one dummy variable for the year that the child was in kindergarten (2007 vs. 2006).

Analytic Procedure

One of the main variables of interest, the child care quality measure (the ECERS-R overall score) is a continuous variable. I dichotomize it with a median split to indicate high quality and low quality. Both the continuous and dichotomous ECERS-R scores are used for demonstration purposes. I first estimate the overall effect of Head Start as compared with other center-based care using conventional OLS regression and propensity score methods (subclassification and weighting), without including the care quality measure and its interaction term with care type. When child care quality is treated as a binary variable, the following approaches are used to examine the differential effect of Head Start between two child care quality groups:

Conventional OLS regression analysis with interaction of care quality (low quality vs. high quality) and care type (Head Start vs. other center-based care).

Inverse of propensity score weighting analysis based on one multinomial propensity score model (Imbens, 2000).

Subclassification based on two binary propensity score models (Imai & van Dyk, 2004). I used a 3 × 3 subclassification based on two propensity score functions estimated by two binary logistic regression models.

When care quality is analyzed as a continuous variable, I use Approach 1 with a modification that care quality is a median-centered continuous variable and Approach 3 with a modification in which subclassification is based on one binary propensity score model and one Gaussian linear regression model.

Findings

As presented in this section, I first checked covariate balance for the analyses before and after applying the propensity score methods. In the end, I report the results of examining whether the effect of the Head Start program varies by the care quality (moderator).

Because my main interest is analyzing both the main and moderator effects, I present the results of covariate balance checking for four comparison groups. I present the results of an F-test of a comparison on covariate means among four comparison groups and the maximum standardized mean difference among four groups. Because an F-test is used to test population parameters and is sensitive to sample size, it may not be appropriate for checking the sample characteristics (Austin, 2007; Imai, King, & Stuart, 2008). Thus, I focus on an interpretation of the maximum standardized mean difference among the four groups.

Tables 12 and 13 present covariate balance checking results for child care type by care quality groups before and after applying the inverse of propensity score weighting, respectively. Stuart (2007) suggests that the standardized mean difference should ideally be less than 0.25 and that a value greater than 0.50 is “particularly problematic.” As seen in Table 12, of the 22 covariates, 9 have a maximum standardized mean difference of greater than 0.50 among four groups. After weighting by the inverse of propensity score, 17 of the 22 covariates have a maximum standardized mean difference of less than 0.25 among four groups, and the other 5 covariates have a maximum standardized mean difference of less than 0.50 among four groups (Table 13)⁴. This suggests that the inverse of propensity score weighting can greatly reduce selection bias. Covariance balance also was examined (but not presented) for the subclassification method, which resulted in the covariates being much more balanced than without the use of any propensity score methods.

Table 10.

Percentage Bias and MSE for ${\hat{β}}_{A}$ in a Nonrandomization Simulation.

Propensity score model	Data	Analysis model	%Bias						MSE

			Proportion of Level 1 of two variables					Average^a	Proportion of Level 1 of two variables					Average^b

			P_a = 0.5, P_b = 0.5	P_a = 0.5, P_b = 0.3	P_a = 0.3, P_b = 0.3	P_a = 0.7, P_b = 0.3	P_a = 0.7, P_b = 0.7		P_a = 0.5, P_b = 0.5	P_a = 0.5, P_b = 0.3	P_a = 0.3, P_b = 0.3	P_a = 0.7, P_b = 0.3	P_a = 0.7, P_b = 0.7
NA	1. Full sample	OLS	16.1	27.6	8.6	34.9	−6.8		18.8	2.93	7.90	1.06	12.50		1.10	5.10
One multinomial propensity score model	2. Common support sample	ANOVA (subclassification)	−25.0	−8.4	−32.2	4.1	−19.1	17.7	7.56	1.77	10.98	1.96	6.21	5.69
		OLS (subclassification)	−27.3	−18.5	−33.7	−6.6	−16.3	20.5	8.53	4.04	11.84	1.79	5.76	6.39
		ANOVA (weighting)	2.1	4.4	−5.0	2.0	5.7	3.8	6.81	12.55	7.96	7.39	7.82	8.50
		OLS (weighting)	−1.9	−2.5	−7.9	2.6	3.9	3.8	2.44	2.49	2.82	2.04	3.31	2.62
	3. Overlap sample of matching one group with all other three groups	ANOVA	−23.7	−14.6	−39.0	1.7	−9.0	17.6	6.33	3.02	16.31	1.40	1.77	5.77
		OLS	−28.0	−28.6	−44.4	−13.4	−14.0	25.7	8.53	9.07	20.80	3.17	2.95	8.90
Two binary propensity score models	4. Common support sample (overlap among four groups based on two propensity scores)	ANOVA (subclassification)	−22.2	−20.9	−29.4	−12.1	−15.8	20.1	5.63	4.79	9.19	2.23	4.54	5.28
Two binary propensity score models		OLS (subclassification)	−3.0	−1.5	−9.3	6.7	3.4	4.8	0.68	0.34	1.31	1.08	2.21	1.12

Note. MSE = mean square error; OLS = ordinary least square; ANOVA = analysis of variance. ^aEntries are the means of the absolute value of percentage bias among five combinations of sample allocations. ^bEntries are the means of the MSE among five combinations of sample allocations. A total of 1,000 replications for each parameter combination. %Bias is the percentage bias $((\overset{\land}{B i a s} (\hat{θ}) / θ) * 100)$ .

Table 11.

Percentage Bias and MSE for ${\hat{β}}_{B}$ in a Nonrandomization Simulation.

Propensity score model	Data	Analysis model	%Bias						MSE

			Proportion of Level 1 of two variables					Average^a	Proportion of Level 1 of two variables					Average^b

			P_a = 0.5, P_b = 0.5	P_a = 0.5, P_b = 0.3	P_a = 0.3, P_b = 0.3	P_a = 0.7, P_b = 0.3	P_a = 0.7, P_b = 0.7		P_a = 0.5, P_b = 0.5	P_a = 0.5, P_b = 0.3	P_a = 0.3, P_b = 0.3	P_a = 0.7, P_b = 0.3	P_a = 0.7, P_b = 0.7
NA	1. Full sample	OLS	15.9	0.0	8.4	−20.5	−6.9		10.4	2.87	0.39	1.01	4.89		1.16	2.07
One multinomial propensity score model	2. Common support sample	ANOVA (subclassification)	−24.9	−56.2	−32.2	−84.0	−18.9	43.3	7.51	32.73	10.94	72.99	6.34	26.10
		OLS (subclassification)	−27.3	−43.7	−33.8	−54.6	−16.4	35.2	8.55	19.82	11.84	31.73	5.78	15.54
		ANOVA (weighting)	1.9	−7.6	−5.3	−19.1	5.6	7.9	6.95	16.00	7.90	18.24	8.56	11.53
		OLS (weighting)	−2.3	−5.4	−7.9	−8.5	3.6	5.5	2.40	3.71	2.80	5.52	3.45	3.57
	3. Overlap sample of matching one group with all other three groups	ANOVA	−23.9	−59.2	−39.2	−81.4	−8.7	42.5	6.35	35.94	16.48	67.73	1.81	25.66
		OLS	−28.4	−44.4	−44.6	−44.3	−14.0	35.1	8.75	20.79	21.03	21.27	2.96	14.96
Two binary propensity score models	4. Common support sample (overlap among four groups based on two propensity scores)	ANOVA (subclassification)	−22.5	−29.0	−29.7	−25.0	−16.0	24.4	5.77	9.39	9.27	8.51	4.61	7.51
Two binary propensity score models		OLS (subclassification)	−3.3	−10.5	−9.5	−9.1	3.4	7.1	0.65	1.97	1.31	2.84	2.17	1.79

Table 12.

Covariate Balance Check Among Child Care Type by Care Quality Groups Before Applying Propensity Score Methods.

Variable	Other center-based care^a		Head Start^a		F ^b	Maximum difference^c

	Low ECERS	High ECERS	Low ECERS	High ECERS
Reading at prekindergarten	27.37 (9.96)	29.85 (12.76)	20.86 (7.53)	22.49 (8.16)	33.24***	0.88
Math at prekindergarten	31.35 (9.75)	32.86 (10.49)	25.20 (8.46)	26.59 (8.30)	29.76***	0.80
Mental ability at 2 years old	127.53 (10.38)	128.79 (11.64)	125.55 (8.66)	123.86 (10.07)	9.44***	0.47
Birth weight	2969.64 (846.14)	3094.14 (833.14)	2996.35 (765.66)	2888.05 (832.67)	2.46	0.25
Age (month)	65.69 (3.70)	65.46 (3.69)	65.18 (3.53)	65.80 (3.36)	0.97	0.17
Mother’s education	14.37 (2.53)	14.77 (2.71)	12.17 (1.82)	12.24 (1.89)	69.00***	1.09
Mother’s age at birth (years)	28.11 (6.49)	29.73 (6.36)	24.78 (5.63)	24.43 (5.44)	37.50***	0.86
Income (thousands)	67.87 (54.10)	82.72 (60.14)	22.54 (17.03)	27.20 (26.16)	79.40***	1.26
SES	0.16 (0.81)	0.34 (0.82)	−0.65 (0.61)	−0.60 (0.65)	94.15***	1.31
Health problem	0.26 (0.54)	0.26 (0.56)	0.24 (0.57)	0.22 (0.46)	0.34	0.08
Violence	0.07 (0.33)	0.07 (0.32)	0.10 (0.37)	0.09 (0.37)	0.42	0.09
Disability	0.82 (1.63)	1.19 (1.85)	1.17 (1.39)	1.42 (1.51)	5.88***	0.36
Kindergarten in 2007 (%)	18 (39)	18 (39)	22 (42)	13 (34)	1.57	0.24
Girl (%)	52 (5)	49 (5)	55 (5)	49 (5)	0.62	0.13
Black (%)	18 (38)	15 (36)	5 (5)	29 (45)	24.58***	0.84
Hispanic (%)	13 (33)	14 (35)	19 (39)	27 (45)	6.83***	0.39
Other race (%)	22 (42)	19 (39)	11 (32)	14 (35)	3.35*	0.28
Single parent (%)	26 (44)	15 (36)	43 (5)	36 (48)	15.49***	0.65
Welfare recipient (%)	7 (25)	7 (26)	18 (39)	22 (42)	13.83***	0.50
Non-English (%)	15 (35)	16 (37)	17 (38)	21 (41)	1.38	0.18
Special education (%)	7 (25)	6 (24)	8 (27)	7 (26)	0.11	0.05
Poverty below 185% (%)	42 (49)	32 (47)	86 (35)	83 (38)	77.88***	1.20
Sample size^d	350	250	150	200

Note. ECERS = Early Childhood Environmental Rating Scale; NCES = National Center for Education Statistics. ^aEntries are the means and standard deviation (in parentheses). ^bEntries are the F-statistics of comparing the means across four groups (*p < .05. **p < .01. ***p < .001). ^cEntries are the maximum standardized mean differences among four groups, which were calculated based on the pooled standard deviations across all four groups. ^dSample size is rounded to the nearest 50 per NCES regulation.

Table 13.

Covariate Balance Check Among Child Care Type by Care Quality Groups Weighted by the Inverse of Propensity Score.

Variable	Other center-based care^a		Head Start^a		F ^b	Maximum difference^c

	Low ECERS	High ECERS	Low ECERS	High ECERS
Reading at prekindergarten	25.82 (16.94)	25.59 (20.84)	23.94 (22.13)	25.55 (21.4)	1.45	0.18
Math at prekindergarten	29.60 (16.79)	29.44 (19.11)	27.75 (24.00)	29.37 (20.88)	1.51	0.19
Mental ability at 2 years old	126.27 (18.75)	126.52 (21.41)	124.76 (21.41)	125.44 (20.26)	1.29	0.17
Birth weight	2982.38 (1437.14)	3014.68 (1573.4)	2964.93 (1935.03)	2876.49 (1738.71)	1.16	0.17
Age (month)	65.42 (6.34)	65.43 (6.98)	65.28 (8.46)	65.69 (6.79)	0.49	0.12
Mother’s education	13.64 (4.40)	13.56 (5.12)	12.61 (4.37)	13.15 (4.60)	8.28***	0.43
Mother’s age at birth (years)	27.14 (11.18)	26.90 (12.76)	25.28 (12.94)	26.42 (13.41)	3.53*	0.30
Income (thousands)	55.24 (88.69)	54.49 (101.5)	31.69 (54.51)	46.87 (98.1)	11.44***	0.49
SES	−0.10 (1.42)	−0.12 (1.65)	−0.45 (1.44)	−0.23 (1.55)	8.83***	0.47
Health problem	0.24 (0.85)	0.25 (1.03)	0.23 (1.32)	0.21 (0.93)	0.16	0.06
Violence	0.09 (0.60)	0.09 (0.67)	0.09 (0.79)	0.08 (0.72)	0.03	0.02
Disability	1.31 (3.92)	1.08 (3.42)	1.22 (3.73)	1.30 (3.37)	0.82	0.14
Kindergarten in 2007 (%)	19 (66)	16 (70)	23 (98)	17 (78)	1.08	0.17
Girl (%)	51 (84)	51 (95)	49 (117)	52 (104)	0.09	0.05
Black (%)	24 (72)	26 (83)	29 (1.07)	26 (92)	0.61	0.14
Hispanic (%)	17 (64)	18 (73)	21 (96)	21 (85)	0.64	0.11
Other race (%)	18 (65)	19 (75)	15 (84)	14 (73)	0.85	0.12
Single parent (%)	27 (75)	28 (85)	31 (109)	36 (99)	1.68	0.20
Welfare recipient (%)	13 (56)	12 (62)	19 (92)	15 (73)	1.79	0.22
Non-English (%)	16 (62)	18 (73)	16 (85)	17 (78)	0.25	0.07
Special education (%)	9 (48)	7 (48)	6 (57)	9 (58)	0.58	0.11
Poverty below 185% (%)	56 (84)	57 (94)	75 (102)	62 (101)	6.46***	0.42
Sample size^d	350	250	150	200

Note. ECERS = Early Childhood Environmental Rating Scale; SES = socioeconomic status. NCES = National Center for Education Statistics. ^aEntries are the means and standard deviation (in parentheses) weighted by the inverse of propensity score. ^bEntries are the F-statistics of comparing the weighted means across four groups (*p < .05. **p < .01. ***p < .001). ^cEntries are the maximum standardized mean differences among four groups, which were calculated based on the same pooled standard deviations as in Table 12 from the unweighted analysis. ^dSample size is rounded to the nearest 50 per NCES regulation.

Table 14 presents the differential effect of Head Start by a binary care quality measure. Neither the conventional OLS nor propensity score analyses detected a statistically significant difference between high- and low-quality care for the effect of Head Start, as compared to other center-based care.

Table 14.

The Differential Effect of Head Start as Compared With Other Center-Based Care on Kindergarten Reading by Binary Care Quality Measure.

Analysis	Estimate	Standard Error	p-Value	N ^a
Conventional OLS	−0.13	1.40	.92	950
Inverse of propensity score weighting	−1.26	1.30	.33	950
Subclassification^b	−1.53	1.93	.43	600

Note. OLS = ordinary least square. ECLS-B = Early Childhood Longitudinal Study–Birth cohort; SES = socioeconomic status; IRT = Item Response Theory. Adapted from ECLS-B. The covariates included in the propensity score models and adjusted in the outcome analysis include the mental IRT scale score measured at 2 years old, reading and math IRT scales measured at prekindergarten, gender, race (White, Black, Hispanic, and other race), birth weight, age at outcome assessment, English speaking at home, disability status, special education status, health problem, family’s SES, mother’s education, income, poverty level, welfare receipt, home violence, household structure, and one dummy variable that indicates the year that the child was in kindergarten (2007 vs. 2006). ^aSample size is rounded to the nearest 50 per National Center for Education Statistics regulations. ^bThe number of the original groups was nine, while the final groups used for weighted impact estimate was six. Three groups were excluded from analysis because at least one of the four groups (cells) contained fewer than five children.

Table 15 presents the differential effect of Head Start by a continuous care quality measure. Neither the conventional OLS nor propensity score analysis detected a statistically significant difference on the coefficient (slope) of care quality measure for Head Start as compared with other center-based care.

Table 15.

The Differential Effect of Head Start as Compared With Other Center-Based Care on Kindergarten Reading by Continuous Care Quality Measure.

Analysis	Estimate	Standard Error	p Value	N ^a
Conventional OLS	−0.19	0.72	.79	950
Subclassification^b	−0.32	0.98	.74	600

Note. OLS = ordinary least square; IRT = Item Response Theory; ECLS-B = Early Childhood Longitudinal Study–Birth cohort. Adapted from ECLS-B. The covariates included in the propensity score models and adjusted in the outcome analysis include the mental IRT scale score measured at 2 years old, reading and math IRT scales measured at prekindergarten, gender, race (White, Black, Hispanic, and other race), birth weight, age at outcome assessment, English speaking at home, disability status, special education status, health problem, family’s SES, mother’s education, income, poverty level, welfare receipt, home violence, household structure, and one dummy variable that indicates the year that the child was in kindergarten (2007 vs. 2006). ^aSample size was rounded to the nearest 50 per NCES regulations. ^bThe number of the original groups was nine, while the number of the final groups used for weighted impact estimate was six. Three groups were excluded from analysis because each contained fewer than five children in Head Start.

In sum, the reading achievement at kindergarten for children in Head Start was significantly lower than that of their peers in other center-based care. The effect of Head Start as compared with other center-based care on kindergarten reading, however, did not differ by care quality.

Conclusions

The Monte Carlo simulation suggested two appropriate applications of propensity score methods for analyzing two binary treatment variables. The inverse of the propensity score weighting method based on one multinomial propensity score model (Imbens, 2000) and subclassification based on two propensity score models (Imai & van Dyk, 2004) are both effective in reducing bias and MSE and can be easily applied. These two methods can be used for analyzing two binary variables. When two variables are multicategorical, for example, the total numbers of levels for two variables are J and K, Imbens’ (2000) method is easier to apply by converting a two-dimensional J × K design to a one-dimensional design with J × K levels. When one or two variables are continuous, Imai and van Dyk’s (2004) subclassification method is a good choice.

Note that this study assumes that the treatment variables and moderator are available for all sample members. However, in some circumstances, the interest variable is not measured for some particular members (e.g., in some experimental studies, the fidelity of implementation is measured only for the treatment group). The propensity score methods can be applied to this type of situation (e.g., Hill et al., 2003; Peck, 2003, 2013; Peck & Bell, 2012; Schochet & Burghardt, 2007); however, this is not the main interest of this article. Further, although this study did not directly address the issues that pertain to clustering, which is common in education evaluation, the methods discussed in this study can be extended to multilevel analysis, for example, the multilevel logistic regression models can be used to estimate the generalized propensity score and the hierarchical linear model for outcome analysis to take account of clustering effects in educational evaluation.

In sum, it is important to eliminate selection bias across variable-A-by-variable B groups in the analysis of the relationship between two variables and an outcome. This article identified and illustrated appropriate applications of propensity score methods for this purpose. Further, although this article focused on the analysis of two independent variables, the method can be extended to analyze mediation effects in future research.

Footnotes

Appendix A

Acknowledgment

The author thanks Mark Lipsey and Peter Schochet, two anonymous reviewers, and the former associate editor Laura Peck for their valuable comments on previous versions; however, all errors remain the author’s.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Notes

References

Austin

P. C.

(2007). The performance of different propensity score methods for estimating marginal odds ratios. Statistics in Medicine, 26, 3078–3094.

Baron

R. M.

Kenny

D. A.

(1986). The moderator-mediator variable distinction in social psychological research: Conceptual, strategic and statistical considerations. Journal of Personality and Social Psychology, 51, 1173–1182.

Burton

Altman

D. G.

Royston

Holder

R. L.

(2006). The design of simulation studies in medical statistics. Statistics in Medicine, 25, 4279–4292. doi:10.1002/sim.2673

D’Agostino

R. B.

Jr (1998). Propensity score methods for bias reduction in the comparison of a treatment to a non-randomized control group. Statistics in Medicine, 17, 2265–2281.

Frazier

P. A.

Tix

A. P.

Barron

K. E.

(2004). Testing moderator and mediator effects in counseling psychology research. Journal of Counseling Psychology, 51, 115–134.

Harms

Clifford

R. M.

Cryer

(1998). Early childhood environment rating scale-revised. New York, NY: Teachers College Press.

Hill

J. L.

Brooks-Gunn

Waldfogel

(2003). Sustained effects of high participation in an early intervention for low-birth-weight premature infants. Developmental Psychology, 39, 730–744.

Holland

P. W.

(1986). Statistics and causal inference. Journal of the American Statistical Association, 81, 945–960.

Hong

Nomi

(2012). Weighting methods for assessing policy effects mediated by peer change. Journal of Research on Educational Effectiveness, 5, 261–289.

10.

Imai

Keele

Tingley

(2010). A general approach to causal mediation analysis. Psychological Methods, 15, 309–334.

11.

Imai

King

Stuart

E. A.

(2008). Misunderstandings among experimentalists and observationalists in causal inference. Journal of the Royal Statistical Society, Series A, 171, 481–502.

12.

Imai

van Dyk

D. A.

(2004). Causal inference with general treatment regimes: Generalizing the propensity score. Journal of the American Statistical Association, 99, 854–866.

13.

Imbens

G. W.

(2000). The role of the propensity score in estimating dose-response functions. Biometrika, 87, 706–710.

14.

Imbens

G. W.

Rubin

(2009). Causal inference in statistics, and in the social and biomedical sciences. New York, NY: Cambridge University Press.

15.

Kraemer

H. C.

Kiernan

Essex

Kupfer

D. J.

(2008). How and why criteria defining moderators and mediators differ between the Baron & Kenny and MacArthur approaches. Health Psychology, 27. S101–S108.

16.

Kraemer

H. C.

Stice

Kazdin

Offord

Kupfer

(2001). How do risk factors work together to produce an outcome? Mediators, moderators, and independent, overlapping and proxy risk factors. American Journal of Psychiatry, 158, 848–856.

17.

Lochman

J. E.

Boxmeyer

C. L.

Powell

N. P.

Roth

Windle

(2006). Masked intervention effects: Analytic methods for addressing low dosage of intervention. New Directions for Evaluation, 110, 19–32.

18.

Greevy

Beck

(2011). Optimal nonbipartite matching and its statistical applications. The American Statistician, 65, 21–29.

19.

Rosenbaum

(2004). Optimal pair matching with two control groups. Journal of Computational and Graphical Statistics, 13, 422–434.

20.

Lunceford

J. K.

Davidian

(2004). Stratification and weighting via the propensity score in estimation of causal treatment effects: A comparative study. Statistics in Medicine, 23, 2937–2960.

21.

Peck

L. R.

(2003). Subgroup analysis in social experiments: Measuring program impacts based on post-treatment choice. American Journal of Evaluation, 24, 157–187.

22.

Peck

L. R.

(2005). Using cluster analysis in program evaluation. Evaluation Review, 29, 178–196.

23.

Peck

L. R.

(2013). On analysis of symmetrically predicted endogenous subgroups: Part one of a method note in three parts. American Journal of Evaluation, 34, 225–236. doi:10.1177/1098214013481666

24.

Peck

L. R.

Bell

S. H.

(2012). Estimating the influence of Head Start quality on child development. Paper presented at the 2012 Annual Fall Research Conference of the Association for Public Policy Analysis and Management, Baltimore, MD.

25.

Peck

L. R.

Camillo

D’Attoma

(2009). A promising new approach to eliminating selection bias. Canadian Journal of Program Evaluation, 24, 31–56.

26.

Potter

F. J.

(1993). The effect of weight trimming on nonlinear survey estimates. In Proceedings of the Section on Survey Research Methods of American Statistical Association. San Francisco, CA: American Statistical Association, San Francisco.

27.

Raudenbush

S. W.

(2011). Modeling mediation: Causes, markers, and mechanisms. Opening Address at the Spring 2011 Society for Research on Educational Effectiveness (SREE) Conference, Washington, DC. Retrieved from http://www.sree.org/conferences/2011/program/downloads/slides/raudenbush.pdf

28.

Rosenbaum

P. R.

(2002). Observational studies (2nd ed.). New York, NY: Springer-Verlag.

29.

Rosenbaum

P. R.

Rubin

D. B.

(1983). The central role of the propensity score in observational studies for causal effects. Biometrika, 70, 41–55.

30.

Rosenbaum

P. R.

Rubin

D. B.

(1984). Reducing bias in observational studies using subclassification on the propensity score. Journal of the American Statistical Association, 79, 516–524.

31.

Rubin

D. B.

(1974). Estimating causal effects of treatments in randomized and nonrandomized studies. Journal of Educational Psychology, 66, 688–701.

32.

Rubin

D. B.

(1980). Comments on “Randomization analysis of experimental data: The fisher randomization test,” by D. Basu. Journal of the American Statistical Association, 75, 591–593.

33.

Rubin

D. B.

(1990). Comments on “On the application of probability theory to agricultural experiments: Essay on principles. Section 9,” by J. Splawa-Neyman (D. M. Dabrowska & T. P. Speed, Trans.). Statistical Science, 5, 472–480.

34.

Scharfstein

D. O.

Rotnitzky

Robins

J. M.

(1999). Adjusting for non-ignorable drop-out using semiparametric non-response models. Journal of the American Statistical Association, 94, 1096–1120.

35.

Schochet

P. Z.

Burghardt

(2007). Using propensity scoring to estimate program-related subgroup impacts in experimental program evaluations. Evaluation Review, 31, 95–120.

36.

Shadish

W. R.

Cook

T. D.

Campbell

D. T.

(2002). Experimental and quasi-experimental designs for generalized causal inference. Boston, MA: Houghton-Mifflin.

37.

Stuart

E. A.

(2007). Estimating causal effects using school-level data sets. Educational Researcher, 36, 187–198.

38.

Stuart

E. A.

(2010). Matching methods for causal inference: A review and a look forward. Statistical Science, 25, 1–21.