library(dplyr)
#> Warning: package 'dplyr' was built under R version 4.5.2
set.seed(1)
N <- 60; TT <- 25
df <- expand.grid(id = 1:N, time = 1:TT)
treat_start <- sample(c(NA, 8:18), N, replace = TRUE)
df$D <- ifelse(is.na(treat_start[df$id]) | df$time < treat_start[df$id],
0, 1)
## Higher intercept (2.0) and smaller noise (sd = 0.1) keep Y safely
## positive throughout, so the v2.4.2+ cell-drop hard-error in
## log.att / aptt does not fire on benign bootstrap noise.
df$Y <- exp(2.0 + 0.05 * df$time + 0.3 * df$D + rnorm(nrow(df), sd = 0.1))
fit <- fect(Y ~ D, data = df, index = c("id", "time"),
method = "fe", force = "two-way",
se = TRUE, nboots = 1000, parallel = FALSE,
keep.sims = TRUE)3 Alternative Estimands
The default output of fect() is the per-event-time average treatment effect on the treated (ATT). This chapter introduces the estimand interface added in v2.4.0. It has two functions: estimand(), which computes alternative summaries directly from a fitted model, and imputed_outcomes(), which returns the cell-level imputed potential outcomes so you can build summaries that estimand() does not include.
estimand() covers four common requests in one line each: cumulative ATT, average proportional treatment effect on the treated (APTT, Chen and Roth 2024), log-scale ATT, and ATT restricted to a window of event times. Each call returns a data frame ready for plotting with fect’s esplot(). For anything else, imputed_outcomes() returns the cell-level data in long format and the rest is standard data manipulation.
R script used in this chapter can be downloaded here.
3.1 Beyond the default ATT
Use estimand() when you want any of the following: the per-event-time level ATT as a data frame (instead of the legacy matrix); the cumulative ATT through each event time, or summed over a fixed window; APTT (Chen and Roth 2024); log-scale ATT; or a level ATT restricted to a window of event times — say, the first five post-treatment periods.
For summaries that estimand() does not cover, imputed_outcomes() returns the cell-level data in long format and you do the aggregation yourself.
3.2 Setup
The examples below use a small synthetic panel with positive outcomes so that log-scale calculations are well defined. The fit sets keep.sims = TRUE so that the cell-level imputed counterfactuals from each bootstrap replicate are saved on the fit object. APTT and log-scale ATT need these saved values to recompute their non-linear formulas inside each replicate. The level ATT and cumulative ATT use the pre-aggregated bootstrap that fect always keeps, so the default keep.sims is fine for those.
3.3 Per-event-time level ATT
estimand(fit, "att", "event.time") |>
head(8)
#> event.time estimate se ci.lo ci.hi n_cells vartype
#> 1 -16 0.34239381 0.3361301 -0.3164091 1.0011967 4 bootstrap
#> 2 -15 0.36167880 0.2798964 -0.1869080 0.9102656 6 bootstrap
#> 3 -14 -0.05907245 0.2214799 -0.4931650 0.3750201 14 bootstrap
#> 4 -13 0.04607633 0.1546846 -0.2570999 0.3492526 22 bootstrap
#> 5 -12 -0.13107862 0.1375386 -0.4006494 0.1384921 25 bootstrap
#> 6 -11 0.05389936 0.1371076 -0.2148265 0.3226252 31 bootstrap
#> 7 -10 -0.09926449 0.1211865 -0.3367856 0.1382566 39 bootstrap
#> 8 -9 -0.10154200 0.1265032 -0.3494837 0.1463997 45 bootstrapEach row is one event time. The columns hold the point estimate, the standard error, the bounds of the confidence interval, the number of treated cells that contributed at that event time, and the variance type used for inference. The numbers match the legacy fit$est.att exactly. The only difference is the data-frame format, which works directly with downstream plotting and aggregation tools — including fect’s own esplot():
est <- estimand(fit, "att", "event.time")
esplot(est, Period = "event.time",
Estimate = "estimate", CI.lower = "ci.lo", CI.upper = "ci.hi",
Count = "n_cells",
main = "Per-event-time ATT")
3.4 Cumulative ATT
estimand(fit, "att.cumu", "event.time") |>
head(8)
#> event.time estimate se ci.lo ci.hi n_cells vartype
#> 1 1 4.610344 0.3146759 4.007585 5.261076 57 bootstrap
#> 2 2 9.904908 0.5322089 8.899812 11.007187 57 bootstrap
#> 3 3 15.444577 0.7849516 13.983251 17.016717 57 bootstrap
#> 4 4 21.857445 1.1807527 19.766302 24.304471 57 bootstrap
#> 5 5 27.611710 1.5039713 24.843169 30.685348 57 bootstrap
#> 6 6 34.241599 1.9335062 30.533162 38.168534 57 bootstrap
#> 7 7 41.299255 2.2068935 37.096000 45.662127 57 bootstrap
#> 8 8 48.151263 2.4887683 43.529078 53.168360 57 bootstrapThe cumulative ATT through each event time. Each row reports the running sum of the per-period average effect, computed cell by cell. This call replaces the legacy effect() function and produces identical numbers.
cumu <- estimand(fit, "att.cumu", "event.time")
esplot(cumu, Period = "event.time",
Estimate = "estimate", CI.lower = "ci.lo", CI.upper = "ci.hi",
Count = "n_cells",
main = "Cumulative ATT",
ylab = "Cumulative Effect")
For a single overall window, request the overall grouping axis instead:
estimand(fit, "att.cumu", "overall", window = c(1, 5))
#> estimate se ci.lo ci.hi n_cells vartype
#> 1 27.61171 1.503971 24.84317 30.68535 285 bootstrapThe result is the cumulative effect over the first five post-treatment periods. It matches the final row of the legacy att.cumu(fit, period = c(1, 5)) output.
3.5 APTT
estimand(fit, "aptt", "event.time") |>
head(8)
#> event.time estimate se ci.lo ci.hi n_cells vartype
#> 1 1 0.3199724 0.02636380 0.2683002 0.3716445 57 bootstrap
#> 2 2 0.3495009 0.02915542 0.2923574 0.4066445 57 bootstrap
#> 3 3 0.3476220 0.02910081 0.2905855 0.4046585 57 bootstrap
#> 4 4 0.3863353 0.03372756 0.3202305 0.4524401 57 bootstrap
#> 5 5 0.3269021 0.03299593 0.3084222 0.3743267 57 bootstrap
#> 6 6 0.3598550 0.03513667 0.3370717 0.4322788 57 bootstrap
#> 7 7 0.3669511 0.02855493 0.3558770 0.4165444 57 bootstrap
#> 8 8 0.3384248 0.02982154 0.2907143 0.4130177 57 bootstrapThe average proportional treatment effect on the treated (Chen and Roth 2024):
\[ \text{APTT}_t \;=\; \frac{\mathbb{E}[\,Y - \widehat{Y}(0) \mid D = 1, t\,]}{\mathbb{E}[\,\widehat{Y}(0) \mid D = 1, t\,]}. \]
Inside each bootstrap replicate, both the numerator and the denominator are recomputed and then divided. The bootstrap distribution is therefore a distribution of ratios, not a ratio of two distributions of means. This requires the cell-level imputed counterfactuals, which is why APTT needs keep.sims = TRUE at fit time.
aptt <- estimand(fit, "aptt", "event.time")
esplot(aptt, Period = "event.time",
Estimate = "estimate", CI.lower = "ci.lo", CI.upper = "ci.hi",
Count = "n_cells",
main = "APTT by event time",
ylab = "Average Proportional TE on the Treated")
3.6 Log-scale ATT
estimand(fit, "log.att", "event.time") |>
head(8)
#> event.time estimate se ci.lo ci.hi n_cells vartype
#> 1 1 0.2707030 0.02347225 0.2246983 0.3167078 57 bootstrap
#> 2 2 0.2921348 0.02577990 0.2416071 0.3426624 57 bootstrap
#> 3 3 0.2943029 0.02479570 0.2457042 0.3429015 57 bootstrap
#> 4 4 0.3163869 0.02750669 0.2624748 0.3702990 57 bootstrap
#> 5 5 0.2784879 0.02844118 0.2398725 0.3471718 57 bootstrap
#> 6 6 0.2995789 0.02932657 0.2606731 0.3679785 57 bootstrap
#> 7 7 0.3066786 0.02433966 0.2747105 0.3727985 57 bootstrap
#> 8 8 0.2850468 0.02517073 0.2281264 0.3215687 57 bootstrapThe mean log-scale treatment effect: at each event time, the average over treated cells of \(\log Y_{it} - \log \widehat{Y}_{it}(0)\). Cells where either the observed outcome or the imputed counterfactual is non-positive are dropped from the point estimate, since the logarithm is undefined there.
When a cell used in the point estimate has \(\widehat{Y}_{it}(0)_b \le 0\) in more than 5% of bootstrap replicates, estimand("log.att", ...) errors with guidance on what to do (filter the unstable cells, transform the outcome, or use a different estimand). The example data above uses an exponential-link outcome with an intercept high enough that no cell’s bootstrap distribution crosses zero.
log_att <- estimand(fit, "log.att", "event.time")
esplot(log_att, Period = "event.time",
Estimate = "estimate", CI.lower = "ci.lo", CI.upper = "ci.hi",
Count = "n_cells",
main = "Log-scale ATT",
ylab = "Mean log(Y) − log(Y0)")
3.7 Placebo tests
Starting in v2.4.2, estimand() takes a test = c("none", "placebo", "carryover") argument that lets you compute placebo and carryover versions of any alternative estimand, not just the default att.placebo / att.carryover scalars (see issue #131). The placebo case requires a fit produced with placeboTest = TRUE and a placebo period.
## Re-fit with placeboTest = TRUE so that estimand() can recompute
## APTT / log-ATT at the placebo cells.
fit_placebo <- fect(Y ~ D, data = df, index = c("id", "time"),
method = "fe", force = "two-way",
se = TRUE, nboots = 200, parallel = FALSE,
keep.sims = TRUE,
placeboTest = TRUE, placebo.period = c(-2, 0))The placebo APTT series should be statistically indistinguishable from zero across the placebo window if the parallel-counterfactual assumption holds; non-zero placebo estimates are evidence against identification.
aptt_placebo <- estimand(fit_placebo, "aptt", "event.time",
test = "placebo")
#> Warning: estimand() with ci.method = "bca" on a fit with nboots = 200 (< 1000):
#> tail-quantile-based CIs may have erratic endpoints at this replicate count
#> (Efron 1987 Section 3; DiCiccio & Efron 1996 Section 4 recommend B >= 1000).
#> The point estimate and SE are unaffected. For publication-grade CIs, refit with
#> `fect(..., nboots = 1000)` and re-call estimand().
esplot(aptt_placebo, Period = "event.time",
Estimate = "estimate", CI.lower = "ci.lo", CI.upper = "ci.hi",
Count = "n_cells",
main = "APTT placebo (pre-treatment)",
ylab = "Average Proportional TE on the Treated")
The call gives a warning that BCa confidence intervals with nboots = 200 may be noisy at the endpoints. This is a caution, not an error. BCa uses the 2.5% and 97.5% quantiles of the bootstrap distribution to form the CI, and these tail quantiles can vary with only 200 bootstrap runs. Efron (1987) and DiCiccio and Efron (1996) recommend at least 1000 runs for publication-quality tail CIs. The point estimate and standard error are not affected by nboots. This example uses 200 to keep the chapter fast. For more stable CIs in your own analysis, refit with fect(..., nboots = 1000).
The same logic carries over to log-scale ATT under placebo:
log_att_placebo <- estimand(fit_placebo, "log.att", "event.time",
test = "placebo")
#> Warning: estimand() with ci.method = "bca" on a fit with nboots = 200 (< 1000):
#> tail-quantile-based CIs may have erratic endpoints at this replicate count
#> (Efron 1987 Section 3; DiCiccio & Efron 1996 Section 4 recommend B >= 1000).
#> The point estimate and SE are unaffected. For publication-grade CIs, refit with
#> `fect(..., nboots = 1000)` and re-call estimand().
esplot(log_att_placebo, Period = "event.time",
Estimate = "estimate", CI.lower = "ci.lo", CI.upper = "ci.hi",
Count = "n_cells",
main = "Log-ATT placebo (pre-treatment)",
ylab = "Mean log(Y) − log(Y0)")
3.8 Carryover tests
The carryover case requires carryoverTest = TRUE and a panel with treatment reversals. The carryover series uses cells in fit$carryover.period—the early post-reversal window where, if the treatment effect dissipates immediately, estimates should be zero. Non-zero carryover estimates suggest the effect persists past treatment removal. We first build a panel with treatment reversals (each treated unit’s spell ends after a uniformly-drawn 5–10 periods) and visualize the on/off pattern with panelview():
## Build a panel with treatment reversals for the carryover demo.
set.seed(2)
df_rev <- df
treat_end <- pmin(treat_start[df_rev$id] + sample(5:10, N, replace = TRUE),
TT + 1L)
df_rev$D <- ifelse(is.na(treat_start[df_rev$id]) |
df_rev$time < treat_start[df_rev$id] |
df_rev$time >= treat_end[df_rev$id],
0, 1)
df_rev$Y <- exp(2.0 + 0.05 * df_rev$time + 0.3 * df_rev$D +
rnorm(nrow(df_rev), sd = 0.1))
aptt_carry <- estimand(fit_carry, "aptt", "event.time",
test = "carryover")
esplot(aptt_carry, Period = "event.time",
Estimate = "estimate", CI.lower = "ci.lo", CI.upper = "ci.hi",
Count = "n_cells",
main = "APTT carryover (post-reversal)",
ylab = "Average Proportional TE on the Treated")
log_att_carry <- estimand(fit_carry, "log.att", "event.time",
test = "carryover")
esplot(log_att_carry, Period = "event.time",
Estimate = "estimate", CI.lower = "ci.lo", CI.upper = "ci.hi",
Count = "n_cells",
main = "Log-ATT carryover (post-reversal)",
ylab = "Mean log(Y) − log(Y0)")
type = "att.cumu" is rejected with a clear error when test != "none". Cumulative ATT is defined relative to treatment onset, so it has no meaning at placebo or post-reversal cells.
3.9 Window-restricted ATT
Pass a window argument to focus on a subset of event times. The function keeps cells whose event time falls inside the closed interval and computes the requested estimand over them.
estimand(fit, "att", "overall", window = c(1, 5))
#> estimate se ci.lo ci.hi n_cells vartype
#> 1 5.522342 0.3049967 4.924559 6.120125 285 bootstrapThe result is a single overall ATT averaged over treated cells with event time between one and five. Use this when the question is “what is the average ATT for the first five post-treatment years?” and you do not need a per-event-time series. The answer is one number with one confidence interval, so a series plot is not relevant; if you want to compare multiple windows, build a small data frame of (window_label, estimate, ci.lo, ci.hi) and pass it to esplot() with Period = "window_label".
3.10 Custom estimands
If the summary you need is not one of the four built-in types, you can compute it yourself from the cell-level data:
po <- imputed_outcomes(fit)
head(po)
#> id time event.time cohort treated Y_obs Y0_hat eff eff_debias W.agg
#> 1 1 15 1 15 TRUE 21.49116 16.03611 5.455053 0 1
#> 2 1 16 2 15 TRUE 22.32793 16.80319 5.524740 0 1
#> 3 1 17 3 15 TRUE 21.36882 17.99862 3.370200 0 1
#> 4 1 18 4 15 TRUE 27.69603 17.98163 9.714395 0 1
#> 5 1 19 5 15 TRUE 23.83178 18.71806 5.113716 0 1
#> 6 1 20 6 15 TRUE 26.23287 20.94672 5.286154 0 1Each row is one treated cell. The columns are the unit and time coordinates, the observed outcome, the imputed counterfactual, the cell-level effect (the cell’s contribution to the ATT estimator), and the aggregation weight. Setting replicates = TRUE expands the data so there is one row per cell per bootstrap replicate.
An example: the per-event-time standard deviation of the cell-level treatment effect, a simple heterogeneity diagnostic.
po |>
group_by(event.time) |>
summarise(sd_eff = sd(eff, na.rm = TRUE),
n = dplyr::n(),
.groups = "drop") |>
head(8)
#> # A tibble: 8 × 3
#> event.time sd_eff n
#> <int> <dbl> <int>
#> 1 1 2.24 57
#> 2 2 2.29 57
#> 3 3 2.27 57
#> 4 4 2.90 57
#> 5 5 2.47 57
#> 6 6 3.09 57
#> 7 7 3.02 57
#> 8 8 3.09 57To get bootstrap confidence intervals for a custom summary, group the replicate-expanded data by event time and replicate, compute the summary inside each group, then take the appropriate quantiles across replicates. ggplot2 is the most flexible option for plotting; the data-frame format keeps the plot code as short as the aggregation code.
3.11 Save bootstrap runs
Setting keep.sims = TRUE at fit time saves the cell-level imputed counterfactuals from every bootstrap replicate on the fit object — a three-dimensional array of dimension \(T \times N \times \text{nboots}\). Some estimands need it, others do not.
The level ATT and the cumulative ATT (over either the event-time axis or an overall window) work without it, because they read the pre-aggregated bootstrap that fect always keeps. The other types — APTT, log-scale ATT, and any call that filters cells or groups by something other than event time — recompute their aggregations over the cell-level array and therefore need the saved runs.
A back-of-envelope memory cost:
| Panel size (T × N) | nboots | Memory cost |
|---|---|---|
| 1,000 | 200 | 1.6 MB |
| 10,000 | 200 | 16 MB |
| 50,000 | 200 | 80 MB |
| 50,000 | 500 | 200 MB |
| 200,000 | 500 | 800 MB |
If you fit without keep.sims = TRUE and then call an estimand that needs it, the function errors with a message pointing back to the option.
# No bootstrap/jackknife results available. Choose keep.sims = TRUE in fect().3.12 Parametric variance
Parametric bootstrap (vartype = "parametric") is one of fect’s fit-time variance options alongside "bootstrap" (the default) and "jackknife". The parametric option resamples residuals from a fitted factor model rather than resampling units; it requires time.component.from = "nevertreated" (see Section 7.2.4).
estimand() treats parametric replicates the same as bootstrap or jackknife replicates: they are stored in fit$eff.boot and all four type values use them through the same code. In v2.4.1, the vartype argument of estimand() accepts "parametric" explicitly so callers can name the method they expect.
For details on how vartype, ci.method, and the parametric para.error sub-option combine to form the bootstrap distribution, see Chapter 7.
fit_para <- fect(Y ~ D, data = sim_linear, index = c("id", "time"),
method = "ife", force = "two-way",
r = 2, CV = FALSE, se = TRUE, nboots = 200,
keep.sims = TRUE,
vartype = "parametric",
time.component.from = "nevertreated",
parallel = FALSE)
est <- estimand(fit_para, "att", "event.time")
head(est)
#> event.time estimate se ci.lo ci.hi n_cells vartype
#> 1 -39 -0.019958112 0.2427538 -0.4957469 0.4558307 80 parametric
#> 2 -38 -0.006695018 0.1693391 -0.3385935 0.3252035 80 parametric
#> ...The vartype column reports "parametric" — the method actually used at fit time — regardless of what you pass as the vartype argument. The argument is informational; it does not re-aggregate replicates. The output matches fit_para$est.att exactly for the per-event-time path.
For the other types ("att.cumu", "aptt", "log.att", and "att" with by = "overall"), estimand() aggregates quantiles over fit$eff.boot. Under parametric, the simulated distribution is centered on the parametric DGP’s mean rather than on the point estimate, so quantile-based confidence intervals reflect the parametric sampling distribution rather than a reflection around the point estimate. This matches what fit$est.att reports under parametric.
3.13 Migration
The legacy functions effect() and att.cumu() still work and produce identical results, but now emit a one-time-per-session deprecation message pointing here. Direct replacements:
| Today | Replacement |
|---|---|
effect(fit, cumu = TRUE) |
estimand(fit, "att.cumu", "event.time") |
effect(fit, cumu = FALSE) |
estimand(fit, "att", "event.time") |
att.cumu(fit, period = c(L, R)) |
estimand(fit, "att.cumu", "overall", window = c(L, R)) |
Outputs from estimand() and the legacy functions match to all decimals (verified by package tests). The only user-facing difference is that estimand() returns a data frame instead of the legacy matrix.