Aggregation

1. Aggregate causal estimands

For event time \(e\), let \(\mathcal{D}(e) = \{d : d + e < D_{\max}\}\) be the treatment groups that have an estimate. Since the control group is \(d' = d + e + 1\), the condition is \(d' \le D_{\max}\), the bound set by max_control_group. With \(\mathbb{E}_D\) the expectation over the treatment distribution, Section 5.1 defines

\[ \rho_{\mathrm{Agg}}(e) = \mathbb{E}_D\big[\rho(D, D, D + e) - \rho(D, \infty, D + e) \mid D + e < D_{\max}\big] \qquad \text{(eq. (10))}, \] \[ \theta_{\mathrm{Agg},1}(g, e) = \mathbb{E}_D\big[\theta(g, D, D + e) \mid G = g,\, D + e < D_{\max}\big], \] \[ \theta_{\mathrm{Agg},2}(g, e) = \frac{\mathbb{E}_D\big[\mathrm{ATE}(g, D, D + e) \mid G = g,\, D + e < D_{\max}\big]}{\mathbb{E}_D\big[\mathrm{APO}(g, D, \infty, D + e) \mid G = g,\, D + e < D_{\max}\big]}, \] \[ \Delta\theta_{\mathrm{Agg}}(e) = \mathbb{E}_D\big[\theta(f, D, D + e) - \theta(m, D, D + e) \mid D + e < D_{\max}\big] \qquad \text{(eq. (11))}. \]

The package replaces each treatment-group causal estimand by its 2x2 descriptive counterpart \(B_d\) and the treatment distribution by weights \(w_d\) (next section), so that \(\hat A(e) = \sum_{d \in \mathcal{D}(e)} w_d \hat B_d\). The rows of aggregate_estimands() are identified by method, estimand and agg_type:

The output has one row per event time, control_offset, estimand, method and agg_type, with est, se, n_groups (treatment groups with positive weight), ci_l, ci_h, t and p.

agg <- aggregate_estimands(res)
agg |>
  filter(event_time == 2, estimand %in% c("theta", "gap_theta", "rho_effect"),
         method %in% c("DID_Female", "Conventional", "NTD")) |>
  select(method, estimand, agg_type, est, se, n_groups)
#>         method   estimand      agg_type     est       se n_groups
#> 1 Conventional  gap_theta  weighted_avg -0.2273 0.004128        7
#> 2   DID_Female      theta ratio_of_avgs -0.2990 0.002595        7
#> 3   DID_Female      theta  weighted_avg -0.3033 0.002533        7
#> 4          NTD rho_effect  weighted_avg -0.2222 0.004226        7

ratio_of_avgs divides the weighted average of the ATE rows by that of the APO rows. Section 5.1 shows that \(\theta_{\mathrm{Agg},2}\) is also an average of the \(\theta(g, d, d + e)\), with weights proportional to the treatment-group share times counterfactual earnings \(\mathrm{APO}(g, d, \infty, d + e)\).

2. Weights

Gender gaps are aggregated as in eq. (11): the within-group gap is averaged with weights pooled over genders, rather than differencing two within-gender aggregates, which would also reflect gender differences in the treatment distribution (Section 5.1). Only DID_Female and DID_Male use within-gender weights.

w_named <- setNames(dnorm(24:30, mean = 27, sd = 3), 24:30)
bind_rows(
  sample  = aggregate_estimands(res, weights = "sample",  methods = "NTD"),
  uniform = aggregate_estimands(res, weights = "uniform", methods = "NTD"),
  named   = aggregate_estimands(res, weights = w_named,   methods = "NTD"),
  .id = "weights"
) |>
  filter(estimand == "rho_effect", event_time %in% c(0, 3)) |>
  select(weights, event_time, est, se, n_groups)
#>   weights event_time     est       se n_groups
#> 1  sample          0 -0.2767 0.005244        7
#> 2  sample          3 -0.1970 0.004096        7
#> 3 uniform          0 -0.2843 0.005427        7
#> 4 uniform          3 -0.2044 0.004060        7
#> 5   named          0 -0.2840 0.005471        7
#> 6   named          3 -0.2045 0.004017        7

