Introduction to etwfe

Background

A canonical research design for social scientists is the so-called “difference-in-differences” (DiD) design. The essential idea behind DiD is that we can estimate the causal impact of an intervention (“treatment”) by comparing the pre- and post-intervention outcomes of indivduals that received treatment (“treated” group) against the outcomes of comparable indivduals that didn’t receive treatment (“control” group).¹

In the classic 2x2 DiD case (two units, two periods), a simple interaction effect between two dummy variables suffices to recover the treatment effect. In base R this might look something like:

lm(y ~ Dtreated_unit * Dpost_treatment, data = somedata)

where the resulting coefficient on the Dtreated_unitTRUE:Dpost_treatmentTRUE interaction term represents the treatment effect.

Rather than manually specify the interaction term, in practice researchers often use an equivalent formulation known as two-way fixed effects (TWFE). The core idea of TWFE is that we can subsume the interaction term from the previous code chunk by adding unit and time fixed effects. A single treatment dummy can then be used to capture the effect of treatment directly. A TWFE regression in base R might look as follows:

lm(y ~ Dtreat + factor(id) + factor(period), data = somedata)

where the treatment effect is now captured by the coefficient on the Dtreat dummy.

The TWFE shortcut is especially nice for more complicated panel data settings with multiple units and multiple times periods. Speaking of which, if you prefer to use a dedicated fixed effects / panel data package like fixest, you could also estimate the previous regression like so:

library(fixest)
feols(y ~ Dtreat | id + period, data = somedata)

These TWFE regressions are easy to run and intuitive, and for a long time everyone was happy. But it was too good to last. A cottage industry of clever research now demonstrates that things are not quite so simple. Among other things, the standard TWFE formulation can impose strange (negative) weighting conditions on key parts of the estimation procedure. One implication is that you risk a high probability of estimate bias in the presence of staggered treatment rollouts, which are very common in real-life applications.

Fortunately, just as econometricians were taking away one of our favourite tools, they were kind enough to replace it with some new ones. Among these, the proposed approach by Wooldridge (2021, 2022) is noteworthy. His idea might be paraphrased as stating that the problem with TWFE is not that we were doing it in the first place. Rather, it’s that we weren’t doing it enough. Instead of only including a single treatment × time interaction, Wooldridge recommends that we saturate our model with all possible interactions between treatment and time variables, including treatment cohorts, as well as other covariates. He goes on to show that this approach actually draws an equivalence between different types of estimators (pooled OLS, twoway Mundlak regression, etc.) So it’s not entirely clear what to call it. But Wooldridge refers to the general idea as as extended TWFE—or, ETWFE—which I rather like and is where this package takes its name.

The Wooldridge ETWFE solution is intuitive and elegant. But it is also rather tedious and error prone to code up manually. You have to correctly specify all the possible interactions, demean control variables within groups, and then recover the treatment effects of interest via an appropriate marginal effect aggregation. The etwfe package aims to simplify the process by providing convenience functions that do all this work for you.

Dataset

To demonstrate the core functionality of etwfe, we’ll use the mpdta dataset on US teen employment from the did package (which you’ll need to install separately).

# install.packages("did")
data("mpdta", package = "did")
head(mpdta)
#>     year countyreal     lpop     lemp first.treat treat
#> 866 2003       8001 5.896761 8.461469        2007     1
#> 841 2004       8001 5.896761 8.336870        2007     1
#> 842 2005       8001 5.896761 8.340217        2007     1
#> 819 2006       8001 5.896761 8.378161        2007     1
#> 827 2007       8001 5.896761 8.487352        2007     1
#> 937 2003       8019 2.232377 4.997212        2007     1

“Treatment” in this dataset refers to an increase in the minimum wage rate. In the examples that follow, our goal is to estimate the effect of this minimum wage treatment (treat) on the log of teen employment (lemp). Notice that the panel ID is at the county level (countyreal), but treatment was staggered across cohorts (first.treat) so that a group of counties were treated at the same time. In addition to these staggered treatment effects, we also observe log population (lpop) as a potential control variable.

Basic usage

Let’s load etwfe and work through its basic functionality. As we’ll see, the core workflow of the package involves two consecutive function calls: 1) etwfe() and 2) emfx().

library(etwfe)

