Robust Inferences from a Before-and-After Study with Multiple Unaffected Control Groups

1.1 The before-and-after study with unaffected control groups for making causal inferences from an observational study: the Mariel Boatlift example

An observational study is an attempt to make causal inferences in a setting in which it is not feasible to randomly assign treatment to units. A major challenge for making causal inferences from observational study is that, because treatment is not randomized, the treatment and control groups may differ in ways other than the treatment. A “natural experiment” is an attempt to find a setting in the world where a treatment was handed to some people and denied to others for no particularly good reason that is haphazard [1, 2]. Good sources of natural experiments may be abrupt changes in government policies. As an example of a natural experiment from an abrupt change in government policy, we will consider Card’s [3] Mariel Boatlift study.

An important concern for immigration policy-makers is to what extent (if at all) do immigrants depress the labor market opportunities of native workers? A regression of native workers labor market opportunities (e.g., proportion employed or wages) on immigrant density in a city may yield a biased answer because immigrants tend to move to cities where the growth in demand for labor can accommodate their supply. Consequently, natural experiments, in which there are close to exogenous increases in the supply of immigrants to a particular labor market, are valuable for studying the effect of immigrants. Card [3] studied one such natural experiment, the experiences of the Miami labor market in the aftermath of the Mariel Boatlift.

Between April 15 and October 31, 1980, an influx of Cuban immigrants arrived in Miami, FL by boat from Cuba’s Mariel Harbor. Before this event there was a downturn in the Cuban economy which led to internal tensions, and subsequently the Cuban government announced that Cubans who wanted to leave could do so. By the end of the October, about 125,000 Cubans had arrived, at which point the Boatlift was ended by mutual agreement of the Cuban and U.S. governments. About half of the Cubans who arrived settled in Miami, thus significantly increasing the available labor force and number of Cuban workers in Miami.

We now discuss the before-and-after study with unaffected control groups study design for making causal inferences from the Mariel Boatlift study. In a causal inference study, a unit is an opportunity to apply or withhold treatment [4]. For the Mariel Boatlift study, we can think of the units as cities at a particular time, for example, Miami in 1979, Miami in 1982, Atlanta in 1979, Atlanta in 1982, etc. The Mariel Boatlift study can be thought of as a natural experiment that assigned some units to have a high amount of immigrants and some to have a low amount of immigrants, and where the outcomes of interest are labor market outcomes for native workers such as the unemployment rate. One way to study the effect of a high amount of immigrants would be to compare Miami after the Mariel Boatlift, say in 1982, to Miami before the Mariel Boatlift, say in 1979; this would be a before-and-after study. However, this study design may lead to biased results because the effect of the increase in immigration caused by the Mariel Boatlift may be confounded with macroeconomic changes between 1979 and 1982. In the before-and-after study with unaffected control groups, the change in outcomes in the place that got the treatment after a policy change or other event is compared to the change in outcomes in places that did not receive the treatment in either the before or after period, which are the unaffected control groups [2]; this is also sometimes called the difference-in-difference study design. The changes in the unaffected control groups control for changes in time, such as changes in the macroeconomy that would have occurred regardless of the treatment. Card [3] considered as unaffected control groups four other major cities in the United States – Atlanta, Houston, Los Angeles, and Tampa Bay–St. Petersburg. According to Card [3], these four cities did not experience a large increase in immigrants between 1979 and 1982 but were otherwise similar to Miami in demographics and pattern of economic growth. Consider the black unemployment rate as an outcome. The differences of blacks’ unemployment between 1979 and 1982 in Miami and the four unaffected control cities are shown in Figure 1.

The key assumption for the before-and-after study design with unaffected control groups to produce unbiased causal estimates is that the natural changes in the outcome that would have occurred in the group that received the treatment in the after period had that group not received treatment would have been the same as the change in the unaffected control groups [2]. This assumption is well discussed in the literature on treatments delivered to a group or cluster [5–7]. In the context of the Mariel Boatlift study, the key assumption is that the change in control cities like Atlanta or Houston serves as a counterfactual substitute and are exchangeable with the treated city, Miami.

