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.
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.571Input 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.
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}\).
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.
| 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).