2  Fect Main Program

In this chapter, we illustrate how to use the fect package to implement counterfactual estimators (or imputation estimators) and conduct diagnostic tests proposed by Liu, Wang, and Xu (2024)–[Paper]. Download the R code used in this chapter here.

2.1 Simulated data

In this demo, we will primarily be using simdata, in which the treatment is allowed to switch on and off. There are 200 units and 35 time periods.

library(fect)
data(fect)
ls()
#> [1] "gs2020"    "hh2019"    "simdata"   "simgsynth" "turnout"

Before conducting any statistical analysis, we use the panelView package to visualize the treatment and outcome variables in simdata:

library(panelView)
panelview(Y ~ D, data = simdata, index = c("id","time"), 
  axis.lab = "time", xlab = "Time", ylab = "Unit", 
  gridOff = TRUE, by.timing = TRUE,
  background = "white", main = "Simulated Data: Treatment Status")

We then take a look at the outcome variable. In the figure below, blue and gray represent treatment and control conditions.

panelview(Y ~ D, data = simdata, index = c("id","time"), 
  axis.lab = "time", xlab = "Time", ylab = "Unit", 
  theme.bw = TRUE, type = "outcome", 
  main = "Simulated Data: Outcome")
#> Treatment has reversals.

2.2 The imputation estimator (FEct)

Using the current version of fect, we can apply several different methods to make counterfactual predictions and estimate treatment effects by specifying the method option: "fe" (two-way fixed effects, default), "ife" (interactive fixed effects), "mc" (matrix completion method), "cfe" (complex fixed effects), and "polynomial" (fixed effects with group-specific time trends). First, we illustrate the main syntax of fect using the "fe" method.

The two-way fixed effects counterfactual estimator (FEct) is also independently proposed by Borusyak, Jaravel, and Spiess (2024) and Gardner (2021), who refer to it as the “imputation method” and “two-stage DID,” respectively.

2.2.1 Estimation

We estimate the average treatment effect on the treated (ATT) using the following information: the outcome variable \(Y\), binary treatment variable \(D\), two observed covariates \(X_{1}\) and \(X_{2}\), and the unit and time indicators \(id\) and \(time\), respectively.

The first variable on the right hand side of the formula is the treatment indicator \(D\); the rest of the right-hand-side variables serve as controls. The index option specifies the unit and time indicators. The force option (“none”, “unit”, “time”, and “two-way”) specifies the additive component(s) of the fixed effects included in the model. The default option is “two-way” (including both unit and time fixed effects).

out.fect <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id","time"),
                 method = "fe", force = "two-way")

2.2.2 Visualization

We can use the plot function to visualize the estimation results. By default, the plot function produces a “gap” plot – as if we type plot(out.fect, type = "gap") — which visualizes the estimated period-wise ATT (dynamic treatment effects). For your reference, the true population average effects in simdata go from 1 to 3 from the 1st to the 10th post-treatment period.

The bar plot at the bottom of the plot shows the number of treated units for each time period. The options cex.main, cex.lab, cex.axis, and cex.text adjust the font sizes of the title, axis labels, axis numbers, and in-graph text, respectively.

Users can choose to plot only those periods whose number of treated observations exceeds a threshold, which is set as a proportion of the largest number of treated observations in a period (the default is proportion = 0.3).

plot(out.fect, main = "Estimated ATT (FEct)", ylab = "Effect of D on Y", 
  cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8)
#> Uncertainty estimates not available.

The uncertainty estimates are unavailable in the plot above because, by default, se = FALSE to save computational power.

The graph is a ggplot2 object; user can conveniently use the ggsave (preferred) function to export the resulting plot. See Chapter 3 for more visualization options.

2.2.3 Uncertainty estimates

The package can produce uncertainty estimates when se = TRUE. One can use the non-parametric bootstrap procedure by setting vartype = "bootstrap". Note that it only works well when the number of units is relatively large and many units in the data set experience the treatment condition. The number of bootstrap runs can be set by nboots.

Parallel computing will speed up both cross-validation and uncertainty estimation significantly. We recommend that users manually set the number of cores using the cores option. If this is not supplied or is NULL, we will automatically select the smaller of 8 and the number of usable system cores minus 2, to prevent excessive use of system resources.

out.fect <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id","time"), 
  method = "fe", force = "two-way", se = TRUE, 
  cores = 8, parallel = TRUE, nboots = 1000)

Clustering.

  • By default, fect() uses cluster-bootstrap at the unit level when se = TRUE, that is vartype = "bootstrap". The uncertainty estimates thus account for arbitrary serial correlation within a unit over time—commonly understood as being “clustered” at the unit level. In the example above, the unit is id, hence, by default, cl = "id".
  • Alternatively, users can obtain uncertainty estimates using the cluster-jackknife method by specifying vartype = "jackknife", also at the unit level. In this case, the algorithm calculates standard errors by iteratively dropping one unit (i.e., the entire time series) from the dataset.
  • To cluster standard errors at a different, usually higher, level, users can specify the clustering variable using the cl option.

The plot() function can visualize the estimated period-wise ATTs as well as their uncertainty estimates. stats = "F.p" shows the p-value for the F test of no-pretrend (more details below).

plot(out.fect, main = "Estimated ATT (FEct)", ylab = "Effect of D on Y", 
  cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, stats = "F.p")

2.2.4 Save estiamtes

