Marginal Effect Confidence Intervals and Standard Errors after Fixed Effects Poisson/PPML

[mikedenly]

Dear Vincent et al, I would like to kindly request your assistance with discerning why marginal effect confidence intervals and standard errors are different in Stata and R after estimating a fixed effects poisson pseudo-maximum likelihood (PMML) model. In R I am using fixest::fepois() and marginaleffects::avg_slopes(). In Stata I am using ppmlhdfe and margins, dydx(*). As one can see from my reproducible example below using the Fatalities dataset from R's AER package, everything is the same in the original models, and the marginal effect coefficients are the same, too. However, I can't obtain the same standard errors or confidence intervals for the marginal effects in R, even after trying different HC types and clustering options. This is odd because, in my reading of package documentation, both Stata and the marginaleffects package in R use the delta method to calculate standard errors.

Here is the Stata output:

. * import the data
. import delimited using  "fatalities.csv"
(34 vars, 336 obs)

* run the model
ppmlhdfe fatal beertax, absorb(state year) vce(cl state) d
Iteration 1:   deviance = 7.5532e+03  eps = .         iters = 3    tol = 1.0e-04  min(eta) =  -1.43  P   
Iteration 2:   deviance = 1.6438e+03  eps = 3.59e+00  iters = 2    tol = 1.0e-04  min(eta) =  -1.97      
Iteration 3:   deviance = 1.4156e+03  eps = 1.61e-01  iters = 2    tol = 1.0e-04  min(eta) =  -2.18      
Iteration 4:   deviance = 1.4138e+03  eps = 1.25e-03  iters = 2    tol = 1.0e-04  min(eta) =  -2.21      
Iteration 5:   deviance = 1.4138e+03  eps = 1.64e-07  iters = 2    tol = 1.0e-04  min(eta) =  -2.21      
Iteration 6:   deviance = 1.4138e+03  eps = 3.08e-15  iters = 1    tol = 1.0e-05  min(eta) =  -2.21   S  
Iteration 7:   deviance = 1.4138e+03  eps = 1.47e-16  iters = 1    tol = 1.0e-08  min(eta) =  -2.21   S O
------------------------------------------------------------------------------------------------------------
(legend: p: exact partial-out   s: exact solver   h: step-halving   o: epsilon below tolerance)
Converged in 7 iterations and 13 HDFE sub-iterations (tol = 1.0e-08)

HDFE PPML regression                              No. of obs      =        336
Absorbing 2 HDFE groups                           Residual df     =         47
Statistics robust to heteroskedasticity           Wald chi2(1)    =       4.04
Deviance             =  1413.784422               Prob > chi2     =     0.0445
Log pseudolikelihood = -2095.542592               Pseudo R2       =     0.9825

Number of clusters (state)  =         48
                                 (Std. Err. adjusted for 48 clusters in state)
------------------------------------------------------------------------------
             |               Robust
       fatal |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     beertax |  -.3472736   .1728049    -2.01   0.044    -.6859649   -.0085822
       _cons |   7.396313   .0928461    79.66   0.000     7.214338    7.578288
------------------------------------------------------------------------------

Absorbed degrees of freedom:
-----------------------------------------------------+
 Absorbed FE | Categories  - Redundant  = Num. Coefs |
-------------+---------------------------------------|
       state |        48          48           0    *|
        year |         7           1           6     |
-----------------------------------------------------+
* = FE nested within cluster; treated as redundant for DoF computation

. margins, dydx(*)

Average marginal effects                        Number of obs     =        336
Model VCE    : Robust

Expression   : Predicted mean of fatal, predict()
dy/dx w.r.t. : beertax

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     beertax |  -322.5004   160.4776    -2.01   0.044    -637.0308   -7.969972
------------------------------------------------------------------------------

Here are my multiple attempts in R using marginaleffects::avg_slopes():

# libraries
library(marginaleffects)
library(fixest)
library(AER)
library(sandwich)

# load fatalities data from AER package
data(Fatalities)

# run regular model with state & year fixed effects -- matches Stata
regular_model <- fixest::fepois(fatal ~ beertax | state + year, data=Fatalities, cluster=~state)
summary(regular_model)

Poisson estimation, Dep. Var.: fatal
Observations: 336
Fixed-effects: state: 48,  year: 7
Standard-errors: Clustered (state) 
         Estimate Std. Error  z value Pr(>|z|)    
beertax -0.347274   0.174639 -1.98852 0.046754 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Log-Likelihood: -2,095.5   Adj. Pseudo R2: 0.981999
           BIC:  4,511.0     Squared Cor.: 0.992513

# Compute the clustered standard errors by state
vcov_clusterDefault <- sandwich::vcovCL(regular_model, cluster = ~state)
vcov_clusterHC0 <- sandwich::vcovCL(regular_model, cluster = ~state, type = "HC0")
vcov_clusterHC1 <- sandwich::vcovCL(regular_model, cluster = ~state, type = "HC1")

# Calculate the marginal effects with clustered standard errors
> marginaleffects::avg_slopes(regular_model, vcov = vcov_clusterDefault)

    Term Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
 beertax     -322        131 -2.47   0.0135 6.2  -578  -66.6

Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
Type:  response 

> marginaleffects::avg_slopes(regular_model, vcov = vcov_clusterHC0)

    Term Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
 beertax     -322        131 -2.47   0.0135 6.2  -578  -66.6

Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
Type:  response 

> marginaleffects::avg_slopes(regular_model, vcov = vcov_clusterHC1)

    Term Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
 beertax     -322        131 -2.47   0.0135 6.2  -578  -66.6

Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
Type:  response 

