I have a dataset with two levels: the city level and person level, at the city level i have the categorical variable "poverty level", in my research estimating the effect of the "poverty level" is of paramount importance.
I am not sure how to include the poverty level in my model.
Here is mini example of my data:
library(tidyverse)
library(lme4)
# dataset with the information of the cities
city_df <- data.frame(
city = LETTERS[1:15],
poverty_level = rep(c("C1", "C2", "C3"), 5)
)
set.seed(123)
# dataset with the information of the persons
persons_df <- data.frame(
age = rpois(100, 20),
earnings = rnorm(100, 100, 10) |> abs(),
city = sample(LETTERS[1:15], 100, replace = T)
)
full_df = persons_df |>
left_join(city_df)
The data looks like this:
age earnings city poverty_level
1 17 113.76610 A C1
2 23 97.58023 A C1
3 19 87.22549 A C1
4 18 76.00802 A C1
5 14 103.62370 A C1
6 22 88.22033 A C1
7 21 117.97140 A C1
8 21 119.00900 A C1
9 27 101.92442 A C1
10 12 103.39957 B C2
11 25 107.90190 B C2
12 18 109.02506 B C2
13 14 99.11315 B C2
14 29 99.81366 B C2
15 17 103.14153 B C2
Including poverty_level as a random and fixed effect:
model_1 <- lmer(earnings ~ age + poverty_level + (1 + poverty_level|city), data=full_df)
summary(model_1)
Linear mixed model fit by REML ['lmerMod']
Formula: earnings ~ age + poverty_level + (1 + poverty_level | city)
Data: full_df
REML criterion at convergence: 747.4
Scaled residuals:
Min 1Q Median 3Q Max
-2.38084 -0.59015 -0.02984 0.64791 2.18420
Random effects:
Groups Name Variance Std.Dev. Corr
city (Intercept) 0.000e+00 0.000e+00
poverty_levelC2 4.633e-11 6.806e-06 NaN
poverty_levelC3 1.246e+01 3.530e+00 NaN -0.92
Residual 1.143e+02 1.069e+01
Number of obs: 100, groups: city, 15
Fixed effects:
Estimate Std. Error t value
(Intercept) 98.4492 5.7722 17.056
age 0.1084 0.2800 0.387
poverty_levelC2 -0.9915 2.6323 -0.377
poverty_levelC3 -1.0670 3.0834 -0.346
Correlation of Fixed Effects:
(Intr) age pvr_C2
age -0.947
pvrty_lvlC2 -0.231 0.003
pvrty_lvlC3 -0.116 -0.083 0.427
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
Including poverty_level only as a fixed effect:
model_2 <- lmer(earnings ~ age + poverty_level + (1|city), data=full_df)
summary(model_2)
Linear mixed model fit by REML ['lmerMod']
Formula: earnings ~ age + poverty_level + (1 | city)
Data: full_df
REML criterion at convergence: 748
Scaled residuals:
Min 1Q Median 3Q Max
-2.34529 -0.64463 -0.03629 0.70258 2.15640
Random effects:
Groups Name Variance Std.Dev.
city (Intercept) 0.0 0.00
Residual 117.7 10.85
Number of obs: 100, groups: city, 15
Fixed effects:
Estimate Std. Error t value
(Intercept) 97.9364 5.7683 16.978
age 0.1347 0.2793 0.482
poverty_levelC2 -0.9907 2.6708 -0.371
poverty_levelC3 -1.3722 2.6670 -0.515
Correlation of Fixed Effects:
(Intr) age pvr_C2
age -0.945
pvrty_lvlC2 -0.234 0.003
pvrty_lvlC3 -0.129 -0.109 0.500
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
I can not include the "poverty level" only as a random effect because i need the estimate of its effect.
I am including the "poverty level" correctly in the model_1 or in the model_2?, if not is there a correct way the include the "poverty level" as fixed effect?
