`fixest`: Incorrect standard errors and error

[vincentarelbundock]

Two problems with the fixest package, initially reported by @passtar

1. Incorrect standard errors
2. Error

Incorrect standard errors

The parameters() function does not return the same summary() function:

library(did)
library(parameters)
library(fixest)

data(mpdta) # from the did package

model_att <- feols(
    lemp ~ sunab(first.treat, year, ref.p = -1:-4, att = TRUE) | countyreal + year,
    data = mpdta, cluster = ~countyreal)

summary(model_att)
> OLS estimation, Dep. Var.: lemp
> Observations: 2,500
> Fixed-effects: countyreal: 500,  year: 5
> Standard-errors: Clustered (countyreal) 
>     Estimate Std. Error  t value   Pr(>|t|)    
> ATT -0.04771   0.013265 -3.59669 0.00035438 ***
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> RMSE: 0.123851     Adj. R2: 0.991531
>                  Within R2: 0.01012

parameters(model_att, digits = 5)
> # Fixed Effects
> 
> Parameter | Coefficient |      SE |               95% CI |   t(499) |      p
> ----------------------------------------------------------------------------
> ATT       |    -0.04771 | 0.02238 | [-0.09168, -0.00374] | -3.59669 | < .001
> 
> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
>   using a Wald t-distribution approximation.

Error

model_es <- feols(
    lemp ~ sunab(first.treat, year, ref.p = -1:-4) | countyreal + year,
    data = mpdta, cluster = ~countyreal)

summary(model_es)
> OLS estimation, Dep. Var.: lemp
> Observations: 2,500
> Fixed-effects: countyreal: 500,  year: 5
> Standard-errors: Clustered (countyreal) 
>          Estimate Std. Error  t value   Pr(>|t|)    
> year::0 -0.031067   0.013621 -2.28083 0.02297891 *  
> year::1 -0.052235   0.018873 -2.76772 0.00585468 ** 
> year::2 -0.136078   0.035455 -3.83800 0.00013997 ***
> year::3 -0.104707   0.033874 -3.09106 0.00210591 ** 
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> RMSE: 0.123851     Adj. R2: 0.991531
>                  Within R2: 0.01012

parameters(model_es, digits = 5)
> Error: Something went wrong... :-/

[strengejacke]

Will look into it. A bit strange, because we directly retrieve the SE from the summary (https://github.com/easystats/parameters/blob/main/R/methods_fixest.R).

[strengejacke]

Ok, that's wird. Sometimes, SEs returned by summary() differ from those in model$coeftable. Sometimns, both return SEs of different length. But I think I have it now, hopefully fixed in #1040

[strengejacke]

data(mpdta, package = "did")
m <- fixest::feols(
  lemp ~ sunab(first.treat, year, ref.p = -1:-4, att = TRUE) | countyreal + year,
  data = mpdta,
  cluster = ~countyreal
)
as.vector(m$coeftable[, "Std. Error"]) # same as summary(m)
as.vector(summary(m)$se) # *not* same as summary(m)

data(base_stagg, package = "fixest")
m <- fixest::feols(y ~ x1 + sunab(year_treated, year) | id + year, base_stagg)
as.vector(m$coeftable[, "Std. Error"]) # same as summary(m)
as.vector(summary(m)$se) # *not* same as summary(m)


data(base_stagg, package = "fixest")
m <- fixest::feols(y ~ x1 + year_treated | id + year, base_stagg)
as.vector(m$coeftable[, "Std. Error"]) # *not* same as summary(m)
as.vector(summary(m)$se) # same as summary(m)

[vincentarelbundock]

@lrberge sorry to bug you, but our goal here is to retrieve the "main" standard errors that the user has requested when fitting a fixest model. Is this logic adequate to cover all cases? Is the unequal length the only case where we should switch from coeftable to summary()$se?

    stats <- summary(model)
    SE <- as.vector(stats$se)
    # sometimes, the correct SEs are in the `coeftable`. In this case,
    # SE from summary have different length
    if (length(SE) != nrow(params)) {
      SE <- NULL
      SE_ct <- .safe(as.vector(model$coeftable[, "Std. Error"]))
      if (!is.null(SE_ct) && length(SE_ct) == nrow(params)) {
        SE <- SE_ct
      }
    }

[lrberge]

Hi Vincernt! Ahh, this is a very special case... totally specific to the sunab function.

When using the sunab function in a fixest call:

• all the treatment_cohort x relative_period are estimated, leading to about N_cohort x N_period coefficients
• the user usually does not want that, but only aggregated statictics based on these coefficients, namely either one coefficient for each relative period, or just one coefficient for the total ATT
• the aggregation is (and must) be a post-estimation procedure, such that the user can request several aggregation types (ATT, or ATT x relative period) without having to re-estimate a costly model
• the aggregation is performed during summary() where all the user requests are resolved. Summary leads to the "right" coeftable
• given that the aggregation of coefficients is an artificial construct out of the initial estimation with N_cohort x N_period coefficients, I do not change the VCOV, and hence the SEs (in the $se slot), that's why the number of rows in summary()$coeftable is different from the number of observations in summary()$se (which refers to the SEs of the original estimation)

tl;dr using summary(model)$coeftable[, "Std. Error"] is the right way, it always works.

Note that is should be summary(model)$coeftable and not model$coeftable, because summary needs to be proc'ed.

On your specific question: the only case I can think of where len(se) > nrow(coeftable) is when there is coefficient aggregation (which happens in sunab). Using summary(model)$se is otherwise safe (or at least I can't think of a counter example).

BTW, there is work under way to compute and return the aggregated VCOV, but I still haven't settled on the design.

Let me know if more information is needed. Sorry to cause you this trouble!

[vincentarelbundock]

Thanks a lot for the info @lrberge I appreciate it!

@strengejacke so this means we need to remove the first part where you take stats$se and always go through summary(model)$coeftable. Do you need help or are you OK?

[strengejacke]

Great, thanks for clairifaction, @lrberge! @vincentarelbundock I'm ok, will revise the PR accordingly.