Add Wild Cluster Bootstrap Support

DESCRIPTION

@@ -32,7 +32,11 @@ Suggests:
   ggplot2,
   knitr,
   rmarkdown,
-  tinytest
+  tinytest, 
+  fwildclusterboot
+Remotes: 
+  kylebutts/fixest/tree/sparse-matrix
+  s3alfisc/fwildclusterboot/tree/etwfe-support
 Encoding: UTF-8
 RoxygenNote: 7.2.3
 URL: https://grantmcdermott.com/etwfe/

R/emfx.R

@@ -25,14 +25,29 @@
 ##'   dataset. But you may want to keep them for presentation reasons (e.g.,
 ##'   plotting an event-study); though be warned that this is strictly
 ##'   performative. This argument will only be evaluated if `type = "event"`.
+##' @param bootstrap Logical. FALSE by default. Should inference be conducted 
+##'        via analytical standard errors or a  wild bootstrap? To run the 
+##'        bootstrap, the `fwildclusterboot` package needs to be installed. 
+##'        If you want to run a bootstrap, you need to 
+##'        pass a number of bootstrap iterations via the `...` through
+##'        `emfx()`, e.g. `B = 9999`. The bootstrap is currently only supported 
+##'        for `type == 'simple'` and clustered errors. 
 ##' @param ... Additional arguments passed to
-##'   [`marginaleffects::marginaleffects`]. For example, you can pass `vcov =
-##'   FALSE` to dramatically speed up estimation times of the main marginal
-##'   effects (but at the cost of not getting any information about standard
-##'   errors; see Performance tips below). Another potentially useful
-##'   application is testing whether heterogeneous treatment effects (i.e. the
-##'   levels of any `xvar` covariate) are equal by invoking the `hypothesis`
-##'   argument, e.g. `hypothesis = "b1 = b2"`.
+##'   [`marginaleffects::marginaleffects`] or 
+##'   [`fwildclusterboot::boot_aggregate`] (the ladder is only relevant when 
+##'   `bootstrap = TRUE`). For example, you can pass `vcov = FALSE` 
+##'   to `marginaleffects` to dramatically speed up estimation times of the
+##'   main marginal effects (but at the cost of not getting any information 
+##'   about standard errors; see Performance tips below). 
+##'   Another potentially useful application is testing whether 
+##'   heterogeneous treatment effects (i.e. the levels of any `xvar` covariate)
+##'   are equal by invoking the `hypothesis` argument, e.g. 
+##'   `hypothesis = "b1 = b2"`. For the bootstrap, you can e.g. pass along the 
+##'   number of bootstrap iterations, the 
+##'   "bootcluster" variable (relevant for the subcluster bootstrap, via the 
+##'   `bootcluster` argument) or the number of threads to use 
+##'   (via the `nthreads` argument). For a comprehensive list 
+##'   of arguments, check `?fwildclusterboot::boot_aggregate()`.
 ##' @return A `slopes` object from the `marginaleffects` package.
 ##' @seealso [marginaleffects::slopes()]
 ##' @inherit etwfe return examples 