# Only robust standard errors
> marginaleffects::avg_slopes(regular_model, vcov = "HC1")

    Term Estimate Std. Error     z Pr(>|z|)    S 2.5 % 97.5 %
 beertax     -322       81.8 -3.94   <0.001 13.6  -483   -162

Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
Type:  response 

> 
> # Default standard errors
> marginaleffects::avg_slopes(regular_model)

    Term Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
 beertax     -322        132 -2.44   0.0145 6.1  -581  -63.9

Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
Type:  response 

Does anyone know why this is happening and how Stata is obtaining the constant? I thought it should be absorbed the fixed effects. By the way, Vincent, your packages are incredible -- I teach them to all of my students. Thank you for everything and apologies if this should have gone on StackOverflow -- I wasn't sure about the correct forum for my question.

And here is sessionInfo():

``` > sessionInfo() R version 4.3.2 (2023-10-31 ucrt) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 11 x64 (build 22631) Matrix products: default locale: [1] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8 LC_MONETARY=English_United States.utf8 [4] LC_NUMERIC=C LC_TIME=English_United States.utf8 time zone: America/Chicago tzcode source: internal attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] rio_1.0.1 marginaleffects_0.16.0 AER_1.2-12 survival_3.5-7 sandwich_3.1-0 [6] lmtest_0.9-40 zoo_1.8-12 car_3.1-2 carData_3.0-5 broom_1.0.5 [11] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2 [16] readr_2.1.4 tidyr_1.3.0 tibble_3.2.1 ggplot2_3.5.0 tidyverse_2.0.0 [21] fixest_0.12.0 loaded via a namespace (and not attached): [1] gtable_0.3.5 insight_0.19.11 lattice_0.21-9 tzdb_0.4.0 numDeriv_2016.8-1.1 [6] vctrs_0.6.5 tools_4.3.2 generics_0.1.3 fansi_1.0.5 pkgconfig_2.0.3 [11] R.oo_1.25.0 Matrix_1.6-5 data.table_1.14.8 checkmate_2.3.1 stringmagic_1.1.2 [16] lifecycle_1.0.4 compiler_4.3.2 munsell_0.5.1 Formula_1.2-5 pillar_1.9.0 [21] R.utils_2.12.3 abind_1.4-5 nlme_3.1-163 tidyselect_1.2.1 mvtnorm_1.2-4 [26] stringi_1.8.2 splines_4.3.2 grid_4.3.2 colorspace_2.1-0 cli_3.6.1 [31] magrittr_2.0.3 utf8_1.2.4 withr_3.0.0 dreamerr_1.4.0 scales_1.3.0 [36] backports_1.4.1 timechange_0.2.0 estimability_1.4.1 emmeans_1.10.0 R.methodsS3_1.8.2 [41] hms_1.1.3 coda_0.19-4.1 rlang_1.1.2 Rcpp_1.0.11 xtable_1.8-4 [46] glue_1.6.2 rstudioapi_0.15.0 R6_2.5.1 ```

[vincentarelbundock]

Hi @mikedenly ,

Glad you find the packages useful!

With the info you give, we can't rule out the possibility that the difference has nothing at all to do with marginaleffects. To make sure this is a marginaleffects issue, you would first need to ensure that standard errors around the simple coefficients are exactly identiical in R and Stata. That doesn't seem to be the case here:

Stata:

beertax |  -.3472736   .1728049    -2.01 

         Estimate Std. Error  z value
beertax -0.347274   0.174639 -1.98852

A small change in the variance-covariance matrix for the coefficients can have a big effect on the marginal effects standard errors, so it would be really important to match the basic estimates first.

I don't have access to Stata and I don't know what ppmlhdfe is, so I can't really help figure out how it computes standard errors.

[mikedenly]

Thank you so much for your fast and helpful response, Vincent! @lrberge 's fixest::fepois() replicates ppmlhdfe in Stata (see confirmation here), so I perhaps should have directed the question to him given what you point out. ppmlhdfe in Stata was put together by Sergio Correia (@sergiocorreia), Paulo Guimares (@pguimaraes99), and Tom Zylkin (@tomzylkin) (see their Stata Journal Article).

[vincentarelbundock]

I'll be closing this now to keep the repo clean, because it does not appear (at first glance) to be a marginaleffects bug or feature request.

If/when it becomes clear that this is a marginaleffects problem, we can re-open the issue.

In any case, I am interested in the resolution, so feel free to keep discussing here if you want.

[mikedenly]

Hi Vincent (@vincentarelbundock),

As you can see from this other open issue that I started, both Laurent Bergé (@lrberge), the maintainer of fixest::fepois(), and Alexander Fischer (@s3alfisc), the maintainer of pyfixest, very kindly provided their expert advice on the standard error difference issue. I thus wanted to kindly ask if your conclusion remains the same regarding the massive differences in resulting marginal effect confidence intervals between the Stata and R after fixest::fepois()? On my side, I ended up submitting my paper to the journal with exponentiated coefficients, even though your work has clearly shown that they are less interpretable, because (i) exponentiated coefficients are consistent across software programs; and (ii) it was hard for me as a non-expert to take a stand one way or the other.

In any case, thank you so much again for your incredible work. It has really made a massive impact for the better on the social sciences and beyond.

[vincentarelbundock]

Hi @mikedenly ,

Thanks for the kind words!

I must admit that I don't have a super strong intuition about what's going on here. If you want to diagnose, I'd probably start with:

1. Exporting the variance-covariance matrices produced in Stata and Python (with sufficient precision), and feeding those to the vcov argument of the avg_slopes() function.
2. Play around with the step size in the numeric differentiation. See: https://marginaleffects.com/vignettes/uncertainty.html#numerical-derivatives-sensitivity-to-step-size

I don't know if any of those paths would pan out, but that's where I would personally start.