In addition to needing this exchangeability assumption to hold, Figure 1 reveals an additional challenge: there is considerable variability among the changes in control groups, which suggests the possibility that the seemingly significant high difference of blacks’ unemployment rates in Miami may be just the result of random variation. To obtain valid causal inferences, we need to control for the natural variability among groups.

Figure 1

Differences of unemployment rates of blacks of 1979 and 1982 in five cities.

1.2 Current inference methods and the challenges

One way of making inferences that account for the variation in groups’ change in outcomes is a permutation approach [8]. The null hypothesis is that the unit to which the treatment was applied was randomly chosen and that the treatment has no effect. Consider the test statistic of the rank of the change between the before and after periods of the unit which was actually treated among all the units. Then, the p-value for testing the null hypothesis is the test statistic divided by the number of units. This permutation approach can work well if there are a large number of control groups, but it cannot find a significant effect, no matter how large the effect is, if there are only a few control groups. For example, if there is one treatment group and four control groups and the treatment group has the most extreme observation among all the groups, as in Figure 1, the permutation p-value is 1/5 0.2 no matter how large the treatment effect is.

Another nonparametric inference approach is the nonparametric bootstrap. Cameron et al. [9] proposed a series of nonparametric bootstrap approaches. However, we found in our simulation study (see Section 3) that when there is only one treatment group and four control groups, the bootstrap methods are not able to yield hypothesis tests with nominal significance levels.

As an alternative to non-parametric inference, we shall consider parametric models for the distribution of the control groups’ differences between the before and after periods, i.e. the group time effects. The most commonly used model is to assume that the groups’ differences between the before and after periods are i.i.d. normal random variables and to test whether the treatment group’s difference between the before and after periods could plausibly have come from the same distribution as the control groups [10]. This means that we conduct a t-test comparing the treatment group’s difference to the control groups’ differences. In the case where there are few treatment and control groups, the t-test is not protected by the central limit theorem and may be sensitive to the normal distribution assumption. From the simulation study in Section 3, the t-test actually performs well under non-normal distributed group variation but it is not totally robust, which means the rejection rate under null hypothesis may exceed its nominal significance level.

To make inferences more robust to violations of the normality assumption, we propose the robust-t method, a parametric method which considers a much wider range of distributions for the control groups’ differences. Specifically, we consider a set of skew-t distributions as in Jones and Faddy [11]. The family of skew-t distributions covers a wide range of shapes of distributions, and the normal distribution is a limiting case of the skew-t distribution. Our approach provides a researcher a method robust to a wide range of misspecification of the distribution of control groups’ differences between the before and after periods. We also propose a two-stage approach, which, compared to the robust-t method, has lower rejection rates under null hypothesis and lower power under alternative hypothesis, but provides an opportunity to conduct sensitivity analysis, and hence provides the researcher an idea about what kind of assumptions are necessary to reject the null hypothesis. Our overall proposal is to use the robust-t method to test for the overall range of shapes of group variation. If the test fails to reject, use the two-stage method to conduct a sensitivity analysis to see if there is a subset range of parameters for which we can be confident that there is a treatment effect.

The rest of the paper is organized as follows. In Section 2, we formulate the problem. In Section 3, we review some traditional methods and show the challenges with simulation. In Section 4, we develop the robust-t method, which is robust to a wide range of misspecification of the distribution of control groups’ differences between the before and after periods. In Section 5, we develop the two-stage method, which provides an opportunity to conduct sensitivity analysis. In Section 6, we compare the above two methods and propose to combine them when studying a real problems. In Sections 7 and Section 8, we apply the proposed methods to two examples. We conclude with discussion in Section 9.

2.1 Notation and model

In this paper, we consider groups of subjects and the outcomes are represented by the mean outcome of subjects of each group. Let be the average outcome of group k in period t assuming treatment was never applied to that group, and we model as

