# Fixed Effects in Linear Regression

Fixed effects is a statistical regression model in which the intercept of the regression model is allowed to vary freely across individuals or groups. It is often applied to panel data in order to control for any individual-specific attributes that do not vary across time.

## Keep in Mind

• To use individual-level fixed effects, you must observe the same person multiple times (panel data).
• In a linear regression context, fixed effects regression is relatively straightforward, and can be thought of as effectively adding a binary control variable for each individual, or subtracting the within-individual mean of each variable (the “within” estimator). However, you may want to apply fixed effects to other models like logit or probit. This is usually possible (depending on the model), but if you just add a set of binary controls or subtract within-individual means, it won’t work very well. Instead, look for a command specifically designed to implement fixed effects for that model.
• If you are using fixed effects to estimate the causal effect of a variable $$X$$, individuals with more variance in $$X$$ will be weighted more heavily (Gibbons, Serrano, & Urbancic 2019, ungated copy here). You may want to consider weighting your regression by the inverse within-individual variance of $$X$$.

## Also Consider

• Instead of fixed effects you may want to use random effects, which requires additional assumptions but is statistically more efficient and also allows the individual effect to be modeled using covariates. See Linear Mixed-Effects Regression
• You may want to consider clustering your standard errors at the same level as (some or more of) your fixed effects.

# Implementations

## Julia

Julia provides support for estimating high-dimensional fixed effect models through the FixedEffectModels.jl package (link). Similarly to felm (R) and reghdfe (Stata), the package uses the method of alternating projections to sweep out fixed effects. However, the Julia implementation is typically quite a bit faster than these other two methods. It also offers further performance gains via GPU computation for users with a working CUDA installation (up to an order of magnitude faster for complicated problems).

# If necessary, install JuliaFixedEffects.jl and some ancilliary packages for reading in the data
# ] add JuliaFixedEffects, CSVFiles, DataFrames

# Read in the example CSV and convert to a data frame
using CSVFiles, DataFrames
# Calculate proportion of graduates working
df[!, :prop_working] = df[!, :count_working] ./ (df[!, :count_working ] .+ df[!, :count_not_working])

using JuliaFixedEffects

# Regress median earnings on the proportion of working graduates.
# We'll control for institution name and year as our fixed effects.
# We'll also cluster our standard errors by institution name.
reg(df, @formula(earnings_med ~ prop_working + fe(inst_name) + fe(year)), Vcov.cluster(:inst_name))

reg(df, @formula(earnings_med ~ prop_working + fe(inst_name) + fe(year)), Vcov.cluster(:inst_name), method = :lsmr_threads)

# GPU example (requires working CUDA installation)
reg(df, @formula(earnings_med ~ prop_working + fe(inst_name) + fe(year)), Vcov.cluster(:inst_name), method = :lsmr_gpu)


## Python

There are a few packages for doing the same task in Python, however, there is a well-known issue with these packages.That is, the calculation of standard deviation might be a little different.

We are going to use linearmodels in python. Installation can be done through pip install linearmodels and the documentation is here

# Import the packages
import pandas as pd
from linearmodels import PanelOLS
import numpy as np

# Set the index for fixed effects
data = data.set_index(['inst_name', 'year'])

# Calculate and drop the NA Values
data['prop_working'] = data['count_working']/(data['count_working'] + data['count_not_working'])
#data = data.dropna(subset=['earnings_med', 'prop_working'])

# Regression
FE = PanelOLS(data.earnings_med, data['prop_working'],
entity_effects = True,
time_effects=True
)

# Result
result = FE.fit(cov_type = 'clustered',
cluster_entity=True,
# cluster_time=True
)


There are also other packages for fixed effect models, such as econtools (link), FixedEffectModelPyHDFE (link), regpyhdfe(link) and econtools (link).

## R

There are numerous packages for estimating fixed effect models in R. We will limit our examples here to the two fastest implementations — lfe::felm and fixest::feols — both of which support high-dimensional fixed effects and standard error correction (multiway clustering, etc.).

We first demonstrate fixed effects in R using felm from the lfe package (link). lfe::felm uses the Method of Alternating Projections to “sweep out” the fixed effects and avoid estimating them directly. By default, this is automatically done in parallel, using all available cores on a user’s machine to maximize performance. (It is also possible to change this behaviour.)

library(lfe)

# Read in data from the College Scorecard
# Calculate proportion of graduates working
df$prop_working <- df$count_working/(df$count_working + df$count_not_working)

# A felm formula is constructed as:
# outcome ~
#   covariates |
#   fixed effects |
#   instrumental variables specification |
#   cluster variables for standard errors

# Here let's regress earnings_med on prop_working
# with institution name and year as our fixed effects
# And clusters for institution name
felm_model <- felm(earnings_med ~ prop_working | inst_name + year | 0 | inst_name, data = df)

# Look at our results
summary(felm_model)


Next, we consider feols from the fixest package (link). The syntax is very similar to lfe::felm and again the estimation will be done in parallel by default. However, rather than the method of alternating projections, fixest::feols uses a concentrated maximum likelihood method to efficiently estimate models with an arbitrary number of fixed effects. Current benchmarks suggest that this can yield significant speed gains, especially for large problems. For the below example, we’ll continue with the same College Scorecard dataset already loaded into memory.

# If necessary, install fixest
# install.packages('fixest')
library(fixest)

# Run the same regression as before
feols_model <- feols(earnings_med ~ prop_working | inst_name + year, data = df)

# Look at our results
# Standard errors are automatically clustered at the inst_name level
summary(feols_model)
# It is also possible to specify additional or different clustering of errors
summary(feols_model, se = "twoway")
summary(feols_model, cluster = c("inst_name", "year")) ## same as the above


As noted above, there are numerous other ways to implement fixed effect models in R. Users may also wish to look at the plm, lme4, and estimatr packages among others. For example, the latter’s estimatr::lm_robust function provides syntax that may be more familar syntax to new R users who are coming over from Stata. Note, however, that it will be less efficient for complicated models.

## Stata

We will estimate fixed effects using Stata in two ways. First, using the built in xtreg command. Second, using the reghdfe package (link), which is more efficient and better handles multiple levels of fixed effects (as well as multiway clustering), but must be downloaded from SSC first.

* Load in College Scorecard data
import delimited "https://raw.githubusercontent.com/LOST-STATS/LOST-STATS.github.io/master/Model_Estimation/Data/Fixed_Effects_in_Linear_Regression/Scorecard.csv", clear

* The missings are written as "NA", let's turn this numeric
destring count_not_working count_working earnings_med, replace force
* Calculate the proportion working
g prop_working = count_working/(count_working + count_not_working)

* xtset requires that the individual identifier be a numeric variable
encode inst_name, g(name_number)

* Set the data as panel data with xtset
xtset name_number

* Use xtreg with the "fe" option to run fixed effects
* Regressing earnings_med on prop_working
* with fixed effects for name_number (implied by fe)
* and also year (which we'll add manually with i.year)
* and standard errors clustered by name_number
xtreg earnings_med prop_working i.year, fe vce(cluster name_number)

* Now, let's demonstrate the same regression with reghdfe.
* Install the package first if necessary.
* ssc install reghdfe

* For reghdfe we don't need to xtset the data. Let's undo that
xtset, clear

* We specify both sets of fixed effects in absorb()
reghdfe earnings_med prop_working, absorb(name_number year) vce(cluster inst_name)