Users can use the print function to take a look at a summary of the estimation results or retrieve relevant statistics by directly accessing the fect object. Specifically, est.avg and est.avg.unit show the ATT averaged over all periods – the former weights each treated observation equally while the latter weights each treated unit equally. est.beta reports the coefficients of the time-varying covariates. est.att reports the average treatment effect on the treated (ATT) by period. Treatment effect estimates from each bootstrap run is stored in eff.boot, an array whose dimension = (#time periods * #treated * #bootstrap runs).

print(out.fect)
#> Call:
#> fect.formula(formula = Y ~ D + X1 + X2, data = simdata, index = c("id", 
#>     "time"), force = "two-way", method = "fe", se = TRUE, nboots = 1000, 
#>     parallel = TRUE, cores = 8)
#> 
#> ATT:
#>                             ATT   S.E. CI.lower CI.upper p.value
#> Tr obs equally weighted   5.121 0.3070    4.519    5.722       0
#> Tr units equally weighted 3.622 0.2446    3.142    4.101       0
#> 
#> Covariates:
#>     Coef    S.E. CI.lower CI.upper p.value
#> X1 1.015 0.03421   0.9481    1.082       0
#> X2 2.937 0.03297   2.8723    3.002       0

To save space, results are not shown here.

out.fect$est.att
out.fect$est.avg
out.fect$est.beta

2.2.5 Leave-on-out approach

Li and Strezhnev (2025) show that in some applications, pre-trend estimates based on in-sample model fit can lead to the mistaken belief that no pre-trend exists, even when a non-parallel pre-trend is present. A simple fix is to use a leave-one-out method by setting loo = TRUE to obtain these estimates, although it is significantly more time-consuming.

We recommend setting loo = TRUE when (i) the event-study plot is intended as a critical piece of evidence to support the parallel trends assumption, which is often the case, or (ii) when implementing an equivalence test for the pre-trend estimates. For more discussion on the LOO pre-trend test, see here

Our most preferred test, however, is the sensitivity analysis discussed in Chapter 6, which combines out-of-sample placebo estimates with post-treatment ATT estimates.

We can implement the leave-one-out pre-trend test by setting loo = TRUE.

out.fect.loo <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id","time"), 
  method = "fe", force = "two-way", se = TRUE, loo = TRUE,
  cores = 8, parallel = TRUE, nboots = 200)

The event study plot utilizing leave-one-out for pretreatment estimates is shown below. This graph is fairly similar to the graphics we presented earlier without using leave-one-out. However, this is not always true.

plot(out.fect.loo,main = "Estimated ATT (FEct) -- LOO", 
  cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8)

2.3 Interactive fixed effects & matrix completion

In addition to FEct, fect also supports the interactive fixed effects counterfactual (IFEct) method proposed by Gobillon and Magnac (2016) and Xu (2017) and the matrix completion (MC) method proposed by Athey et al. (2021)—method = "ife" and method = "mc", respectively. We use EM algorithm to impute the counterfactuals of treated observations.

2.3.1 IFE

For the IFE approach, we need to specify the number of factors using option r. For the MC method, we need to specify the tuning parameter in the penalty term using option lambda. By default, the algorithm will select an optimal hyper-parameter via a built-in cross-validation procedure.

Choosing the number of factors. We provide a cross-validation procedure (by setting CV = TRUE) to help determine the tuning parameter in IFE, MC, and Gsynth (new addition, see next chapter for further details) methods. By default, the cross-validation procedure is run for k rounds (k = 10) and the candidate tuning parameter corresponding to the minimal mean squared prediction error is selected (criterion = "mspe").

In each round, some untreated observations are removed as the testing set to evaluate the prediction performance of the model with a tuning parameter. The option cv.prop specifies the size of testing set comparing to the set of observations under control (default: cv.prop = 0.1). If we want to restrict the testing set to untreated observations only from treated units (those whose treatment statuses have changed), set cv.treat = TRUE.

An additional issue is the serial correlation within a unit. We remove a consecutive number of observations from a unit as elements in the testing set in order to avoid over fitting caused by serial correlation. The consecutive number is specified in option cv.nobs (e.g. when cv.nobs = 3, the test set is a number of triplets).

We can also remove triplets in the fitting procedure but only include the middle observation of each triplet in the test set using the option cv.donut (e.g. when cv.donut = 1, the first and the last observation in each removed triplet will not be included in the test set).

Hyper-parameter tuning The package offers several criteria when tuning hyper-parameters. For the IFE method, we can set criterion = "pc" to select the hyper-parameter based on the information criterion. If we want to select the hyper-parameter based on mean-squared prediction errors from cross-validation to get a better prediction ability, set criterion = "mspe" (default), and to alleviate the impact of some outlier prediction errors, we allow the criterion of geometric-mean squared prediction errors (criterion = "gmspe"). If one wants to select the hyper-parameter that yields a better pre-trend fitting on test sets rather than a better prediction ability, set criterion = "moment" (we average the residuals in test sets by their relative periods to treatments and then average the squares of these period-wise deviations weighted by the number of observations at each period) .

For the IFE method, we need to specify an interval of candidate number of unobserved factors in option r like r=c(0,5). When cross-validation is switched off, the first element in r will be set as the number of factors. Below we use the MSPE criterion and search the number of factors from 0 to 5.

out.ife <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id","time"), 
          force = "two-way", method = "ife", CV = TRUE, r = c(0, 5), 
          se = TRUE, cores = 8, nboots = 1000, parallel = TRUE)
