Running 'fixest_tests.R' [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
  > stvec = stringmagic::string_vec_alias()
  > 
  > 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)
  > 
  > # 1 obs (after FE removal) estimation
  > base_1obs = data.frame(y = c(1, 0), fe = c(1, 2), x = c(1, 0))
  > test(fepois(y ~ x | fe, base_1obs), "err")
  > # no error
  > res = fepois(y ~ 1 | fe, base_1obs)
  > 
  > # warning when demeaning algo reaches max iterations 
  > data(trade)
  > test(feols(Euros ~ log(dist_km) | Destination + Origin + Product, 
  +            trade, fixef.iter = 1), "warn")
  > 
  > 
  > ####
  > #### ... 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))
  > 
  > ####
  > #### global variables ####
  > ####
  > 
  > chunk("globals")
  GLOBALS
  
  > 
  > est_reg = function(df, yvar, xvar, refgrp) {
  +   feols(.[yvar] ~ i(.[xvar], ref = refgrp), data = df)
  + }
  > 
  > (est = est_reg(iris, "Sepal.Length", "Species", ref = "setosa"))
  OLS estimation, Dep. Var.: Sepal.Length
  Observations: 150
  Standard-errors: IID 
                      Estimate Std. Error  t value   Pr(>|t|)    
  (Intercept)            5.006   0.072802 68.76164  < 2.2e-16 ***
  Species::versicolor    0.930   0.102958  9.03282 8.7702e-16 ***
  Species::virginica     1.582   0.102958 15.36551  < 2.2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  RMSE: 0.509616   Adj. R2: 0.613518
  > 
  > # checking when it should not work
  > base = setNames(iris, c("y", "x1", "x2", "x3", "species"))
  > 
  > z = base$x1
  > test(feols(y ~ z, base), "err")
  > 
  > 
  > 
  > ####
  > #### ... Collinearity ####
  > ####
  > 
  > chunk("COLLINEARITY")
  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)]
  > 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)]
  > 
  > 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("|")
  + }
  .....|.....|.....|.....|.....|.....|> cat("\n")
  
  > 
  > 
  > ####
  > #### ... split ####
  > ####
  > 
  > chunk("split")
  SPLIT
  
  > 
  > base = setNames(iris, c("y", "x1", "x2", "x3", "species"))
  > 
  > # simple: formula
  > est = feols(y ~ x.[1:3], base, split = ~species %keep% "@^v")
  > test(length(est), 2)
  > 
  > est = feols(y ~ x.[1:3], base, fsplit = ~species %keep% c("set", "vers"))
  > test(length(est), 3)
  > 
  > est = feols(y ~ x.[1:3], base, split = ~species %drop% "set")
  > test(length(est), 2)
  > 
  > # simple: vector
  > est = feols(y ~ x.[1:3], base, split = base$species %keep% "@^v")
  > test(length(est), 2)
  > 
  > est = feols(y ~ x.[1:3], base, split = base$species %keep% c("set", "vers"))
  > test(length(est), 2)
  > 
  > est = feols(y ~ x.[1:3], base, split = base$species %drop% "set")
  > test(length(est), 2)
  > 
  > # with bin
  > est = feols(y ~ x.[1:2], base,
  +       split = ~bin(x3, c("cut::5", "saint emilion", "pessac leognan",
  +                  "margaux", "saint julien", "entre deux mers")) %keep% c("saint e", "pe"))
  > test(length(est), 2)
  > 
  > est = feols(y ~ x.[1:2], base,
  +       split = ~bin(x3, c("cut::5", "saint emilion", "pessac leognan", NA)) %drop% "@\\d")
  > test(length(est), 2)
  > 
  > # with argument
  > est = feols(y ~ x.[1:3], base, split = ~species, split.keep = "@^v")
  > test(length(est), 2)
  > 
  > est = feols(y ~ x.[1:3], base, fsplit = ~species, split.keep = c("set", "vers"))
  > test(length(est), 3)
  > 
  > est = feols(y ~ x.[1:3], base, split = ~species, split.drop = "set")
  > test(length(est), 2)
  > 
  > 
  > ####
  > #### ... Multiple estimations ####
  > ####
  > 
  > chunk("Multiple")
  MULTIPLE
  
  > 
  > set.seed(2)
  > base = iris
  > names(base) = c("y1", "x1", "x2", "x3", "species")
  > base$y2 = 10 + rnorm(150) + 0.5 * base$x1
  > base$x4 = rnorm(150) + 0.5 * base$y1
  > base$fe2 = rep(letters[1:15], 10)
  > base$fe2[50:51] = NA
  > base$y2[base$fe2 == "a" & !is.na(base$fe2)] = 0
  > base$x2[1:5] = NA
  > base$x3[6] = NA
  > base$x5 = rnorm(150)
  > base$x6 = rnorm(150) + base$y1 * 0.25
  > base$fe3 = rep(letters[1:10], 15)
  > 
  > 
  > for(id_fun in 1:5){
  +   estfun = switch(as.character(id_fun),
  +                   "1" = feols,
  +                   "2" = feglm,
  +                   "3" = fepois,
  +                   "4" = femlm,
  +                   "5" = feNmlm)
  + 
  +   # Following weird bug ASAN on CRAN I cannot replicate, check 4/5 not performed on non Windows
  +   if(Sys.info()["sysname"] != "Windows"){
  +     if(id_fun %in% 4:5) next
  +   }
  + 
  + 
  +   est_multi = estfun(c(y1, y2) ~ x1 + sw(x2, x3), base, split = ~species)
  + 
  +   k = 1
  +   for(s in c("setosa", "versicolor", "virginica")){
  +     for(lhs in c("y1", "y2")){
  +       for(rhs in c("x2", "x3")){
  +         res = estfun(.[lhs] ~ x1 + .[rhs], base[base$species == s, ], notes = FALSE)
  + 
  +         test(coef(est_multi[[k]]), coef(res))
  +         test(se(est_multi[[k]], cluster = "fe3"), se(res, cluster = "fe3"))
  +         k = k + 1
  +       }
  +     }
  +   }
  + 
  +   cat("__")
  + 
  +   est_multi = estfun(c(y1, y2) ~ x1 + csw0(x2, x3) + x4 | species + fe2, base, fsplit = ~species)
  +   k = 1
  +   all_rhs = c("", "x2", "x3")
  +   for(s in c("all", "setosa", "versicolor", "virginica")){
  +     for(lhs in c("y1", "y2")){
  +       for(n_rhs in 1:3){
  +         if(s == "all"){
  +           res = estfun(xpd(..lhs ~ x1 + ..rhs + x4 | species + fe2, ..lhs = lhs, ..rhs = all_rhs[1:n_rhs]), base, notes = FALSE)
  +         } else {
  +           res = estfun(xpd(..lhs ~ x1 + ..rhs + x4 | species + fe2, ..lhs = lhs, ..rhs = all_rhs[1:n_rhs]), base[base$species == s, ], notes = FALSE)
  +         }
  + 
  +         vname = names(coef(res))
  +         test(coef(est_multi[[k]])[vname], coef(res), "~" , 1e-6)
  +         test(se(est_multi[[k]], cluster = "fe3")[vname], se(res, cluster = "fe3"), "~" , 1e-6)
  +         k = k + 1
  +       }
  +     }
  +   }
  + 
  +   cat("|")
  + }