---
title: "Aggregation"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Aggregation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 8, fig.height = 5,
                      message = FALSE, warning = FALSE)
library(dplyr)
data.table::setDTthreads(1L)
options(scipen = 999, digits = 4)
```

`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*](https://arxiv.org/abs/2602.07486)). Symbols and row labels are in [Notation and output](notation.html), the estimators and influence functions in [Estimation](estimators.html). 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.

```{r estimate}
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`.

```{r agg-rows}
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)
```

`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.

```{r 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)
```

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

**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.

```{r dmax}
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)
```

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](estimators.html).