#> Parallel computing ...
#> Cross-validating ...
#> Criterion: Mean Squared Prediction Error
#> Interactive fixed effects model...
#> r = 0; sigma2 = 6.35460; IC = 2.21891; PC = 6.08178; MSPE = 6.78393
#> r = 1; sigma2 = 4.52698; IC = 2.24325; PC = 5.26760; MSPE = 4.93004
#> r = 2; sigma2 = 3.89603; IC = 2.45349; PC = 5.33953; MSPE = 4.50855
#> *
#> r = 3; sigma2 = 3.79056; IC = 2.78325; PC = 5.98062; MSPE = 5.13365
#> r = 4; sigma2 = 3.67967; IC = 3.10762; PC = 6.56967; MSPE = 5.72471
#> r = 5; sigma2 = 3.57625; IC = 3.43005; PC = 7.12886; MSPE = 6.49848
#> 
#>  r* = 2
print(out.ife)
#> Call:
#> fect.formula(formula = Y ~ D + X1 + X2, data = simdata, index = c("id", 
#>     "time"), force = "two-way", r = c(0, 5), CV = TRUE, method = "ife", 
#>     se = TRUE, nboots = 1000, parallel = TRUE, cores = 8)
#> 
#> ATT:
#>                             ATT   S.E. CI.lower CI.upper p.value
#> Tr obs equally weighted   3.073 0.2985    2.488    3.658       0
#> Tr units equally weighted 1.935 0.2048    1.533    2.336       0
#> 
#> Covariates:
#>     Coef    S.E. CI.lower CI.upper p.value
#> X1 0.998 0.02979   0.9396    1.056       0
#> X2 2.977 0.02743   2.9228    3.030       0

The figure below shows the estimated ATT using the IFE method. The cross-validation procedure selects the correct number of factors (\(r=2\)).

plot(out.ife, main = "Estimated ATT (IFEct)")

2.3.2 MC

For the MC method, we also need to specify a sequence of candidate tuning parameters. For example, we can specify lambda = c(1, 0.8, 0.6, 0.4, 0.2, 0.05). If users don’t have any prior knowledge to set candidate tuning parameters, a number of candidate tuning parameters can be generated automatically based on the information from the outcome variable. We specify the number in option nlambda, e.g. nlambda = 10.

out.mc <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id","time"), 
          force = "two-way", method = "mc", CV = TRUE, 
          se = TRUE, cores = 8, nboots = 1000, parallel = TRUE)
#> Parallel computing ...
#> Cross-validating ...
#> Criterion: Mean Squared Prediction Error
#> Matrix completion method...
#> lambda.norm = 1.00000; MSPE = 7.02065; MSPTATT = 0.28561; MSE = 5.80999
#> lambda.norm = 0.42170; MSPE = 5.45943; MSPTATT = 0.12218; MSE = 4.15356
#> lambda.norm = 0.17783; MSPE = 5.17432; MSPTATT = 0.04007; MSE = 1.57548
#> *
#> lambda.norm = 0.07499; MSPE = 5.37711; MSPTATT = 0.00815; MSE = 0.28590
#> lambda.norm = 0.03162; MSPE = 5.63317; MSPTATT = 0.00170; MSE = 0.05133
#> lambda.norm = 0.01334; MSPE = 6.55505; MSPTATT = 0.00034; MSE = 0.00914
#> 
#>  lambda.norm* = 0.177827941003892
#> 

print(out.mc)
#> Call:
#> fect.formula(formula = Y ~ D + X1 + X2, data = simdata, index = c("id", 
#>     "time"), force = "two-way", CV = TRUE, method = "mc", se = TRUE, 
#>     nboots = 1000, parallel = TRUE, cores = 8)
#> 
#> ATT:
#>                             ATT   S.E. CI.lower CI.upper p.value
#> Tr obs equally weighted   4.262 0.2978    3.678    4.846       0
#> Tr units equally weighted 2.826 0.2283    2.378    3.273       0
#> 
#> Covariates:
#>     Coef    S.E. CI.lower CI.upper p.value
#> X1 1.017 0.02825   0.9618    1.073       0
#> X2 2.959 0.02907   2.9022    3.016       0
plot(out.mc, main = "Estimated ATT (MC)")

2.4 Weighitng treatment effects

After obtaining the individual treatment effects using one of the counterfactual estimators, we can weight these estimates using a constructed balanced treated sample or other user-supplied weighting schemes.

2.4.1 Balanced treated sample

fect also provides the option balance.period, which allows the calculation of the average treatment effects only for treated units that exhibit complete data in specified pre- and post-treatment periods. For instance, if the option is set to balance.period = c(-3,4), the algorithm will calculate the average treatment effects for units that have at least four consecutive non-missing observations in the pre-treatment periods (-3, -2, -1, 0) and at least four consecutive non-missing observations in the post-treatment periods (1, 2, 3, 4). Note that this option does not affect whether a never-treated unit enters estimation.

out.bal <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id","time"), 
  balance.period = c(-3, 4), force = "two-way", method = "ife", 
  CV = FALSE, r = 2, se = TRUE, nboots = 200, parallel = TRUE)
#> Parallel computing ...

We can then visualize the dynamic treatment effects using the inbuilt function plot. By default, it displays the dynamic treatment effects of the “balanced” sample.

plot(out.bal, main = "Estimated ATT (Balanced Sample)")

The usual plotting options can be used to adjust the balanced plot as well.

plot(out.bal, main = "Estimated ATT (Balanced Sample)", 
  color = "red", count.color = "blue")

2.4.2 Weighted treatment effect

The package offers the option W to calculate the weighted average treatment effects. The weighting variable does not affect the estimation of fixed effects or factors. Only the weighted average treatment effects or weighted dynamic treatment effects are obtained by aggregating the treatment effects using the weight W.

