Randomization and matching

Program details

In this hypothetical situation, an NGO is planning on launching a training program designed to boost incomes. Based on their experiences in running pilot programs in other countries, they’ve found that older, richer men tend to self-select into the training program. The NGO’s evaluation consultant (you!) drew this causal model explaining the effect of the program on participant incomes, given the confounding caused by age, sex, and prior income:

income_dag <- dagify(post_income ~ program + age + sex + pre_income,
                     program ~ age + sex + pre_income,
                     exposure = "program",
                     outcome = "post_income",
                     labels = c(post_income = "Post income",
                                program = "Program",
                                age = "Age",
                                sex = "Sex",
                                pre_income = "Pre income"),
                     coords = list(x = c(program = 1, post_income = 5, age = 2, 
                                         sex = 4, pre_income = 3),
                                   y = c(program = 2, post_income = 2, age = 1, 
                                         sex = 1, pre_income = 3)))

ggdag_status(income_dag, use_labels = "label", text = FALSE, seed = 1234) + 
  guides(color = FALSE) +
  theme_dag()

The NGO just received funding to run a randomized controlled trial (RCT) in a village, and you’re excited because you can finally manipulate access to the program—you can calculate \(E(\text{Post-income} | do(\text{Program})\). Following the rules of causal diagrams, you get to delete all the arrows going into the program node:

income_dag_rct <- dagify(post_income ~ program + age + sex + pre_income,
                         exposure = "program",
                         outcome = "post_income",
                         labels = c(post_income = "Post income",
                                    program = "Program",
                                    age = "Age",
                                    sex = "Sex",
                                    pre_income = "Pre income"),
                         coords = list(x = c(program = 1, post_income = 5, age = 2, 
                                             sex = 4, pre_income = 3),
                                       y = c(program = 2, post_income = 2, age = 1, 
                                             sex = 1, pre_income = 3)))

ggdag_status(income_dag_rct, use_labels = "label", text = FALSE, seed = 1234) +
  guides(color = FALSE) +
  theme_dag()

1. Check balance

You ran the study on 1,000 participants over the course of 6 months and you just got your data back.

Before calculating the effect of the program, you first check to see how well balanced assignment was, and you find that assignment to the program was pretty much split 50/50, which is good:

village_randomized %>%
  count(program) %>% 
  mutate(prop = n / sum(n))

## # A tibble: 2 x 3
##   program        n  prop
##   <chr>      <int> <dbl>
## 1 No program   503 0.503
## 2 Program      497 0.497

You then check to see how well balanced the treatment and control groups were in participants’ pre-treatment characteristics:

village_randomized %>% 
  group_by(program) %>% 
  summarize(prop_male = mean(sex_num),
            avg_age = mean(age),
            avg_pre_income = mean(pre_income))

## # A tibble: 2 x 4
##   program    prop_male avg_age avg_pre_income
##   <chr>          <dbl>   <dbl>          <dbl>
## 1 No program     0.584    34.9           803.
## 2 Program        0.604    34.9           801.

These variables appear fairly well balanced. To check that there aren’t any statistically significant differences between the groups, you make some graphs (you could run t-tests too, but graphs are easier for your non-statistical audience to read later).

There were more men in both the treatment and control groups, but the proportion is the same in both, and there’s no substantial difference in sex proportion:

# Here we save each plot as an object so that we can combine the two plots with
# + (which comes from the patchwork package). If you want to see what an
# individual plot looks like, you can run `plot_diff_sex`, or whatever the plot
# object is named.
#
# stat_summary() here is a little different from the geom_*() layers you've seen
# in the past. stat_summary() takes a function (here mean_se()) and runs it on
# each of the program groups to get the average and standard error. It then
# plots those with geom_pointrange. The fun.args part of this lets us pass an
# argument to mean_se() so that we can multiply the standard error by 1.96,
# giving us the 95% confidence interval.
plot_diff_sex <- ggplot(village_randomized, aes(x = program, y = sex_num, color = program)) +
  stat_summary(geom = "pointrange", fun.data = "mean_se", fun.args = list(mult = 1.96)) +
  guides(color = FALSE) +
  labs(x = NULL, y = "Proportion male")
# plot_diff_sex  # Uncomment this if you want to see this plot by itself

plot_prop_sex <- ggplot(village_randomized, aes(x = program, fill = sex)) +
  # Using position = "fill" makes the bars range from 0-1 and show the proportion
  geom_bar(position = "fill") +
  labs(x = NULL, y = "Proportion", fill = NULL) +
  scale_fill_manual(values = c("darkblue", "darkred"))

# Show the plots side-by-side
plot_diff_sex + plot_prop_sex

The distribution of ages looks basically the same in the treatment and control groups, and there’s no substantial difference in the average age across groups:

plot_diff_age <- ggplot(village_randomized, aes(x = program, y = age, color = program)) +
  stat_summary(geom = "pointrange", fun.data = "mean_se", fun.args = list(mult = 1.96)) +
  guides(color = FALSE) +
  labs(x = NULL, y = "Age")

plot_hist_age <- ggplot(village_randomized, aes(x = age, fill = program)) +
  geom_histogram(binwidth = 1, color = "white") +
  guides(fill = FALSE) +
  labs(x = "Age", y = "Count") +
  facet_wrap(vars(program), ncol = 1)

plot_diff_age + plot_hist_age

Pre-program income is also distributed the same—and has no substantial difference in averages—across treatment and control groups:

plot_diff_income <- ggplot(village_randomized, aes(x = program, y = pre_income, color = program)) +
  stat_summary(geom = "pointrange", fun.data = "mean_se", fun.args = list(mult = 1.96)) +
  guides(color = FALSE) +
  labs(x = NULL, y = "Pre income")

plot_hist_income <- ggplot(village_randomized, aes(x = pre_income, fill = program)) +
  geom_histogram(binwidth = 20, color = "white") +
  guides(fill = FALSE) +
  labs(x = "Pre income", y = "Count") +
  facet_wrap(vars(program), ncol = 1)

plot_diff_income + plot_hist_income

All our pre-treatment covariates look good and balanced! You can now estimate the causal effect of the program.

2. Estimate difference

You can find the causal effect (or average treatment effect) of the program with this formula, based on potential outcomes:

\[ \text{ATE} = E(\overline{\text{Post income }} | \text{ Program} = 1) - E(\overline{\text{Post income }} | \text{ Program} = 0) \]

This is simply the average outcome for people in the program minus the average outcome for people not in the program. You calculate the group averages:

village_randomized %>% 
  group_by(program) %>% 
  summarize(avg_post = mean(post_income))

## # A tibble: 2 x 2
##   program    avg_post
##   <chr>         <dbl>
## 1 No program    1180.
## 2 Program       1279.

That’s 1279 − 1180, or 99, which means that the program caused an increase of $99 in incomes, on average.

Finding that difference required some manual math, so as a shortcut, you run a regression model with post-program income as the outcome variable and the program indicator variable as the explanatory variable. The coefficient for program is the causal effect (and it includes information about standard errors and significance). You find the same result:

model_rct <- lm(post_income ~ program, data = village_randomized)
tidy(model_rct)

## # A tibble: 2 x 5
##   term           estimate std.error statistic  p.value
##   <chr>             <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)      1180.       4.27     276.  0.      
## 2 programProgram     99.2      6.06      16.4 1.23e-53

Based on your RCT, you conclude that the program causes an average increase of $99.25 in incomes.

Program details

A consortium of MPA and MPP programs are interested in improving the quantitative skills of their students before they begin their programs. Some schools have developed a two-week math camp that reviews basic algebra, probability theory, and microeconomics as a way to jumpstart students’ quantitative skills for classes like statistics, microeconomics, and program evaluation (👋 y’all!).

These schools have collected data on student outcomes, measuring final degree outcomes with a (totally fake) 120-160 point scale. You’re curious about whether math camps actually have a causal effect on final degree scores.

Unfortunately for you, these schools did not use an RCT to provide access to math camp—students self selected into the program, and all you have is observational data. However, armed with the knowledge of DAGs, confounders, and do-calculus, you think you can still estimate a causal effect!

(For reference, the true causal effect of this (totally fake) program is 10)

Based on your extensive knowledge of MPA/MPP grades and math classes, you draw the following causal model:

math_camp_dag <- dagify(
  final_grade ~ math_camp + gre_quant + gre_verbal + 
    undergraduate_gpa + background,
  math_camp ~ needs_camp, 
  needs_camp ~ background + undergraduate_gpa + gre_quant,
  gre_quant ~ background + undergraduate_gpa,
  gre_verbal ~ background + undergraduate_gpa,
  undergraduate_gpa ~ background,
  exposure = "math_camp",
  outcome = "final_grade",
  latent = c("background", "needs_camp"),
  coords = list(x = c(math_camp = 2, final_grade = 4, needs_camp = 1, gre_quant = 2.5, 
                      gre_verbal = 5, background = 2, undergraduate_gpa = 4), 
                y = c(math_camp = 1, final_grade = 1, needs_camp = 2, gre_quant = 2, 
                      gre_verbal = 2, background = 3, undergraduate_gpa = 3)),
  labels = c(math_camp = "Math camp", final_grade = "Final grade", 
             needs_camp = "Needs camp", gre_quant = "GRE quantitative", 
             gre_verbal = "GRE verbal", background = "Background",
             undergraduate_gpa = "Undergraduate GPA")
)

ggdag_status(math_camp_dag, use_labels = "label", text = FALSE, seed = 1234) + 
  guides(color = guide_legend(title = NULL)) +
  theme_dag() + 
  theme(legend.position = "bottom")

Your final grade in the program is caused by a host of things, including your quantitative and verbal GRE scores (PROBABLY DEFINITELY NOT in real life, but go with it), your undergraduate GPA, and your unmeasured background factors (age, parental income, math anxiety, level of interest in the program, etc.). Your undergraduate GPA is determined by your background, and your GRE scores are determined by both your undergraduate GPA and your background. Because this math camp program is open to anyone, there is self-selection in who chooses to do it. We can pretend that this is decided by your undergraduate GPA, your quantitative GRE score, and your background. If the program was need-based and only offered to people with low GRE scores, we could draw an arrow from GRE quantitative to math camp, but we don’t. Finally, needing the math camp causes people to do it.

There is a direct path between our treatment and outcome (math camp → final grade), but there is also some possible backdoor confounding. Both GRE quantitative and undergraduate GPA have arrows pointing to final grade and math camp (through “needs camp”), which means they’re a common cause, and background is both a confounder and unmeasurable. But you don’t need to give up! If you adjust or control for “needs camp,” you can block the association between background, GRE quantitative, and undergraduate GPA. With this backdoor closed, you’ve isolated the math camp → final grade relationship and can find the causal effect.

However, you don’t really have a measure for needing math camp—we can’t read peoples’ minds and see if they need the program—so while it’d be great to just include a needs_camp variable in a regression model, you’ll have to use other techniques to close the backdoor.

Since you don’t have a variable to indicate needing math camp, you can draw a slightly simpler DAG. Note how the “needs camp” node is an intermediate node on the path between GPA and GRE scores and actually participating in the math camp program. If you can guess what causes people to enroll in the program, that’s roughly the same as predicting their need for the camp. That means you can remove that node (and get rid of background too, just for the sake of extra simplicity; technically it’s still unmeasured too).

math_camp_dag_simpler <- dagify(
  final_grade ~ math_camp + gre_quant + gre_verbal + undergraduate_gpa,
  math_camp ~ undergraduate_gpa + gre_quant,
  gre_quant ~ undergraduate_gpa,
  gre_verbal ~ undergraduate_gpa,
  exposure = "math_camp",
  outcome = "final_grade",
  coords = list(x = c(math_camp = 2, final_grade = 4, gre_quant = 2.5, 
                      gre_verbal = 5, undergraduate_gpa = 4), 
                y = c(math_camp = 1, final_grade = 1, gre_quant = 2, 
                      gre_verbal = 2, undergraduate_gpa = 3)),
  labels = c(math_camp = "Math camp", final_grade = "Final grade", gre_quant = "GRE quantitative", 
             gre_verbal = "GRE verbal", undergraduate_gpa = "Undergraduate GPA")
)

ggdag_status(math_camp_dag_simpler, use_labels = "label", text = FALSE, seed = 1234) + 
  guides(color = guide_legend(title = NULL)) +
  theme_dag() + 
  theme(legend.position = "bottom")

Naive difference in means

For fun, you can calculate the difference in average grades for those who did/didn’t participate in math camp. This is most definitely not the actual causal effect—this is the “correlation is not causation” effect that doesn’t account for any of the backdoors in the DAG.

You can do this with a table (but then you have to do manual math):

math_camp %>% 
  group_by(math_camp) %>% 
  summarize(number = n(),
            avg = mean(final_grade))

## # A tibble: 2 x 3
##   math_camp    number   avg
##   <fct>         <int> <dbl>
## 1 No math camp   1182  138.
## 2 Math camp       818  144.

Or you can do it with regression:

model_wrong <- lm(final_grade ~ math_camp, data = math_camp)
tidy(model_wrong)

## # A tibble: 2 x 5
##   term               estimate std.error statistic   p.value
##   <chr>                 <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)          138.       0.185     744.  0.       
## 2 math_campMath camp     6.56     0.290      22.6 3.70e-101

According to this estimate, participating in math camp is associated with a 6.56 point increase in final grades. We can’t legally talk about casual effects though.

Adjustment using educated-guess-based matching

One big issue with trying to isolate or identify the causal effect between math camp and final grade is that the “needs camp” node is unmeasurable and unobserved (or “latent” in DAG terminology). There’s something that makes people need to participate in math camp, and based on the DAG it’s related to GPA and quantitative GRE scores, but we don’t know what it is exactly. We can attempt to match observations by their need for camp using a kind of manual coarsened exact matching (CEM). If we make an arbitrary decision based on GRE scores or GPA and assume that people below some grade or test score need math camp, that can substitute for the missing “needs camp” node.

You can plot a histogram of GRE scores to see if there’s any possible systematic reason for people to participate in the camp:

ggplot(math_camp, aes(x = gre_quant, fill = math_camp)) +
  geom_histogram(binwidth = 2, color = "white") + 
  guides(fill = FALSE) +
  facet_wrap(vars(math_camp), ncol = 1)

There’s a visible break in the bottom panel. There are a lot more people in math camp who scored under 145 than those who scored above 145. You can guess that people who scored under 145 need the camp and make a new variable indicating that:

math_camp_guessed_need <- math_camp %>% 
  mutate(maybe_needs_camp = gre_quant < 145)

We can now adjust for “needs camp” in a regression model, which closes that backdoor and gets us a more accurate estimate of the causal effect:

model_adj_needs_camp_guess <- lm(final_grade ~ math_camp + maybe_needs_camp, 
                                 data = math_camp_guessed_need)
tidy(model_adj_needs_camp_guess)

## # A tibble: 3 x 5
##   term                 estimate std.error statistic   p.value
##   <chr>                   <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)            139.       0.193     724.  0.       
## 2 math_campMath camp       8.57     0.290      29.6 7.41e-160
## 3 maybe_needs_campTRUE    -5.24     0.286     -18.3 1.31e- 69

According to this rough estimate of needing math camp, participating in the program causes a 8.57 point increase in final grade.

Adjustment with Mahalanobis nearest-neighbor matching

Instead of trying to figure out the hidden “needs camp” node, we can try modeling participation in math camp directly and essentially get rid of the need for camp. We can use matching techniques to pair up similar observations and make the unconfoundedness assumption—that if we see two observations that are pretty much identical, and one went to math camp and one didn’t, that choice was random.

Because we know from the DAG that undergraduate GPA and quantitative GRE scores help cause participation in math camp, we’ll try to find observations with similar values of GPA and test scores that both went to math camp and didn’t go to math camp.

You can use the matchit() function from the MatchIt R package to match points based on Mahalanobis distance. There are lots of other options available—see the online documentation for details.

You can include the replace = TRUE option to make it so that points that have been matched already can be matched again (that is, we’re not forcing a one-to-one matching; we have one-to-many matching instead).

# For whatever reason, matchit() doesn't work with categorical variables, so we
# have to use math_camp_num instead of math_camp here
matched <- matchit(math_camp_num ~ undergrad_gpa + gre_quant, data = math_camp,
                   method = "nearest", distance = "mahalanobis", replace = TRUE)
summary(matched)

## 
## Call:
## matchit(formula = math_camp_num ~ undergrad_gpa + gre_quant, 
##     data = math_camp, method = "nearest", distance = "mahalanobis", 
##     replace = TRUE)
## 
## Summary of balance for all data:
##               Means Treated Means Control SD Control Mean Diff eQQ Med eQQ Mean
## undergrad_gpa           2.8           3.2       0.44     -0.38    0.41     0.38
## gre_quant             142.9         146.5       4.56     -3.61    4.00     3.60
##               eQQ Max
## undergrad_gpa    0.45
## gre_quant        5.00
## 
## 
## Summary of balance for matched data:
##               Means Treated Means Control SD Control Mean Diff eQQ Med eQQ Mean
## undergrad_gpa           2.8           2.8       0.47    -0.007    0.22     0.22
## gre_quant             142.9         142.9       4.81    -0.056    2.00     2.29
##               eQQ Max
## undergrad_gpa    0.37
## gre_quant        4.00
## 
## Percent Balance Improvement:
##               Mean Diff. eQQ Med eQQ Mean eQQ Max
## undergrad_gpa         98      47       41      18
## gre_quant             98      50       36      20
## 
## Sample sizes:
##           Control Treated
## All          1182     818
## Matched       368     818
## Unmatched     814       0
## Discarded       0       0

Here you can see that all 818 of the math camp participants were paired with similar-looking non-participants (368 of them). 814 people weren’t matched and will get discarded. If you’re curious, you can see which treated rows got matched to which control rows by running matched$match.matrix.

You can create a new data frame of those matches with match.data():

math_camp_matched <- match.data(matched)

Now that the data has been matched, it should work better for modeling. Also, because we used undergraduate GPA and quantitative GRE scores in the matching process, we’ve adjusted for those DAG nodes and have closed those backdoors, so our model can be pretty simple here:

model_matched <- lm(final_grade ~ math_camp, data = math_camp_matched)
tidy(model_matched)

## # A tibble: 2 x 5
##   term               estimate std.error statistic  p.value
##   <chr>                 <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)          136.       0.343     398.  0.      
## 2 math_campMath camp     8.10     0.413      19.6 1.41e-74

The 8.1 point increase here is better than the naive estimate, but worse than the educated guess model. Perhaps that’s because the matches aren’t great, or maybe we threw away too much data. There are a host of diagnostics you can look at to see how well things are matched (check the documentation for MatchIt for examples.)

One nice thing about using matchit() is that it also generates a kind of weight based on the distance between points. You can incorporate those weights into the model and get a more accurate estimate:

model_matched_weighted <- lm(final_grade ~ math_camp, data = math_camp_matched, weights = weights)
tidy(model_matched_weighted)

## # A tibble: 2 x 5
##   term               estimate std.error statistic   p.value
##   <chr>                 <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)          135.       0.336     400.  0.       
## 2 math_campMath camp     9.83     0.405      24.3 4.16e-106

After using the weights, we find a 9.83 point causal effect. That’s much better than any of the other estimates we’ve tried.

Adjustment with inverse probability weighting

One potential downside to matching is that you generally have to throw away a sizable chunk of your data—anything that’s unmatched doesn’t get included in the final matched data.

An alternative approach to matching is to assign every observation some probability of receiving treatment, and then weighting each observation by its inverse probability—observations that are predicted to get treatment and then don’t, or observations that are predicted to not get treatment and then do will receive more weight than the observations that get/don’t get treatment as predicted.

Generating these inverse probability weights requires a two step process: (1) you first generate propensity scores, or the probability of receiving treatment, and then (2) you use a special formula to convert those propensity scores into weights. Once you have weights, you can incorporate them into your regression model like we did above with the matched and weighted data.

ggplot(math_camp, aes(x = gre_quant, y = math_camp_num)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = "lm", color = "red", size = 0.5) +
  geom_smooth(method = "glm", method.args = list(family = binomial(link = "logit"))) +
  labs(x = "Quantitative GRE score", y = "Probability of attending math camp")

# Notice how we use glm() instead of lm(). The "g" stands for "generalized"
# linear model. We have to specify a family in any glm() model. You can
# technically run a regular OLS model (like you do with lm()) if you use 
# glm(y ~ x1 + x2, family = gaussian(link = "identity")), but people rarely do that.
#
# To use logistic regression, you have to specify a binomial/logit family like so:
# family = binomial(link = "logit")
model_logit <- glm(math_camp ~ gre_quant, data = math_camp,
                   family = binomial(link = "logit"))

tidy(model_logit)

## # A tibble: 2 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   23.4      1.59        14.7 5.27e-49
## 2 gre_quant     -0.165    0.0110     -14.9 2.34e-50

tidy(model_logit, exponentiate = TRUE)

## # A tibble: 2 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept) 1.51e+10    1.59        14.7 5.27e-49
## 2 gre_quant   8.48e- 1    0.0110     -14.9 2.34e-50

needs_camp_model <- glm(math_camp ~ undergrad_gpa + gre_quant, data = math_camp, 
                        family = binomial(link = "logit"))

# We could look at these results if we wanted, but we don't need to for this class
# tidy(needs_camp_model, exponentiate = TRUE)

# augment_columns() handles the plugging in of values. You need to feed it the
# name of the model and the name of the dataset you want to add the predictions
# to. The type.predict = "response" argument makes it so the predictions are in
# the 0-1 scale. If you don't include that, you'll get predictions in an
# uninterpretable log odds scale.
math_camp_propensities <- augment_columns(needs_camp_model, math_camp, 
                                          type.predict = "response") %>% 
  rename(p_camp = .fitted)

# Look at the first few rows of a few columns
math_camp_propensities %>% 
  select(id, math_camp, undergrad_gpa, gre_quant, p_camp) %>% 
  head()

## # A tibble: 6 x 5
##      id math_camp    undergrad_gpa gre_quant p_camp
##   <dbl> <fct>                <dbl>     <dbl>  <dbl>
## 1     1 No math camp          3.9        151 0.0996
## 2     2 No math camp          3.2        143 0.373 
## 3     3 Math camp             2.83       140 0.563 
## 4     4 No math camp          2.63       154 0.301 
## 5     5 No math camp          3.24       148 0.258 
## 6     6 No math camp          2.95       146 0.381

math_camp_ipw <- math_camp_propensities %>% 
  mutate(ipw = (math_camp_num / p_camp) + ((1 - math_camp_num) / (1 - p_camp)))

# Look at the first few rows of a few columns
math_camp_ipw %>% 
  select(id, math_camp, undergrad_gpa, gre_quant, p_camp, ipw) %>% 
  head()

## # A tibble: 6 x 6
##      id math_camp    undergrad_gpa gre_quant p_camp   ipw
##   <dbl> <fct>                <dbl>     <dbl>  <dbl> <dbl>
## 1     1 No math camp          3.9        151 0.0996  1.11
## 2     2 No math camp          3.2        143 0.373   1.60
## 3     3 Math camp             2.83       140 0.563   1.78
## 4     4 No math camp          2.63       154 0.301   1.43
## 5     5 No math camp          3.24       148 0.258   1.35
## 6     6 No math camp          2.95       146 0.381   1.62

model_ipw <- lm(final_grade ~ math_camp, 
                data = math_camp_ipw, weights = ipw)
tidy(model_ipw)

## # A tibble: 2 x 5
##   term               estimate std.error statistic   p.value
##   <chr>                 <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)           137.      0.223     613.  0.       
## 2 math_campMath camp     10.9     0.308      35.3 8.50e-213

math_camp_ipw <- math_camp_ipw %>% 
  mutate(ipw_truncated = ifelse(ipw > 10, 10, ipw))

model_ipw_truncated <- lm(final_grade ~ math_camp, 
                          data = math_camp_ipw, weights = ipw_truncated)
tidy(model_ipw_truncated)

## # A tibble: 2 x 5
##   term               estimate std.error statistic   p.value
##   <chr>                 <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)           137.      0.219     622.  0.       
## 2 math_campMath camp     10.6     0.306      34.7 4.30e-207

Comparison of all results

Let’s compare all the ATEs that we just calculated: