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