simdata$Weight <- abs(rnorm(n = dim(simdata)[1]))
out.w <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id","time"), 
  force = "two-way", method = "ife", W = 'Weight',
  CV = FALSE, r = 2, se = TRUE, nboots = 200, parallel = TRUE)
#> Parallel computing ...

We can then visualize the weighted dynamic treatment effects using the inbuilt function plot, it by default shows the weighted dynamic treatment effects.

plot(out.w, main = "Estimated Weighted ATT")

2.5 Effect heterogeneity

We provide several methods for researchers to explore heterogeneous treatment effects (HTE).

2.5.1 Box plots

One way to understand HTE is to use a series of box plots to visualize the estimated individualistic treatment effects of observations under the treatment condition (by setting type = "box"). Although these effects are not identified at the individual observation level, their level of dispersion is informative of treatment effects heterogeneity at different (relative) time periods, as well as model performance.

plot(out.ife, type = "box", xlim = c(-15, 10))

2.5.2 CATT by calendar time

Another way to explore HTE is to investigate how the treatment effect evolves over time. In the plot below, the point estimates represents the ATTs by calendar time; the blue curve and band represent a lowess fit of the estimates and its 95% confidence interval, respectively; and the red horizontal dashed line represents the ATT (averaged over all time periods).

plot(out.ife, type = "calendar", xlim = c(1, 35))

2.5.3 CATT by a covariate

By setting type = "hte" or type = "heterogeneous, we can also plot the HTE by arbitrary covariates that are unaffected by the treatment. As before, the blue curve and band represent a lowess fit of the estimates and its 95% confidence interval, respectively. The red dashed line represents the ATT. The histogram at the bottom of the figure illustrates the distribution of the covariates, and can be turned off using show.count = FALSE. In our simulated case, the effect size is unrelated to the values of covariate X1.

plot(out.ife, type = "hte", covariate = "X1")

We can also plot the CATT when a covariate is discrete. To demonstrate this, we artificially create a moderating variable X3, which must be included in the outcome model and then specified in the heterogeneous treatment effect plot.

simdata$X3 <- sample(1:3, size = nrow(simdata), replace = TRUE)
out.ife.X3 <- fect(Y ~ D + X1 + X2 + X3, data = simdata, index = c("id","time"), 
                   method = "ife", r = 2, se = TRUE, seed = 123,
                   cores = 8, nboots = 1000, parallel = TRUE)
#> Parallel computing ...

As expected, there is not much effect heterogeneity along X3. In the resulting figure, we can also assign labels to the discrete values in the moderator.

plot(out.ife.X3, type="hte", covariate = "X3", 
     xlab = "", ylab = "Effet of D on Y",
     covariate.labels = c("USA", "China", "UK"),
     ylim = c(-2, 6))

Our next update will accommodate time-invariant covariates and allow users to explore effect heterogeneity around them.

2.6 Diagnostic tests

We provide three types of diagnostic tests: (1) a placebo test, (2) a test for (no) pretrend, and (3) a test for (no) carry-over effects. For each test, we support both the difference-in-means (DIM) approach and the equivalence approach. The details are provided in the paper.

2.6.1 Placebo tests

We provide a placebo test for a settled model—hence, cross-validation is not allowed—by setting placeboTest = TRUE. We specify a range of pretreatment periods as “placebo periods” in option placebo.period to remove observations in the specified range for model fitting, and then test whether the estimated ATT in this range is significantly different from zero. Below, we set c(-2, 0) as the placebo periods.

out.fect.p <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id", "time"),
  force = "two-way", parallel = TRUE, se = TRUE, CV = 0,
  nboots = 200, placeboTest = TRUE, placebo.period = c(-2, 0))

out.ife.p <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id", "time"),
  force = "two-way", method = "ife",  r = 2, CV = 0,
  parallel = TRUE, se = TRUE,
  nboots = 200, placeboTest = TRUE, placebo.period = c(-2, 0))

out.mc.p <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id", "time"),
  force = "two-way", method = "mc",  lambda = out.mc$lambda.cv, 
  CV = 0, parallel = TRUE, se = TRUE,
  nboots = 200, placeboTest = TRUE, placebo.period = c(-2, 0))

The placebo test conducts two types of tests:

t test. If t-test p-value is smaller than a pre-specified threshold (e.g. 5%), we reject the null of no-differences. Hence, the placebo test is deemed failed.

TOST. The TOST checks whether the 90% confidence intervals for estimated ATTs in the placebo period exceed a pre-specified range (defined by a threshold), or the equivalence range. A TOST p-value smaller than a pre-specified threshold suggests that the null of difference bigger than the threshold is rejected; hence, the placebo test is passed.

By default, the plot will display the p-value of the \(t\)-test (stats = "placebo.p"). Users can also add the p-value of a corresponding TOST test by setting stats = c("placebo.p","equiv.p"). A larger placebo p-value from a t-test and a smaller placebo TOST p-value are preferred.

plot(out.fect.p, cex.text = 0.8, stats = c("placebo.p","equiv.p"), 
     main = "Estimated ATT (TWFE)")

plot(out.ife.p, ylab = "Effect of D on Y", main = "Estimated ATT (IFE)", 
     cex.text = 0.8, stats = c("placebo.p","equiv.p"))

plot(out.mc.p, cex.text = 0.8, stats = c("placebo.p","equiv.p"),
     main = "Estimated ATT (MC)")

The results in the placebo test confirm that IFEct is a better model than MC for this particular DGP.

2.6.2 Tests for (no) pre-trend

