I am trying to run a staggered difference-in-difference model in R using the did2s package (Kyle Butts 2023). My base model specification is the following, set up with fake mtcars data and using the fixest package, which is the base for the did2s approach.

library(tidyverse)
library(fixest)
set.seed(123)

mtcars$time_til <- sample(-10:10, nrow(mtcars), replace = TRUE)

mtcars = mtcars |> 
  mutate(treat = case_when(time_til >= 0 ~ TRUE,
                           time_til < 0 ~ FALSE)) |> 
  mutate(am_factor = factor(am))

feols(fml = mpg ~ am_factor + am_factor:(time_til) + disp | vs, vcov = ~cyl,data = mtcars)

)

Where time_til is the time until treatment variable, and treat is the treatment status. am_factor is a factor variable interacted with treatment, as I want to calculate the treatment effect by groups, to assess variation in the treatment effect. I therefore adapt this model to the did2s model like so:

library(did2s)

es = did2s(
  data = mtcars, 
  yname = "mpg", 
  treatment = "treat",
  first_stage = ~ 0 + disp + am | vs, 
    second_stage = ~ am_factor + am_factor:i(time_til, ref = c(-1, Inf)) + disp | vs,
    cluster_var = "cyl", 
  verbose = FALSE

However, this draws the following error:

Error: in vcov.fixest(object, vcov = vcov, ssc = ...:
 Value `vcov` must be a square matrix with, in this context, 22 rows. Problem: it has 3 rows instead of 22.

with the traceback:

Error: in vcov.fixest(object, vcov = vcov, ssc = ...:
 Value `vcov` must be a square matrix with, in this context, 22 rows. Problem: it has 3 rows instead of 22.
10.
stop("in ", my_call, ":n", fit_screen(full_msg), call. = FALSE)
9.
send_error(all_reasons, x_name = x_names[qui], type, .message,
.choices, .up, .value = .value, .data = .data)
8.
check_arg_core(.x = .x, .type = .type, .message = .message, .arg_name = .arg_name,
.prefix = .prefix, .choices = .choices, .data = .data, .value = .value,
.env = .env, .up = .up, .mc = mc, .is_plus = FALSE, .is_value = TRUE)
7.
check_value(vcov, "square matrix nrow(value)", .value = n_coef)
6.
vcov.fixest(object, vcov = vcov, ssc = ssc, forceCovariance = forceCovariance,
vcov_fix = vcov_fix, keepBounded = keepBounded, nthreads = nthreads,
attr = TRUE, se = se, cluster = cluster, ...)
5.
summary.fixest(est$second_stage, .vcov = cov)
4.
summary(est$second_stage, .vcov = cov)
3.
withCallingHandlers(expr, warning = function(w) if (inherits(w,
classes)) tryInvokeRestart("muffleWarning"))
2.
base::suppressWarnings(summary(est$second_stage, .vcov = cov))
1.
did2s(data = mtcars, yname = "mpg", treatment = "treat", first_stage = ~0 +
disp + am | vs, second_stage = ~am_factor + am_factor:i(time_til,
ref = c(-1, Inf)) + disp | vs, cluster_var = "cyl", verbose = FALSE)

Why would this error be generated? Is it because there is no interaction term in the first stage? If so, how can I resolve it, given that there is no time_til term in the first stage? Is there any way to assess treatment effects by group with the did2s, if not this way?