Modern Estimators: fixest and alpaca

stargazer2 supports fixest and alpaca model objects natively, with two features that matter for modern applied work:

  1. Fixed effects are reported as indicator rows (“Yes / No”) at the bottom of the table rather than as coefficient rows. This follows the convention in the trade and IO literatures where high-dimensional FEs are controls, not objects of interest.

  2. Standard error types are auto-detected from the model object and reported in the table note. When different columns use different SE types, the note breaks them out by column group automatically.

Standard error interface

stargazer2 accepts SEs through three mechanisms, applied in order of precedence:

Priority Mechanism When to use
1 (highest) vcov = list(V1, V2, ...) Any vcov matrix — most flexible
2 se = list(se1, se2, ...) Numeric SE vectors (drop-in for original stargazer scripts)
3 (lowest) Auto-extraction fixest and alpaca models — SE type read from the model object

Passing NULL for a specific entry in the vcov list tells stargazer2 to fall back to auto-extraction for that column. This is useful when mixing lm models (with externally supplied vcov matrices) and fixest models (which already carry their SE type).

Fixed effects with fixest

library(wooldridge)
data(wage1)

wage1$region <- factor(
  ifelse(wage1$northcen == 1, "northcen",
  ifelse(wage1$south    == 1, "south",
  ifelse(wage1$west     == 1, "west", "northeast"))),
  levels = c("northeast", "northcen", "south", "west")
)

wage1$occupation <- factor(
  ifelse(wage1$profocc == 1, "professional",
  ifelse(wage1$clerocc == 1, "clerical",
  ifelse(wage1$servocc == 1, "service", "other"))),
  levels = c("other", "professional", "clerical", "service")
)

wage1$industry <- factor(
  ifelse(wage1$construc == 1, "construction",
  ifelse(wage1$ndurman  == 1, "nondurable_manuf",
  ifelse(wage1$trcommpu == 1, "transport",
  ifelse(wage1$trade    == 1, "trade",
  ifelse(wage1$services == 1, "services",
  ifelse(wage1$profserv == 1, "prof_services", "other")))))),
  levels = c("other", "construction", "nondurable_manuf",
             "transport", "trade", "services", "prof_services")
)

We estimate five log-wage regressions with feols, varying the set of fixed effects and using two-way clustering on region × industry throughout. The interacted FE specification region^industry absorbs a full set of region-by-industry cells.

library(fixest)

f1 <- feols(lwage ~ educ + exper + tenure + female + married |
              region,               wage1, vcov = ~region^industry)
f2 <- feols(lwage ~ educ + exper + tenure + female + married |
              occupation,           wage1, vcov = ~region^industry)
f3 <- feols(lwage ~ educ + exper + tenure + female + married |
              region + occupation,  wage1, vcov = ~region^industry)
f4 <- feols(lwage ~ educ + exper + tenure + female + married |
              region + occupation + industry,
                                    wage1, vcov = ~region^industry)
f5 <- feols(lwage ~ educ + exper + tenure + female + married |
              region^industry,      wage1, vcov = ~region^industry)

A bare call — no labels specified — already produces a complete table. stargazer2 extracts the dependent variable name (lwage) from the model formula and, since all five columns share the same estimator, omits the redundant model-type row:

stargazer(f1, f2, f3, f4, f5, type = "text")
======================================================================
                                    Dependent variable:               
                     -------------------------------------------------
                                           lwage                      
                        (1)       (2)       (3)       (4)       (5)   
----------------------------------------------------------------------
educ                 0.083***  0.061***  0.060***  0.057***  0.081*** 
                      (0.007)   (0.007)   (0.007)   (0.008)   (0.008) 
                                                                      
exper                 0.003*     0.002     0.002     0.002    0.004** 
                      (0.002)   (0.002)   (0.002)   (0.001)   (0.002) 
                                                                      
tenure               0.017***  0.016***  0.016***  0.014***  0.014*** 
                      (0.004)   (0.004)   (0.004)   (0.003)   (0.004) 
                                                                      
female               -0.291*** -0.256*** -0.263*** -0.264*** -0.285***
                      (0.039)   (0.042)   (0.043)   (0.045)   (0.039) 
                                                                      