We introduce two statistical tests for the presence of a pre-trend (or the lack thereof). The first test is an \(F\) test for zero residual averages in the pretreatment periods. The second test is a two-one-sided \(t\) (TOST) test, a type of equivalence tests.

F test. We offer a goodness-of-fit test (a variant of the \(F\) test) and to gauge the presence of pretreatment (differential) trends. A larger F-test p-value suggests a better pre-trend fitting. Users can specify a test range in option pre.periods. For example, pre.periods = c(-4,0) means that we test pretreatment trend of the last 5 periods prior to the treatment (from period -4 to period 0). If pre.period = NULL (default), all pretreatment periods in which the number of treated units exceeds the total number of treated units * proportion will be included in the test.

TOST. The TOST checks whether the 90% confidence intervals for estimated ATTs in the pretreatment periods (again, subject to the proportion option) exceed a pre-specified range, or the equivalence range. A smaller TOST p-value suggests a better pre-trend fitting. While users can check the values of confidence intervals, we give a visualization of the equivalence test. We can plot the pretreatment residual average with the equivalence confidence intervals by setting type = "equiv". Option tost.threshold sets the equivalence range (the default is \(0.36\sigma_{\epsilon}\) in which \(\sigma_{\epsilon}\) is the standard deviation of the outcome variable after two-way fixed effects are partialed out). By setting range = "both", both the minimum range (in gray) and the equivalence range (in red) are drawn.

On the topleft corner of the graph, we show several statistics of the user’s choice. User can choose which statistics to show by setting stats = c("none", "F.stat", "F.p", "F.equiv.p", "equiv.p") which corresponds to not showing any, the \(F\) statistic, the p-value for the \(F\) test, the p-value for the equivalence \(F\) test, the (maximum) p-value for the the TOST tests, respectively. For the gap plot, the default is stats = "none". For the equivalence plot, the default is stats = c("equiv.p, F.p"). Users can also change the labels of statistics using the stats.labs options. Users can adjust its position using the stats.pos option, for example stats.pos = c(-30, 4). To turn off the statistics, set stats = "none".

Below, we visualize the result of the equivalence test for each of the three estimators using our simulated data. These figures show that both the IFE and MC methods pass the equivalence test while the FE method does not.

plot(out.fect, type = "equiv", ylim = c(-4,4), 
     cex.legend = 0.6, main = "Testing Pre-Trend (FEct)", cex.text = 0.8)

plot(out.ife, type = "equiv", ylim = c(-4,4), 
     cex.legend = 0.6, main = "Testing Pre-Trend (IFEct)", cex.text = 0.8)

plot(out.mc, type = "equiv", ylim = c(-4,4),
     cex.legend = 0.6, main = "Testing Pre-Trend (MC)", cex.text = 0.8)

From the above plots, we see that FEct fails both tests; IFEct passes both tests using a conventional test size (5%); and MC fails the F tests, but passes the TOST (equivalence) test. Hence, we may conclude that IFEct is a more suitable model.

2.6.3 LOO pre-trend test

Instead of using estimated ATTs for periods prior to the treatment to test for pre-trends, we recommend users employ a leave-one-out (LOO) approach (loo = TRUE) to consecutively hide one pretreatment period (relative to the timing of the treatment) and repeatedly estimate the pseudo treatment effects for that pretreatment period. The LOO approach can be understood as an extension of the placebo test. It has the benefit of providing users with a more holistic view of whether the identifying assumptions likely hold. However, as the program needs to conduct uncertainty estimates for each turn, it is much more time-consuming than the original one.

out.fect.loo <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id","time"), 
  method = "fe", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200, loo = TRUE)
out.ife.loo <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id","time"), 
  method = "ife", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200, loo = TRUE)
out.mc.loo <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id","time"), 
  method = "mc", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200, loo = TRUE)

After the LOO estimation, one can plot these LOO pre-trends in the gap plot or the equivalence plot by setting loo = TRUE in the plot function.

plot(out.fect.loo, type = "equiv", ylim = c(-4,4), loo = TRUE,
     cex.legend = 0.6, main = "Testing Pre-Trend LOO (FEct)", cex.text = 0.8)

plot(out.ife.loo, type = "equiv", ylim = c(-4,4), loo = TRUE,
     cex.legend = 0.6, main = "Testing Pre-Trend LOO (IFEct)", cex.text = 0.8)

plot(out.mc.loo, type = "equiv", ylim = c(-4,4), loo = TRUE, 
     cex.legend = 0.6, main = "Testing Pre-Trend LOO (MC)", cex.text = 0.8)

Note that the LOO test usually takes lots of computational power. For our example, we find that the IFE estimator still passes both the F test and the equivalence test based on its LOO pre-trends, while the MC estimator fails both tests.

2.6.4 Exiting the treatment

fect allows the treatment to switch back and forth and provides diagnostic tools for this setting. After the estimation, we can visualize the period-wise ATTs relative to the exit of treatments by setting type = "exit" (one can still draw the classic gap plot by setting type = "gap"). The x-axis is then realigned based on the timing of the treatment’s exit, not onset, e.g., 1 represents 1 period after the treatment ends.

plot(out.fect, type = "exit", ylim = c(-2.5,4.5), main = "What Happens after the Treatment Switches Off?")

plot(out.ife, type = "exit", ylim = c(-2.5,4.5), main = "Exit Plot (IFE)")

plot(out.mc, type = "exit", ylim = c(-2.5,4.5), main = "Exit Plot (MC)")

2.6.5 Tests for (no) carryover effects

