Three way fixed effects vs combining two of the effects

I have panel data on employment that varies by year, sector and location and thus would like to run a fixed effects regression considering these 3 dimensions.

The issue is that I use R and the plm package and it only allows me to have two fixed effects. If I use them and add dummy variables for the third dimension, I get an error since each observation is not uniquely determined by the two fixed effects.

As a solution, I have tried combining my location and sector dimensions into a single one, and then running twoways fixed effects.

Is this right? I am suspicious about that, since I also ran the same specification for each of my 13 sectors using year and location fixed effects. The results are really different: in the aggregate regression, the coefficient is largely positive whereas by sector they are are in general much smaller.

Am I doing right? Can someone help on using threeway fixed effects in R?

• Any particular reason you do not accept the answer given? Commented Mar 11, 2020 at 8:07
• Sorry! It was an awesome answer. I just took some time to test it. Commented Mar 11, 2020 at 8:09
• @StopClosingQuestionsFast It would be nice to get an answer on this follow-up: stats.stackexchange.com/questions/453585/… Commented Mar 11, 2020 at 16:56
• Yes, ok I have looked at it but I cannot really provide an answer. Commented Mar 13, 2020 at 13:05

Consider the model $$(1) \ \ w_{it} = \mathbf x_{it}^\top \beta + \delta_t +\psi_{a(i,t)} + \eta_{k(i,t)} + \epsilon_{it},$$ with the area effect $$\psi_a$$ and sector effect $$\eta_k$$ unobserved. Assuming that $$\mathbf x_{it}$$ is correlated with the area and sector effect the OLS estimator $$\hat \beta_{OLS}:=(\sum_i \sum_t\mathbf x_{it}\mathbf x_{it}^\top)^{-1}(\sum_i \sum_t\mathbf x_{it}y_{it})$$ associated with the estimation equation $$w_{it} = \mathbf x_{it}^\top \beta + u_{it}$$ is inconsistent because $$\mathbb E[\mathbf x_{it}u_{it}]=\mathbb E[\mathbf x_{it}(\delta_t +\psi_{a(i,t)} + \eta_{k(i,t)} + \epsilon_{it})]\not=0$$. Doing the estimation with fixed effects for area $$\psi_a$$, sector $$\eta_k$$ and time $$\delta_t$$ will give you consistent estimates assuming $$\mathbb E[x_{it}\epsilon_{it}]$$.

But then again so will doing the estimation with area-sector fixed $$\phi_{ak}$$ hence using the estimation equation

$$(2) \ \ w_{it} = \mathbf x_{it}^\top \beta + \delta_t + \phi_{a(i,t),k(i,t)} + \epsilon_{it},$$

where the area-sector specific fixed effect is $$\phi_{a(i,t),k(i,t)}$$. This is perhaps most simply seen by just recognizing that the first model is an instance of the latter by the restriction that $$\phi_{a(i,t),k(i,t)} = \psi_{a(it)} + \eta_{k(i,t)},$$ however the two estimators are not the same and so estimates might differ. Also model (2) cannot alwys be estimated consistently using model (1) as estimation equation.

In R you should use the lfe-package by Simen Gaure and you can find documentation here. Remember to cluster standard errors on id for panel data.

Here is a simulation example (I leave it to you to figure out how to add the fixed effect for time):

library(lfe)
library(data.table)

# Simulate a population of N workers observed over T timeperiods.
# Balanced panel
# Workers are assigned to A different areas
# Workers are assigned to K different sectors
N <- 1000
T <- 10
NT <- T*N
A <- 30
K <- 10
vA <- 10 # strength of area effect
vK <- 10 # strength of sector effect

# Initialize vectors for area and sector assignment
area <- rep(NA,NT)
sector <- rep(NA,NT)

# Choose probabilities for assigning individual to sector and area
# Probabilities are increasing in index pA[j+1]>pA[j] this is used
# to assign certain individual according to observed skill x to certain
# sectors and areas.
# Economic literature suggest that workers sort over sectors and areas
# according to skill (see for example Glaeser and Mare (2001) Cities and Skills and
# Combes (2008) Spatial wage disparities: Sorting matters!)
pA <- (1:A)^4/sum((1:A)^4)
pK <- (1:K)^0.7/sum((1:K)^0.7)