A named vector aggregates the same treatment-group estimates under another treatment distribution, as in Figure 9 of the paper (Replicating the paper).

Groups that drop out. A treatment group has a cell at event time \(e\) only if its control group \(d + e + 1\) is at most max_control_group (\(D_{\max}\) in the paper). Groups beyond the bound drop out at later event times, and n_groups records how many remain.

res_cap <- multiple_treatment_group_analysis(sim, 24:30, periods_post = 3, periods_pre = NULL,
                                             max_control_group = 32, verbose = FALSE)
bind_rows(
  `Inf` = aggregate_estimands(res,     methods = "NTD"),
  `32`  = aggregate_estimands(res_cap, methods = "NTD"),
  .id = "max_control_group"
) |>
  filter(estimand == "rho_effect") |>
  select(max_control_group, event_time, est, n_groups)
#>   max_control_group event_time     est n_groups
#> 1               Inf          0 -0.2767        7
#> 2               Inf          1 -0.2516        7
#> 3               Inf          2 -0.2222        7
#> 4               Inf          3 -0.1970        7
#> 5                32          0 -0.2767        7
#> 6                32          1 -0.2516        7
#> 7                32          2 -0.2360        6
#> 8                32          3 -0.2148        5

With max_control_group = 32, group 30 has cells at \(e = 0, 1\) only and group 29 at \(e \le 2\), so \(\mathcal{D}(e)\) shrinks and the weights change with \(e\). As the footnote to eq. (10) notes (following Callaway and Sant’Anna 2021), a comparison of aggregates across \(e\) then combines variation in causal effects with a shift in the composition of treatment groups. The composition is the same at every event time when each group in the aggregate has a cell at each of them, as with max_control_group = Inf above (n_groups is 7 throughout).

3. Standard errors

Let \(B_d\) be the estimate of treatment group \(d\), \(p_d = \Pr(D = d)\) its population share, \(P_e = \sum_{d \in \mathcal{D}(e)} p_d\) and \(w_d = p_d / P_e\). The influence function of \(A(e) = \sum_{d \in \mathcal{D}(e)} w_d B_d\) is (the aggregate estimator section of the appendix on estimation and inference) \[ \psi_{A(e)} = \frac{1}{P_e} \sum_{d \in \mathcal{D}(e)} \Big[ \big(B_d - A(e)\big)\,\psi_{p_d} + p_d\,\psi_{B_d} \Big], \qquad \psi_{p_d} = 1\{D = d\} - p_d . \] For the within-gender rows the share is the joint share of gender \(g\) and group \(d\), with \(1\{G = g, D = d\} - p_{g,d}\) in place of \(\psi_{p_d}\). With "uniform" or fixed weights the first term is absent; with "sample" it accounts for the estimated shares. For ratio_of_avgs, \(\hat\theta_{\mathrm{Agg},2} = \hat N / \hat M\) with \(\hat N = \sum_d w_d\,\hat\delta_{\mathrm{ATE},d}\) and \(\hat M = \sum_d w_d\,\hat\delta_{\mathrm{APO},d}\), and by the chain rule (not in the paper) \[ \psi_{\theta_{\mathrm{Agg},2}} = \frac{\psi_N}{M} - \frac{N\,\psi_M}{M^2}, \] where \(\psi_N\) and \(\psi_M\) are the influence functions of the two weighted averages, each with its share term. The influence function of share_other is minus that of share_parenthood.

The group-level influence functions are summed per person across treatment groups before they are squared, with clusters at the person level, \[ \widehat{\mathrm{Var}}(\hat A) = \frac{1}{n^2} \sum_{c} \Big( \sum_{i \in c} \hat\psi_i \Big)^2 \frac{G}{G - 1}, \] for \(n\) person-age rows and \(G\) persons. A person who serves as a control group member for several treatment groups therefore enters several terms of the same sum. The share term is put on the person-age scale with the average number of ages per person, which is exact for a balanced panel. The estimators and the influence functions of the group-level rows are in Estimation.