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?
 A: 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)
    {
        # 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
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