married              0.131***  0.115***  0.120***   0.094**  0.097*** 
                      (0.034)   (0.036)   (0.037)   (0.035)   (0.033) 
                                                                      
----------------------------------------------------------------------
Region FE               Yes       No        Yes       Yes       No    
Occupation FE           No        Yes       Yes       Yes       No    
Industry FE             No        No        No        Yes       No    
Region x Industry FE    No        No        No        No        Yes   
----------------------------------------------------------------------
Observations            526       526       526       526       526   
R2                     0.414     0.450     0.460     0.510     0.473  
Adjusted R2            0.405     0.442     0.449     0.494     0.439  
Within R2              0.407     0.282     0.287     0.267     0.378  
Adjusted Within R2     0.401     0.275     0.280     0.260     0.371  
======================================================================
Note:                standard errors clustered by region-industry;
                      *p<0.1; **p<0.05; ***p<0.01 

Labels can be overridden when needed. The LaTeX source below also renames the covariates for presentation:

stargazer(f1, f2, f3, f4, f5,
          type             = "latex",
          title            = "Log Wages: Varying Fixed Effects",
          label            = "tab:fe-wages",
          dep.var.labels   = "log(Wage)",
          covariate.labels = c("Education", "Experience", "Tenure",
                               "Female", "Married"))
% Table produced by stargazer2 v.0.1.0 by Tom Zylkin, University of Richmond (tzylkin@richmond.edu)
% Original stargazer package by: Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
%   R package version 5.2.3. https://CRAN.R-project.org/package=stargazer

\begin{table}[!htbp] \centering 
  \caption{Log Wages: Varying Fixed Effects} 
  \label{tab:fe-wages} 
\begin{tabular}{@{\extracolsep{5pt}}lccccc} 
\hline 
 & \multicolumn{5}{c}{\textit{Dependent variable:}} \\ 
 & \multicolumn{5}{c}{log(Wage)} \\ 
 & (1) & (2) & (3) & (4) & (5)\\ 
\hline 
 Education & 0.083$^{***}$ & 0.061$^{***}$ & 0.060$^{***}$ & 0.057$^{***}$ & 0.081$^{***}$ \\ 
  & (0.007) & (0.007) & (0.007) & (0.008) & (0.008) \\ 
 Experience & 0.003$^{*}$ & 0.002 & 0.002 & 0.002 & 0.004$^{**}$ \\ 
  & (0.002) & (0.002) & (0.002) & (0.001) & (0.002) \\ 
 Tenure & 0.017$^{***}$ & 0.016$^{***}$ & 0.016$^{***}$ & 0.014$^{***}$ & 0.014$^{***}$ \\ 
  & (0.004) & (0.004) & (0.004) & (0.003) & (0.004) \\ 
 Female & $-$0.291$^{***}$ & $-$0.256$^{***}$ & $-$0.263$^{***}$ & $-$0.264$^{***}$ & $-$0.285$^{***}$ \\ 
  & (0.039) & (0.042) & (0.043) & (0.045) & (0.039) \\ 
 Married & 0.131$^{***}$ & 0.115$^{***}$ & 0.120$^{***}$ & 0.094$^{**}$ & 0.097$^{***}$ \\ 
  & (0.034) & (0.036) & (0.037) & (0.035) & (0.033) \\ 
\hline 
Region FE & Yes & No & Yes & Yes & No \\ 
Occupation FE & No & Yes & Yes & Yes & No \\ 
Industry FE & No & No & No & Yes & No \\ 
Region x Industry FE & No & No & No & No & Yes \\ 
\hline 
Observations & 526 & 526 & 526 & 526 & 526 \\ 
R$^{2}$ & 0.414 & 0.450 & 0.460 & 0.510 & 0.473 \\ 
Adjusted R$^{2}$ & 0.405 & 0.442 & 0.449 & 0.494 & 0.439 \\ 
Within R$^{2}$ & 0.407 & 0.282 & 0.287 & 0.267 & 0.378 \\ 
Adjusted Within R$^{2}$ & 0.401 & 0.275 & 0.280 & 0.260 & 0.371 \\ 
\hline 
\hline 
\multicolumn{6}{p{\linewidth}}{\textit{Note:} standard errors clustered by region-industry; $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.} \\ 
\end{tabular} 
\end{table}  

