Running ‘fixest_tests.R’ [11s/17s]
Running the tests in ‘tests/fixest_tests.R’ failed.
Complete output:
  > #----------------------------------------------#
  > # Author: Laurent Berge
  > # Date creation: Fri Jul 10 09:03:06 2020
  > # ~: package sniff tests
  > #----------------------------------------------#
  > 
  > # Not everything is currently covered, but I'll improve it over time
  > 
  > # Some functions are not trivial to test properly though
  > 
  > library(fixest)
  > 
  > test = fixest:::test ; chunk = fixest:::chunk
  > vcovClust = fixest:::vcovClust
  > 
  > setFixest_notes(FALSE)
  > 
  > if(fixest:::is_r_check()){
  +     if(requireNamespace("data.table", quietly = TRUE)){
  +         library(data.table)
  +         data.table::setDTthreads(1)
  +     }
  +     setFixest_nthreads(1)
  + }
  > 
  > ####
  > #### ESTIMATIONS ####
  > ####
  > 
  > ####
  > #### ... Main ####
  > ####
  > 
  > 
  > chunk("ESTIMATION")
  ESTIMATION 
  
  > 
  > set.seed(0)
  > 
  > base = iris
  > names(base) = c("y", "x1", "x2", "x3", "species")
  > base$fe_2 = rep(1:5, 30)
  > base$fe_3 = sample(15, 150, TRUE)
  > base$constant = 5
  > base$y_int = as.integer(base$y)
  > base$w = as.vector(unclass(base$species) - 0.95)
  > base$offset_value = unclass(base$species) - 0.95
  > base$y_01 = 1 * ((scale(base$x1) + rnorm(150)) > 0)
  > # what follows to avoid removal of fixed-effects (logit is pain in the neck)
  > base$y_01[1:5 + rep(c(0, 50, 100), each = 5)] = 1
  > base$y_01[6:10 + rep(c(0, 50, 100), each = 5)] = 0
  > # We enforce the removal of observations
  > base$y_int_null = base$y_int
  > base$y_int_null[base$fe_3 %in% 1:5] = 0
  > 
  > for(model in c("ols", "pois", "logit", "negbin", "Gamma")){
  +     cat("Model: ", format(model, width = 6), sep = "")
  +     for(use_weights in c(FALSE, TRUE)){
  +         my_weight = NULL
  +         if(use_weights) my_weight = base$w
  + 
  +         for(use_offset in c(FALSE, TRUE)){
  +             my_offset = NULL
  +             if(use_offset) my_offset = base$offset_value
  + 
  +             for(id_fe in 0:9){
  + 
  +                 cat(".")
  + 
  +                 tol = switch(model, "negbin" = 1e-2, "logit" = 3e-5, 1e-5)
  + 
  +                 # Setting up the formula to accommodate FEs
  +                 if(id_fe == 0){
  +                     fml_fixest = fml_stats = y ~ x1
  +                 } else if(id_fe == 1){
  +                     fml_fixest = y ~ x1 | species
  +                     fml_stats = y ~ x1 + factor(species)
  +                 } else if(id_fe == 2){
  +                     fml_fixest = y ~ x1 | species + fe_2
  +                     fml_stats = y ~ x1 + factor(species) + factor(fe_2)
  +                 } else if(id_fe == 3){
  +                     # varying slope
  +                     fml_fixest = y ~ x1 | species[[x2]]
  +                     fml_stats = y ~ x1 + x2:species
  +                 } else if(id_fe == 4){
  +                     # varying slope -- 1 VS, 1 FE
  +                     fml_fixest = y ~ x1 | species[[x2]] + fe_2
  +                     fml_stats = y ~ x1 + x2:species + factor(fe_2)
  +                 } else if(id_fe == 5){
  +                     # varying slope -- 2 VS
  +                     fml_fixest = y ~ x1 | species[x2]
  +                     fml_stats = y ~ x1 + x2:species + species
  +                 } else if(id_fe == 6){
  +                     # varying slope -- 2 VS bis
  +                     fml_fixest = y ~ x1 | species[[x2]] + fe_2[[x3]]
  +                     fml_stats = y ~ x1 + x2:species + x3:factor(fe_2)
  +                 } else if(id_fe == 7){
  +                     # Combined clusters
  +                     fml_fixest = y ~ x1 + x2 | species^fe_2
  +                     fml_stats = y ~ x1 + x2 + paste(species, fe_2)
  +                 } else if(id_fe == 8){
  +                     fml_fixest = y ~ x1 | species[x2] + fe_2[x3] + fe_3
  +                     fml_stats = y ~ x1 + species + i(species, x2) + factor(fe_2) + i(fe_2, x3) + factor(fe_3)
  +                 } else if(id_fe == 9){
  +                     fml_fixest = y ~ x1 | species + fe_2[x2,x3] + fe_3
  +                     fml_stats = y ~ x1 + species + factor(fe_2) + i(fe_2, x2) + i(fe_2, x3) + factor(fe_3)
  +                 }
  + 
  +                 # ad hoc modifications of the formula
  +                 if(model == "logit"){
  +                     fml_fixest = xpd(y_01 ~ ..rhs, ..rhs = fml_fixest[[3]])
  +                     fml_stats = xpd(y_01 ~ ..rhs, ..rhs = fml_stats[[3]])
  + 
  +                     # The estimations are OK, conv differences out of my control
  +                     if(id_fe %in% 8:9) tol = 0.5
  + 
  +                 } else if(model == "pois"){
  +                     fml_fixest = xpd(y_int_null ~ ..rhs, ..rhs = fml_fixest[[3]])
  +                     fml_stats = xpd(y_int_null ~ ..rhs, ..rhs = fml_stats[[3]])
  + 
  +                 } else if(model %in% c("negbin", "Gamma")){
  +                     fml_fixest = xpd(y_int ~ ..rhs, ..rhs = fml_fixest[[3]])
  +                     fml_stats = xpd(y_int ~ ..rhs, ..rhs = fml_stats[[3]])
  +                 }
  + 
  +                 adj = 1
  +                 if(model == "ols"){
  +                     res = feols(fml_fixest, base, weights = my_weight, offset = my_offset)
  +                     res_bis = lm(fml_stats, base, weights = my_weight, offset = my_offset)
  + 
  +                 } else if(model %in% c("pois", "logit", "Gamma")){
  +                     adj = 0
  +                     if(model == "Gamma" && use_offset) next
  + 
  +                     my_family = switch(model, pois = poisson(), logit = binomial(), Gamma = Gamma())
  + 
  +                     res = feglm(fml_fixest, base, family = my_family, weights = my_weight, offset = my_offset)
  + 
  +                     if(!is.null(res$obs_selection$obsRemoved)){
  +                         qui = res$obs_selection$obsRemoved
  + 
  +                         # I MUST do that.... => subset does not work...
  +                         base_tmp = base[qui, ]
  +                         base_tmp$my_offset = my_offset[qui]
  +                         base_tmp$my_weight = my_weight[qui]
  +                         res_bis = glm(fml_stats, base_tmp, family = my_family, weights = my_weight, offset = my_offset)
  +                     } else {
  +                         res_bis = glm(fml_stats, data = base, family = my_family, weights = my_weight, offset = my_offset)
  +                     }
  + 
  +                 } else if(model == "negbin"){
  +                     # no offset in glm.nb + no VS in fenegbin + no weights in fenegbin
  +                     if(use_weights || use_offset || id_fe > 2) next
  + 
  +                     res = fenegbin(fml_fixest, base, notes = FALSE)
  +                     res_bis = MASS::glm.nb(fml_stats, base)
  + 
  +                 }
  + 
  +                 test(coef(res)["x1"], coef(res_bis)["x1"], "~", tol)
  +                 test(se(res, se = "st", ssc = ssc(adj = adj))["x1"], se(res_bis)["x1"], "~", tol)
  +                 test(pvalue(res, se = "st", ssc = ssc(adj = adj))["x1"], pvalue(res_bis)["x1"], "~", tol*10**(model == "negbin"))
  +                 # cat("Model: ", model, ", FE: ", id_fe, ", weight: ", use_weights,  ", offset: ", use_offset, "\n", sep="")
  + 
  +             }
  +             cat("|")
  +         }
  +     }
  +     cat("\n")
  + }
  Model: ols   ..........|..........|..........|..........|
  Model: pois  ..........|..........|..........|..........|
  Model: logit ..........|..........|..........|..........|
  Model: negbin..........|..........|..........|..........|
  Model: Gamma ..........|..........|..........|..........|
  There were 36 warnings (use warnings() to see them)
  > 
  > ####
  > #### ... Corner cases ####
  > ####
  > 
  > chunk("Corner cases")
  CORNER CASES 
  
  > 
  > 
  > # We test the absence of bugs
  > 
  > base = iris
  > names(base) = c("y", "x1", "x2", "x3", "fe1")
  > base$fe2 = rep(1:5, 30)
  > base$y[1:5] = NA
  > base$x1[4:8] = NA
  > base$x2[4:21] = NA
  > base$x3[110:111] = NA
  > base$fe1[110:118] = NA
  > base$fe2[base$fe2 == 1] = 0
  > base$fe3 = sample(letters[1:5], 150, TRUE)
  > base$period = rep(1:50, 3)
  > base$x_cst = 1
  > 
  > res = feols(y ~ 1 | csw(fe1, fe1^fe2), base)
  > 
  > res = feols(y ~ 1 + csw(x1, i(fe1)) | fe2, base)
  > 
  > res = feols(y ~ csw(f(x1, 1:2), x2) | sw0(fe2, fe2^fe3), base, panel.id = ~ fe1 + period)
  > 
  > res = feols(d(y) ~ -1 + d(x2), base, panel.id = ~ fe1 + period)
  > test(length(coef(res)), 1)
  > 
  > res = feols(c(y, x1) ~ 1 | fe1 | x2 ~ x3, base)
  > 
  > res = feols(y ~ x1 | fe1[x2] + fe2[x2], base)
  > 
  > #
  > # NA models (ie all variables are collinear with the FEs)
  > #
  > 
  > # Should work when warn = FALSE or multiple est
  > for(i in 1:2){
  +     fun = switch(i, "1" = feols, "2" = feglm)
  + 
  +     res = feols(y ~ x_cst | fe1, base, warn = FALSE)
  +     res         # => no error
  +     etable(res) # => no error
  + 
  +     # error when warn = TRUE
  +     test(feols(y ~ x_cst | fe1, base), "err")
  + 
  +     # multiple est => no error
  +     res = feols(c(y, x1) ~ x_cst | fe1, base)
  +     res         # => no error
  +     etable(res) # => no error
  + }
  > 
  > 
  > # Removing the intercept!!!
  > 
  > res = feols(y ~ -1 + x1 + i(fe1), base)
  > test("(Intercept)" %in% names(res$coefficients), FALSE)
  > 
  > res = feols(y ~ -1 + x1 + factor(fe1), base)
  > test("(Intercept)" %in% names(res$coefficients), FALSE)
  > 
  > res = feols(y ~ -1 + x1 + i(fe1) + i(fe2), base)
  > test("(Intercept)" %in% names(res$coefficients), FALSE)
  > test(is.null(res$collin.var), TRUE)
  > 
  > 
  > # IV + interacted FEs
  > res = feols(y ~ x1 | fe1^fe2 | x2 ~ x3, base)
  > 
  > # IVs no exo var
  > res = feols(y ~ 0 | x2 ~ x3, base)
  > # Same in stepwise
  > res = feols(y ~ 0 | sw0(fe1) | x2 ~ x3, base)
  > 
  > # IVs + lags
  > res = feols(y ~ x1 | fe1^fe2 | l(x2, -1:1) ~ l(x3, -1:1), base, panel.id = ~ fe1 + period)
  > 
  > # functions in interactions
  > res = feols(y ~ x1 | factor(fe1)^factor(fe2), base)
  > res = feols(y ~ x1 | round(x2^2), base)
  > test(feols(y ~ x1 | factor(fe1^fe2), base), "err")
  > 
  > res = feols(y ~ x1 | bin(x2, "bin::1")^fe1 + fe1^fe2, base)
  > 
  > 
  > 
  > 
  > 
  > ####
  > #### ... Fit methods ####
  > ####
  > 
  > chunk("Fit methods")
  FIT METHODS 
  
  > 
  > base = iris
  > names(base) = c("y", "x1", "x2", "x3", "species")
  > base$y_int = as.integer(base$y)
  > base$y_log = sample(c(TRUE, FALSE), 150, TRUE)
  > 
  > res = feglm.fit(base$y, base[, 2:4])
  > res_bis = feglm(y ~ -1 + x1 + x2 + x3, base)
  > test(coef(res), coef(res_bis))
  > 
  > res = feglm.fit(base$y_int, base[, 2:4])
  > res_bis = feglm(y_int ~ -1 + x1 + x2 + x3, base)
  > test(coef(res), coef(res_bis))
  > 
  > res = feglm.fit(base$y_log, base[, 2:4])
  > res_bis = feglm(y_log ~ -1 + x1 + x2 + x3, base)
  > test(coef(res), coef(res_bis))
  > 
  > 
  > 
  > res = feglm.fit(base$y, base[, 2:4], family = "poisson")
  > res_bis = feglm(y ~ -1 + x1 + x2 + x3, base, family = "poisson")
  > test(coef(res), coef(res_bis))
  > 
  > res = feglm.fit(base$y_int, base[, 2:4], family = "poisson")
  > res_bis = feglm(y_int ~ -1 + x1 + x2 + x3, base, family = "poisson")
  > test(coef(res), coef(res_bis))
  > 
  > res = feglm.fit(base$y_log, base[, 2:4], family = "poisson")
  > res_bis = feglm(y_log ~ -1 + x1 + x2 + x3, base, family = "poisson")
  > test(coef(res), coef(res_bis))
  > 
  > 
  > 
  > ####
  > #### ... Collinearity ####
  > ####
  > 
  > cat("COLLINEARITY\n\n")
  COLLINEARITY
  
  > 
  > base = iris
  > names(base) = c("y", "x1", "x2", "x3", "species")
  > base$constant = 5
  > base$y_int = as.integer(base$y)
  > base$w = as.vector(unclass(base$species) - 0.95)
  > 
  > for(useWeights in c(FALSE, TRUE)){
  +     for(model in c("ols", "pois")){
  +         for(use_fe in c(FALSE, TRUE)){
  +             cat(".")
  + 
  +             my_weight = NULL
  +             if(useWeights) my_weight = base$w
  + 
  +             adj = 1
  +             if(model == "ols"){
  +                 if(!use_fe){
  +                     res = feols(y ~ x1 + constant, base, weights = my_weight)
  +                     res_bis = lm(y ~ x1 + constant, base, weights = my_weight)
  +                 } else {
  +                     res = feols(y ~ x1 + constant | species, base, weights = my_weight)
  +                     res_bis = lm(y ~ x1 + constant + species, base, weights = my_weight)
  +                 }
  +             } else {
  +                 if(!use_fe){
  +                     res = fepois(y_int ~ x1 + constant, base, weights = my_weight)
  +                     res_bis = glm(y_int ~ x1 + constant, base, weights = my_weight, family = poisson)
  +                 } else {
  +                     res = fepois(y_int ~ x1 + constant | species, base, weights = my_weight)
  +                     res_bis = glm(y_int ~ x1 + constant + species, base, weights = my_weight, family = poisson)
  +                 }
  +                 adj = 0
  +             }
  + 
  +             test(coef(res)["x1"], coef(res_bis)["x1"], "~")
  +             test(se(res, se = "st", ssc = ssc(adj=adj))["x1"], se(res_bis)["x1"], "~")
  +             # cat("Weight: ", useWeights, ", model: ", model, ", FE: ", use_fe, "\n", sep="")
  + 
  +         }
  +     }
  + }
  ........> cat("\n")
  
  > 
  > 
  > ####
  > #### ... Non linear tests ####
  > ####
  > 
  > chunk("NON LINEAR")
  NON LINEAR 
  
  > 
  > base = iris
  > names(base) = c("y", "x1", "x2", "x3", "species")
  > 
  > tab = c("versicolor" = 5, "setosa" = 0, "virginica" = -5)
  > 
  > fun_nl = function(a, b, spec){
  +     res = as.numeric(tab[spec])
  +     a*res + b*res^2
  + }
  > 
  > est_nl = feNmlm(y ~ x1, base, NL.fml = ~fun_nl(a, b, species), NL.start = 1, family = "gaussian")
  > 
  > base$var_spec = as.numeric(tab[base$species])
  > 
  > est_lin = feols(y ~ x1 + var_spec + I(var_spec^2), base)
  > 
  > test(coef(est_nl), coef(est_lin)[c(3, 4, 1, 2)], "~")
  > 
  > ####
  > #### ... Lagging ####
  > ####
  > 
  > # Different types of lag
  > # 1) check no error in wide variety of situations
  > # 2) check consistency
  > 
  > chunk("LAGGING")
  LAGGING 
  
  > 
  > data(base_did)
  > base = base_did
  > 
  > n = nrow(base)
  > 
  > set.seed(0)
  > base$y_na = base$y ; base$y_na[sample(n, 50)] = NA
  > base$period_txt = letters[base$period]
  > ten_dates = c("1960-01-15", "1960-01-16", "1960-03-31", "1960-04-05", "1960-05-12", "1960-05-25", "1960-06-20", "1960-07-30", "1965-01-02", "2002-12-05")
  > base$period_date = as.Date(ten_dates, "%Y-%m-%d")[base$period]
  > base$y_0 = base$y**2 ; base$y_0[base$id == 1] = 0
  > 
  > # We compute the lags "by hand"
  > base = base[order(base$id, base$period), ]
  > base$x1_lag = c(NA, base$x1[-n]) ; base$x1_lag[base$period == 1] = NA
  > base$x1_lead = c(base$x1[-1], NA) ; base$x1_lead[base$period == 10] = NA
  > base$x1_diff = base$x1 - base$x1_lag
  > 
  > # we create holes
  > base$period_bis = base$period ; base$period_bis[base$period_bis == 5] = 50
  > base$x1_lag_hole = base$x1_lag ; base$x1_lag_hole[base$period %in% c(5, 6)] = NA
  > base$x1_lead_hole = base$x1_lead ; base$x1_lead_hole[base$period %in% c(4, 5)] = NA
  > 
  > # we reshuffle the base
  > base = base[sample(n), ]
  > 
  > #
  > # Checks consistency
  > #
  > 
  > cat("consistentcy...")
  consistentcy...> 
  > test(lag(x1 ~ id + period, data = base), base$x1_lag)
  > test(lag(x1 ~ id + period, -1, data = base), base$x1_lead)
  > 
  > test(lag(x1 ~ id + period_bis, data = base), base$x1_lag_hole)
  > test(lag(x1 ~ id + period_bis, -1, data = base), base$x1_lead_hole)
  > 
  > test(lag(x1 ~ id + period_txt, data = base), base$x1_lag)
  > test(lag(x1 ~ id + period_txt, -1, data = base), base$x1_lead)
  > 
  > test(lag(x1 ~ id + period_date, data = base), base$x1_lag)
  > test(lag(x1 ~ id + period_date, -1, data = base), base$x1_lead)
  > 
  > cat("done.\nEstimations...")
  done.
  Estimations...> 
  > #
  > # Estimations
  > #
  > 
  > # Poisson
  > 
  > for(depvar in c("y", "y_na", "y_0")){
  +     for(p in c("period", "period_txt", "period_date")){
  + 
  +         base$per = base[[p]]
  + 
  +         cat(".")
  + 
  +         base$y_dep = base[[depvar]]
  +         pdat = panel(base, ~ id + period)
  + 
  +         if(depvar == "y_0"){
  +             estfun = fepois
  +         } else {
  +             estfun = feols
  +         }
  + 
  +         est_raw  = estfun(y_dep ~ x1 + x1_lag + x1_lead, base)
  +         est      = estfun(y_dep ~ x1 + l(x1) + f(x1), base, panel.id = "id,per")
  +         est_pdat = estfun(y_dep ~ x1 + l(x1, 1) + f(x1, 1), pdat)
  +         test(coef(est_raw), coef(est))
  +         test(coef(est_raw), coef(est_pdat))
  + 
  +         # Now diff
  +         est_raw  = estfun(y_dep ~ x1 + x1_diff, base)
  +         est      = estfun(y_dep ~ x1 + d(x1), base, panel.id = "id,per")
  +         est_pdat = estfun(y_dep ~ x1 + d(x1, 1), pdat)
  +         test(coef(est_raw), coef(est))
  +         test(coef(est_raw), coef(est_pdat))
  + 
  +         # Now we just check that calls to l/f works without checking coefs
  + 
  +         est = estfun(y_dep ~ x1 + l(x1) + f(x1), base, panel.id = "id,per")
  +         est = estfun(y_dep ~ l(x1, -1:1) + f(x1, 2), base, panel.id = c("id", "per"))
  +         est = estfun(y_dep ~ l(x1, -1:1, fill = 1), base, panel.id = ~ id + per)
  +         if(depvar == "y") test(est$nobs, n)
  +         est = estfun(f(y_dep) ~ f(x1, -1:1), base, panel.id = ~ id + per)
  +     }
  + }
  .........> 
  > cat("done.\n\n")
  done.
  
  > 
  > #
  > # Data table
  > #
  > 
  > cat("data.table...")
  data.table...> # We just check there is no bug (consistency should be OK)
  > 
  > library(data.table)
  > 
  > base_dt = data.table(id = c("A", "A", "B", "B"),
  +                      time = c(1, 2, 1, 3),
  +                      x = c(5, 6, 7, 8))
  > 
  > base_dt = panel(base_dt, ~id + time)
  > 
  > base_dt[, x_l := l(x)]
  [1] TRUE
  > test(base_dt$x_l, c(NA, 5, NA, NA))
  > 
  > lag_creator = function(dt) {
  +     dt2 = panel(dt, ~id + time)
  +     dt2[, x_l := l(x)]
  +     return(dt2)
  + }
  > 
  > base_bis = lag_creator(base_dt)
  > 
  > base_bis[, x_d := d(x)]
  [1] TRUE
  > 
  > cat("done.\n\n")
  done.
  
  > 
  > #
  > # Panel
  > #
  > 
  > # We ensure we get the right SEs whether we use the panel() or the panel.id method
  > data(base_did)
  > 
  > # Setting a data set as a panel...
  > pdat = panel(base_did, ~id+period)
  > pdat$fe = sample(15, nrow(pdat), replace = TRUE)
  > 
  > base_panel = unpanel(pdat)
  > 
  > est_pdat = feols(y ~ x1 | fe, pdat)
  > est_panel = feols(y ~ x1 | fe, base_panel, panel.id = ~id+period)
  > 
  > test(attr(vcov(est_pdat, attr = TRUE), "type"),
  +      attr(vcov(est_panel, attr = TRUE), "type"))
  > 
  > ####
  > #### ... subset ####
  > ####
  > 
  > chunk("SUBSET")
  SUBSET 
  
  > 
  > set.seed(5)
  > base = iris
  > names(base) = c("y", "x1", "x2", "x3", "species")
  > base$fe_bis = sample(letters, 150, TRUE)
  > base$x4 = rnorm(150)
  > base$x1[sample(150, 5)] = NA
  > 
  > fml = y ~ x1 + x2
  > 
  > # Errors
  > test(feols(fml, base, subset = ~species), "err")
  > test(feols(fml, base, subset = -1:15), "err")
  > test(feols(fml, base, subset = integer(0)), "err")
  > test(feols(fml, base, subset = c(TRUE, TRUE, FALSE)), "err")
  > 
  > # Valid use
  > for(id_fun in 1:6){
  +     estfun = switch(as.character(id_fun),
  +                     "1" = feols,
  +                     "2" = feglm,
  +                     "3" = fepois,
  +                     "4" = femlm,
  +                     "5" = fenegbin,
  +                     "6" = feNmlm)
  + 
  +     for(id_fe in 1:5){
  + 
  +         cat(".")
  + 
  +         fml = switch(as.character(id_fe),
  +                      "1" = y ~ x1 + x2,
  +                      "2" = y ~ x1 + x2 | species,
  +                      "3" = y ~ x1 + x2 | fe_bis,
  +                      "4" = y ~ x1 + x2 + i(fe_bis),
  +                      "5" = y ~ x1 + x2 | fe_bis[x3])
  + 
  +         if(id_fe == 5 && id_fun %in% 4:6) next
  + 
  +         if(id_fun == 6){
  +             res_sub_a = estfun(fml, base, subset = ~species == "setosa", NL.fml = ~ a*x4, NL.start = 0)
  +             res_sub_b = estfun(fml, base, subset = base$species == "setosa", NL.fml = ~ a*x4, NL.start = 0)
  +             res_sub_c = estfun(fml, base, subset = which(base$species == "setosa"), NL.fml = ~ a*x4, NL.start = 0)
  +             res       = estfun(fml, base[base$species == "setosa", ], NL.fml = ~ a*x4, NL.start = 0)
  +         } else {
  +             res_sub_a = estfun(fml, base, subset = ~species == "setosa")
  +             res_sub_b = estfun(fml, base, subset = base$species == "setosa")
  +             res_sub_c = estfun(fml, base, subset = which(base$species == "setosa"))
  +             res       = estfun(fml, base[base$species == "setosa", ])
  +         }
  + 
  +         test(coef(res_sub_a), coef(res))
  +         test(coef(res_sub_b), coef(res))
  +         test(coef(res_sub_c), coef(res))
  +         test(se(res_sub_c, cluster = "fe_bis"), se(res, cluster = "fe_bis"))
  +     }
  +     cat("|")
  + }
  .....|.....|.....|..
   *** caught segfault ***
  address 0x555f1ac34aa8, cause 'memory not mapped'
  
  Traceback:
   1: cpp_quf_gnl(x)
   2: quickUnclassFactor(cluster)
   3: vcovClust(index, bread, scores, adj = ssc$cluster.adj && ssc$cluster.df ==     "conventional", sandwich = sandwich, nthreads = nthreads)
   4: vcov_cluster_internal(bread = c(0.0343333522550164, -0.0116081067476734, -0.0116081067476734, 0.138400979709878), scores = c(0.00528619655138787, -0.083038964096984, 0.0234407877065867, 0.0636098350465548, -0.0169291270411196, 0.000980504335863097, 0.00914764235448007, -0.000215745957475467, 0.127292629322581, -0.033957622026217, 0.0562447753614339, 0.00555999678779484, -0.0405164784060599, 0.137662706663015, 0.279100079490618, 0.0596250527972676, 0.00528619655138787, 0.14002074950487, -0.0618727152442642, -0.00883990460774831, -0.0261877675657947, -0.0664549403395713, -0.0146981498144392, 0.00777785780819418, -0.10227158196829, 0.000515025405946973, 0.0105810231597044, -0.00598727179106318, 0.0423900407922186, 0.00761701623776041, -0.010305688721171, -0.189023505522232, -0.0168071601778402, -0.033957622026217, -0.0503703299859697, 0.0373663463024261, -0.0344066413501957, 0.118029228757164, -0.00273823164839931, -2.12252429544547e-05, -0.269542250543915, 0.0910082447793581, -0.00656936144751319, -0.108083385205665, -0.0405164784060599, -0.0733252905051341, 0.0522431721345435, 0.0287672610523576, -0.0129149292612105, -0.00455208286956575, -0.012574403273988, 0.017109349676563, -0.00696469721251414, 0.00623704680462295, 0.000486579392713255, 0.0233509509048727, 0.00030456195941838, 0.0156056582223915, 0.00371805013008616, 0.00728897475670426, -0.0298906148219057, -0.00613531905485473, 0.117968243869575, -0.128383257277603, -0.0206451681131329, -0.00455208286956575, 0.088026571331577, -0.00587881913142302, 0.0825678610282907, -0.00339377258589604, 0.176575828612207, 0.0276544033237395, -0.134316166849775, 0.0326155982490932, -0.00276878326017211, 0.00503885793991393, -0.0152835544099267, -0.0255232828658178, -0.00317607075558158, 0.0145482250828854, -0.00997512713177012, 0.00139682746451892, 0.00371805013008616, -0.0591295923070008, -0.082148758054704, 0.0126761310237562, 0.0456297979840936, 0.00386547774028514, 4.66630409621984e-05, -0.0393790425984478, 0.0664266023339626, -0.0119139266929476, -0.125627628380752, -0.00613531905485473, -0.0265321248615174, 0.0149361038504095, 0.0037280589758375, -0.00664087939484378), vars = list(cl1 = c(19L, 21L, 11L, 18L, 21L, 11L, 8L, 22L, 24L, 6L, 17L, 1L, 12L, 10L, 18L, 23L, 7L, 10L, 2L, 8L, 16L, 20L, 22L, 10L, 10L, 17L, 5L, 10L, 16L, 4L, 13L, 18L, 6L, 21L, 16L, 3L, 14L, 7L, 11L, 20L, 20L, 9L, 4L, 15L, 9L, 9L, 13L, 21L, 15L)), ssc = list(adj = TRUE, fixef.K = "nested",     cluster.adj = TRUE, cluster.df = "min", t.df = "min", fixef.force_exact = FALSE),     sandwich = TRUE, nthreads = 1, var_names_all = c(cl1 = "fe_bis"))
   5: do.call(fun_name, args)
   6: vcov.fixest(object, vcov = vcov, ssc = ssc, forceCovariance = forceCovariance,     keepBounded = keepBounded, nthreads = nthreads, attr = TRUE,     se = se, cluster = cluster, ...)
   7: summary.fixest(object, vcov = vcov, ssc = ssc, cluster = cluster,     ...)
   8: summary(object, vcov = vcov, ssc = ssc, cluster = cluster, ...)
   9: coeftable.fixest(object, vcov = vcov, ssc = ssc, cluster = cluster,     keep = keep, drop = drop, order = order, ...)
  10: coeftable(object, vcov = vcov, ssc = ssc, cluster = cluster,     keep = keep, drop = drop, order = order, ...)
  11: se.fixest(res, cluster = "fe_bis")
  12: se(res, cluster = "fe_bis")
  13: test(se(res_sub_c, cluster = "fe_bis"), se(res, cluster = "fe_bis"))
  An irrecoverable exception occurred. R is aborting now ...
  Segmentation fault