three-level logit regression with random slopes

Question

I am trying to find the correct code for a three-level multivariate logit regression. I have protest participation, education level, and income data from a number of countries. Data are also divided by regions within each country (each region across all countries has a numeric code different from all others so they can be aggregated without risk of overlapping labels). Protest participation is my binary dependent variable (have you participated in the last year? Yes-No) while income (divided in quintiles) and education level (7 levels from illiteracy to Phd) are my independent variables.

I want to find the link between participation, income and education at the aggregated level (all countries together) and I know that slopes of income and education change significantly both at the country and regional levels. So I guess I need to run a 3-level regression (citizens-regions-countries) using income and education both as fixed and random effects. However, since income is already divided in quintiles (based on the national income distribution), I guess that income should only be nested within regions and not countries.

I wrote the following R code and I really cannot understand if the results I obtain are correct or if I am doing something wrong:

model <- glmer(participation ~ Education + income + (Education|region/country) + (income|region), data = NorthAfrica, family = 'binomial')

Thanks for your help!

Answer 1

When you include a term like (1|country/region) in your model, that is equivalent with including the following two terms: (1|country) and (1|region:country). See the Cross Validated post Crossed vs nested random effects: how do they differ and how are they specified correctly in lme4? for details.

This means that you can replace a model formulation like:

model1a <- glmer(participation ~ Education + income + (1|country/region), 
               data = NorthAfrica, 
               family = 'binomial')

with

model1b <- glmer(participation ~ Education + income + (1|country) + (1|region:country),
               data = NorthAfrica, 
               family = 'binomial')

The two model formulations are equivalent and allow for a random intercept for country and a random intercept for region nested within country. The second formulation is however more flexible as it allows you to include random slopes for predictor variables at the appropriate level of your data hierarchy. For instance, if you would like the effect of education to vary at the region level, you would use a model formulation like this:

model2 <- glmer(participation ~ Education + income + (1|country) + (1 + Education|region:country),
               data = NorthAfrica, 
               family = 'binomial')

This model formulation assumes that each study subject within a region comes with his/her own value for Education. If each study subject also came with his/her own income value, you could expand the above model to also include a random slope for income:

model3 <- glmer(participation ~ Education + income + (1|country) + (1 + Education + income|region:country),
               data = NorthAfrica, 
               family = 'binomial')

To sum up, your data hierarchy includes 3 levels:

countries

regions

subjects 

The region and the country are considered as random grouping factors, with region nested within country. To reflect this data hierarchy, your model should include, at a minimum, random intercepts for country and region nested within country:

(1|country) + (1|region:country)

If you have a predictor variable measured for each subject (e.g., income), you can include a random slope for it in your model like so:

(1|country) + (1 + income|region:country)

This is telling the model that the effect of that predictor (e.g., income) varies across regions (nested within countries).

If you have a predictor variable measured for each region (say, proportion of people with graduate and post-graduate education in the region), then you can include a random slope for that predictor in your model as follows:

 (1 + proportion|country) + (1|region:country)

This is telling the model that the effect of the proportion predictor varies across countries.

If your predictor is measured for the entire country (e.g., proportion of people with graduate and post-graduate degrees in the country), then you cannot include a random slope for it in your model - to do that, you would need an additional level for your data hierarchy determined by a random grouping factor (e.g., continent).