Two things to notice in both outputs:

  • FE indicator rows appear below the coefficient block — one row per unique fixed effect across all columns, with “Yes” or “No” per model. The interacted specification region^industry is rendered as a single “Region x Industry FE” row.
  • SE type is auto-detected from the fixest model and reported uniformly in the note, since all five models use the same clustering strategy.

Multiple estimators and clustering: the gravity model

The fixest package’s built-in trade dataset provides a natural setting for showcasing different estimators. We estimate a standard gravity equation for trade flows using OLS, Poisson pseudo-maximum likelihood (PPML), and negative binomial — all with the same four-way fixed effects.

data(trade, package = "fixest")

gravity_ols    <- feols(log(Euros) ~ log(dist_km) |
                          Origin + Destination + Product + Year, trade)
gravity_pois   <- fepois(Euros ~ log(dist_km) |
                           Origin + Destination + Product + Year, trade)
gravity_negbin <- fenegbin(Euros ~ log(dist_km) |
                             Origin + Destination + Product + Year, trade)

Two further Poisson columns demonstrate the clustering-label machinery: ~Origin^Destination clusters by the interaction (one cluster per origin-destination pair); ~Origin+Destination is two-way clustering by origin and by destination separately.

gravity_pois1  <- fepois(Euros ~ log(dist_km) |
                           Origin + Destination + Product + Year,
                         trade, vcov = ~Origin^Destination)
gravity_pois2  <- fepois(Euros ~ log(dist_km) |
                           Origin + Destination + Product + Year,
                         trade, vcov = ~Origin + Destination)

With no arguments beyond the model list, stargazer2 automatically extracts the dependent variable names from each formula (log(Euros) for the OLS model, Euros for the count models), detects the estimator type for each column (OLS, Poisson, Neg. Binomial), and reads the SE method from each model object:

stargazer(gravity_ols, gravity_pois, gravity_negbin,
          gravity_pois1, gravity_pois2,
          type = "text")
=========================================================================
                                    Dependent variable:                  
                   ------------------------------------------------------
                   log(Euros)   Euros       Euros       Euros     Euros  
                      OLS      Poisson  Neg. Binomial  Poisson   Poisson 
                      (1)        (2)         (3)         (4)       (5)   
-------------------------------------------------------------------------
log(dist_km)       -2.170***  -1.528***   -1.711***   -1.528*** -1.528***
                    (0.021)    (0.000)     (0.017)     (0.077)   (0.131) 
                                                                         
-------------------------------------------------------------------------
Origin FE             Yes        Yes         Yes         Yes       Yes   
Destination FE        Yes        Yes         Yes         Yes       Yes   
Product FE            Yes        Yes         Yes         Yes       Yes   
Year FE               Yes        Yes         Yes         Yes       Yes   
-------------------------------------------------------------------------
Observations         38,325    38,325      38,325      38,325    38,325  
R2                   0.706      0.612                   0.612     0.612  
Adjusted R2          0.705                                               
Within R2            0.219      0.307       0.016       0.307     0.307  
Adjusted Within R2   0.219                                               
Theta                                       0.549                        
=========================================================================
Note:              (1) OLS standard errors; (2)-(3) MLE standard errors;
                    (4) standard errors clustered by Origin-Destination;
                    (5) standard errors clustered by Origin and
                    Destination; *p<0.1; **p<0.05; ***p<0.01 

The same call in LaTeX produces submission-ready output:

stargazer(gravity_ols, gravity_pois, gravity_negbin,
          gravity_pois1, gravity_pois2,
          type  = "latex",
          title = "Gravity Equation for Trade Flows",
          label = "tab:gravity")
% Table produced by stargazer2 v.0.1.0 by Tom Zylkin, University of Richmond (tzylkin@richmond.edu)
% Original stargazer package by: Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
%   R package version 5.2.3. https://CRAN.R-project.org/package=stargazer

\begin{table}[!htbp] \centering 
  \caption{Gravity Equation for Trade Flows} 
  \label{tab:gravity} 