mod =
  etwfe(
    fml  = lemp ~ lpop, # outcome ~ controls
    tvar = year,        # time variable
    gvar = first.treat, # group variable
    data = mpdta,       # dataset
    vcov = ~countyreal  # vcov adjustment (here: clustered)
    )

mod
#> OLS estimation, Dep. Var.: lemp
#> Observations: 2,500
#> Fixed-effects: first.treat: 4,  year: 5
#> Varying slopes: lpop (first.treat): 4,  lpop (year): 5
#> Standard-errors: Clustered (countyreal) 
#>                                               Estimate Std. Error   t value
#> .Dtreat:first.treat::2004:year::2004         -0.021248   0.021728 -0.977890
#> .Dtreat:first.treat::2004:year::2005         -0.081850   0.027375 -2.989963
#> .Dtreat:first.treat::2004:year::2006         -0.137870   0.030795 -4.477097
#> .Dtreat:first.treat::2004:year::2007         -0.109539   0.032322 -3.389024
#> .Dtreat:first.treat::2006:year::2006          0.002537   0.018883  0.134344
#> .Dtreat:first.treat::2006:year::2007         -0.045093   0.021987 -2.050907
#> .Dtreat:first.treat::2007:year::2007         -0.045955   0.017975 -2.556568
#> .Dtreat:first.treat::2004:year::2004:lpop_dm  0.004628   0.017584  0.263184
#> .Dtreat:first.treat::2004:year::2005:lpop_dm  0.025113   0.017904  1.402661
#> .Dtreat:first.treat::2004:year::2006:lpop_dm  0.050735   0.021070  2.407884
#> .Dtreat:first.treat::2004:year::2007:lpop_dm  0.011250   0.026617  0.422648
#> .Dtreat:first.treat::2006:year::2006:lpop_dm  0.038935   0.016472  2.363731
#> .Dtreat:first.treat::2006:year::2007:lpop_dm  0.038060   0.022477  1.693276
#> .Dtreat:first.treat::2007:year::2007:lpop_dm -0.019835   0.016198 -1.224528
#>                                                Pr(>|t|)    
#> .Dtreat:first.treat::2004:year::2004         3.2860e-01    
#> .Dtreat:first.treat::2004:year::2005         2.9279e-03 ** 
#> .Dtreat:first.treat::2004:year::2006         9.3851e-06 ***
#> .Dtreat:first.treat::2004:year::2007         7.5694e-04 ***
#> .Dtreat:first.treat::2006:year::2006         8.9318e-01    
#> .Dtreat:first.treat::2006:year::2007         4.0798e-02 *  
#> .Dtreat:first.treat::2007:year::2007         1.0866e-02 *  
#> .Dtreat:first.treat::2004:year::2004:lpop_dm 7.9252e-01    
#> .Dtreat:first.treat::2004:year::2005:lpop_dm 1.6134e-01    
#> .Dtreat:first.treat::2004:year::2006:lpop_dm 1.6407e-02 *  
#> .Dtreat:first.treat::2004:year::2007:lpop_dm 6.7273e-01    
#> .Dtreat:first.treat::2006:year::2006:lpop_dm 1.8474e-02 *  
#> .Dtreat:first.treat::2006:year::2007:lpop_dm 9.1027e-02 .  
#> .Dtreat:first.treat::2007:year::2007:lpop_dm 2.2133e-01    
#> ... 10 variables were removed because of collinearity (.Dtreat:first.treat::2006:year::2004, .Dtreat:first.treat::2006:year::2005 and 8 others [full set in $collin.var])
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> RMSE: 0.537131     Adj. R2: 0.87167 
#>                  Within R2: 8.449e-4

emfx(mod)
#> 
#>     Term                 Contrast .Dtreat Estimate Std. Error     z Pr(>|z|)
#>  .Dtreat mean(TRUE) - mean(FALSE)    TRUE  -0.0506     0.0125 -4.05   <0.001
#>     S   2.5 %  97.5 %
#>  14.3 -0.0751 -0.0261
#> 
#> Columns: term, contrast, .Dtreat, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted 
#> Type:  response

