Running 'fixest_tests.R' [11s]
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(dreamerr) ; library(fixest)
  > 
  > 
  > 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: ", dreamerr::sfill(model, 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(x2, species) + factor(fe_2) + i(x3, fe_2) + 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(x2, fe_2) + i(x3, fe_2) + 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$obsRemoved)){
  +                         qui = -res$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", dof = dof(adj = adj))["x1"], se(res_bis)["x1"], "~", tol)
  +                 test(pvalue(res, se = "st", dof = dof(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)
  > 
  > 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(c(y, x1) ~ 1 | fe1 | x2 ~ x3, 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))
  > 
  > 
  > ####
  > #### 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, dof = dof(fixef.K = "nested"), attr = TRUE), "dof.K"), 6)
  > 
  > # fixef.K = FALSE
  > #  => adjustment K = 1 (i.e. only x)
  > test(attr(vcov(est, dof = dof(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, dof = dof(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, dof = dof(fixef.K = "full", fixef.force_exact = TRUE), attr = TRUE), "dof.K"), 7)
  > 
  > #
  > # Manual checks of the SEs
  > #
  > 
  > n = est$nobs
  > VCOV_raw = est$cov.unscaled / ((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", dof = dof(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.unscaled / est$sigma2
  > H = vcovClust(est$fixef_id$fe1, VCOV_raw, scores = est$scores, dof = FALSE)
  > 
  > 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 - 8)
  +             my_adj = ifelse(adj, (n - 1) / (n - K), 1)
  + 
  +             V = H * cluster_factor
  + 
  +             # test SE
  +             test(vcov(est, se = "cluster", dof = dof(adj = adj, fixef.K = k_val, cluster.adj = c_adj)), V * my_adj)
  + 
  +             # test pvalue
  +             my_tstat = tstat(est, se = "cluster", dof = dof(adj = adj, fixef.K = k_val, cluster.adj = c_adj))
  +             test(pvalue(est, se = "cluster", dof = dof(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.unscaled / est$sigma2
  > M_i  = vcovClust(est$fixef_id$fe1, VCOV_raw, scores = est$scores, dof = FALSE)
  > M_t  = vcovClust(est$fixef_id$fe2, VCOV_raw, scores = est$scores, dof = FALSE)
  > M_it = vcovClust(paste(base$fe1, base$fe2), VCOV_raw, scores = est$scores, dof = FALSE, do.unclass = TRUE)
  > 
  > M_i + M_t - M_it
              [,1]
  [1,] 0.005391594
  > vcov(est, se = "two", dof = dof(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 - 8)
  +                     my_adj = ifelse(adj, (n - 1) / (n - K), 1)
  + 
  +                     # test SE
  +                     test(vcov(est, se = "two", dof = dof(adj = adj, fixef.K = k_val, cluster.adj = c_adj, cluster.df = cdf)),
  +                          V * my_adj)
  + 
  +                     # test pvalue
  +                     my_tstat = tstat(est, se = "two", dof = dof(adj = adj, fixef.K = k_val, cluster.adj = c_adj, cluster.df = cdf))
  +                     test(pvalue(est, se = "two", dof = dof(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"])
  > se_CL_grp_stata  = 0.165385 # vce(cluster grp)
  > 
  > # How to get the lm
  > test(se(est_feols, dof = dof(fixef.K = "full")), se_CL_grp_lm_HC1)
  > test(se(est_feols, dof = dof(adj = FALSE, fixef.K = "full")), se_CL_grp_lm_HC0)
  > 
  > # How to get the Stata
  > test(se(est_feols), se_CL_grp_stata, "~")
  > 
  > #
  > # 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", dof = dof(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", dof = dof(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")
  > 
  > test(se(est_pois, se = "twoway", cluster = c("Origin^Destination", "Product", "error")), "err")
  > test(se(est_pois, se = "twoway", cluster = base[, 1:4]), "err")
  > test(se(est_pois, se = "twoway", cluster = ~Origin+Destination+Product), "err")
  > test(se(est_pois, se = "fourway", cluster = ~Origin+Destination+Product), "err")
  > 
  > ####
  > #### Residuals ####
  > ####
  > 
  > chunk("RESIDUALS")
  RESIDUALS 
  
  > 
  > base = iris
  > names(base) = c("y", "x1", "x2", "x3", "species")
  > base$y_int = as.integer(base$y) + 1
  > 
  > # OLS + GLM + FENMLM
  > 
  > for(method in c("ols", "feglm", "femlm", "fenegbin")){
  +     cat("Method: ", dreamerr::sfill(method, 8))
  +     for(do_weight in c(FALSE, TRUE)){
  +         cat(".")
  + 
  +         if(do_weight){
  +             w = unclass(as.factor(base$species))
  +         } else {
  +             w = NULL
  +         }
  + 
  +         if(method == "ols"){
  +             m = feols(y_int ~ x1 | species, base, weights = w)
  +             mm = lm(y_int ~ x1 + species, base, weights = w)
  + 
  +         } else if(method == "feglm"){
  +             m = feglm(y_int ~ x1 | species, base, weights = w)
  +             mm = glm(y_int ~ x1 + species, base, weights = w, family = poisson())
  + 
  +         } else if(method == "femlm"){
  +             if(!is.null(w)) next
  +             m = femlm(y_int ~ x1 | species, base)
  +             mm = glm(y_int ~ x1 + species, base, family = poisson())
  + 
  +         } else if(method == "fenegbin"){
  +             if(!is.null(w)) next
  +             m = fenegbin(y_int ~ x1 | species, base, notes = FALSE)
  +             mm = MASS::glm.nb(y_int ~ x1 + species, base)
  +         }
  + 
  +         tol = ifelse(method == "fenegbin", 1e-2, 1e-6)
  + 
  +         test(resid(m, "r"), resid(mm, "resp"), "~", tol = tol)
  +         test(resid(m, "d"), resid(mm, "d"), "~", tol = tol)
  +         test(resid(m, "p"), resid(mm, "pearson"), "~", tol = tol)
  + 
  +         test(deviance(m), deviance(mm), "~", tol = tol)
  +     }
  +     cat("\n")
  + }
  Method:       ols..
  Method:     feglm..
  Method:     femlm..
  Method:  fenegbin..
  Warning messages:
  1: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace >  :
    iteration limit reached
  2: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace >  :
    iteration limit reached
  > cat("\n")
  
  > 
  > 
  > ####
  > #### fixef ####
  > ####
  > 
  > chunk("FIXEF")
  FIXEF 
  
  > 
  > set.seed(0)
  > base = iris
  > names(base) = c("y", "x1", "x2", "x3", "species")
  > base$x4 = rnorm(150) + 0.25*base$y
  > base$fe_bis = sample(10, 150, TRUE)
  > base$fe_ter = sample(15, 150, TRUE)
  > 
  > get_coef = function(all_coef, x){
  +     res = all_coef[grepl(x, names(all_coef))]
  +     names(res) = gsub(x, "", names(res))
  +     res
  + }
  > 
  > #
  > # With 2 x 1 FE
  > #
  > 
  > m = feols(y ~ x1 + x2 | species + fe_bis, base)
  > all_coef = coef(feols(y ~ -1 + x1 + x2 + species + factor(fe_bis), base))
  > 
  > m_fe = fixef(m)
  > c1  = get_coef(all_coef, "species")
  > test(var(c1 - m_fe$species[names(c1)]), 0)
  > 
  > c2  = get_coef(all_coef, "factor\\(fe_bis\\)")
  > test(var(c2 - m_fe$fe_bis[names(c2)]), 0)
  > 
  > 
  > #
  > # With 1 FE + 1 FE 1 VS
  > #
  > 
  > m = feols(y ~ x1 + x2 | species + fe_bis[x3], base)
  > all_coef = coef(feols(y ~ -1 + x1 + x2 + species + factor(fe_bis) + i(x3, fe_bis), base))
  > 
  > m_fe = fixef(m)
  > c1  = get_coef(all_coef, "species")
  > test(var(c1 - m_fe$species[names(c1)]), 0, "~")
  > 
  > c2  = get_coef(all_coef, "factor\\(fe_bis\\)")
  > test(var(c2 - m_fe$fe_bis[names(c2)]), 0, "~")
  > 
  > c3  = get_coef(all_coef, "x3:fe_bis::")
  > test(c3, m_fe[["fe_bis[[x3]]"]][names(c3)], "~", tol = 1e-5)
  > 
  > #
  > # With 2 x (1 FE + 1 VS) + 1 FE
  > #
  > 
  > m = feols(y ~ x1 | species[x2] + fe_bis[x3] + fe_ter, base)
  > all_coef = coef(feols(y ~ -1 + x1 + species + i(x2, species) + factor(fe_bis) + i(x3, fe_bis) + factor(fe_ter), base))
  > 
  > m_fe = fixef(m)
  > c1  = get_coef(all_coef, "^species")
  > test(var(c1 - m_fe$species[names(c1)]), 0, "~")
  > 
  > c2  = get_coef(all_coef, "^factor\\(fe_bis\\)")
  > test(var(c2 - m_fe$fe_bis[names(c2)]), 0, "~")
  > 
  > c3  = get_coef(all_coef, "x3:fe_bis::")
  > test(c3, m_fe[["fe_bis[[x3]]"]][names(c3)], "~", tol = 2e-4)
  > 
  > c4  = get_coef(all_coef, "x2:species::")
  > test(c4, m_fe[["species[[x2]]"]][names(c4)], "~", tol = 2e-4)
  > 
  > #
  > # With 2 x (1 FE) + 1 FE 2 VS
  > #
  > 
  > m = feols(y ~ x1 | species + fe_bis[x2,x3] + fe_ter, base)
  > all_coef = coef(feols(y ~ x1 + species + factor(fe_bis) + i(x2, fe_bis) + i(x3, fe_bis) + factor(fe_ter), base))
  > 
  > m_fe = fixef(m)
  > c1  = get_coef(all_coef, "^species")
  > test(var(c1 - m_fe$species[names(c1)]), 0, "~")
  > 
  > c2  = get_coef(all_coef, "^factor\\(fe_bis\\)")
  > test(var(c2 - m_fe$fe_bis[names(c2)]), 0, "~")
  > 
  > c3  = get_coef(all_coef, "x2:fe_bis::")
  > test(c3, m_fe[["fe_bis[[x2]]"]][names(c3)], "~", tol = 2e-4)
  > 
  > c4  = get_coef(all_coef, "x3:fe_bis::")
  > test(c4, m_fe[["fe_bis[[x3]]"]][names(c4)], "~", tol = 2e-4)
  > 
  > 
  > #
  > # With weights
  > #
  > 
  > w = 3 * (as.integer(base$species) - 0.95)
  > m = feols(y ~ x1 | species + fe_bis[x2,x3] + fe_ter, base, weights = w)
  > all_coef = coef(feols(y ~ x1 + species + factor(fe_bis) + i(x2, fe_bis) + i(x3, fe_bis) + factor(fe_ter), base, weights = w))
  > 
  > m_fe = fixef(m)
  > c1  = get_coef(all_coef, "^species")
  > test(var(c1 - m_fe$species[names(c1)]), 0, "~")
  > 
  > c2  = get_coef(all_coef, "^factor\\(fe_bis\\)")
  > test(var(c2 - m_fe$fe_bis[names(c2)]), 0, "~")
  > 
  > c3  = get_coef(all_coef, "x2:fe_bis::")
  > test(c3, m_fe[["fe_bis[[x2]]"]][names(c3)], "~", tol = 2e-4)
  > 
  > c4  = get_coef(all_coef, "x3:fe_bis::")
  > test(c4, m_fe[["fe_bis[[x3]]"]][names(c4)], "~", tol = 2e-4)
  > 
  > 
  > ####
  > #### To Integer ####
  > ####
  > 
  > chunk("TO_INTEGER")
  TO_INTEGER 
  
  > 
  > base = iris
  > names(base) = c("y", "x1", "x2", "x3", "species")
  > base$z = sample(5, 150, TRUE)
  > 
  > # Normal
  > m = to_integer(base$species)
  > test(length(unique(m)), 3)
  > 
  > m = to_integer(base$species, base$z)
  > test(length(unique(m)), 15)
  > 
  > # with NA
  > base$species_na = base$species
  > base$species_na[base$species == "setosa"] = NA
  > 
  > m = to_integer(base$species_na, base$z)
  > test(length(unique(m)), 11)
  > 
  > m = to_integer(base$species_na, base$z, add_items = TRUE, items.list = TRUE)
  > test(length(m$items), 10)
  > 
  > 
  > ####
  > #### 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", dof = dof(adj=adj))["x1"], se(res_bis)["x1"], "~")
  +             # cat("Weight: ", useWeights, ", model: ", model, ", FE: ", use_fe, "\n", sep="")
  + 
  +         }
  +     }
  + }
  ........> cat("\n")
  
  > 
  > 
  > ####
  > #### Interact ####
  > ####
  > 
  > 
  > chunk("Interact")
  INTERACT 
  
  > 
  > base = iris
  > names(base) = c("y", "x1", "x2", "x3", "species")
  > base$fe_2 = round(rnorm(150))
  > 
  > #
  > # We just ensure it works without error
  > #
  > 
  > m = feols(y ~ x1 + i(fe_2), base)
  > coefplot(m)
  > etable(m, dict = c("0" = "zero"))
                                  m
  Dependent Var.:                 y
                                   
  (Intercept)     6.267*** (0.8997)
  x1               -0.2669. (0.155)
  fe_2 = -2          0.4708 (0.855)
  fe_2 = -1           0.518 (0.827)
  fe_2 = zero       0.3394 (0.8238)
  fe_2 = 1          0.1872 (0.8272)
  fe_2 = 2           1.166 (0.8753)
  fe_2 = 3          -0.3331 (1.159)
  _______________ _________________
  S.E. type                Standard
  Observations                  150
  R2                       0.083689
  Adj. R2                  0.038519
  > 
  > m = feols(y ~ x1 + i(fe_2) + i(x2, fe_2), base)
  > coefplot(m, only.inter = FALSE)
  > etable(m, dict = c("0" = "zero"))
                                     m
  Dependent Var.:                    y
                                      
  (Intercept)        4.196*** (0.3714)
  x1               0.5616*** (0.07032)
  fe_2 = -2         -1.658*** (0.4188)
  fe_2 = -1         -1.992*** (0.3579)
  fe_2 = zero       -1.825*** (0.3509)
  fe_2 = 1           -1.722*** (0.356)
  fe_2 = 2           -1.544*** (0.452)
  fe_2 = 3             -1.162* (0.468)
  x2 x fe_2 = -2   0.4329*** (0.05832)
  x2 x fe_2 = -1    0.5021*** (0.0295)
  x2 x fe_2 = zero 0.4568*** (0.02712)
  x2 x fe_2 = 1    0.4358*** (0.03367)
  x2 x fe_2 = 2    0.4999*** (0.06698)
  ________________ ___________________
  S.E. type                   Standard
  Observations                     150
  R2                          0.856442
  Adj. R2                     0.843868
  > 
  > a = i(base$fe_2)
  > b = i(base$fe_2, drop = 0:1)
  > d = i(base$fe_2, keep = 0:1)
  > 
  > test(ncol(a), ncol(b) + 2)
  > test(ncol(d), 2)
  > 
  > 
  > 
  > ####
  > #### demean ####
  > ####
  > 
  > chunk("DEMEAN")
  DEMEAN 
  
  > 
  > data(trade)
  > 
  > base = trade
  > base$ln_euros = log(base$Euros)
  > base$ln_dist = log(base$dist_km)
  > 
  > X = base[, c("ln_euros", "ln_dist")]
  > fe = base[, c("Origin", "Destination")]
  > 
  > base_new = demean(X, fe)
  Error in anyNA.data.frame(f) : could not find function "anyNA.data.frame"
  Calls: demean
  Execution halted