\begin{tabular}{@{\extracolsep{5pt}}lccccc} 
\hline 
 & \multicolumn{5}{c}{\textit{Dependent variable:}} \\ 
 & log(Euros) & Euros & Euros & Euros & Euros \\ 
 & \textit{OLS} & \textit{Poisson} & \textit{Neg. Binomial} & \textit{Poisson} & \textit{Poisson} \\ 
 & (1) & (2) & (3) & (4) & (5)\\ 
\hline 
 log(dist\_km) & $-$2.170$^{***}$ & $-$1.528$^{***}$ & $-$1.711$^{***}$ & $-$1.528$^{***}$ & $-$1.528$^{***}$ \\ 
  & (0.021) & (0.000) & (0.017) & (0.077) & (0.131) \\ 
\hline 
Origin FE & Yes & Yes & Yes & Yes & Yes \\ 
Destination FE & Yes & Yes & Yes & Yes & Yes \\ 
Product FE & Yes & Yes & Yes & Yes & Yes \\ 
Year FE & Yes & Yes & Yes & Yes & Yes \\ 
\hline 
Observations & 38,325 & 38,325 & 38,325 & 38,325 & 38,325 \\ 
R$^{2}$ & 0.706 & 0.612 &  & 0.612 & 0.612 \\ 
Adjusted R$^{2}$ & 0.705 &  &  &  &  \\ 
Within R$^{2}$ & 0.219 & 0.307 & 0.016 & 0.307 & 0.307 \\ 
Adjusted Within R$^{2}$ & 0.219 &  &  &  &  \\ 
Theta &  &  & 0.549 &  &  \\ 
\hline 
\hline 
\multicolumn{6}{p{\linewidth}}{\textit{Note:} (1) OLS standard errors; (2)-(3) MLE standard errors; (4) standard errors clustered by Origin-Destination; (5) standard errors clustered by Origin and Destination; $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.} \\ 
\end{tabular} 
\end{table}  

The SE note shows per-column types: OLS standard errors for column (1), MLE standard errors for (2)–(3) (auto-detected from the fixest objects), standard errors clustered by Origin-Destination for column (4), and two-way clustered by Origin and Destination for column (5).

Non-linear models with alpaca

The alpaca package offers an alternative implementation of fixed-effects GLMs (logit, probit, Poisson). stargazer2 supports alpaca::feglm objects through a companion pair of vcov helpers:

  • alpaca_vcovSandwich() — heteroskedasticity-robust (sandwich) SEs
  • alpaca_vcovCL() — clustered SEs, with a formula interface matching sandwich::vcovCL
library(alpaca)

# Logit model: P(married) as a function of wages and human capital,
# with occupation and industry fixed effects.
# industry must be in the FE specification for clustering by industry.
m_alp <- feglm(married ~ lwage + educ + exper | occupation + industry,
               wage1, binomial("logit"))

V_robust    <- alpaca_vcovSandwich(m_alp)
V_clustered <- alpaca_vcovCL(m_alp, cluster = ~industry)

stargazer(m_alp, m_alp,
          type             = "text",
          dep.var.labels   = "Married (0/1)",
          covariate.labels = c("log(Wage)", "Education", "Experience"),
          column.labels    = c("Sandwich-robust", "Industry-clustered"),
          vcov             = list(V_robust, V_clustered))
=========================================
                  Dependent variable:    
              ---------------------------
                     Married (0/1)       
              Sandwich-robust Industry-clustered
                   (1)           (2)     
-----------------------------------------
log(Wage)       0.890***      0.890***   
                 (0.263)       (0.344)   
                                         
Education         0.078       0.078***   
                 (0.047)       (0.026)   
                                         
Experience      0.058***      0.058***   
                 (0.010)       (0.007)   
                                         
-----------------------------------------
Occupation FE      Yes           Yes     
Industry FE        Yes           Yes     
-----------------------------------------
Observations       526           526     
R2                0.190         0.190    
=========================================
Note:         (1)
               heteroskedasticity-robust
               standard errors; (2)
               standard errors clustered
               by industry; *p<0.1;
               **p<0.05; ***p<0.01 

The SE note names the type for each column exactly as it does for fixest and sandwich models, confirming that the reporting machinery is consistent across all supported packages.