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

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
options(scipen = 999)
data.table::setDTthreads(1L)
```

This page describes how the package runs the paper's 2x2 comparisons, then writes out every
estimator it returns and the influence function its standard error is built from. Notation
follows [Notation and output](notation.html). Section, theorem and appendix numbers refer to Leventer (2025),
*Identification of Child Penalties*; appendices are referred to by name, because their letters
differ between the arXiv and journal versions: *Estimation appendix* is the appendix on
estimation and inference, *Proportional appendix* the online appendix on proportional DID and
TD. Formulas that are not in the paper are marked **not in the paper**.

## Running the comparisons

`multiple_treatment_group_analysis()` runs the grid of 2x2 comparisons that the paper's figures
use. For each treatment group $d$ in `treatment_groups` and event time $e = 0, \dots, H$
(`periods_post = H`), the target age is $a = d + e$ and the control group is the closest
not-yet-treated group, $d' = a + 1$ (Section 1.3). Pre-treatment event times
$e = -(\text{pre} + K), \dots, -(\text{pre} + 1)$ (`periods_pre = K`) are estimated against every
control group used after treatment, $d' = d + 1, \dots, d + H + 1$: the identifying assumptions are
about each treatment-control pair, so each pair has its own validation test (Section 2.3). The
column `control_offset` $= d' - d$ identifies the pair. `max_control_group` is the largest age at
first birth used as a control (the paper's $D_{max}$). `single_treatment_group_analysis()` runs
one comparison.

```{r grid, eval = FALSE}
res <- multiple_treatment_group_analysis(data, treatment_groups = 24:34,
                                         periods_post = 5, periods_pre = 3,
                                         max_control_group = 40)
```

With the paper's sample this grid is the one behind its figures: treatment groups 24 to 34,
event times 0 to 5, pre-treatment event times $-4$ to $-2$ (baseline $d - 1$), and controls up to
40 ([Replicating the paper](replication.html)).

**One comparison.** `single_treatment_group_analysis()` returns the rows of one cell
$(d, d', a)$. Treatment group $d = 26$ against control group $d' = 28$, after childbirth at age
$a = 27$ (event time 1) and before childbirth at age $a = 23$ (event time $-3$), both with baseline
age 25:

```{r one-cell}
library(childpen)
sim  <- simulate_data(5000, 24:30)
fmt  <- function(v) format(signif(v, 3), scientific = FALSE, drop0trailing = TRUE, trim = TRUE)
show <- function(x) transform(x[c("event_time", "control_offset", "method", "estimand", "est", "se")],
                               est = fmt(est), se = fmt(se))

post <- single_treatment_group_analysis(sim, d = 26, dp = 28, a = 27)
show(post)