where indicates the period before and after the treatment, is the overall baseline level, is the time effect of group k, and ’s are independent normally distributed random errors with known standard deviation which signify the random measurement errors in our observations.

Let and represent the treatment and control groups, respectively. Letting be the average observed response for group k in period t, it follows

[1]

where is the causal effect of applying treatment in the groups k belonging to in the after period.

We assume that has the same distribution as . This means that we assume that the group(s) to which the treatment was applied was effectively randomly chosen. Presence of control groups allows us to separate the treatment effect from the time effect . If there is only one treatment and one control group, then it is common to assume time effects to be constant across groups, while if there are several treatment and/or control groups, the variation in the time effects can be taken into account by assuming that, for example, ’s follow a normal distribution with unknown variance. However, there is usually no strong a priori evidence for the normal distribution assumption.

In order to consider flexible group time effects, , we assume that are independently and identically distributed as a skew- distribution [11] with parameters controlling the skewness and scale, where is the nuisance scaling parameter.

The skew-t distribution is given by the following density function

where , denotes the beta function, and a and b are positive real numbers. When reduces to the density of the t-distribution with 2a degrees of freedom, and hence when and they are large, this distribution is similar to the normal distribution. When or is skewed to the left or right, respectively. The density functions of skew-t distributions with a range of parameters as well as the normal distribution are drawn in Figure 2, which shows that the skew-t distributions allow for flexibility in skewness and tail shapes.

Figure 2

Density functions of skew-t and normal distributions.

The differences in observed outcomes between the before and after treatment periods are

[2]

2.2 Miami Boatlift effect study on Miami employment rates

The above structure could be applied to the Miami unemployment rates data to study the effect of the influx of immigrants on the unemployment rates. The dataset is from the Current Population Survey (CPS), which contains the employment status of people sampled between 1979 and 1982 from Miami and four other major cities which were considered to have similar employment situation to Miami. For each city in each time period, an unemployment rate was obtained from a complex survey, along with standard error of the estimator . Considering the size of the dataset, it is appropriate to assume that the estimated unemployment rate is the true unemployment rate with an error following distribution. We regard the unemployment rates as observations. Hence, there is one treatment group (Miami) and four control groups (the other four cities), and each group has one observation before and after the treatment period. Consider to be the unemployment rate in kth city in year 1979 (if ) or 1982 (if ) estimated from the survey, along with be the standard error of the estimator (see Section 7.3 for details). Let be indicator of Miami and be the indicators of the four unaffected cities. Letting and the Mariel Boatlift effect be denoted by the unemployment rates can be modeled by eq. [2] exactly.

We would like to view the massive immigration in the Mariel Boatlift as the “treatment” and know whether it increased the unemployment rates in Miami. The hypothesis test of interest is given by

[3]

Because there are only a small number of groups, it is not possible to estimate the parameters a and b of the skew-t distribution with any accuracy, so instead we consider a range of a and b’s which imply different skewnesses and tail shapes. We propose a robust-t method which guarantees nominal significance level across a wide range of plausible skew-t distributions and hence is robust to misspecification of the distribution of group time effects (Section 4). We also consider a two-stage method that has lower power than the robust-t method but provides an opportunity to infer the sensitivity with respect to specification of , and hence provides the researcher an idea about what kind of assumptions are necessary to reject the null hypothesis (Section 5). In the study of a real problem, we propose to use the robust-t method to test for the overall range of skew-t distributions of interest. If the test fails to reject, use the two-stage method to conduct a sensitivity analysis to see if there is a subset range of parameters for which we can be confident that there is a treatment effect.

We apply the proposed method to two datasets. One dataset is from CPS to study the impact of the Mariel Boatlift on Miami unemployment rates between 1979 and 1982. The other dataset contains the student enrollment and grade repeating data in West Germany in the 1960s to study the impact of the short school year 1966–1967 on the grade repeating rates. Both datasets have few control and treatment groups, and each group has only one observation before and after the treatment period.

5.1 Two-stage inference

For a given pair of we may use a simulation-based two-stage method to test hypothesis (3).