The idea of the placebo test can be extended to testing the presence of carryover effects. Instead of hiding a few periods right before the treatment starts, we hide a few periods right after the treatment ends. If carryover effects do not exist, we would expect the average prediction error in those periods to be close to zero. To perform the carryover test, we set the option carryoverTest = TRUE. We can treat a range of exit-treatment periods in option carryover.period to remove observations in the specified range for model fitting, and then test whether the estimated ATT in this range is significantly different from zero.

Below, we set carryover.period = c(1, 3). As we deduct the treatment effect from the outcome in simdata, we expect the average prediction error for these removed periods to be close to zero.

out.fect.c <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id", "time"),
  force = "two-way", parallel = TRUE, se = TRUE, CV = 0, 
  nboots = 200, carryoverTest = TRUE, carryover.period = c(1, 3))

out.ife.c <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id", "time"),
  force = "two-way", method = "ife", r = 2, CV = 0, 
  parallel = TRUE, se = TRUE,
  nboots = 200, carryoverTest = TRUE, carryover.period = c(1, 3))

out.mc.c <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id", "time"),
  force = "two-way", method = "mc",  lambda = out.mc$lambda.cv, 
  CV = 0, parallel = TRUE, se = TRUE, 
  nboots = 200, carryoverTest = TRUE, carryover.period = c(1, 3))

Like the placebo test, the plot will display the p-value of the carryover effect test (stats = "carryover.p"). Users can also add the p-value of a corresponding TOST test by setting stats = c("carryover.p","equiv.p").

plot(out.fect.c, type = "exit", ylim = c(-2.5,4.5), 
          cex.text = 0.8, main = "Carryover Effects (FE)")

plot(out.ife.c, type = "exit", ylim = c(-2.5,4.5), 
          cex.text = 0.8, main = "Carryover Effects (IFE)")

plot(out.mc.c, type = "exit", ylim = c(-2.5,4.5), 
          cex.text = 0.8, main = "Carryover Effects (MC)")

Once again, the IFE estimator outperforms the other two.

Using real-world data, researchers will likely find that carryover effects exist. If such effects are limited, researchers can consider removing a few periods after the treatment ended for the treated units from the first-stage estimation (using the carryover.period option) and re-estimated the model (and re-conduct the test). We provide such an example in the paper. Here, we illustrate the option using simdata.

out.ife.rm.test <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id", "time"),
  force = "two-way", method = "ife", r = 2, CV = 0, 
  parallel = TRUE, se = TRUE,  carryover.rm = 3,
  nboots = 200, carryoverTest = TRUE, carryover.period = c(1, 3))# remove three periods

plot(out.ife.rm.test, cex.text = 0.8, stats.pos = c(5, 2))

In the above plot, the three periods in blue are droppred from the first-stage estimation of the factor model while the periods in red are reserved for the (no) carryover effects test.

2.7 Cumulative effects

Users can use effect() to calculate cumulative treatment effects. The behavior of effect() is similar to the function of the same name in gsynth. - Calculation of cumulative effects will need unit-time level bootstrap results. Choose the option keep.sims=TRUE to record them.

out <- fect(Y ~ D + X1 + X2, data = simgsynth, index = c("id","time"),
                        method = "gsynth", force = "two-way", CV = TRUE, r = c(0, 5),
                        se = TRUE, nboots = 200, vartype = 'bootstrap',
                        parallel = FALSE, keep.sims=TRUE)
#> Cross-validating ...
#> Criterion: Mean Squared Prediction Error
#> Interactive fixed effects model...
#> Cross-validating ...
#> r = 0; sigma2 = 1.84865; IC = 1.02023; PC = 1.74458; MSPE = 2.37280
#> r = 1; sigma2 = 1.51541; IC = 1.20588; PC = 1.99818; MSPE = 1.71743
#> r = 2; sigma2 = 0.99737; IC = 1.16130; PC = 1.69046; MSPE = 1.14540
#> *
#> r = 3; sigma2 = 0.94664; IC = 1.47216; PC = 1.96215; MSPE = 1.15032
#> r = 4; sigma2 = 0.89411; IC = 1.76745; PC = 2.19241; MSPE = 1.21397
#> r = 5; sigma2 = 0.85060; IC = 2.05928; PC = 2.40964; MSPE = 1.23876
#> 
#>  r* = 2
#> 
#> Can't calculate the F statistic because of insufficient treated units.
cumu.out <- effect(out)

Print and plot cumulative effects

print(cumu.out)
#> 
#> Overall cumulative effect:
#>  [1]  1.236  2.866  5.578  9.045 14.785 20.065 28.502 36.342 45.797 55.435
#> 
#> Period-by-period cumulative effect:
#>       ATT   S.E. CI.lower CI.upper p.value
#> 1   1.236 0.8291  -0.2039    2.848    0.16
#> 2   2.866 1.2167   0.1849    4.877    0.04
#> 3   5.578 1.1270   3.2334    7.608    0.00
#> 4   9.045 1.2726   6.3574   11.097    0.00
#> 5  14.785 1.4454  11.5459   17.154    0.00
#> 6  20.065 1.3558  17.2273   22.298    0.00
#> 7  28.502 1.6676  25.3248   31.361    0.00
#> 8  36.342 1.7104  32.4730   39.082    0.00
#> 9  45.797 1.6984  41.8843   48.439    0.00
#> 10 55.435 2.0912  50.6100   59.217    0.00
#> NULL
plot(cumu.out)

Users can choose to calculate by-period average effects by setting cumu=FALSE.