mod_es = emfx(mod, type = "event")
mod_es
#> 
#>     Term                 Contrast event Estimate Std. Error     z Pr(>|z|)    S
#>  .Dtreat mean(TRUE) - mean(FALSE)     0  -0.0332     0.0134 -2.48    0.013  6.3
#>  .Dtreat mean(TRUE) - mean(FALSE)     1  -0.0573     0.0172 -3.34   <0.001 10.2
#>  .Dtreat mean(TRUE) - mean(FALSE)     2  -0.1379     0.0308 -4.48   <0.001 17.0
#>  .Dtreat mean(TRUE) - mean(FALSE)     3  -0.1095     0.0323 -3.39   <0.001 10.5
#>    2.5 %   97.5 %
#>  -0.0594 -0.00701
#>  -0.0910 -0.02373
#>  -0.1982 -0.07751
#>  -0.1729 -0.04619
#> 
#> Columns: term, contrast, event, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted 
#> Type:  response

library(modelsummary)
#> Version 2.0.0 of `modelsummary`, to be released soon, will introduce a
#>   breaking change: The default table-drawing package will be `tinytable`
#>   instead of `kableExtra`. All currently supported table-drawing packages
#>   will continue to be supported for the foreseeable future, including
#>   `kableExtra`, `gt`, `huxtable`, `flextable, and `DT`.
#>   
#>   You can always call the `config_modelsummary()` function to change the
#>   default table-drawing package in persistent fashion. To try `tinytable`
#>   now:
#>   
#>   config_modelsummary(factory_default = 'tinytable')
#>   
#>   To set the default back to `kableExtra`:
#>   
#>   config_modelsummary(factory_default = 'kableExtra')

# Quick renaming function to replace ".Dtreat" with something more meaningful
rename_fn = function(old_names) {
  new_names = gsub(".Dtreat", "Years post treatment =", old_names)
  setNames(new_names, old_names)
}

modelsummary(
  mod_es,
  shape       = term:event:statistic ~ model,
  coef_rename = rename_fn,
  gof_omit    = "Adj|Within|IC|RMSE",
  title       = "Event study",
  notes       = "Std. errors are clustered at the county level"
)

library(ggplot2)
theme_set(
  theme_minimal() + theme(panel.grid.minor = element_blank())
)