The two-stage inference is usually applied to hypothesis test with nuisance parameter(s) [12–14]. Generally, in a hypothesis test about with significance level with an unknown nuisance parameter the basic two-stage inference is conducted as follows.

1. Specify where .

2. Derive a confidence interval of with coverage .

3. For conduct the hypothesis test with significance level .

4. If we accept the null hypothesis for we accept the null hypothesis. Then under null hypothesis, the false rejection rate is upper-bounded by .

In order to test the hypothesis (3), we follow the above two-stage inference procedure, where is the interval with the left end to be 0 and the right end to be upper bound of denoted by . One can prove that given larger , the null hypothesis is more likely to be rejected in the second stage. Thus in the second stage, we only need to test the null hypothesis with . The whole test procedure is summarized as follows.

1. Specify a range of plausible values.

2. Stage 1: For each pair of estimate an upper bound of with coverage with Monte Carlo simulation. The procedure to find the upper bound of is:

   1. Set a tolerance error for the upper bound, for example, 0.001. Set sufficiently large number to be and 0 to be . Set k0 and .

   2. Use Monte Carlo simulation to estimate the distribution of the standard deviation of the control groups’ time effects with given and , and determine the empirical quantile (denoted by ) of the observed standard deviation of the control groups’ time effects.

   3. If set . If set k to be and return to step (b). If exit and use as the upper bound of .

   4. If set . If set k to be and return to step (b). If exit and use as the upper bound of .

3. Stage 2: Denote the upper bound of to be . Let the test statistic to be Compute the value of the test statistic T from the data. With and eq. [2] do the following:

   1. Use Monte Carlo simulation to estimate the distribution of the test statistic T with .

   2. Either reject or accept depending on the empirical quantile of the observed test statistic value among the estimated distribution computed on the previous step, with significance level .

With the above two stage procedure, under null hypothesis the rejection rate would be upper-bounded by . Usually, the two-stage procedure is conservative, which means the rejection rate under the null hypothesis could be much smaller than .

The rejection rate under the null hypothesis depends on specific choice of and . We conducted a simulation study, and found that the choice of and has only slight influence on the rejection rate under the null hypothesis. In the rest of the paper, we set .

In Section 6, we compare the rejection rates under the null hypothesis and the powers under alternative hypothesis of the robust-t method and the two-stage method. We find that the power of this method is lower than the robust-t method proposed in Section 4.

5.2 Sensitivity analysis

When dealing with practical problems, usually some information about is available. For given assuming is known, one can conduct the hypothesis test as the second stage in the two-stage method. The test result may be heavily influenced by , and a sensitivity analysis can be conducted to find a subset of with which the null hypothesis is rejected.

Specifically, we set to range from low to high and investigate how the p-value of the hypothesis test changes along with . The procedure is summarized as follows:

1. Specify a range of plausible values.

2. For each pair of specify a range of ’s.

3. For in the range considered, compute the value of the test statistic T from the data, calculated the p-value of the test, and see how the p-value would change with .

From the plot of ’s and the corresponding p-values, the method provides a researcher with an idea of what kinds of assumptions about the variation, skewness, and heavy-tailedness of the control groups’ differences would be necessary to infer that the treatment has an effect. See the applications in Sections 7 and 8.

7.1 Dataset and problem

Between April 15 and October 31, 1980, an influx of Cuban immigrants arrived in Miami, Florida by boat from Cuba’s Mariel Harbor. In order to study the effect of mass immigration on unemployment rates, we obtained the CPS data of 1979 and 1982 in Miami and four other cities: Atlanta, Houston, Los Angeles, and Tampa Bay–St. Petersburg, which did not experience a large increase in immigrants between 1979 and 1982 but, according to Card [3], were otherwise similar to Miami in demographics and pattern of economic growth. The dataset contained a “race” variable which divided subjects in “white,” “black,” and “other.” The dataset also contained an “ethnic” variable which divided subjects into “Mexican American,” “Chicano,” “Mexicano,” “Puerto Rican,” “Cuban,” “Central or South American,” “Other Spanish,” “All other,” and “Don’t know.” All the non-hispanics except unknowns were labeled as “All other.” Within the dataset we get, all people labeled as “Cuban” were labeled as “white” in terms of race.

