PR for #42 enable estimation using fixest

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: eventstudyr
 Title: Estimation and Visualization of Linear Panel Event Studies
-Version: 1.1.5
+Version: 1.2.0
 Authors@R: 
     c(person(given = "Simon",
            family = "Freyaldenhoven",
@@ -43,6 +43,7 @@ Imports:
     data.table,
     dplyr,
     estimatr,
+    fixest,
     ggplot2,
     MASS,
     rlang,
@@ -53,7 +54,7 @@ VignetteBuilder:
     knitr
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 Suggests: 
     rmarkdown,
     knitr,

NAMESPACE

@@ -19,12 +19,20 @@ importFrom(data.table,setnames)
 importFrom(data.table,setorder)
 importFrom(data.table,setorderv)
 importFrom(data.table,shift)
+importFrom(dplyr,rename)
+importFrom(dplyr,select)
+importFrom(fixest,coeftable)
+importFrom(fixest,feols)
+importFrom(fixest,fitstat)
+importFrom(fixest,ssc)
 importFrom(pracma,inv)
 importFrom(pracma,pinv)
 importFrom(rlang,.data)
 importFrom(stats,as.formula)
+importFrom(stats,coef)
 importFrom(stats,pnorm)
 importFrom(stats,qchisq)
 importFrom(stats,qnorm)
 importFrom(stats,reformulate)
 importFrom(stats,setNames)
+importFrom(stats,vcov)

R/AddSuptBand.R

@@ -17,6 +17,8 @@
 #' for each event-study coefficient.
 #' @import estimatr
 #' @importFrom MASS mvrnorm
+#' @importFrom fixest coeftable
+#' @importFrom dplyr rename
 #' @keywords internal
 #' @noRd
 #'
@@ -47,16 +49,18 @@
 #'   eventstudy_coefficients = eventstudy_estimates$arguments$eventstudy_coefficients
 #')
 
-AddSuptBand <- function(estimates, num_sim = 1000, conf_level = .95, eventstudy_coefficients) {
+AddSuptBand <- function(model_estimates, num_sim = 1000, conf_level = .95, eventstudy_coefficients) {
 
-    if (! class(estimates) %in% c("lm_robust", "iv_robust")) {
+    if (! class(model_estimates) %in% c("lm_robust", "iv_robust", "fixest")) {
         stop("estimates is not a data frame with coefficient estimates and standard errors")
     }
     if (! is.numeric(num_sim) | num_sim %% 1 != 0 | num_sim <= 0) {stop("num_sim should be a natural number.")}
     if (! is.numeric(conf_level) | conf_level < 0 | conf_level > 1) {stop("conf_level should be a real number between 0 and 1, inclusive.")}
     if (! is.character(eventstudy_coefficients)) {stop("eventstudy_coefficients should be a character.")}
 
-    vcov_matrix_all <- estimates$vcov
+    fixest = class(model_estimates) == "fixest"
+
+    vcov_matrix_all <- if(fixest){vcov(model_estimates)} else {model_estimates$vcov}
     v_terms_to_keep <- colnames(vcov_matrix_all) %in% eventstudy_coefficients
     vcov_matrix <- vcov_matrix_all[v_terms_to_keep, v_terms_to_keep]
 
@@ -75,12 +79,18 @@ AddSuptBand <- function(estimates, num_sim = 1000, conf_level = .95, eventstudy_
         critical_value = t[floor(conf_level_num_sim) + 1]
     }
 
-    df_estimates_tidy <- estimatr::tidy(estimates)
+    df_estimates_tidy <- if(fixest){
+        coef_table <- model_estimates |> 
+        fixest::coeftable() |> 
+        as.data.frame()
+        coef_table$term <- rownames(coef_table)
+        coef_table |> 
+        dplyr::rename(estimate = Estimate, std.error = `Std. Error`) |> 
+        dplyr::select(term, estimate, std.error)
+    } else {estimatr::tidy(model_estimates)}
 
     df_estimates_tidy["suptband_lower"] <- df_estimates_tidy$estimate - (critical_value * df_estimates_tidy$std.error)
     df_estimates_tidy["suptband_upper"] <- df_estimates_tidy$estimate + (critical_value * df_estimates_tidy$std.error)
 
-
     return(df_estimates_tidy)
-
 }

R/EventStudy.R

@@ -34,6 +34,7 @@
 #' when there are anticipation effects. If set to FALSE, does not make the switch. Defaults to TRUE.
 #' @param allow_duplicate_id If TRUE, the function estimates a regression where duplicated ID-time rows are weighted by their duplication count. If FALSE, the function raises an error if duplicate unit-time keys exist in the input data. Default is FALSE.
 #' @param avoid_internal_copy If TRUE, the function avoids making an internal deep copy of the input data, and instead directly modifies the input data.table. Default is FALSE.
+#' @param kernel Accepts one of "estimatr" or "fixest". If "estimatr" is specified, uses the estimatr package for estimation. If "fixest" is specified, uses the fixest package for estimation. Defaults to "estimatr" (deprecated - will change to "fixest" in a future release).
 #'
 #' @return A list that contains, under "output", the estimation output as an lm_robust object, and under "arguments", the arguments passed to the function.
 #' @import dplyr
@@ -54,7 +55,8 @@
 #'     idvar = "id",
 #'     timevar = "t",
 #'     pre = 0, post = 3,
-#'     normalize = -1
+#'     normalize = -1,
+#'     kernel = "fixest"
 #'   )
 #'
 #' ### Access estimated model
@@ -63,7 +65,8 @@
 #' summary(eventstudy_model$output)
 #'
 #' ### data.frame of estimates
-#' estimatr::tidy(eventstudy_model$output)
+#' fixest::coeftable(eventstudy_model$output) # for kernel='fixest'
+#' # estimatr::tidy(eventstudy_model$output) # for kernel='estimatr'
 #'
 #' ### Access arguments
 #' eventstudy_model$arguments
@@ -83,6 +86,7 @@
 #'     pre  = 2, overidpre  = 4,
 #'     normalize = - 3,
 #'     cluster = TRUE,
+#'     kernel = "fixest",
 #'     anticipation_effects_normalization = TRUE
 #'   )
 #'
@@ -100,7 +104,8 @@
 #'     FE = TRUE, TFE = TRUE,
 #'     post = 0, overidpost = 0,
 #'     pre  = 0, overidpre  = 0,
-#'     cluster = TRUE
+#'     cluster = TRUE,
+#'     kernel = "fixest"
 #'   )
 #'
 #' summary(eventstudy_model_static$output)
@@ -117,7 +122,8 @@
 #'     idvar = "id",
 #'     timevar = "t",
 #'     pre = 0, post = 3,
-#'     normalize = -1
+#'     normalize = -1,
+#'     kernel = "fixest"
 #'   )
 #'
 #' summary(eventstudy_model_unbal$output)
@@ -136,7 +142,8 @@
 #'     post = 2, overidpost = 1,
 #'     pre  = 0, overidpre  = 3,
 #'     normalize = -1,
-#'     cluster = TRUE
+#'     cluster = TRUE,
+#'     kernel = "fixest"
 #'   )
 #'
 #' summary(eventstudy_model_iv$output)
@@ -145,10 +152,13 @@
 EventStudy <- function(estimator, data, outcomevar, policyvar, idvar, timevar, controls = NULL,
                        proxy = NULL, proxyIV = NULL, FE = TRUE, TFE = TRUE, post, overidpost = 1, pre, overidpre = post + pre,
                        normalize = -1 * (pre + 1), cluster = TRUE, anticipation_effects_normalization = TRUE,
-                       allow_duplicate_id = FALSE, avoid_internal_copy = FALSE) {
+                       allow_duplicate_id = FALSE, avoid_internal_copy = FALSE, kernel = "estimatr") {
 
     # Check for errors in arguments
     if (! estimator %in% c("OLS", "FHS")) {stop("estimator should be either 'OLS' or 'FHS'.")}
+    if (! kernel %in% c("estimatr", "fixest")) {stop("kernel should be either 'estimatr' or 'fixest'.")}
+    if (missing(kernel)) {warning("Argument 'kernel' was not specified; using 'estimatr' as default; we strongly recommend explicitly specifying a kernel because the default is scheduled to change.")}
+    if (kernel == "estimatr") {warning("'estimatr' selected as kernel. We no longer maintain it and will depreciate it in a future release. We recommend using 'fixest' instead.")}
     if (! is.data.frame(data))            {stop("data should be a data frame.")}
     for (var in c(idvar, timevar, outcomevar, policyvar)) {
         if ((! is.character(var))) {
@@ -336,14 +346,17 @@ EventStudy <- function(estimator, data, outcomevar, policyvar, idvar, timevar, c
     }
 
     if (estimator == "OLS") {
-        event_study_formula <- PrepareModelFormula(estimator, outcomevar, str_policy_vars,
-                                                   static, controls, proxy, proxyIV)
-
-        output       <- EventStudyOLS(event_study_formula, data, idvar, timevar, FE, TFE, cluster)
+        formula <- PrepareModelFormula(estimator, outcomevar, str_policy_vars,
+                                       static, controls, proxy, proxyIV,
+                                       kernel, idvar, timevar, FE, TFE)
+
+        if (kernel == "estimatr") {
+            output <- EventStudyOLS(formula, data, idvar, timevar, FE, TFE, cluster)
+        } else if (kernel == "fixest") {
+            output <- EventStudyFEOLS(formula, data, idvar, timevar, FE, TFE, cluster)
+        }
         coefficients <- str_policy_vars
-    }
-    if (estimator == "FHS") {
-
+    } else if (estimator == "FHS") {
         if (is.null(proxyIV)) {
             Fstart <- 0
             str_fd_leads <- str_policy_vars[grepl("^z_fd_lead", str_policy_vars)]
@@ -360,10 +373,15 @@ EventStudy <- function(estimator, data, outcomevar, policyvar, idvar, timevar, c
                            ". To specify a different proxyIV use the proxyIV argument."))
         }
 
-        event_study_formula <- PrepareModelFormula(estimator, outcomevar, str_policy_vars,
-                                                   static, controls, proxy, proxyIV)
-
-        output       <- EventStudyFHS(event_study_formula, data, idvar, timevar, FE, TFE, cluster)
+        formula <- PrepareModelFormula(estimator, outcomevar, str_policy_vars,
+                                       static, controls, proxy, proxyIV,
+                                       kernel, idvar, timevar, FE, TFE)
+
+        if (kernel == "estimatr") {
+            output <- EventStudyFHS(formula, data, idvar, timevar, FE, TFE, cluster)
+        } else if (kernel == "fixest") {
+            output <- EventStudyFEOLS_FHS(formula, data, idvar, timevar, FE, TFE, cluster)
+        }
         coefficients <- dplyr::setdiff(str_policy_vars, proxyIV)
     }
 

R/EventStudyFHS.R

@@ -12,7 +12,8 @@
 #'
 #' @return A data.frame that contains the estimates for the event study coefficients.
 #' @import estimatr
-#' @importFrom stats qnorm pnorm
+#' @importFrom fixest feols
+#' @importFrom stats qnorm pnorm coef
 #' @keywords internal
 #' @noRd
 #'
@@ -127,3 +128,51 @@ EventStudyFHS <- function(prepared_model_formula, prepared_data,
   return(fhs_output)
 }
 
+
+EventStudyFEOLS_FHS <- function(formula, prepared_data,
+                                idvar, timevar, FE, TFE, cluster) {
+
+    if (! inherits(formula, "formula")) {stop("formula should be a formula")}
+    if (! is.data.frame(prepared_data)) {stop("data should be a data frame.")}
+    if (! is.character(idvar)) {stop("idvar should be a character.")}
+    if (! is.character(timevar)) {stop("timevar should be a character.")}
+    if (! is.logical(FE)) {stop("FE should be either TRUE or FALSE.")}
+    if (! is.logical(TFE)) {stop("TFE should be either TRUE or FALSE.")}
+    if (! is.logical(cluster)) {stop("cluster should be either TRUE or FALSE.")}
+    if (FE & !cluster) {stop("cluster=TRUE required when FE=TRUE.")}
+
+    if (cluster) {
+        vcov_arg <- as.formula(paste0("~", idvar))
+    } else {
+        vcov_arg <- "HC1"
+    }
+
+    fhs_output <- fixest::feols(
+        fml = formula,
+        data = prepared_data,
+        vcov = vcov_arg
+    )
+
+    if (FE & cluster) {
+        coefs <- coef(fhs_output)
+        N <- fhs_output$nobs
+        n <- length(unique(prepared_data[[idvar]]))
+
+        if (TFE) {
+            K <- fhs_output$fixef_sizes[[timevar]] + length(coefs)
+        } else {
+            K <- 1 + length(coefs)
+        }
+
+        adjustment_factor <- (N - K) / (N - n - K + 1)
+        fhs_output$se <- fhs_output$se / sqrt(adjustment_factor)
+        fhs_output$cov.scaled <- fhs_output$cov.scaled / adjustment_factor
+
+        fhs_output$tstat <- coefs / fhs_output$se
+        fhs_output$pvalue <- 2 * stats::pnorm(abs(fhs_output$tstat), lower.tail = FALSE)
+        fhs_output$conf.low <- coefs - stats::qnorm(0.975) * fhs_output$se
+        fhs_output$conf.high <- coefs + stats::qnorm(0.975) * fhs_output$se
+    }
+
+    return(fhs_output)
+}

R/EventStudyOLS.R

@@ -1,6 +1,6 @@
 #' Runs Ordinary Least Squares (OLS) with optional fixed effects and clustering
 #'
-#' @param prepared_model_formula A formula object created in [PrepareModelFormula()] that is passed to [EventStudy()].
+#' @param prepared_model_formula A formula object created in `PrepareModelFormula()` that is passed to [EventStudy()].
 #' @param prepared_data Data frame containing all of the parameters required for [EventStudy()] plus leads and
 #' lags of the first differenced policy variable and leads and lags of the policy variable.
 #' @param idvar Character indicating column of units.
@@ -12,6 +12,7 @@
 #'
 #' @return A data.frame that contains the estimates for the event study coefficients.
 #' @import estimatr
+#' @importFrom fixest feols ssc
 #' @keywords internal
 #' @noRd
 #'
@@ -123,3 +124,23 @@ EventStudyOLS <- function(prepared_model_formula, prepared_data,
 
     return(ols_output)
 }
+
+EventStudyFEOLS <- function(formula, prepared_data,
+                          idvar, timevar, FE, TFE, cluster) {
+
+    if (cluster) {
+        vcov_fixest <- as.formula(paste0("~", idvar))
+        small_sample_correction <- fixest::ssc(K.fixef = "full")
+    } else {
+        vcov_fixest <- "iid"
+        small_sample_correction <- fixest::ssc()
+    }
+
+    ols_output <- fixest::feols(
+        fml = formula,
+        data = prepared_data,
+        vcov = vcov_fixest,
+        ssc = small_sample_correction
+    )
+    return(ols_output)
+}