@@ -46,6 +61,7 @@ emfx = function(
     by_xvar = "auto",
     collapse = "auto",
     post_only = TRUE,
+    bootstrap = FALSE, 
     ...
 ) {
 
@@ -140,14 +156,71 @@ emfx = function(
 
   if (by_xvar) by_var = c(by_var, xvar)
 
-  mfx = marginaleffects::slopes(
-    object,
-    newdata = dat,   
-    wts = "N",
-    variables = ".Dtreat",
-    by = by_var,
-    ...
-  )
+  if(!bootstrap){
+
+    mfx = marginaleffects::slopes(
+      object,
+      newdata = dat,   
+      wts = "N",
+      variables = ".Dtreat",
+      by = by_var,
+      ...
+    )
+
+  } else {
+
+    if(!requireNamespace("fwildclusterboot")){
+      stop(
+        "To run the bootstrap, the `fwildclusterboot` package",
+        "needs to be installed. ", 
+        "However, the package cannot be found. Please install it by",
+        "running `install.packages('fwildclusterboot')`\n."
+      ) 
+    }
+
+    if(mod$method != "feols"){
+      stop(
+        "Bootstrapping is only supported for models estimated", 
+        "via OLS / `feols()`."
+      )
+    }
+
+    if(type == "simple"){
+      agg = c("ATT"="^.Dtreat:first.treat::[0-9]{4}:year::[0-9]{4}$")
+    } else {
+      stop("Only type = 'auto' is currently supported.")
+    }
+
+    clustid_long <- attr(object$cov.scaled, "type")
+    clustid <- sub(".*\\((.*?)\\).*", "\\1", clustid_long)
+    ssc <- attr(object$cov.scaled, "ssc")
+    boot_ssc <- fwildclusterboot::boot_ssc(
+      adj = ssc$adj, 
+      fixef.K = "none", 
+      cluster.adj = ssc$cluster.adj, 
+      cluster.df = ssc$cluster.df
+    )
+
+    if(ssc$fixef.K != "none"){
+      warning(
+        paste0("The bootstrap does not support the ssc() argument", 
+        "`fixef.K='", ssc$fixef.K, "'`."), 
+        "Using `fixef.K='none' instead.", 
+        "This will lead to a slightly different non-bootstrapped t-statistic`",
+        "but will not affect bootstrapped p-values and CIs.\n"
+      )
+    }
+
+    mfx <- fwildclusterboot::boot_aggregate(
+      x = mod, 
+      agg = agg, 
+      ssc = boot_ssc, 
+      clustid = clustid, 
+      ...
+    )
+
+  }
+
 
 
   return(mfx)

R/etwfe.R

@@ -234,7 +234,7 @@ etwfe = function(
   if (is.numeric(ivar)) ivar = names(data)[ivar]
   xvar = eval(substitute(xvar), nl, parent.frame())
   if (is.numeric(xvar)) xvar = names(data)[xvar]
-
+  
   if (is.null(gvar)) stop("A non-NULL `gvar` argument is required.\n")
   if (is.null(tvar)) stop("A non-NULL `tvar` argument is required.\n")
   if (!is.null(family)) ivar = NULL
@@ -278,13 +278,13 @@ etwfe = function(
   ref_string = paste0(", ref = ", gref)
 
   if (is.null(tref)) {
-      tref = min(data[[tvar]], na.rm = TRUE)
+    tref = min(data[[tvar]], na.rm = TRUE)
   } else if (!(tref %in% unique(data[[tvar]]))) {
-      stop("Proposed reference level ", tref, " not found in ", tvar, ".\n")
+    stop("Proposed reference level ", tref, " not found in ", tvar, ".\n")
   }
   if (length(tref) > 1) {
-      tref = min(tref, na.rm = TRUE) ## placeholder. could do something a bit smarter here like bin post periods.
-      ## also: what about NA vals?
+    tref = min(tref, na.rm = TRUE) ## placeholder. could do something a bit smarter here like bin post periods.
+    ## also: what about NA vals?
   }
   ref_string = paste0(ref_string, ", ref2 = ", tref)
 
@@ -324,9 +324,9 @@ etwfe = function(
           ctrls,
           paste(paste0("i(", gvar, ", ", ictrls, ", ref = ", gref, ")"), collapse = " + "),
           paste(paste0("i(", tvar, ", ", ictrls, ", ref = ", tref, ")"), collapse = " + ")
-          ),
+        ),
         collapse = " + "
-        )
+      )
       rhs = paste(rhs, "+", ictrls) 
     }
   }
@@ -344,14 +344,14 @@ etwfe = function(
       xvar_fml_vars = paste0("(",paste(names(xvar_dm_df), collapse = "+"), ")")
     }
     data = cbind(data, xvar_dm_df)
-
+    
     if (is.null(ctrls)) {
       rhs = paste0(
         rhs, 
         " / ", xvar_fml_vars,
         " +  i(", tvar, ", ", xvar_fml_vars, ", ref = ", tref, ")"
       )
-    # splice together with ctrl vars if necessary
+      # splice together with ctrl vars if necessary
     } else {
       rhs = paste0(
         gsub(
@@ -362,7 +362,7 @@ etwfe = function(
         " +  i(", tvar, ", ", xvar_fml_vars, ", ref = ", tref, ")"
       )
     }
-
+    
   }
 
   ## Fixed effects ----
@@ -378,7 +378,7 @@ etwfe = function(
     rhs = paste0(
       rhs, 
       "+ i(", gvar, ", ref = ", gref, ") + i(", tvar, ", ref = ", tref, ")"
-      )
+    )
   }
 
   ## Estimation ----
@@ -408,16 +408,16 @@ etwfe = function(
 
   ## Overload class and new attributes (for post-estimation) ----
   class(est) = c("etwfe", class(est))
-    attr(est, "etwfe") = list(
-      gvar = gvar,
-      tvar = tvar,
-      xvar = xvar,
-      ivar = ivar,
-      gref = gref,
-      tref = tref
-      )
-
+  attr(est, "etwfe") = list(
+    gvar = gvar,
+    tvar = tvar,
+    xvar = xvar,
+    ivar = ivar,
+    gref = gref,
+    tref = tref
+  )
+  
   ## Return ----
   return(est)
-
+  
 }