Table 2 and Figure 3 show the unemployment rates in 1979 (before the Boatlift) and 1982 (after the Boatlift) in Miami and four control cities. Figure 4 shows the increment of the unemployment rates from 1979 to 1982. All of the cities experienced increase in unemployment between 1979 and 1982, but Miami experienced a larger increase, particularly among blacks. Is there strong evidence that Miami experienced a larger increase in unemployment than could be expected if the immigration had no effect? The effect of the immigration to Miami is confounded with the downturn in the economy between 1979 and 1982. Therefore, we control for time by comparing the difference of unemployment rates between the before and after periods in Miami to those of the unaffected cities.

Figure 3

Unemployment rates of 1979 and 1982 in the five cities.

Figure 4

Differences of unemployment rates of 1979 and 1982 in the five cities.

Since the treatment effect is the same for all the individuals sampled in Miami and the city time effect is the same for all individuals sampled in the corresponding city, by aggregating the individual employment status to unemployment rates we would not lose the information about the city time effects and treatment effect. Also because of the complex sampling method in the CPS, the usual assumption that the observations are independent is not plausible for this dataset. Hence, for this dataset, for each city in each time period, an unemployment rate was obtained from the complex survey, along with the standard error of the estimators (see Section 7.2 for details). We regard the unemployment rates as the observations, and hence there is one treatment group (Miami) and four control groups (the other four cities), and each group has one observation before and after the treatment period. The years 1979 and 1982 correspond to the time period before and after the migration, respectively.

We would like to know whether the massive immigration in the Mariel Boatlift had an increasing effect on the unemployment rates in Miami. Let be the estimated unemployment rate in the kth city in year 1979 (if ) or 1982 (if ). Let be indicator of Miami and be the indicators of the four unaffected cities. The standard error of the estimation is estimated with binomial approximation corrected by a sampling strategy factor provided by CPS. Considering the size of the dataset, it is appropriate to assume that the estimated is the true unemployment rate with an error following distribution, where is the estimated standard error of estimated . Letting and the Mariel Boatlift effect denoted by the unemployment rate data can be models by eq. [2]. The hypothesis test of interest is given by

[4]

We test the hypothesis with the unemployment rates of four different groups: the unemployment rate of the population, the unemployment rate of whites without Hispanics, the unemployment rate of blacks without Hispanics, and the unemployment rate of Hispanics.

7.2 Data and setup

We obtain and as follows.

We estimated the unemployment rate by the ratio of number of unemployed sampled individuals and available labor force in the dataset. We use the variable “ftpt79 (Full-time or part-time labor force status)” as the indicator of employment condition. The value of this variable includes: not in labor force (0), employed full-time (1), part-time for economic reasons (2), unemployed full-time (3), employed part-time (4), and unemployed part-time (5). We estimate the unemployment rate by the ratio of the number of individuals with value (3) and (5) to the number of individuals with value values (1), (2), (3), (4) or (5).

is estimated as the standard error of estimated . For each city k at time t, let x be the number of unemployed people and y be the labor force, we have . The CPS calculate the standard error of by

where r is the correlation between x and y. According to the U.S. Bureau of Census (2000), the analyst should assume r is equal to zero [15]. For each city k at time t, let n be the total number of observations, we calculate

and

where are the design effect coefficients on standard errors. We use the design effects for total monthly variances after composing on Unemployment and Civilian Labor Force as Table 14–5 in CPS [15, Chapter 14]. For each city in each time period, we calculate , and obtain as the standard error of unemployment rate. The estimated standard error of city k at time t are treated as the standard deviation of the error terms , for and

7.3 Robust-t test of the mass migration effect

7.4 Sensitivity analysis

