FixedEffectModels.jl

By default, the FixedEffectModel object is lightweight and does not store the original dataset or residuals. To access these, you must use the save option during estimation.

Options for save:

• :none (default): Saves nothing.
• :residuals: Saves residuals.
• :fe: Saves fixed effects.
• :all: Saves both residuals and fixed effects.

Accessing Data:

Once saved, use the following functions. Note that predict and residuals require the original DataFrame as a second argument:

• residuals(m, df)
• fe(m) (returns fixed effects)
• predict(m, df)

# To enable retrieval of residuals and fixed effects
m = reg(df, @formula(Sales ~ NDI + fe(State)), save = :all)

# Accessing them
res = residuals(m, df)
fixed_effects = fe(m)

For faster performance on large problems, you can use experimental GPU support by setting the method option.

NVIDIA GPUs (CUDA)

1. Run using CUDA before using FixedEffectModels.
2. Set method = :CUDA in reg().
3. (Optional) Use double_precision = true to use Float64 (defaults to Float32).

Apple Silicon GPUs (Metal)

1. Run using Metal before using FixedEffectModels.
2. Set method = :Metal in reg().
3. Note: Metal only supports Float32.

CPU

• method = :cpu (default).

# NVIDIA CUDA Example
using CUDA, FixedEffectModels, RDatasets
@assert CUDA.functional()
df = dataset("plm", "Cigar")
reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)), method = :CUDA)

# Apple Metal Example
using Metal, FixedEffectModels, RDatasets
@assert Metal.functional()
reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)), method = :Metal)

The package is registered in the Julia General registry. You can install it via the Julia REPL using the add command.

] add FixedEffectModels

To reproduce the performance benchmarks for the R package fixest, use the following code. This compares feols with various fixed effects and clustering options.

library(fixest)
N = 10000000
K = 100
df = data.frame(
  id1 =  as.factor(sample(N/K, N, replace = TRUE)),
  id2 =  as.factor(sample(K, N, replace = TRUE)),
  x1 =  runif(N),
  x2 =  runif(N)
)
df[, "y"] =  3 * df[, "x1"] + 2 * df[, "x2"] + sin(as.numeric(df[, "id1"])) + cos(as.numeric(df[, "id2"])) + runif(N)

system.time(feols(y ~ x1 + x2, df))
system.time(feols(y ~ x1 + x2|id1, df))
system.time(feols(y ~ x1 + x2|id1 + id2, df))
system.time(feols(y ~ x1 + x2, cluster = "id1", df))
system.time(feols(y ~ x1 + x2, cluster = c("id1", "id2"), df))

To reproduce the performance benchmarks for FixedEffectModels.jl, use the following Julia code. Note that the timings in the documentation are based on warm calls using Julia 1.12.7 on an Apple M4 (16GB RAM).

using DataFrames, CategoricalArrays, FixedEffectModels
N = 10_000_000
K = 100
id1 = rand(1:(N/K), N)
id2 = rand(1:K, N)
x1 =  randn(N)
x2 =  randn(N)
y= 3 .* x1 .+ 2 .* x2 .+ sin.(id1) .+ cos.(id2).^2 .+ randn(N)
df = DataFrame(id1 = categorical(id1), id2 = categorical(id2), x1 = x1, x2 = x2, y = y)

# Benchmark different model specifications
@time reg(df, @formula(y ~ x1 + x2))
@time reg(df, @formula(y ~ x1 + x2 + fe(id1)))
@time reg(df, @formula(y ~ x1 + x2 + fe(id1) + fe(id2)))
@time reg(df, @formula(y ~ x1 + x2), Vcov.cluster(:id1))
@time reg(df, @formula(y ~ x1 + x2), Vcov.cluster(:id1, :id2))

To reproduce the performance benchmarks for the R package lfe, use the following code. This compares felm with various fixed effect specifications.

library(lfe)
N = 10000000
K = 100
df = data.frame(
  id1 =  as.factor(sample(N/K, N, replace = TRUE)),
  id2 =  as.factor(sample(K, N, replace = TRUE)),
  x1 =  runif(N),
  x2 =  runif(N)
)
df[, "y"] =  3 * df[, "x1"] + 2 * df[, "x2"] + sin(as.numeric(df[, "id1"])) + cos(as.numeric(df[, "id2"])) + runif(N)

system.time(felm(y ~ x1 + x2, df))
system.time(felm(y ~ x1 + x2|id1, df))
system.time(felm(y ~ x1 + x2|id1 + id2, df))
system.time(felm(y ~ x1 + x2|0|0|id1, df))
system.time(felm(y ~ x1 + x2|0|0|id1 + id2, df))

To reproduce the performance benchmarks for Stata using reghdfe and ivreg2, use the following Stata commands.

clear all
local N = 10000000
local K = 100
set obs `N'
gen  id1 =  floor(runiform() * (`N'+1)/`K')
gen  id2 =  floor(runiform() * (`K'+1))
gen   x1 =  runiform()
gen   x2 =  runiform()
gen   y =  3 * x1 + 2 * x2 + sin(id1) + cos(id2) + runiform()
timer clear

set rmsg on
reg y x1 x2
reghdfe y x1 x2, a(id1)
reghdfe y x1 x2, a(id1 id2)
reg y x1 x2, cl(id1)
ivreg2 y x1 x2, cluster(id1 id2)

The object returned by reg() implements StatsAPI methods. Use these to extract model statistics and coefficients.

Coefficients and Covariance

• coef(m)
• vcov(m)
• confint(m)
• coefnames(m)
• responsename(m)

Model Statistics

• nobs(m): Number of observations.
• dof(m): Degrees of freedom.
• dof_residual(m): Residual degrees of freedom.
• r2(m): R-squared.
• adjr2(m): Adjusted R-squared.
• loglikelihood(m): Log-likelihood.
• coeftable(m): Coefficient table.
• formula(m): The model formula.

Other Statistics

• islinear(m)
• deviance(m)
• nulldeviance(m)
• rss(m): Residual sum of squares.
• mss(m): Model sum of squares.
• nullloglikelihood(m)

Standard Errors (Vcov)

Standard errors are specified using the Vcov prefix (from the Vcov.jl package):

• Vcov.robust(): Robust standard errors.
• Vcov.cluster(:State): Clustered standard errors by a single variable.
• Vcov.cluster(:State, :Year): Multi-way clustering.

Weights

Use the weights option to specify a column name for weights:

• weights = :Pop

Dummy Coding (Contrasts)

Use the contrasts option to treat columns as dummy variables via DummyCoding():

• contrasts = Dict(:Year => DummyCoding())
• contrasts = Dict(:Year => DummyCoding(base = 80))

# Example with dummy coding
reg(df, @formula(Sales ~ Price + Year); contrasts = Dict(:Year => DummyCoding(base = 80)))