Aggregation

aggregate_estimands() averages the treatment-group-specific estimates of multiple_treatment_group_analysis() across treatment groups, by event time (Section 5 of Leventer 2025, Identification of Child Penalties). Symbols and row labels are in Notation and output, the estimators and influence functions in Estimation. Equation, section and theorem numbers refer to the main text; appendices are referred to by name.

The examples use one simulated panel. Treatment groups 24-30 at event times 0-3 use control groups up to 34, so every treatment group has a cell at every event time.

library(childpen)
sim <- simulate_data(20000, 24:34)
res <- multiple_treatment_group_analysis(sim, 24:30, periods_post = 3,
                                         periods_pre = NULL, verbose = FALSE)

The result carries the person-level influence functions that aggregate_estimands() reads, so it is passed on unmodified.

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:

method	estimand	agg_type	aggregate	$B_d$
DID_Female, DID_Male	theta	weighted_avg	$\theta_{\mathrm{Agg},1}(g, e)$	$\delta_\theta(g)$
DID_Female, DID_Male	theta	ratio_of_avgs	$\theta_{\mathrm{Agg},2}(g, e)$	$\delta_{\mathrm{ATE}}(g)$ and $\delta_{\mathrm{APO}}(g)$, averaged separately and divided
NTD	rho_effect	weighted_avg	$\rho_{\mathrm{Agg}}(e)$	$\delta_{\Delta\rho}$
Conventional	gap_theta	weighted_avg	$\Delta\theta_{\mathrm{Agg}}(e)$	$\delta_{\Delta\theta}$
TD_E_LOGS	rho_effect	weighted_avg	$\rho_{\mathrm{Agg}}(e)$ under TD-E-LOGS	$\delta_{\Delta\rho}$ with $\delta^{\times}_{\mathrm{APO}}$
TD, TD_PR, accounting rows of NTD and TD_E_LOGS	gap_ATE, gap_theta_pre, rho_obs, rho_cf, share_parenthood, share_other	weighted_avg	weighted average of the group-level row; no symbol in the paper	the group-level row

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

weights	weight of treatment group $d$ in $\mathcal{D}(e)$
"sample" (default)	the share of persons in the group: $w_d = n_d / \sum_{k \in \mathcal{D}(e)} n_k$, where $n_d$ counts persons of both genders with $D = d$. For DID_Female and DID_Male, shares within gender: $w_{g,d} = n_{g,d} / \sum_{k \in \mathcal{D}(e)} n_{g,k}$
"uniform" or NULL	$1/m_e$, with $m_e$ the number of groups in $\mathcal{D}(e)$
a named numeric vector	proportional to the entry named $d$, renormalized to sum to one over $\mathcal{D}(e)$. Groups missing from the vector get weight zero, with a warning

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.