Using fixed effects predictors with same cluster structure as random intercept in multilevel model

Question

Say, I have following data from different countries, with following variables:

• for_xs - should people, not being from that country, have access to publicly funded health care? (yes/no, my dependent variable)

• health_expend - expenditure in each country for health care

• forgo - do people not use health services because they cannot afford

• income, age - equivalent income and age of respondent

• country - country where respondent lives

The variables health_expend and forgo have just one value within each country, because they represent the mean amount of money spend for health care (health_expend) and the proportion of people in that country that forgo care. So, these variables do not vary for observations within a country.

First, here is the sample data:

structure(list(country = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 
5L, 5L, 5L), .Label = c("AU-Australia", "BE-Belgium", "BG-Bulgaria", 
"CL-Chile", "CN-China", "TW-Taiwan", "HR-Croatia", "CZ-Czech Republic", 
"DK-Denmark", "FI-Finland", "FR-France", "DE-Germany", "IL-Israel", 
"IT-Italy", "JP-Japan", "KR-Korea (South)", "LT-Lithuania", "NL-Netherlands", 
"NO-Norway", "PH-Philippines", "PL-Poland", "PT-Portugal", "RU-Russia", 
"SK-Slovak Republic", "SI-Slovenia", "ZA-South Africa", "ES-Spain", 
"SE-Sweden", "CH-Switzerland", "TR-Turkey", "GB-Great Britain", 
"US-United States"), class = "factor"), for_xs = structure(c(1L, 
1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 
2L, 2L, 2L, 1L, 1L, 1L, 2L, 1L), .Label = c("Disagree", "Agree"
), class = "factor"), health_expend = c(3970, 3970, 3970, 3970, 
3970, 4079, 4079, 4079, 4079, 4079, 1080, 1080, 1080, 1080, 1080, 
1478, 1478, 1478, 1478, 1478, 0, 0, 0, 0, 0), forgo = c(9.2, 
9.2, 9.2, 9.2, 9.2, 11.5, 11.5, 11.5, 11.5, 11.5, 8.4, 8.4, 8.4, 
8.4, 8.4, 17.7, 17.7, 17.7, 17.7, 17.7, 17.7, 17.7, 17.7, 17.7, 
17.7), income = structure(c(NA, 1L, 2L, 2L, NA, 3L, 1L, 3L, 3L, 
3L, 3L, NA, 1L, 3L, 1L, 3L, 3L, NA, NA, 2L, 2L, 2L, 1L, 3L, 3L
), .Label = c("low", "medium", "high"), class = "factor"), AGE = c(52, 
66, 34, 69, 64, 91, 48, 33, 67, 26, 56, 19, 44, 18, 60, 70, 31, 
48, 48, 44, 35, 62, 59, 70, 54)), .Names = c("country", "for_xs", 
"health_expend", "forgo", "income", "AGE"), row.names = c(NA, 
-25L), class = c("tbl_df", "tbl", "data.frame"))

This is how the data looks like:

# A tibble: 25 × 6
        country   for_xs health_expend forgo income   AGE
         <fctr>   <fctr>         <dbl> <dbl> <fctr> <dbl>
1  AU-Australia Disagree          3970   9.2     NA    52
2  AU-Australia Disagree          3970   9.2    low    66
3  AU-Australia Disagree          3970   9.2 medium    34
4  AU-Australia Disagree          3970   9.2 medium    69
5  AU-Australia Disagree          3970   9.2     NA    64
6    BE-Belgium    Agree          4079  11.5   high    91
7    BE-Belgium Disagree          4079  11.5    low    48
8    BE-Belgium    Agree          4079  11.5   high    33
9    BE-Belgium    Agree          4079  11.5   high    67
10   BE-Belgium    Agree          4079  11.5   high    26
# ... with 15 more rows

When I compute logistic regressions for each country subset, estimates for these variables can't be computed, of course:

glm(for_xs ~ health_expend + forgo + income + AGE, 
    data = filter(tmp, country ==  "BE-Belgium"), family = binomial)

## Call:  glm(formula = for_xs ~ health_expend + forgo + income + AGE, 
##     family = binomial, data = filter(tmp, country == "BE-Belgium"))
##
## Coefficients:
##   (Intercept)  health_expend          forgo     incomehigh            AGE  
##    -2.457e+01             NA             NA      4.913e+01     -1.360e-11  
##
## Degrees of Freedom: 4 Total (i.e. Null);  2 Residual
## Null Deviance:       5.004 
## Residual Deviance: 2.143e-10     AIC: 6

However, I want to run a multilevel model, with country as random intercept:

glmer(for_xs ~ health_expend + forgo + income + AGE + (1 | country), 
    data = tmp, family = binomial)

There are some warnings, but I have coefficients for all predictors.

My question is: Does it make sense to include health_expend and forgo in the multilevel-model, because they virtually are just "copies" of the country-variable? They have no variance within each country. Why would I include health_expend and forgo in a multilevel-model anyway?

My feeling is, I would not. However, some of my colleagues say, it's a useful information. This is true from a theoretical perspective, but from a statistical perspective I would argue that this only makes sense if these variables also vary within countries.

Answer 1

Perhaps to answer this we could consider a hypothetical situation. Suppose you have also got average age for each country as well as individual age. Does it make sense to include either or both of these in your model. I would argue that it does. The effect of average age tells you the effect of living in a country which has an older population while individual age tells you the effect of being old. You could even have an interaction between them: does the effect of being old depend on whether you live in a country of old people.

What would upset the apple cart would be trying to use country as a fixed effect and then adding country level covariates but that is not what you have done.