Modern Estimators: fixest and alpaca
In stargazer2: Beautiful Regression Tables for Modern Econometrics Packages

knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "",
  message  = FALSE,
  warning  = FALSE
)
library(stargazer2)

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

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.
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")

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"))

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")

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")

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))

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.