ggplot(mod_es, aes(x = event, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_pointrange(col = "darkcyan") +
  labs(x = "Years post treatment", y = "Effect on log teen employment")

# Use post_only = FALSE to get the "zero" pre-treatment effects
mod_es2 = emfx(mod, type = "event", post_only = FALSE)

ggplot(mod_es2, aes(x = event, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = -1, lty = 2) +
  geom_pointrange(col = "darkcyan") +
  labs(
    x = "Years post treatment", y = "Effect on log teen employment",
    caption = "Note: Zero pre-treatment effects for illustrative purposes only."
  )
#> Warning: Removed 4 rows containing missing values (`geom_segment()`).

Heterogeneous treatment effects

So far we’ve limited ourselves to homogeneous treatment effects, where the impact of treatment (i.e., minimum wage hike) is averaged across all US counties in our dataset. However, many research problems require us to estimate treatment effects separately across groups and then, potentially, test for differences between them. For example, we might want to test whether the efficacy of a new vaccine differs across age groups, or whether a marketing campaign was equally successful across genders. The ETWFE framework naturally lends itself to these kinds of heterogeneous treatment effects.

Consider the following example, where we first create a logical dummy variable for all US counties in the eight Great Lake states (GLS).

gls_fips = c("IL" = 17, "IN" = 18, "MI" = 26, "MN" = 27,
             "NY" = 36, "OH" = 39, "PA" = 42, "WI" = 55)

mpdta$gls = substr(mpdta$countyreal, 1, 2) %in% gls_fips

Now imagine that we are interested in estimating separate treatment effects for GLS versus non-GLS counties. We do this simply by invoking the optional xvar argument as part of our etwfe() call.³ Any subsequent emfx() call on this object will automatically recognize that we want to recover treatment effects by these two distinct groups.

hmod = etwfe(
   lemp ~ lpop, tvar = year, gvar = first.treat, data = mpdta, 
   vcov = ~countyreal,
   xvar = gls           ## <= het. TEs by gls
   )

# Heterogeneous ATTs (could also specify `type = "event"`, etc.) 
emfx(hmod)
#> 
#>     Term                 Contrast .Dtreat   gls Estimate Std. Error     z
#>  .Dtreat mean(TRUE) - mean(FALSE)    TRUE FALSE  -0.0637     0.0376 -1.69
#>  .Dtreat mean(TRUE) - mean(FALSE)    TRUE  TRUE  -0.0472     0.0271 -1.74
#>  Pr(>|z|)   S  2.5 %  97.5 %
#>    0.0906 3.5 -0.137 0.01007
#>    0.0817 3.6 -0.100 0.00594
#> 
#> Columns: term, contrast, .Dtreat, gls, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted 
#> Type:  response

The above point estimates might tempt us to think that minimum wage hikes caused less teen disemployment in GLS counties than in the rest of the US on average. However, to test this formally we can invoke the powerful hypothesis infrastructure of the underlying marginaleffects package. Probably the easiest way to do this is by using b[i]-style positional arguments, where “[i]” denotes the row of the emfx() return object. Thus, by specifying hypothesis = "b1 = b2", we can test whether the ATTs from row 1 (non-GLS) and row 2 (GLS) are different from one another.

emfx(hmod, hypothesis = "b1 = b2")
#> 
#>   Term Estimate Std. Error      z Pr(>|z|)   S  2.5 % 97.5 %
#>  b1=b2  -0.0164     0.0559 -0.294    0.769 0.4 -0.126  0.093
#> 
#> Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
#> Type:  response

Here we see that there is actually no statistical difference in the average disemployment effect between GLS and non-GLS counties.

One final aside is that you can easily display the results of heterogeneous treatment effects in plot or table form. Here’s an example of the latter, where we again make use of the modelsummary(..., shape = ...) argument.

modelsummary(
    models      = list("GLS county" = emfx(hmod)),
    shape       = term + statistic ~ model + gls, # add xvar variable (here: gls)
    coef_map    = c(".Dtreat" = "ATT"),
    gof_map     = NA,
    title       = "Comparing the ATT on GLS and non-GLS counties"
)

While the simple example here has been limited to a binary comparison of group
ATTs, note the same logic carries over to richer settings. We can use the exact same workflow to estimate heterogeneous treatment effects by different aggregations (e.g., event studies) and across groups with many levels.

Other families

Another key feature of the ETWFE approach—one that sets it apart from other advanced DiD implementations and extensions—is that it supports nonlinear model (distribution / link) families. Users need simply invoke the family argument. Here’s a brief example, where we recast our earlier event-study as a Poisson regression.

mpdta$emp = exp(mpdta$lemp)

etwfe(
  emp ~ lpop, tvar = year, gvar = first.treat, data = mpdta, vcov = ~countyreal,
  family = "poisson"
  ) |>
  emfx("event")
#> 
#>     Term                 Contrast event Estimate Std. Error       z Pr(>|z|)
#>  .Dtreat mean(TRUE) - mean(FALSE)     0   -25.35       15.9 -1.5942    0.111
#>  .Dtreat mean(TRUE) - mean(FALSE)     1     1.09       41.8  0.0261    0.979
#>  .Dtreat mean(TRUE) - mean(FALSE)     2   -75.12       22.3 -3.3696   <0.001
#>  .Dtreat mean(TRUE) - mean(FALSE)     3  -101.82       28.1 -3.6234   <0.001
#>     S  2.5 % 97.5 %
#>   3.2  -56.5   5.82
#>   0.0  -80.9  83.09
#>  10.4 -118.8 -31.43
#>  11.7 -156.9 -46.75
#> 
#> Columns: term, contrast, event, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted 
#> Type:  response

Performance tips

Thinking of the etwfe workflow again as a pair of consecutive functional calls, the first etwfe() stage tends to be very fast. We’re leveraging the incredible performance of fixest and also taking some shortcuts to avoid wasting time on nuisance parameters. See the Regarding fixed effects section below for more details about this.

For its part, the second emfx() stage also tends to be pretty performant. If your data has less than 100k rows, it’s unlikely that you’ll have to wait more than a few seconds to obtain results. However, emfx’s computation time does tend to scale non-linearly with the size of the original data, as well as the number of interactions from the underlying etwfe model object. Without getting too deep into the weeds, we are relying on a numerical delta method of the (excellent) marginaleffects package underneath the hood to recover the ATTs of interest. This method requires estimating two prediction models for each coefficient in the model and then computing their standard errors. So it’s a potentially expensive operation that can push the computation time for large datasets (> 1m rows) up to several minutes or longer.

Fortunately, there are two complementary strategies that you can use to speed things up. The first is to turn off the most expensive part of the whole procedure—standard error calculation—by calling emfx(..., vcov = FALSE). Doing so should bring the estimation time back down to a few seconds or less, even for datasets in excess of a million rows. Of course, the loss of standard errors might not be an acceptable trade-off for projects where statistical inference is critical. But the good news is this first strategy can still be combined our second strategy: it turns out that collapsing the data by groups prior to estimating the marginal effects can yield substantial speed gains on its own. Users can do this by invoking the emfx(..., collapse = TRUE) argument. While the effect here is not as dramatic as the first strategy, collapsing the data does have the virtue of retaining information about the standard errors. The trade-off this time, however, is that collapsing our data does lead to a loss in accuracy for our estimated parameters. On the other hand, testing suggests that this loss in accuracy tends to be relatively minor, with results equivalent up to the 1st or 2nd significant decimal place (or even better).

Summarizing, here is a quick plan of attack for you to try if you are worried about the estimation time for large datasets and models:

0. Estimate mod = etwfe(...) as per usual.
1. Run emfx(mod, vcov = FALSE, ...).
2. Run emfx(mod, vcov = FALSE, collapse = TRUE, ...).
3. Compare the point estimates from steps 1 and 2. If they are are similar enough to your satisfaction, get the approximate standard errors by running emfx(mod, collapse = TRUE, ...).

It’s a bit of performance art, since all of the examples in this vignette complete very quickly anyway. But here is a reworking of our earlier event study example to demonstrate this performance-conscious workflow.

# Step 0 already complete: using the same `mod` object from earlier...

# Step 1
emfx(mod, type = "event", vcov = FALSE)
#> 
#>     Term                 Contrast event Estimate
#>  .Dtreat mean(TRUE) - mean(FALSE)     0  -0.0332
#>  .Dtreat mean(TRUE) - mean(FALSE)     1  -0.0573
#>  .Dtreat mean(TRUE) - mean(FALSE)     2  -0.1379
#>  .Dtreat mean(TRUE) - mean(FALSE)     3  -0.1095
#> 
#> Columns: term, contrast, event, estimate, predicted_lo, predicted_hi, predicted 
#> Type:  response

# Step 2
emfx(mod, type = "event", vcov = FALSE, collapse = TRUE)
#> 
#>     Term                 Contrast event Estimate
#>  .Dtreat mean(TRUE) - mean(FALSE)     0  -0.0332
#>  .Dtreat mean(TRUE) - mean(FALSE)     1  -0.0573
#>  .Dtreat mean(TRUE) - mean(FALSE)     2  -0.1379
#>  .Dtreat mean(TRUE) - mean(FALSE)     3  -0.1095
#> 
#> Columns: term, contrast, event, estimate, predicted_lo, predicted_hi, predicted 
#> Type:  response

# Step 3: Results from 1 and 2 are similar enough, so get approx. SEs
mod_es2 = emfx(mod, type = "event", collapse = TRUE)

To put a fine point on it, we can can compare our original event study with the collapsed estimates and see that the results are indeed very similar.

modelsummary(
    list("Original" = mod_es, "Collapsed" = mod_es2),
    shape       = term:event:statistic ~ model,
    coef_rename = rename_fn,
    gof_omit    = "Adj|Within|IC|RMSE",
    title       = "Event study",
    notes       = "Std. errors are clustered at the county level"
)

Under the hood

Now that you’ve seen etwfe in action, let’s circle back to what the package is doing under the hood. This section isn’t necessary for you to use the package; feel free to skip it. But a review of the internal details should help you to optimize for different scenarios and also give you a better understanding of etwfe’s default choices.

mod$fml_all
#> $linear
#> lemp ~ .Dtreat:i(first.treat, i.year, ref = 0, ref2 = 2003)/lpop_dm
#> <environment: 0x7fde98dc3478>
#> 
#> $fixef
#> ~first.treat + first.treat[[lpop]] + year + year[[lpop]]

# First construct the dataset
mpdta2 = mpdta |>
  transform(
    .Dtreat = year >= first.treat & first.treat != 0,
    lpop_dm = ave(lpop, first.treat, year, FUN = \(x) x - mean(x, na.rm = TRUE))
  )

# Then estimate the manual version of etwfe
mod2 = fixest::feols(
  lemp ~ .Dtreat:i(first.treat, i.year, ref = 0, ref2 = 2003) / lpop_dm |
    first.treat[lpop] + year[lpop],
  data = mpdta2,
  vcov = ~countyreal
)

modelsummary(
  list("etwfe" = mod, "manual" = mod2),
  gof_map = NA # drop all goodness-of-fit info for brevity
)

library(marginaleffects)

# Simple ATT
slopes(
  mod2, 
  newdata   = subset(mpdta2, .Dtreat), # we only want rows where .Dtreat is TRUE
  variables = ".Dtreat", 
  by        = ".Dtreat"
  )
#> 
#>     Term                 Contrast .Dtreat Estimate Std. Error     z Pr(>|z|)
#>  .Dtreat mean(TRUE) - mean(FALSE)    TRUE  -0.0506     0.0125 -4.05   <0.001
#>     S   2.5 %  97.5 %
#>  14.3 -0.0751 -0.0261
#> 
#> Columns: term, contrast, .Dtreat, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted 
#> Type:  response

# Event study
slopes(
  mod2, 
  newdata   = transform(subset(mpdta2, .Dtreat), event = year - first.treat),
  variables = ".Dtreat", 
  by        = "event"
  )
#> 
#>     Term                 Contrast event Estimate Std. Error     z Pr(>|z|)    S
#>  .Dtreat mean(TRUE) - mean(FALSE)     0  -0.0332     0.0134 -2.48    0.013  6.3
#>  .Dtreat mean(TRUE) - mean(FALSE)     1  -0.0573     0.0172 -3.34   <0.001 10.2
#>  .Dtreat mean(TRUE) - mean(FALSE)     2  -0.1379     0.0308 -4.48   <0.001 17.0
#>  .Dtreat mean(TRUE) - mean(FALSE)     3  -0.1095     0.0323 -3.39   <0.001 10.5
#>    2.5 %   97.5 %
#>  -0.0594 -0.00701
#>  -0.0910 -0.02373
#>  -0.1982 -0.07751
#>  -0.1729 -0.04619
#> 
#> Columns: term, contrast, event, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted 
#> Type:  response

# fe = "feo" (fixed effects only)
mod_feo = etwfe(
  lemp ~ lpop, tvar = year, gvar = first.treat, data = mpdta, vcov = ~countyreal,
  fe = "feo"
)
# ... which is equivalent to the manual regression
mod_feo2 = fixest::feols(
  lemp ~ .Dtreat:i(first.treat, i.year, ref = 0, ref2 = 2003) / lpop_dm +
    lpop + i(first.treat, lpop, ref = 0) + i(year, lpop, ref = 2003) |
    first.treat + year,
  data = mpdta2, vcov = ~countyreal
)

# fe = "none"
mod_none = etwfe(
  lemp ~ lpop, tvar = year, gvar = first.treat, data = mpdta, vcov = ~countyreal,
  fe = "none"
)
# ... which is equivalent to the manual regression
mod_none2 = fixest::feols(
  lemp ~ .Dtreat:i(first.treat, i.year, ref = 0, ref2 = 2003) / lpop_dm +
    lpop + i(first.treat, lpop, ref = 0) + i(year, lpop, ref = 2003) +
    i(first.treat, ref = 0) + i(year, ref = 2003),
  data = mpdta2, vcov = ~countyreal
)

mods = list(
  "etwfe"         = mod,
  "manual"        = mod2,
  "etwfe (feo)"   = mod_feo,
  "manual (feo)"  = mod_feo2,
  "etwfe (none)"  = mod_none,
  "manual (none)" = mod_none2
)

modelsummary(mods, gof_map = NA)

mod_es_i = etwfe(
  lemp ~ lpop, tvar = year, gvar = first.treat, data = mpdta,
  ivar = countyreal  # NEW: Use unit-level (county) FEs
  ) |>
  emfx("event")

modelsummary(
  list("Group-level FEs (default)" = mod_es, "Unit-level FEs" = mod_es_i),
  shape       = term:event:statistic ~ model,
  coef_rename = rename_fn,
  gof_omit    = "Adj|Within|IC|RMSE"
)