# Check distribution
layout(matrix(1:2,nrow=1))
barplot(table(sample(1:A,size=1000,prob=pA,replace=TRUE)))
barplot(table(sample(K:1,size=1000,prob=pK,replace=TRUE)))

# Set unobserved individual parameter deciding individual tendency to sort
# Individuals with high mu[i] will be in high sector versus low sector
# and in high area versus low area
mu <- rnorm(N)

# Sart loop to assign individuals to sector and area
ii <- 1
for (i in 1:N)
{
a <- ifelse(mu[i] > 0,sample(1:A,size=1,prob=pA),sample(A:1,size=1,prob=pA))
k <- ifelse(mu[i] > 0,sample(1:K,size=1,prob=pK),sample(K:1,size=1,prob=pK))

# The above assigns individuals with high mu to high index sector and area
# because probabilities are increasing in index
for (t in 1:T)
{
# Individual worker do not frequently change sector and area
# here the probability of change is chosen to 0.2 (very high)
# probably more around 5-10% (but we need variation in panel)
if (runif(1)<0.2)
{
a <- ifelse(mu[i] > 0,sample(1:A,size=1,prob=pA),sample(A:1,size=1,prob=pA))
}

if (runif(1)<0.2)
{
k <- ifelse(mu[i] > 0,sample(1:K,size=1,prob=pK),sample(K:1,size=1,prob=pK))
}

# Assign and note that a and k have changed from last period with probability 0.2
area[ii] <- a
sector[ii] <- k
ii <- ii + 1
}
}

# Specify area and sector effect, vA and vK controls size of effect
# The are sorted so higher index sector is high wage sector and higher
# index area is high wage area (where to individuals of high mu sort)
area_effect <- sort(vA*runif(A))
sector_effect <- sort(vK*runif(K))

# Define id and time period for observation
id <- rep(1:N,each=T)
time <- rep(1:T,N)

# Make some covariate ... here made correlated with area and sector
# mu[i] is used as mean of individual i's time varying observed skill x
x <-  rnorm(NT,mean=rep(mu,each=T)) + area_effect[area] + sector_effect[sector]

# rnorm(NT,mean=rep(mu,each=T))  ... check strength of covariance
# high covariance implies larger bias in OLS estimates
cov(x,area_effect[area])
cov(x,sector_effect[sector])

# Make dependent variable using the Mincer wage equation
y <- beta*x + area_effect[area] + sector_effect[sector] + (rt(NT,10)+abs(rt(NT,7)))
dt <- data.table(id=id,time=time,y=y,x=x,area=area,sector=sector,as=interaction(area,sector))
setkey(dt,id,time)

# Start estimation, first OLS is seen to be inconsistent
lm(y~x,data=dt)

# Must control for sector and area using fixed effects
# Dummy estimators (break down on large number of fixed effects)
# Both consistent but not good with many fixed effects
# Also standard error is not clustered on id ...
model1 <- lm(y ~ x + as.factor(sector) + as.factor(area),data=dt)
model2 <- lm(y ~ x + as,data=dt)

# Use lfe (designed to handle many fixed effects)
# Cluster on id because it is panel (actually whether this is necessary depends on assumption about variance)
#                modelspec   fixed effect    instru     cluster sd
model3 <- felm(    y~x     |  sector+area |     0     |      id      , data=dt)
model4 <- felm(y~x|as|0|id,data=dt)

# check estimates ... all consistent
coef(model1)[2]
coef(model2)[2]
coef(model3)
coef(model4)


Good reads (armed with patience) on the topic:

John M. Abowd, Francis Kramarz and David N. Margolis Source: Econometrica, Vol. 67, No. 2 (Mar., 1999), pp. 251-333

Pierre-Philippe Combes & Laurent Gobillon (2015) The Empirics of Agglomeration Economies in Handbook of Regional and Urban Economics