library(dplyr)
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. Only two calls work without them, because they read the pre-aggregated bootstrap that fect always keeps: the per-event-time level ATT and the cumulative ATT over a window (by = "overall"). See “Save bootstrap runs” below.
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. Before fect 2.4.7 the cumulative ATT was count-weighted: at event time \(k\) it was \(k\) times the average effect over all treated cells in event times 1 to \(k\). att.cumu(fit, weighted = TRUE) still returns that number, which differs from the running sum when the number of treated units changes over event time.
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.02260264 0.2712569 0.3607629 57 bootstrap
#> 2 2 0.3495009 0.02558756 0.2945536 0.3954525 57 bootstrap
#> 3 3 0.3476220 0.02519393 0.2966826 0.3964558 57 bootstrap
#> 4 4 0.3863353 0.03517108 0.3254791 0.4622917 57 bootstrap
#> 5 5 0.3269021 0.03063056 0.2662682 0.3866661 57 bootstrap
#> 6 6 0.3598550 0.03714676 0.2901648 0.4381495 57 bootstrap
#> 7 7 0.3669511 0.02606676 0.3177827 0.4148571 57 bootstrap
#> 8 8 0.3384248 0.02963756 0.2839744 0.3978245 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.01708885 0.2320849 0.3013446 57 bootstrap
#> 2 2 0.2921348 0.01844194 0.2490247 0.3228842 57 bootstrap
#> 3 3 0.2943029 0.01796569 0.2585549 0.3280302 57 bootstrap
#> 4 4 0.3163869 0.02405264 0.2732990 0.3688683 57 bootstrap
#> 5 5 0.2784879 0.02278222 0.2298704 0.3206099 57 bootstrap
#> 6 6 0.2995789 0.02642032 0.2492163 0.3540273 57 bootstrap
#> 7 7 0.3066786 0.01957622 0.2670976 0.3420557 57 bootstrap
#> 8 8 0.2850468 0.02244078 0.2440915 0.3292031 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.3007943 4.932796 6.111888 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".
window and cells work only with by = "overall". For "att", they choose the treated cells to average. For "att.cumu", use window: it sets the range of event times whose per-period ATTs are summed. A cells filter stops there, because the cumulative ATT is a sum over event times, not an average over chosen cells (before v2.4.7 it was ignored without a message). "aptt" and "log.att" take neither. The table lists what each type takes in this release; by = "cohort" and by = "calendar.time" are not available yet, and they stop with a message that names the values that work.
type |
by = "event.time" |
by = "overall" |
cells or window
|
|---|---|---|---|
"att" |
yes | yes | with by = "overall"
|
"att.cumu" |
yes | yes |
window only, with by = "overall"
|
"aptt", "log.att"
|
yes | no | no |
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 returns the cells of every bootstrap replicate, with a replicate column. Each replicate is a refit on a resampled panel, so its rows are the treated cells of the units it drew: a unit drawn twice appears twice, and a unit not drawn does not appear. The number of rows therefore varies by replicate (under the parametric bootstrap, every treated cell appears once per replicate). Averaging eff within a replicate gives that replicate’s overall ATT.
po_rep <- imputed_outcomes(fit, replicates = TRUE)
## rows per replicate: that replicate's own treated cells
summary(as.vector(table(po_rep$replicate)))
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 676.0 723.0 738.0 738.3 753.0 798.0
## the mean of eff within each replicate is its overall ATT draw
all.equal(c(tapply(po_rep$eff, po_rep$replicate, mean)),
c(fit$att.avg.boot), check.attributes = FALSE)
#> [1] TRUEAn 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 treatment effects from every bootstrap replicate on the fit object — fit$eff.boot, a three-dimensional array of dimension \(T \times N \times \text{nboots}\). Slice \(b\) holds the units that replicate \(b\) drew, in the order given by fit$colnames.boot[[b]] (their column positions in the original panel), not the units \(1, \dots, N\) in order. With a cluster bootstrap (cl), replicates differ in size, and the array is padded with NA to the widest one. estimand(), imputed_outcomes() and effect() read each replicate through colnames.boot. Some estimands need these saved runs, others do not.
Two calls work without them, because they read the pre-aggregated bootstrap that fect always keeps: the per-event-time level ATT, estimand(fit, "att", "event.time"), and the cumulative ATT over a window, estimand(fit, "att.cumu", "overall"). Everything else — the overall level ATT, the cumulative ATT by event time, APTT, log-scale ATT, and any call that filters cells — recomputes its aggregation over the cell-level array and therefore needs 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 -1.3251852 0.1476319 -1.614538 -1.0358320 80 parametric
#> 2 -38 -0.8431369 0.1299598 -1.097853 -0.5884204 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 "aptt", "log.att", and "att" with by = "overall", estimand() recomputes the statistic inside each replicate from fit$eff.boot, reading replicate \(b\) through its own units (fit$colnames.boot[[b]]; under the parametric bootstrap the treated units come first). The parametric replicates are simulated under no treatment effect, so they are centered near zero; estimand() shifts them to center on the point estimate before building the interval, which leaves their spread, and so the SE, unchanged. "att.cumu" takes its interval from att.cumu() (by = "overall") or effect() (by = "event.time"), which for parametric fits use the estimate plus or minus the normal critical value times the SE.
3.13 Migration
The legacy functions effect() and att.cumu() still work, 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)) |
For the cumulative ATT the replacements give the same numbers to all decimals: estimand(fit, "att.cumu", "event.time") returns what effect(fit, cumu = TRUE) returns, and the overall form returns the last row of att.cumu(). For effect(fit, cumu = FALSE), the point estimates and standard errors equal those of estimand(fit, "att", "event.time"). For bootstrap fits the intervals differ: effect() takes the 2.5% and 97.5% quantiles of the bootstrap replicates, while estimand() returns the normal intervals of fit$est.att. In the example above they differ by up to 0.57; at event time 1 they are [4.01, 5.26] and [3.99, 5.23]. For parametric and jackknife fits both give the normal intervals (with the fit’s default ci.method = "normal"), so they agree. The other difference is that estimand() returns a data frame instead of the legacy matrix.