I am trying to run a plm panel fixed effects model. The panel index runs at the county-year level. The challenge? I want to run fixed effects at the (higher order) state level. I can easily program this with fixest but I am unsure how to do it with plm.
Here is the sample data
df <- data.frame(year= rep(c(2000, 2001, 2003, 2000, 2002, 2004), each=10),
county= rep(c("c1", "c2", "c3", "c4", "c5", "c6", "c7", "c8", "c9", "c10"), times=6),
state= rep(c("s1", "s2", "s3", "s4", "s5"), each=2, times=6),
temperature= runif(60, min=40, max=80),
trade= runif(60, min=10, max=20),
policy.outcome= runif(60, min=0, max=100))
And this is the fixest model I would like to replicate with plm
m1.fixest <- feols(policy.outcome
~ temperature + trade|state+year,
panel.id= ~county+year,
data=df)
Thank you,
Frist note that your data set has duplicated index for year = 2000. feols does not complain about this but plm does. Assuming this is an error, I slightly changed your data generation to end up with unique index combinations.
With plm, if you want fixed effects other than the onces of the index dimensions, you will need to explicitly specify them as a regressor (LSDV approach). Here, you can estimate a time fixed effect model via effect = "time" to for the year-fixed effect and explicitly specify + factor(state) to account for state-fixed effects as the individual dimension is given by county.
df <- data.frame(year= rep(c(2000, 2001, 2003, 2005, 2002, 2004), each=10),
county= rep(c("c1", "c2", "c3", "c4", "c5", "c6", "c7", "c8", "c9", "c10"), times=6),
state= rep(c("s1", "s2", "s3", "s4", "s5"), each=2, times=6),
temperature= runif(60, min=40, max=80),
trade= runif(60, min=10, max=20),
policy.outcome= runif(60, min=0, max=100))
library(fixest)
m1.fixest <- feols(policy.outcome
~ temperature + trade|state+year,
panel.id= ~county+year,
data=df)
library(plm)
m1.plm <- plm(policy.outcome ~ temperature + trade + factor(state),
index = c("county", "year"), effect = "time",
data=df)