Estimation

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

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

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:

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)
#>    event_time control_offset       method         estimand     est     se
#> 1           1              2   DID_Female              APO   24900    377
#> 2           1              2     DID_Male              APO   26200    388
#> 3           1              2   DID_Female              ATE   -8340    327
#> 4           1              2     DID_Male              ATE   -1160    361
#> 5           1              2   DID_Female            theta  -0.335 0.0095
#> 6           1              2     DID_Male            theta -0.0444 0.0134
#> 7           1              2           TD          gap_ATE   -7170    487
#> 8           1              2 Conventional        gap_theta   -0.29 0.0164
#> 9           1              2          NTD       rho_effect  -0.289 0.0176
#> 10          1              2    TD_E_LOGS       rho_effect  -0.288 0.0202
#> 11          1              2        TD_PR    gap_theta_pre   -0.34 0.0241
#> 12          1              2          NTD          rho_obs   0.662 0.0124
#> 13          1              2          NTD           rho_cf   0.951 0.0202
#> 14          1              2          NTD share_parenthood   0.856 0.0572
#> 15          1              2          NTD      share_other   0.144 0.0572
#> 16          1              2    TD_E_LOGS          rho_obs   0.662 0.0124
#> 17          1              2    TD_E_LOGS           rho_cf    0.95 0.0235
#> 18          1              2    TD_E_LOGS share_parenthood   0.853  0.067
#> 19          1              2    TD_E_LOGS      share_other   0.147  0.067

pre <- single_treatment_group_analysis(sim, d = 26, dp = 28, a = 23)
show(pre)
#>    event_time control_offset       method         estimand     est     se
#> 1          -3              2   DID_Female              APO   17400    346
#> 2          -3              2     DID_Male              APO   17700    361
#> 3          -3              2   DID_Female              ATE     213    289
#> 4          -3              2     DID_Male              ATE     492    324
#> 5          -3              2   DID_Female            theta  0.0122 0.0168
#> 6          -3              2     DID_Male            theta  0.0278 0.0187
#> 7          -3              2           TD          gap_ATE    -280    434
#> 8          -3              2 Conventional        gap_theta -0.0155 0.0251
#> 9          -3              2          NTD       rho_effect -0.0148  0.024
#> 10         -3              2    TD_E_LOGS       rho_effect -0.0156 0.0218
#> 11         -3              2        TD_PR    gap_theta_pre -0.0132 0.0184
#> 12         -3              2          NTD          rho_obs   0.963 0.0185
#> 13         -3              2          NTD           rho_cf   0.978 0.0279
#> 14         -3              2          NTD share_parenthood   0.403  0.656
#> 15         -3              2          NTD      share_other   0.597  0.656
#> 16         -3              2    TD_E_LOGS          rho_obs   0.963 0.0185
#> 17         -3              2    TD_E_LOGS           rho_cf   0.979  0.024
#> 18         -3              2    TD_E_LOGS share_parenthood   0.427  0.571
#> 19         -3              2    TD_E_LOGS      share_other   0.573  0.571

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.

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