Since the robust-t approach fails to reject the null hypothesis, we use the two-stage method to do a sensitivity analysis as in Section 5.2 to see if there is a subset range of parameters for which we can be confident that there is a treatment effect.

We set to be various values and look at how the p-value of the hypothesis test changes along with the . As an example, we summarize the results with Blacks’ unemployment rate with a small collection of ’s in Figure 5. The graph implies that, even with , the treatment effect is not statistically significant because of the variation in the estimated unemployment rates, and hence there is not significant evident to conclude that the migration affects the Blacks’ unemployment rates.

Figure 5

p-value changes with for blacks.

The height of the horizontal line is 0.025 and 0.05 respectively.

8.1 Dataset and problem

The number of students enrolled in each grade and the number of students repeating the corresponding grade was published annually by the Federal Statistical Office in the serial Fachserie A. Bevölkerung und Kultur, Reihe 10, I, Allgemeines Bildungswesen. Our dataset contains the numbers of students and repeating students in grade one to four in the eleven states of West Germany from the year 1961–1971. In 1965, all the West German states had a normal long school year, while in 1966, some of the states experienced a short school year in order to fulfill the requirement of the agreement by the beginning of the 1967 school year. We view this short school year in 1966–1967 as the “treatment” to the states, and would like to study the influence of short school year on grade repeating rates, assuming state time effects following skew-t distribution.

Among the eleven states, Schleswig-Holstein, Bremen, Nordrhein-Westphalen, Hessen, Rheinland-Pfalz, Baden-Württemberg and Saarland experienced short school year in 1966–1967, and hence they were the treatment states. Hamburg, Bayern and Berlin stuck to a regular long school year, and thus they were the control states. Niedersachsen adopted the short school years in 1966, but added additional time in subsequent years for some schools, so it was neither purely treatment nor control group, and we decided not to consider this state in this study.

We summarize the average grade repeating rates of treatment and control states in 1965 and 1966 in Table 5 and draw them in Figure 6. They suggest that the short school year may have increased the grade repeating rates. The differences of the grade repeating rates of year 1966 and 1965 in Figure 7 show the same implication.

With the model (2) described in Section 2, we study the short school year effect on each grade with aggregate grade repeating rates in 1965 and 1966. There are seven treatment groups and three control groups. For each grade let be the grade repeating rate in kth state in year 1965 (if ) or 1966 (if ), which is estimated by the ratio of the number of students repeating a grade and the number of students enrolled in the corresponding grade. is the standard error of the corresponding grade repeating rate, which is estimated by binomial approximation. Let be the indicators of the seven treatment states (Schleswig-Holstein, Bremen, Nordrhein-Westphalen, Hessen, Rheinland-Pfalz, Baden-Württemberg and Saarland) and be the indicators of the three unaffected states (Hamburg, Bayern and Berlin). The short school year effect is represented by and the hypothesis test of interest is

Figure 6

Grade repeating rates of 1965 and 1966 in the ten states.

Figure 7

Differences of grade repeating rates of 1965 and 1966 in the ten states.

8.2 Robust-t test of the short school year effect

8.3 Sensitivity analysis

We use the two-stage method to do a sensitivity analysis as in Section 5.2 to see if there is a subset range of parameters for which we can be confident that there is a treatment effect. We again set to be various value and look at how the p-value of the hypothesis test changes along with the . As an example, we summarize the results with the grade 2 repeating rate with a small collection of ’s in Figure 8. The graph implies that, one has to assume that is small to conclude that the short school year affects the grade 2 repeating rate. Typically needs to be smaller than 0.005 for the collection of ’s we tested.

In the two-stage estimation process, we estimate the 0.975 upper bound of in stage 1. For the collection of ’s we test and the grade 2 repeating rate data, the 0.975 upper bounds of varies from 0.002 to 0.01, which is not large comparing to the 0.005 threshold. If one assumes to be smaller than 0.005, it may be concluded that the short school year affects the grade 2 repeating rate.

Figure 8

p-value changes with for grade 2 repeating rate.

The height of the horizontal line is 0.025 and 0.05 respectively.