Running 'fixest_tests.R' [26s/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)
  (Permanently remove the following message with fixest_startup_msg(FALSE).)
  fixest 0.10.0:
  - vcov: new argument 'vcov' that replaces 'se' and 'cluster' in all functions
  (retro compatibility is ensured).
  - function 'dof()' has been renamed into 'ssc()' (i.e. small sample correction).
  From fixest 0.9.0 onward: BREAKING changes! 
  - In i():
      + the first two arguments have been swapped! Now it's i(factor_var,
  continuous_var) for interactions.
      + argument 'drop' has been removed (put everything in 'ref' now).
  - In feglm(): 
      + the default family becomes 'gaussian' to be in line with glm(). Hence, for
  Poisson estimations, please use fepois() instead.
  > 
  > 
  > test = fixest:::test ; chunk = fixest:::chunk
  > vcovClust = fixest:::vcovClust
  > 
  > setFixest_notes(FALSE)
  > 
  > ####
  > #### Estimations ####
  > ####
  > 
  > 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 50 or more warnings (use warnings() to see the first 50)
  > 
  > ####
  > #### 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
  + }
  Warning messages:
  1: In feols(env = current_env, xwx = xwx[qui, qui, drop = FALSE], xwy = xwy[[ii]][qui],  :
    The only variable 'x_cst' is collinear with the fixed effects. In such circumstances, the estimation is void.
  2: In feols(env = current_env, xwx = xwx[qui, qui, drop = FALSE], xwy = xwy[[ii]][qui],  :
    The only variable 'x_cst' is collinear with the fixed effects. In such circumstances, the estimation is void.
  3: In feols(env = current_env, xwx = xwx[qui, qui, drop = FALSE], xwy = xwy[[ii]][qui],  :
    The only variable 'x_cst' is collinear with the fixed effects. In such circumstances, the estimation is void.
  4: In feols(env = current_env, xwx = xwx[qui, qui, drop = FALSE], xwy = xwy[[ii]][qui],  :
    The only variable 'x_cst' is collinear with the fixed effects. In such circumstances, the estimation is void.
  > 
  > 
  > # 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 + lags
  > res = feols(y ~ x1 | fe1^fe2 | l(x2, -1:1) ~ l(x3, -1:1), base, panel.id = ~ fe1 + period)
  > 
  > 
  > ####
  > #### 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))
  > 
  > ####
  > #### Standard-errors ####
  > ####
  > 
  > chunk("STANDARD ERRORS")
  STANDARD ERRORS 
  
  > 
  > #
  > # Fixed-effects corrections
  > #
  > 
  > # We create "irregular" FEs
  > set.seed(0)
  > base = data.frame(x = rnorm(20))
  > base$y = base$x + rnorm(20)
  > base$fe1 = rep(rep(1:3, c(4, 3, 3)), 2)
  > base$fe2 = rep(rep(1:5, each = 2), 2)
  > est = feols(y ~ x | fe1 + fe2, base)
  > 
  > # fe1: 3 FEs
  > # fe2: 5 FEs
  > 
  > #
  > # Clustered standard-errors: by fe1
  > #
  > 
  > # Default: fixef.K = "nested"
  > #  => adjustment K = 1 + 5 (i.e. x + fe2)
  > test(attr(vcov(est, ssc = ssc(fixef.K = "nested"), attr = TRUE), "dof.K"), 6)
  > 
  > # fixef.K = FALSE
  > #  => adjustment K = 1 (i.e. only x)
  > test(attr(vcov(est, ssc = ssc(fixef.K = "none"), attr = TRUE), "dof.K"), 1)
  > 
  > # fixef.K = TRUE
  > #  => adjustment K = 1 + 3 + 5 - 1 (i.e. x + fe1 + fe2 - 1 restriction)
  > test(attr(vcov(est, ssc = ssc(fixef.K = "full"), attr = TRUE), "dof.K"), 8)
  > 
  > # fixef.K = TRUE & fixef.exact = TRUE
  > #  => adjustment K = 1 + 3 + 5 - 2 (i.e. x + fe1 + fe2 - 2 restrictions)
  > test(attr(vcov(est, ssc = ssc(fixef.K = "full", fixef.force_exact = TRUE), attr = TRUE), "dof.K"), 7)
  > 
  > #
  > # Manual checks of the SEs
  > #
  > 
  > n = est$nobs
  > VCOV_raw = est$cov.iid / ((n - 1) / (n - est$nparams))
  > 
  > # standard
  > for(k_val in c("none", "nested", "full")){
  +     for(adj in c(FALSE, TRUE)){
  + 
  +         K = switch(k_val, none = 1, nested = 8, full = 8)
  +         my_adj = ifelse(adj, (n - 1) / (n - K), 1)
  + 
  +         test(vcov(est, se = "standard", ssc = ssc(adj = adj, fixef.K = k_val)), VCOV_raw * my_adj)
  + 
  +         # cat("adj = ", adj, " ; fixef.K = ", k_val, "\n", sep = "")
  +     }
  + }
  > 
  > # Clustered, fe1
  > VCOV_raw = est$cov.iid / est$sigma2
  > H = vcovClust(est$fixef_id$fe1, VCOV_raw, scores = est$scores, adj = FALSE)
  > n = nobs(est)
  > 
  > for(tdf in c("conventional", "min")){
  +     for(k_val in c("none", "nested", "full")){
  +         for(c_adj in c(FALSE, TRUE)){
  +             for(adj in c(FALSE, TRUE)){
  + 
  +                 K = switch(k_val, none = 1, nested = 6, full = 8)
  +                 cluster_factor = ifelse(c_adj, 3/2, 1)
  +                 df = ifelse(tdf == "min", 2, 20 - K)
  +                 my_adj = ifelse(adj, (n - 1) / (n - K), 1)
  + 
  +                 V = H * cluster_factor
  + 
  +                 # test SE
  +                 test(vcov(est, se = "cluster", ssc = ssc(adj = adj, fixef.K = k_val, cluster.adj = c_adj)), V * my_adj)
  + 
  +                 # test pvalue
  +                 my_tstat = tstat(est, se = "cluster", ssc = ssc(adj = adj, fixef.K = k_val, cluster.adj = c_adj))
  +                 test(pvalue(est, se = "cluster", ssc = ssc(adj = adj, fixef.K = k_val, cluster.adj = c_adj, t.df = tdf)), 2*pt(-abs(my_tstat), df))
  + 
  +                 # cat("adj = ", adj, " ; fixef.K = ", k_val, " ; cluster.adj = ", c_adj, " t.df = ", tdf, "\n", sep = "")
  +             }
  +         }
  +     }
  + }
  > 
  > 
  > # 2-way Clustered, fe1 fe2
  > VCOV_raw = est$cov.iid / est$sigma2
  > M_i  = vcovClust(est$fixef_id$fe1, VCOV_raw, scores = est$scores, adj = FALSE)
  > M_t  = vcovClust(est$fixef_id$fe2, VCOV_raw, scores = est$scores, adj = FALSE)
  > M_it = vcovClust(paste(base$fe1, base$fe2), VCOV_raw, scores = est$scores, adj = FALSE, do.unclass = TRUE)
  > 
  > M_i + M_t - M_it
              [,1]
  [1,] 0.005391594
  > vcov(est, se = "two", ssc = ssc(adj = FALSE, cluster.adj = FALSE))
              x
  x 0.005391594
  > 
  > for(cdf in c("conventional", "min")){
  +     for(tdf in c("conventional", "min")){
  +         for(k_val in c("none", "nested", "full")){
  +             for(c_adj in c(FALSE, TRUE)){
  +                 for(adj in c(FALSE, TRUE)){
  + 
  +                     K = switch(k_val, none = 1, nested = 2, full = 8)
  + 
  +                     if(c_adj){
  +                         if(cdf == "min"){
  +                             V = (M_i + M_t - M_it) * 3/2
  +                         } else {
  +                             V = M_i  * 3/2 + M_t * 5/4 - M_it * 6/5
  +                         }
  +                     } else {
  +                         V = M_i + M_t - M_it
  +                     }
  + 
  +                     df = ifelse(tdf == "min", 2, 20 - K)
  +                     my_adj = ifelse(adj, (n - 1) / (n - K), 1)
  + 
  +                     # test SE
  +                     test(vcov(est, se = "two", ssc = ssc(adj = adj, fixef.K = k_val, cluster.adj = c_adj, cluster.df = cdf)),
  +                          V * my_adj)
  + 
  +                     # test pvalue
  +                     my_tstat = tstat(est, se = "two", ssc = ssc(adj = adj, fixef.K = k_val, cluster.adj = c_adj, cluster.df = cdf))
  +                     test(pvalue(est, se = "two", ssc = ssc(adj = adj, fixef.K = k_val, cluster.adj = c_adj, cluster.df = cdf, t.df = tdf)),
  +                          2*pt(-abs(my_tstat), df))
  + 
  +                     # cat("adj = ", adj, " ; fixef.K = ", k_val, " ; cluster.adj = ", c_adj, " t.df = ", tdf, "\n", sep = "")
  +                 }
  +             }
  +         }
  +     }
  + }
  > 
  > 
  > #
  > # Comparison with sandwich and plm
  > #
  > 
  > library(sandwich)
  > 
  > # Data generation
  > set.seed(0)
  > N = 20; G = N/5; T = N/G
  > d = data.frame( y=rnorm(N), x=rnorm(N), grp=rep(1:G,T), tm=rep(1:T,each=G) )
  > 
  > # Estimations
  > est_lm    = lm(y ~ x + as.factor(grp) + as.factor(tm), data=d)
  > est_feols = feols(y ~ x | grp + tm, data=d)
  > 
  > #
  > # Standard
  > #
  > 
  > test(se(est_feols, se = "st")["x"], se(est_lm)["x"])
  > 
  > #
  > # Clustered
  > #
  > 
  > # Clustered by grp
  > se_CL_grp_lm_HC1 = sqrt(vcovCL(est_lm, cluster = d$grp, type = "HC1")["x", "x"])
  > se_CL_grp_lm_HC0 = sqrt(vcovCL(est_lm, cluster = d$grp, type = "HC0")["x", "x"])
  > 
  > # How to get the lm
  > test(se(est_feols, ssc = ssc(fixef.K = "full")), se_CL_grp_lm_HC1)
  > test(se(est_feols, ssc = ssc(adj = FALSE, fixef.K = "full")), se_CL_grp_lm_HC0)
  > 
  > #
  > # Heteroskedasticity-robust
  > #
  > 
  > se_white_lm_HC1 = sqrt(vcovHC(est_lm, type = "HC1")["x", "x"])
  > se_white_lm_HC0 = sqrt(vcovHC(est_lm, type = "HC0")["x", "x"])
  > 
  > test(se(est_feols, se = "hetero"), se_white_lm_HC1)
  > test(se(est_feols, se = "hetero", ssc = ssc(adj = FALSE, cluster.adj = FALSE)), se_white_lm_HC0)
  > 
  > #
  > # Two way
  > #
  > 
  > # Clustered by grp & tm
  > se_CL_2w_lm    = sqrt(vcovCL(est_lm, cluster = ~ grp + tm, type = "HC1")["x", "x"])
  > se_CL_2w_feols = se(est_feols, se = "twoway")
  > 
  > test(se(est_feols, se = "twoway", ssc = ssc(fixef.K = "full", cluster.df = "conv")), se_CL_2w_lm)
  > 
  > #
  > # Checking the calls work properly
  > #
  > 
  > data(trade)
  > 
  > est_pois = femlm(Euros ~ log(dist_km)|Origin+Destination, trade)
  > 
  > se_clust = se(est_pois, se = "cluster", cluster = "Product")
  > test(se(est_pois, cluster = trade$Product), se_clust)
  > test(se(est_pois, cluster = ~Product), se_clust)
  > 
  > se_two = se(est_pois, se = "twoway", cluster = trade[, c("Product", "Destination")])
  > test(se_two, se(est_pois, cluster = c("Product", "Destination")))
  > test(se_two, se(est_pois, cluster = ~Product+Destination))
  > 
  > se_clu_comb = se(est_pois, cluster = "Product^Destination")
  > test(se_clu_comb, se(est_pois, cluster = paste(trade$Product, trade$Destination)))
  > test(se_clu_comb, se(est_pois, cluster = ~Product^Destination))
  > 
  > se_two_comb = se(est_pois, cluster = c("Origin^Destination", "Product"))
  > test(se_two_comb, se(est_pois, cluster = list(paste(trade$Origin, trade$Destination), trade$Product)))
  > test(se_two_comb, se(est_pois, cluster = ~Origin^Destination + Product))
  > 
  > # With cluster removed
  > base = trade
  > base$Euros[base$Origin == "FR"] = 0
  > est_pois = femlm(Euros ~ log(dist_km)|Origin+Destination, base)
  > 
  > se_clust = se(est_pois, se = "cluster", cluster = "Product")
  > test(se(est_pois, cluster = base$Product), se_clust)
  > test(se(est_pois, cluster = ~Product), se_clust)
  > 
  > se_two = se(est_pois, se = "twoway", cluster = base[, c("Product", "Destination")])
  > test(se_two, se(est_pois, cluster = c("Product", "Destination")))
  > test(se_two, se(est_pois, cluster = ~Product+Destination))
  > 
  > se_clu_comb = se(est_pois, cluster = "Product^Destination")
  > test(se_clu_comb, se(est_pois, cluster = paste(base$Product, base$Destination)))
  > test(se_clu_comb, se(est_pois, cluster = ~Product^Destination))
  > 
  > se_two_comb = se(est_pois, cluster = c("Origin^Destination", "Product"))
  > test(se_two_comb, se(est_pois, cluster = list(paste(base$Origin, base$Destination), base$Product)))
  > test(se_two_comb, se(est_pois, cluster = ~Origin^Destination + Product))
  > 
  > # With cluster removed and NAs
  > base = trade
  > base$Euros[base$Origin == "FR"] = 0
  > base$Euros_na = base$Euros ; base$Euros_na[sample(nrow(base), 50)] = NA
  > base$Destination_na = base$Destination ; base$Destination_na[sample(nrow(base), 50)] = NA
  > base$Origin_na = base$Origin ; base$Origin_na[sample(nrow(base), 50)] = NA
  > base$Product_na = base$Product ; base$Product_na[sample(nrow(base), 50)] = NA
  > 
  > est_pois = femlm(Euros ~ log(dist_km)|Origin+Destination_na, base)
  > 
  > se_clust = se(est_pois, se = "cluster", cluster = "Product")
  > test(se(est_pois, cluster = base$Product), se_clust)
  > test(se(est_pois, cluster = ~Product), se_clust)
  > 
  > se_two = se(est_pois, se = "twoway", cluster = base[, c("Product", "Destination")])
  > test(se_two, se(est_pois, cluster = c("Product", "Destination")))
  > test(se_two, se(est_pois, cluster = ~Product+Destination))
  > 
  > se_clu_comb = se(est_pois, cluster = "Product^Destination")
  > test(se_clu_comb, se(est_pois, cluster = paste(base$Product, base$Destination)))
  > test(se_clu_comb, se(est_pois, cluster = ~Product^Destination))
  > 
  > se_two_comb = se(est_pois, cluster = c("Origin^Destination", "Product"))
  > test(se_two_comb, se(est_pois, cluster = list(paste(base$Origin, base$Destination), base$Product)))
  > test(se_two_comb, se(est_pois, cluster = ~Origin^Destination + Product))
  > 
  > #
  > # Checking errors
  > #
  > 
  > # Should report error
  > test(se(est_pois, cluster = "Origin_na"), "err")
  > test(se(est_pois, cluster = base$Origin_na), "err")
  > test(se(est_pois, cluster = list(base$Origin_na)), "err")
  > test(se(est_pois, cluster = ~Origin_na^Destination), "err")
  > 
  > test(se(est_pois, se = "cluster", cluster = ~Origin_na^not_there), "err")
  > 
  > #
  > # Checking that the aliases work fine
  > #
  > 
  > se_hetero = se(est_pois, se = "hetero")
  > se_hc1    = se(est_pois, se = "hc1")
  > se_white  = se(est_pois, se = "white")
  > 
  > test(se_hetero, se_hc1)
  > test(se_hetero, se_white)
  > 
  > #
  > # New argument vcov
  > #
  > 
  > # We mostly check the absence of errors
  > data(base_did)
  > 
  > est_panel = feols(y ~ x1, base_did, panel.id = ~id + period, subset = 1:500)
  > 
  > se_est = se(est_panel)
  > test(se(est_panel, ~id), se_est)
  > 
  > # changing ssc argument
  > test(se(est_panel, ssc = ssc(adj = FALSE)), se(est_panel, ~id + ssc(adj = FALSE)))
  > 
  > # using vcov_cluster
  > test(se_est, se(est_panel, vcov_cluster("id")))
  > test(se_est, se(vcov_cluster(est_panel, "id")))
  > 
  > # NW
  > se_NW = se(est_panel, "NW")
  > test(se_NW, se(est_panel, NW ~ id + period))
  > test(se_NW, se(est_panel, newey ~ id + period))
  > test(se_NW, se(est_panel, vcov_NW("id", "period")))
  > test(se_NW, se(est_panel, vcov_NW(time = "period"))) # here unit is deduced
  > 
  > se_NW2 = se(est_panel, NW(2))
  > test(se_NW2, se(est_panel, NW(2) ~ id + period))
  > test(se_NW2, se(est_panel, vcov_NW(lag = 2)))
  > 
  > # errors
  > est = feols(y ~ x1, base_did)
  > test(se(est, NW ~ period), "err")
  > 
  > # DK
  > se_DK = se(est_panel, "DK")
  > test(se_DK, se(est_panel, DK ~ period))
  > test(se_DK, se(est_panel, dris ~ period))
  > test(se_DK, se(est_panel, vcov_DK("period")))
  > 
  > se_DK2 = se(est_panel, DK(2))
  > test(se_DK2, se(est_panel, DK(2) ~ period))
  > test(se_DK2, se(est_panel, vcov_DK(lag = 2)))
  > 
  > 
  > # Conley
  > data(quakes)
  > 
  > est = feols(depth ~ mag, quakes, "conley")
  > 
  > se_conley = se(est)
  > test(se_conley, se(est, conley(90) ~ 1))
  > test(se_conley, se(est, conley(90) ~ lat + long))
  > 
  > se_conley200 = se(est, conley(200) ~ lat + long)
  > test(se_conley200, se(est, vcov_conley(cutoff = 200)))
  > test(se_conley200, se(est, vcov_conley("lat", "long", cutoff = 200)))
  > 
  > se_conleyExtra = se(est, conley(pixel = 20, distance = "spherical"))
  > test(se_conleyExtra, se(vcov_conley(est, pixel = 20, distance = "spherical")))
  > 
  > 
  > # Checking the value of Conley SEs with equivalences
  > # we generate data that leads to simple values
  > base = iris
  > names(base) = c("y", "x1", "x2", "x3", "species")
  > 
  > # scattered along 111km
  > base$lat = rep(seq(-0.5, 0.5, length.out = 50), 3)
  > 
  > # scattered across very long distances
  > base$lon = rep(c(0, 80, 160), each = 50)
  > 
  > est = feols(y ~ x1, base)
  > 
  > # Equivalence 1 -- clustered SEs
  > se_clu = se(est, ~lon + ssc(adj = FALSE, cluster.adj = FALSE))
  > test(se_clu, se(est, conley(200) ~ ssc(adj = FALSE)))
  malloc(): unsorted double linked list corrupted
  Aborted