effect(out, cumu=FALSE)
#> 
#> Overall cumulative effect:
#>  [1] 1.236 1.630 2.712 3.467 5.740 5.280 8.436 7.840 9.455 9.639
#> 
#> Period-by-period cumulative effect:
#>      ATT   S.E. CI.lower CI.upper p.value
#> 1  1.236 0.8291  -0.2039    2.848    0.16
#> 2  1.630 0.7986  -0.0948    3.112    0.08
#> 3  2.712 0.6297   1.4700    3.936    0.00
#> 4  3.467 0.6770   2.1126    4.714    0.00
#> 5  5.740 0.6849   4.2557    6.948    0.00
#> 6  5.280 0.3867   4.6240    6.026    0.00
#> 7  8.436 0.5291   7.5012    9.410    0.00
#> 8  7.840 0.6340   6.5970    8.937    0.00
#> 9  9.455 0.5291   8.4700   10.554    0.00
#> 10 9.639 0.9232   8.2871   11.765    0.00

Calculate the cumulative effect of certain units at certain periods.

effect(out, cumu=TRUE, id=c(101,102,103), period=c(1,5))
#> 
#> Overall cumulative effect:
#> [1]  2.249  4.642  6.949  9.530 14.275
#> 
#> Period-by-period cumulative effect:
#>      ATT   S.E. CI.lower CI.upper p.value
#> 1  2.249 0.9759   0.2055    4.042 0.03015
#> 2  4.642 0.5503   3.5683    5.498 0.00000
#> 3  6.949 1.1317   4.8012    9.226 0.00000
#> 4  9.530 1.7910   6.1976   12.933 0.00000
#> 5 14.275 2.2008   9.8936   17.864 0.00000

effect() also accepts results from other estimation and inference methods. For example, we can use matrix completion:

out_mc <- fect(Y ~ D + X1 + X2, data = simgsynth, index = c("id","time"),
                        method = "mc", force = "two-way", CV = TRUE, r = c(0, 5),
                        se = TRUE, nboots = 200, vartype = 'bootstrap',
                        parallel = FALSE, keep.sims=TRUE)
#> Cross-validating ...
#> Criterion: Mean Squared Prediction Error
#> Matrix completion method...
#> lambda.norm = 1.00000; MSPE = 2.09509; MSPTATT = 0.69437; MSE = 1.97846
#> lambda.norm = 0.42170; MSPE = 1.69635; MSPTATT = 0.31226; MSE = 1.02780
#> *
#> lambda.norm = 0.17783; MSPE = 1.91180; MSPTATT = 0.09063; MSE = 0.31477
#> lambda.norm = 0.07499; MSPE = 2.07144; MSPTATT = 0.01738; MSE = 0.06379
#> lambda.norm = 0.03162; MSPE = 2.03016; MSPTATT = 0.00308; MSE = 0.01170
#> 
#>  lambda.norm* = 0.421696503428582
#> 
#> Can't calculate the F statistic because of insufficient treated units.
plot(effect(out_mc))

We can also use jackknife instead of bootstrap for inference:

out_jack <- fect(Y ~ D + X1 + X2, data = simgsynth, index = c("id","time"),
                        method = "mc", force = "two-way", CV = TRUE, r = c(0, 5),
                        se = TRUE, nboots = 200, vartype = 'jackknife',
                        parallel = FALSE, keep.sims=TRUE)
#> Cross-validating ...
#> Criterion: Mean Squared Prediction Error
#> Matrix completion method...
#> lambda.norm = 1.00000; MSPE = 2.46863; MSPTATT = 0.69437; MSE = 1.97846
#> lambda.norm = 0.42170; MSPE = 1.93789; MSPTATT = 0.31226; MSE = 1.02780
#> *
#> lambda.norm = 0.17783; MSPE = 2.10134; MSPTATT = 0.09063; MSE = 0.31477
#> lambda.norm = 0.07499; MSPE = 2.23331; MSPTATT = 0.01738; MSE = 0.06379
#> lambda.norm = 0.03162; MSPE = 2.31343; MSPTATT = 0.00308; MSE = 0.01170
#> 
#>  lambda.norm* = 0.421696503428582
#> 
#> Can't calculate the F statistic because of insufficient treated units.
plot(effect(out_jack))

2.8 Other estimators

The counterfacutal/imputation estimator framework can be extended to more settings.

Complex Fixed Effects. When there exists more dimensions of fixed effects in addition to the unit and time fixed effects, we can resort to the “cfe” (complex fixed effects) estimator to impute the counterfactual based on a linear model with multiple levels of fixed effects.

Note, fect allows the method to run without requiring a baseline two-way fixed effects model. That is, users do not need to set force = "two-way" to impute.

It accepts two options: sfe specifies simple (additive) fixed effects in addition to the unit and time fixed effects and cfe receives a list object and each component in the list is a vector of length 2.

The value of the first element of each component is the name of group variable for which fixed effects are to be estimated (e.g. unit names); the value of the second element is the name of a regressor (e.g., a time trend). For example, we can estimate a model with an additional fixed effects FE3 along with a unit-specific time trend.

simdata[,"FE3"] <- sample(c(1,2,3,4,5), size = dim(simdata)[1], replace = TRUE)
out.cfe <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id","time"), 
  method = "cfe", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200,
  sfe = c("FE3"), cfe = list(c("id","time")))
plot(out.cfe)

Polynomial. Sometimes researchers may want to include unit-specific time trends in the model estimated using non-treated data. We can set method = "polynomial" to achieve this.

In addition, By setting degree = 2, we can estimate the ATT based on a linear model with unit and time fixed effects, along with a unit-specific quadratic time trend. Similar to cfe estimator, a two-way fixed effects model, while encouraged, is not required.

