Performance on 3+ fixed effects slow

[wklimowicz]

Performance on some data with compress is very slow, compared to feols. I'm not sure how much this is an "issue" with this package, versus a testament to how quick fixest is. The main issue seems to be the 3rd fixed effect, since up to 2 we can use strategy = within/demean, and its either a bit quicker or much quicker.

I was also getting "matrix too big for memory" type errors on dbreg compared to some regressions which took ~20m on feols, also perhaps related to strategy = "compress"?

strategy = "mundlak" works fast on everything I've tested, but it would be nice to have a drop in replacement for feols.

(This was run on #49 , since I also got matrix inversion errors when running on master: might be related?)

# From PR 49: on main it errors out with matrix inversion.
# pak::pak("grantmcdermott/dbreg#49")

library(duckdb)
library(fixest)
library(dplyr)
library(dbplyr)
library(dbreg)

# From 100k to 500k rows, increasing by 100k
n_rows <- (1:5)*1e5

benchmark_single_regression <- function(n_rows) {
  set.seed(44)

  df <- data.frame(
    fe1 = rep(1:(n_rows / 100), each = 100),
    fe2 = rep(1:(n_rows / 200), each = 200),
    fe3 = rep(1:(n_rows / 500), each = 500),
    x = rnorm(n_rows),
    y = rbinom(n_rows, 1, 0.2)
  )

  con <- dbConnect(duckdb())
  df_db <- copy_to(con, df, overwrite = TRUE)

  t1 <- Sys.time()
  dbreg(
    table = df_db,
    strategy = "auto",
    verbose = TRUE,
    y ~ x | fe1 + fe2 + fe3,
  )
  t2 <- Sys.time()
  time_dbreg <- as.numeric(t2 - t1, units = "secs")

  t1 <- Sys.time()
  feols(
    data = df,
    y ~ x | fe1 + fe2 + fe3,
  )
  t2 <- Sys.time()
  time_feols <- as.numeric(t2 - t1, units = "secs")

  return(data.frame(
    n_rows = n_rows,
    type = c("feols", "dbreg"),
    time = c(time_feols, time_dbreg)
  ))
}

results <- lapply(n_rows, benchmark_single_regression) |> dplyr::bind_rows()

Number of rows	feols (seconds)	dbreg (seconds)
1e+05	0.018	1.451
2e+05	0.035	30.890
3e+05	0.044	105.122
4e+05	0.057	258.883
5e+05	0.077	550.779

[jamesbrandecon]

Thanks for continuing to test! Unless I'm forgetting something, compress is not really meant for cases with continuous x. Compression is most value when (xvars, fe) take on not very many total values. strategy = mundlak, on the other hand, is fine with continuous variables.

[grantmcdermott]

+1 to @jamesbrandecon’s comment. I plan to eventually add a "When should I (not) use dbreg?" section to the vignette and this is a classic case where the compression logic won’t buy you much. (Same story for "true" panels, where each id+time observation represents a unique combination.)

Fwiw, I think we may have a path to proper FWL partialling-out in dbreg, which would then open the door to alternating projections (AP) for high-dimensional FEs. AP is generally slower than fixest’s algo (which uses a clever fixed point solution), but it’s still quite competitive and also works cleanly with unbalanced panels.

All of this is still far in the future (and currently untested). In the meantime, I'd also say there's nothing inherently "wrong" with Mundlak/CRE. In fact, it’s arguably a more defensible estimand than standard FE in a lot of applied settings (FEs aren't truly "fixed" over time, staggered DiD, etc.). YMMV

[wklimowicz]

Thanks both! I've been playing around and got a bit further on "debugging" my specific use case, e.g. why mundlak was so different from feols. I think I've managed to narrow it down to both of these things being true:

1. Correlation between the x's and fixed effects
2. Unbalanced panel

From my (limited) understanding and testing, if I just have correlation (but balanced panel), mundlak is identical. If I have an unbalanced panel but no correlation the results are very similar. If you have both, mundlak can give biased results:

library(tidyverse)
library(fixest)
library(dbreg)

set.seed(123)

# 1. Simulate Data
n_ids <- 2000
n_times <- 5
n_grps <- 50

data_unbalanced <- expand.grid(id = 1:n_ids, time = 1:n_times) %>%
  as_tibble() %>%
  mutate(
    grp = sample(1:n_grps, n(), replace = TRUE),
    id_fe = rnorm(n_ids)[id],
    grp_fe = rnorm(n_grps)[grp],
    time_fe = rnorm(n_times)[time],
    # x1 is correlated with both fe's
    x1 = 4.0 * grp_fe + 0.5 * id_fe + rnorm(n()),
    # True coefficient for x1 is -1.0
    y = -1.0 * x1 + id_fe + 5.0 * grp_fe + time_fe + rnorm(n())
  )

# 2. feols: true coefficient recovered
fit_feols <- feols(y ~ x1 | id_fe + grp_fe + time_fe, data = data_unbalanced)
print(fit_feols)
#    Estimate Std. Error  t value  Pr(>|t|)    
# x1 -1.00534   0.011234 -89.4907 < 2.2e-16 ***

# 3. Mundlak: biased?
fit_dbreg <- dbreg(
  y ~ x1 | id_fe + grp_fe + time_fe,
  strategy = "mundlak",
  data = data_unbalanced
)
print(fit_dbreg)
#     Estimate Std. Error  t value  Pr(>|t|)    
# x1 -0.634544    0.01393 -45.5514 < 2.2e-16 ***

Some prompting got me this "add second-order means by hand to fix the bias", which seems to work:

# Augment data with group means and second-order means
data_aug <- data_unbalanced %>%
  group_by(grp) %>% mutate(x_mean_grp = mean(x1)) %>% ungroup() %>%
  group_by(id)  %>% mutate(x_mean_id  = mean(x1)) %>% ungroup() %>%
  group_by(time)%>% mutate(x_mean_time= mean(x1)) %>% ungroup() %>%
  # Second-Order Means
  group_by(id)  %>% mutate(m_id_grp = mean(x_mean_grp)) %>% ungroup() %>%
  group_by(grp) %>% mutate(m_grp_id = mean(x_mean_id))  %>% ungroup()

fit <- dbreg(
  y ~ x1 + x_mean_id + x_mean_grp + x_mean_time + m_id_grp + m_grp_id,
  strategy = "moments",
  data = data_aug
)
print(fit)
#              Estimate Std. Error    t value   Pr(>|t|)
# (Intercept)  0.291749   0.035963   8.112452 5.5440e-16
# x1          -1.002012   0.015689 -63.868868  < 2.2e-16

I plan to eventually add a "When should I (not) use dbreg?" section to the vignette

This vignette is probably the only thing that could come out of this issue, since the problem isn't with the package. Unless there's econometric theory out there where these second order averages are justified in a formal way where it can be automated (and you guys will be more familiar with the theory), I think we can close this.

I was half-hoping that simply because dbreg runs in duckdb it would be quicker everwhere, but there's already a lot of C++ code and clever optimisations in fixest: I think for the specific use case I'm working on right now I'll just stick with that.

Fwiw, I think we may have a path to proper FWL partialling-out in dbreg, which would then open the door to alternating projections (AP) for high-dimensional FEs. AP is generally slower than fixest’s algo (which uses a clever fixed point solution), but it’s still quite competitive and also works cleanly with unbalanced panels.

This would be great: being able to run arbitrary regressions quickly without leaving tidyverse/dbplyr would be a huge QOL improvement for me: the main problem with feols is that I need to download the data from a database first.

[grantmcdermott]

Your "second-order means" example is actually a good one, and something that I wrestled with whilst writing #49. I ultimately decided to keep it simple, since correctly expanding (nesting) these terms can be a PITA. But I've had it in the back of mind to revisit; especially once/if we expand our coverage to dedicated staggered event study models.

[grantmcdermott]

P.S. As an aside, in newer variants of dplyr, you can use the ephemeral .by = <by> argument instead of group_by and then having to worry about ungrouping. https://tidyverse.org/blog/2023/02/dplyr-1-1-0-per-operation-grouping/