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?


1 Answer 1


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


# 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

# 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)
        # Assign individual to sector
        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

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

# Start estimation, first OLS is seen to be inconsistent

# 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

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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.