Estimation of level two coefficients despite the inclusion of fixed effects

[Felix Kaysers]

Dear Statalist,

I have a basic setup of a cross-sectional sample of individuals nested in regions. I am interested in estimating how, for example, the number of natural disasters in region j is associated with the likelihood of voting against the incumbents (binary outcome).

Problem: In my main model, I had naively run a simple logistic regression including my level two predictor , i.e. the number of natural disasters, and regional dummies, i.e. fixed effects, and Stata presented coefficient estimates for both. I have now found out that the coefficients of my level two variable should be non-estimable when including the fixed effects for the region (Schunk, 2017, p. 8 and McNeish, 2023, p. 8).

Questions: I would like to know if a) I misunderstood something and b) if my understanding is correct, what problem the regional dummies introduce? Since Stata provices estimates of the coefficients, the claim that the effect of level two variables are non-estimable seems wrong. If this were the case, I would expect a . as estimate for my level two coefficient. Is it simply that the estimate of level two predictors, number of natural disasters, is heavily biased with fixed effects? In my original regression, the coefficients for my level two predictor vary substantially between the regression withou (beta = 0.88) and without a fixed effect dummies (beta = 0.07).

Below is an MRE recreating a data-generating process similar to my question. I use the community-contributed esttab for saving the regressions, but I am unsure about the construction of y:

clear all
set seed 56
set obs 1000
gen id = _N
gen rand = runiform()

* IVs
sort rand
gen region = ceil(_n/20)
bysort region: gen leveltwo = rpoisson(1) if _n == 1
bysort region: replace leveltwo = leveltwo[1]
gen age = runiformint(18,60)
gen sex = runiformint(0,1)

* DVs
gen linpred = -0.4+2*leveltwo+0.05*age+1.5*sex
gen p = invlogit(linpred)
gen y = (runiform() < p)

** Regressions
* Run a logistic regression
eststo m1: logistic y leveltwo age sex

* Run a logistic regression with fixed effects
eststo m2: logistic y leveltwo age sex i.region

* Run a mulitlevel model
eststo m3: melogit y leveltwo age sex || region:

* Run a multilevel model with region fixed effects
eststo m4: melogit y leveltwo age sex i.region || region:

** Reporting
esttab m1 m2 m3 m4, eform

References
McNeish, D. (2023). A practical guide to selecting and blending approaches for clustered data: Clustered errors, multilevel models, and fixed-effect models. Psychological Methods. https://doi.org/10.1037/met0000620
Schunck, R., & Perales, F. (2017). Within- and Between-cluster Effects in Generalized Linear Mixed Models: A Discussion of Approaches and the Xthybrid command. The Stata Journal: Promoting Communications on Statistics and Stata, 17(1), 89–115. https://doi.org/10.1177/1536867X1701700106

[Jeff Wooldridge]

There are some problems with the way you generated the data because in several regions you have no variation in y. That can happen with relatively few observations per region in practice, and it happens in your data. But then those regions contribute nothing to the estimation. And I because of the incidental parameters problem, you should only include region dummies in a logit if you have fairly large sample sizes per region.

To see what's going on, you should generate data so that there is some variation in y within each region. Then, you'll see that if leveltwo is kept in the equation, Stata will have dropped one of the region dummies it should've kept. There is no way you can estimate the leveltwo coefficient if you include a full set of region dummies. When Stata sees this situation of perfect collinearity, it has rules about the variable it drops. Often it is to drop the last in the list because that is when collinearity is first detected. But the bottom line is that leveltwo should be the variable that gets dropped, but it often doesn't work out that way. Put it last and see what happens. But only after you solve the other problem of no variation in y within some regions.

[Felix Kaysers]

Thank you for the detailed comments, Prof. Wooldridge! My apologies for the faulty data-generating process (first time). I updated the code by decreasing the number of regions and manually substracted 0.25 when generating y to generate some variation in y.

gen region = ceil(_n/50) // Line 10
gen y = (runiform() < p-0.25) // Line 20

If I understand your comments correctly, the model is non-estimable due to the perfect collinearity. I am convinced (and strongly hope) that this is sufficient to justify not including governorate fixed effects in the model.