out.poly <- fect(Y ~ D + X1 + X2, data = simdata, index = c("id","time"), 
  method = "polynomial", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200,
  degree = 2)
plot(out.poly)

2.9 Other options

We provide a few other options for estimation and visualization.

2.9.1 More visualization options

The plot function shipped in fect offers some options that help to improve the visualization.

We can remove the bar plot at the bottom of the plot by setting show.count = FALSE

plot(out.ife, show.count = FALSE)

By setting the option type = "counterfactual", we visualize the period-wise average treated and counterfactual outcomes with shaded confidence intervals.

plot(out.ife, type = "counterfactual")

By setting the option type = "status", we can visualize the treatment status of all observations. We only present the label of the time by setting axis.lab = "time".

plot(out.fect, type = 'status', axis.lab = "time", cex.axis  = 0.6)

For the placebo test, the manually hided observations are marked in cyan. We can show only a sub-group’s treatment status by specifying the option id to certain units.

plot(out.fect.p, type = 'status', axis.lab = "both", id = c(101:120), cex.axis  = 0.6)

For the carryover test, the manually hidden observations are marked in light red. We can also remove grid lines by setting gridOff = TRUE.

plot(out.fect.c, type = 'status', axis.lab = "off", gridOff = TRUE)

For the carryover test with removed observation, the removed observations are marked in yellow.

plot(out.ife.rm.test, type = 'status', axis.lab = "off", gridOff = TRUE)

2.9.2 Average cohort treatment effect

fect allows us to estimate and visualize the ATTs for sub-groups of treated units. For example, it can draw the gap plot for units that adopt the treatment at the same time under staggered adoption, which is defined as “Cohort” in Sun & Abraham (2021). Our simulated dataset is not ideal to demonstrate this functionality because the treatment switches on and off. To improve feasibility, we define a cohort as a group of treated units that first adopt the treatment at the same time.

panelview(Y ~ D, data = simdata, index = c("id","time"), by.timing = TRUE,
  axis.lab = "time", xlab = "Time", ylab = "Unit", 
  background = "white", main = "Simulated Data: Treatment Status")

The get.cohort() function (originally from the paneltools package) can generate a new variable “Cohort” based on the timing when treated units first get treated. The new version of fect incorporates the feature: users no longer need to install any complementary packages to replicate the tutorial.

#  devtools:: install_github("xuyiqing/paneltools" if not already installed
simdata.cohort <- get.cohort(data = simdata,D = 'D',index = c("id","time"))
print(table(simdata.cohort[,'Cohort']))
#> 
#>  Cohort:1 Cohort:10 Cohort:11 Cohort:12 Cohort:13 Cohort:14 Cohort:15 Cohort:16 
#>        35       245       105        35       175       175       210       245 
#> Cohort:17 Cohort:18 Cohort:19  Cohort:2 Cohort:20 Cohort:21 Cohort:22 Cohort:23 
#>       210       210       140       280       315       105       140       210 
#> Cohort:24 Cohort:25 Cohort:26 Cohort:27 Cohort:28 Cohort:29  Cohort:3 Cohort:30 
#>       105       280       175       210        70       105       175        35 
#> Cohort:31 Cohort:33 Cohort:34 Cohort:35  Cohort:4  Cohort:5  Cohort:6  Cohort:7 
#>       105        70       175       140        70       140       105       210 
#>  Cohort:8  Cohort:9   Control 
#>       105       140      1750

We can also pass a list of intervals for first get-treated time into the entry.time option of get.cohort(). For example, we can categorize all treated units into the group that adopts the treatment between time 21 and 27, and the group that adopts the treatment in time 30 and 33.

simdata.cohort2 <- get.cohort(data = simdata,D = 'D',index = c("id","time"),
                               entry.time = list(c(21,27),c(30,33)))
print(table(simdata.cohort2[,'Cohort']))
#> 
#> Cohort:21-27 Cohort:30-33 Cohort:Other      Control 
#>         1225          210         3815         1750

By setting the option group = "Cohort", fect estimates the ATT for each specified sub-group and saves it for further visualization.

out.ife.g <- fect(Y ~ D + X1 + X2, data = simdata.cohort, index = c("id","time"), 
          force = "two-way", method = "ife", CV = TRUE, r = c(0, 5), 
          se = TRUE, nboots = 200, parallel = TRUE, group = 'Cohort')
out.ife.g.p <- fect(Y ~ D + X1 + X2, data = simdata.cohort, index = c("id","time"), 
          force = "two-way", method = "ife", CV = FALSE, 
          placeboTest = TRUE, placebo.period = c(-2,0), 
          se = TRUE, nboots = 200, parallel = TRUE, group = 'Cohort')

Then one can draw the gap plot, as well as the equivalence plot, for each sub-group. Here we present the gap plot for Cohort 22.

plot(out.ife.g, show.group = "Cohort:22", 
          xlim = c(-15, 10), ylim = c(-10, 10))

2.10 Additional notes

  1. By default, the program will drop the units that have no larger than 5 observations under control, which is the reason why sometimes there are less available units in the placebo test or carryover test than in the original estimation. We can specify a preferred criteria in the option min.T0 (default to 5). As a rule of thumb for the IFE estimator, the minimum number of observations under control for a unit should be larger than the specified number of factor r.

  2. We can get replicable results by setting the option seed to a certain integer, no matter whether the parallel computing is used.

  3. When na.rm = FALSE (default), the program allows observations to have missing outcomes \(Y\) or covariates \(X\) but decided treatment statuses \(D\). Otherwise the program will drop all observations that have missing values in outcomes, treatments, or covariates.