pre <- single_treatment_group_analysis(sim, d = 26, dp = 28, a = 23)
show(pre)
```

**Input data.** `data` is a data.frame or data.table in long format, one row per person and age;
the panel need not be balanced. The columns are located by name:

| argument | default | content |
|---|---|---|
| `id_name` | `"id"` | person identifier (the clustering unit) |
| `age_name` | `"age"` | age, numeric |
| `D_name` | `"D"` | age at first birth, constant within person; a missing (or infinite) value means the person never has a child |
| `female_name` | `"female"` | gender, constant within person; `1` / `TRUE` is female and `0` / `FALSE` male; labels such as `"female"` / `"male"`, `"f"` / `"m"`, `"woman"` / `"man"` are also accepted |
| `Y_name` | `"Y"` | earnings, numeric; zeros are allowed, while the ratio rows need positive cell means |

Rows with a missing outcome are dropped with a message. The functions stop with an error if a
person-age pair is repeated, if `female` or `D` varies within a person, if a gender label is not
recognised, or if only one gender is present. `data` is not modified. Persons with a missing `D`
stay in the data, belong to no treatment or control group, and count in the totals the standard
errors are scaled by.

## Setup

Data are one row per person $i$ and age $t$, with earnings $Y$, gender $g \in \{f, m\}$ and age at
first birth $D$. A 2x2 comparison has a treatment group $d$, a control group $d' > a$, a target
age $a$ and a baseline age $b = d - \text{pre}$ (default $b = d - 1$).

**Cell means.** For a cell $(g, k, t)$, with indicator $S_{it} = 1\{G_i = g, D_i = k, \text{age} = t\}$,
$$
\hat\mu(g,k,t) = \frac{\sum_{i,t'} S_{it'} Y_{it'}}{\sum_{i,t'} S_{it'}},
\qquad
\psi_{\mu(g,k,t)} = \frac{S}{\hat P(S)}\,\big(Y - \hat\mu(g,k,t)\big),
$$
where $\hat P(S)$ is the share of person-age rows in the cell (Estimation appendix). Every estimator below
is a smooth function of eight cell means,
$$
\mu_1 = \mu(g,d,b),\quad \mu_2 = \mu(g,d',a),\quad \mu_3 = \mu(g,d',b),\quad \mu_4 = \mu(g,d,a),
\qquad g \in \{f, m\},
$$
and its influence function $\psi$ follows from theirs by the chain rule.

**Standard errors.** For an estimator with influence function $\psi$, with $n$ person-age rows,
$G$ persons and clusters $c$ = persons,
$$
\widehat{\mathrm{Var}} = \frac{1}{n^2}\sum_{c}\Big(\sum_{i \in c}\psi_i\Big)^2\,\frac{G}{G-1}.
$$
Confidence intervals are $\text{est} \pm 1.96\,\text{se}$.

## Rows of `single_treatment_group_analysis()`

**Imputations of counterfactual earnings.** Additive (DID) and multiplicative:
$$
\hat\delta_{APO}(g) = \mu_1 + \mu_2 - \mu_3,
\qquad
\psi_{\delta_{APO}(g)} = \psi_1 + \psi_2 - \psi_3;
$$
$$
\hat\delta^{\times}_{APO}(g) = \frac{\mu_1\,\mu_2}{\mu_3},
\qquad
\psi_{\delta^{\times}_{APO}(g)} = \hat\delta^{\times}_{APO}(g)\Big(\frac{\psi_1}{\mu_1} + \frac{\psi_2}{\mu_2} - \frac{\psi_3}{\mu_3}\Big).
$$

| method | estimand | estimator | influence function | identifies its target under | in the paper |
|---|---|---|---|---|---|
| `DID_Female`, `DID_Male` | `APO` | $\hat\delta_{APO}(g)$ | $\psi_{\delta_{APO}(g)}$ | NA, DID-PT | Lemma 1; eq. (3), Estimation app. |
| `DID_Female`, `DID_Male` | `ATE` | $\hat\delta_{ATE}(g) = \mu_4 - \hat\delta_{APO}(g)$ | $\psi_4 - \psi_{\delta_{APO}(g)}$ | NA, DID-PT | Lemma 1; eq. (3), Estimation app. |
| `DID_Female`, `DID_Male` | `theta` | $\hat\delta_\theta(g) = \hat\delta_{ATE}(g)/\hat\delta_{APO}(g)$ | $\dfrac{\psi_{\delta_{ATE}}}{\delta_{APO}} - \dfrac{\delta_{ATE}\,\psi_{\delta_{APO}}}{\delta_{APO}^2}$ | NA, DID-PT | Lemma 1; eq. (3), Estimation app. |
| `TD` | `gap_ATE` | $\hat\delta_{ATE}(f) - \hat\delta_{ATE}(m)$ | $\psi_{\delta_{ATE}(f)} - \psi_{\delta_{ATE}(m)}$ | NA, TD-PT | Lemma 1; Estimation app. |
| `Conventional` | `gap_theta` | $\hat\delta_\theta(f) - \hat\delta_\theta(m)$ | $\psi_{\delta_\theta(f)} - \psi_{\delta_\theta(m)}$ | NA, DID-PT; under NA and NTD-PT it equals $\text{Bias} \times [\theta(f) - \theta(m)]$ | Thm. 1; Estimation app. |
| `NTD` | `rho_obs` | $\hat\rho_{obs} = \mu_4(f)/\mu_4(m)$ | $\dfrac{\psi_4(f)}{\mu_4(m)} - \dfrac{\mu_4(f)\,\psi_4(m)}{\mu_4(m)^2}$ | consistency | Thm. 2; Estimation app. |
| `NTD` | `rho_cf` | $\hat\rho_{cf} = \hat\delta_{APO}(f)/\hat\delta_{APO}(m)$ | $\dfrac{\psi_{\delta_{APO}(f)}}{\delta_{APO}(m)} - \dfrac{\delta_{APO}(f)\,\psi_{\delta_{APO}(m)}}{\delta_{APO}(m)^2}$ | NA, NTD-PT | Thm. 2; Estimation app. |
| `NTD` | `rho_effect` | $\hat\rho_{obs} - \hat\rho_{cf}$ | $\psi_{\rho_{obs}} - \psi_{\rho_{cf}}$ | NA, NTD-PT | Thm. 2; Estimation app. |
| `TD_E_LOGS` | `rho_cf` | $\hat\rho^{\times}_{cf} = \hat\delta^{\times}_{APO}(f)/\hat\delta^{\times}_{APO}(m)$ | as `NTD` `rho_cf`, with $\delta^{\times}_{APO}$ | NA, TD-E-LOGS | Proportional app.; IF **not in the paper** |
| `TD_E_LOGS` | `rho_obs`, `rho_effect` | as `NTD`, with $\hat\rho^{\times}_{cf}$ | as `NTD`, with $\psi_{\rho^{\times}_{cf}}$ | `rho_obs`: consistency; `rho_effect`: NA, TD-E-LOGS | Proportional app.; IF **not in the paper** |
| `TD_PR` | `gap_theta_pre` | $\hat\theta_{pre}(f) - \hat\theta_{pre}(m)$, $\hat\theta_{pre}(g) = \dfrac{\mu_4 - \hat\delta^{\times}_{APO}(g)}{\mu_1}$ | $\psi_{\theta_{pre}(f)} - \psi_{\theta_{pre}(m)}$, $\psi_{\theta_{pre}(g)} = \dfrac{\psi_4 - \psi_{\delta^{\times}_{APO}(g)}}{\mu_1} - \dfrac{\theta_{pre}(g)\,\psi_1}{\mu_1}$ | NA, TD-PR | Proportional app.; IF **not in the paper** |
| `NTD`, `TD_E_LOGS` | `share_parenthood` | $\hat s = \dfrac{\hat\rho_{cf} - \hat\rho_{obs}}{1 - \hat\rho_{obs}}$ | $\dfrac{\psi_{\rho_{cf}} - \psi_{\rho_{obs}}}{1 - \rho_{obs}} + \dfrac{(\rho_{cf} - \rho_{obs})\,\psi_{\rho_{obs}}}{(1 - \rho_{obs})^2}$ | as `rho_cf` of the same method | Sec. 4.2 (estimate); **not in the paper** (IF) |
| `NTD`, `TD_E_LOGS` | `share_other` | $1 - \hat s$ | $-\psi_s$ | as `rho_cf` of the same method | **not in the paper** |

NA (no anticipation), DID-PT (parallel trends), TD-PT and NTD-PT are the assumptions of Sections 2
and 3 of the paper; TD-E-LOGS and TD-PR are those of the Proportional appendix. For `TD_E_LOGS`
the counterfactual ratio uses $\hat\delta^{\times}_{APO}$; the `share_parenthood` row of
`TD_E_LOGS` uses $\hat\rho^{\times}_{cf}$ in place of $\hat\rho_{cf}$.

The standard errors of `share_parenthood` and `share_other` are **not in the paper**, which
reports parenthood's share as a point estimate (Section 4.2); the influence function in the table
is the chain rule applied to $s = (\rho_{cf} - \rho_{obs})/(1 - \rho_{obs})$.

## `bias_bounds()`

The bias-corrected conventional gap at an assumed fathers' effect $\theta(m)$ (Proposition 2, in
the form of the Estimation appendix), with $\hat\Delta$ the `NTD` `rho_effect` estimate:
$$
\hat B(\theta(m)) = \hat\Delta\,\frac{\hat\delta_{APO}(m)}{\hat\delta_{APO}(f)}\,\big(1 + \theta(m)\big),
$$
$$
\psi_{B} = \big(1 + \theta(m)\big)\Big[\frac{\delta_{APO}(m)}{\delta_{APO}(f)}\psi_{\Delta}
 + \frac{\Delta}{\delta_{APO}(f)}\psi_{\delta_{APO}(m)}
 - \frac{\Delta\,\delta_{APO}(m)}{\delta_{APO}(f)^2}\psi_{\delta_{APO}(f)}\Big].
$$
The estimator is in the paper; its influence function is **not in the paper**. The standard error
scales with $1 + \theta(m)$.

## `aggregate_estimands()`

The aggregate estimators, their influence functions and standard errors are in
[Aggregation](aggregation.html).

## Where this lives in the code

| quantity | code |
|---|---|
| cell means and their IF; clustered variance | `R/utility.R`: `compute_mean_if()`, `se_cluster()` |
| all rows of the table above | `R/single_treatment_group_analysis.R` |
| aggregation | `R/aggregate_estimands.R` |
| bias-corrected gap | `R/bias_bounds.R` |

Standard errors checked against a person-level (cluster) bootstrap in `tests/testthat/`:
`TD_E_LOGS` and `TD_PR` `rho_effect` / `gap_theta_pre`, `NTD` `share_parenthood`, and
`bias_bounds()`. `TD_E_LOGS` and `TD_PR` are also checked by Monte Carlo
(`dev/sims/td_sim_results.md` in the source repository).