Random Coefficients Modeling
Modelos multinivel para variables binarias
Gonzalo J. Muñoz
2019-07-10
Modelos multinivel para variables binarias
En un modelo multinivel el valor de la persona i en el grupo j se denota como \(y_{ij}\) y una variable explicativa/predictora de nivel 1 como \(x_ij\).
\[y_{ij}=\beta_0+\beta_1x_{ij}+u_j+e_{ij}\]
El modelo correspondiente para outcomes binarios se especifica exactamente del mismo modo pero de acuerdo al modelo generalizado ahora incluyen una función de link (logit) que es el inverso de una distribución acumulativa con una distribución conocida, \(F^-1\). El link logit es una de estas funciones, que es la que hemos usado hasta ahora y se escribiría de esta forma.
\[log\bigg{(}\frac{\pi}{1-\pi}\bigg{)}=\beta_0+\beta_1X_{ij}+u_j\]
El valor de \(\beta_1\) es el efecto de aumentar x en una unidad sobre el log odds de que y = 1, sólo que ahora es el efecto de x controlando el efecto grupal. A menos que \(\sigma^2_{u}\) = 0, en el modelo logístico multinivel el promedio de X no es igual a la gran media.
Intercept-only model
Primero vamos a leer la misma base de datos que analizamos previamente.
mydata <- read.table("6.1.txt", sep = ",", header = TRUE)
La sintaxis para obtener el modelo multinivel para outcomes binarios es muy similar a la anterior, pero ahora usando la función glmer especificando que la función de link es el logit y que hay una variable de agrupación (que en este caso es la comunidad a que pertence la persona o comm)
require(lme4)
## Loading required package: lme4
## Loading required package: Matrix
fit <- glmer(antemed ~ (1 | comm), family = binomial("logit"), data = mydata)
summary(fit)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: antemed ~ (1 | comm)
## Data: mydata
##
## AIC BIC logLik deviance df.resid
## 6639.5 6652.7 -3317.8 6635.5 5364
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7779 -0.7458 0.3423 0.7118 2.6784
##
## Random effects:
## Groups Name Variance Std.Dev.
## comm (Intercept) 1.464 1.21
## Number of obs: 5366, groups: comm, 361
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.14809 0.07178 2.063 0.0391 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
De acuerdo a este modelo, podemos decir que el log odds de recibir atención prenatal de una comunidad promedio (donde \(u_{0j}\)) es \(\hat\beta_0=0.148\). El intercepto para la counidad j es 0.148 + \(u_{0j}\) con varianza \(\hat\sigma^2_{u_{0}}\) = 1.464.
Podemos determinar si \(\hat\sigma^2_{u_{0}}\) es estadísticamente significativo comparando el ajuste del modelo multinivel con el ajuste de un nivelo uninivel.
fita <- glm(antemed ~ 1, data = mydata, family = binomial("logit"))
summary(fita)
##
## Call:
## glm(formula = antemed ~ 1, family = binomial("logit"), data = mydata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.200 -1.200 1.155 1.155 1.155
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.05219 0.02731 1.911 0.056 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 7435.2 on 5365 degrees of freedom
## Residual deviance: 7435.2 on 5365 degrees of freedom
## AIC: 7437.2
##
## Number of Fisher Scoring iterations: 3
La diferencia entre el ajuste de este último modelo (7435.2) y el modelo con el efecto del grupo (6635.524) es de 799.68. La prueba de \(\chi^2\) con gl = 1 es obviamente estadísticamente significativa.
# se pueden extraer los log likelihood de un modelo con la función `logLik`.
# El valor de deviance es -2*loglikelihood
-2*logLik(fita)
## 'log Lik.' 7435.202 (df=1)
-2*logLik(fit)
## 'log Lik.' 6635.524 (df=2)
-2*logLik(fita)-(-2*(logLik(fit)))
## 'log Lik.' 799.6782 (df=1)
u0 <- ranef(fit, condVar = TRUE)
u0se <- sqrt(attr(u0[[1]], "postVar")[1, , ])
commid <- as.numeric(rownames(u0[[1]]))
u0tab <- cbind("commid" = commid, "u0" = u0[[1]], "u0se" = u0se)
colnames(u0tab)[2] <- "u0"
u0tab <- u0tab[order(u0tab$u0), ]
u0tab <- cbind(u0tab, c(1:dim(u0tab)[1]))
u0tab <- u0tab[order(u0tab$commid), ]
colnames(u0tab)[4] <- "u0rank"
plot(u0tab$u0rank, u0tab$u0, type = "n", xlab = "u_rank", ylab = "conditional
modes of r.e. for comm_id:_cons", ylim = c(-4, 4))
segments(u0tab$u0rank, u0tab$u0 - 1.96*u0tab$u0se, u0tab$u0rank, u0tab$u0 +
1.96*u0tab$u0se)
points(u0tab$u0rank, u0tab$u0, col = "blue")
abline(h = 0, col = "red")
Este gráfico simplemente sirve para ilustrar que cada una de las 361 comunidades en la muestra tiene un efecto \(u_0\) y un residuo \(\sigma^2_e\) y que algunas comunidades están por sobre o por debajo del promedio poblacional, como lo demuestra el intervalo de confianza del 95%.
Means-as-outcomes model
A continuación incluiremos la edad de la madre (centrada) como variable explicativa en el modelo.
\[log\bigg{(}\frac{\pi}{1-\pi}\bigg{)}=\beta_0+\beta_1magec_{ij}+u_j\]
mydata$meduc2 <- mydata$meduc == 2
mydata$meduc3 <- mydata$meduc == 3
summary(mydata$mage)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 13.00 19.00 23.00 23.63 28.00 49.00
mydata$magec <- mydata$mage - mean(mydata$mage)
mydata$magecsq <- mydata$magec * mydata$magec
fit2 <- glmer(antemed ~ magec + (1 | comm), family = binomial("logit"), data = mydata)
summary(fit2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: antemed ~ magec + (1 | comm)
## Data: mydata
##
## AIC BIC logLik deviance df.resid
## 6603.4 6623.2 -3298.7 6597.4 5363
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9757 -0.7431 0.3357 0.7190 3.2357
##
## Random effects:
## Groups Name Variance Std.Dev.
## comm (Intercept) 1.462 1.209
## Number of obs: 5366, groups: comm, 361
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.144604 0.071781 2.015 0.044 *
## magec -0.032357 0.005235 -6.181 6.37e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## magec 0.008
La varianza estimada del intercepto aleatorio es casi igual a la obtenida en el paso anterior, lo que sugiere que la comunidades tienen una composción similar en términos de las edades de las madres.
La ecuación para los efectos fijos quedaría como
\[log\bigg{(}\frac{\pi}{1-\pi}\bigg{)}=0.144-0.032magec_{ij}\]
Si graficamos los datos del modelo resultante obtenemos una línea recta que representa cada una de las 361 comunidades representadas en la data. Fijarse que la pendiente es la misma para todos los grupos.
## Loading required package: VGAM
## Loading required package: stats4
## Loading required package: splines
## Loading required package: lattice
Podemos incluir la varianza aleatoria de la pendiente especificando (magec | comm). Sin embargo, como se verá a continuación este modelo anterior tienen un singular fit. La razón en este caso es que la varianza de la pendiente es cercana a cero. De todas maneras se puede parametrizar el modelo y hacer un gráfico, pero claramente las pendientes siguen siendo casi paralelas.
fit2.B <- glmer(antemed ~ magec + (magec | comm), family = binomial("logit"), data = mydata, glmerControl(calc.derivs = FALSE))
summary(fit2.B)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: antemed ~ magec + (magec | comm)
## Data: mydata
## Control: glmerControl(calc.derivs = FALSE)
##
## AIC BIC logLik deviance df.resid
## 6606.8 6639.8 -3298.4 6596.8 5361
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9334 -0.7463 0.3422 0.7179 3.4338
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## comm (Intercept) 1.461e+00 1.208861
## magec 3.266e-05 0.005715 1.00
## Number of obs: 5366, groups: comm, 361
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.143236 0.071355 2.007 0.0447 *
## magec -0.031658 0.005198 -6.091 1.12e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## magec 0.068
Por otra parte, podemos incorporar el hecho de que la relación entre la edad de la madre y la probabilidad de recibir atención prenatal es no-lineal. Para ello incluimos la variable magecsq en el modelo.
fit3 <- glmer(antemed ~ magec + magecsq + (1 | comm), family = binomial("logit"), data = mydata)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00146687
## (tol = 0.001, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
summary(fit3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: antemed ~ magec + magecsq + (1 | comm)
## Data: mydata
##
## AIC BIC logLik deviance df.resid
## 6602.8 6629.2 -3297.4 6594.8 5362
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0432 -0.7434 0.3364 0.7142 3.6900
##
## Random effects:
## Groups Name Variance Std.Dev.
## comm (Intercept) 1.452 1.205
## Number of obs: 5366, groups: comm, 361
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.1847867 0.0758194 2.437 0.0148 *
## magec -0.0275872 0.0060187 -4.584 4.57e-06 ***
## magecsq -0.0010630 0.0006627 -1.604 0.1087
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) magec
## magec 0.168
## magecsq -0.330 -0.485
## convergence code: 0
## Model failed to converge with max|grad| = 0.00146687 (tol = 0.001, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
De nuevo el modelo falla en converger y, además, la variable magecsq no es estadísticamente significativa. De todas maneras, podemos visualizar el resultado para apreciar el modelo que se está proponinendo.
predprob <- fitted(fit3)
predlogit <- logit(predprob)
datapred <- unique(data.frame(cbind(predlogit = predlogit, comm =
mydata$comm, mage = mydata$mage)))
xyplot(predlogit ~ mage, data = datapred, groups = comm, type = c("p", "l",
"g"), col = "blue", xlim = c(9, 51), ylim = c(-4, 4))
Means-as-outcomes model (alternativo)
Evaluemos un modelo para ilustrar cómo funciona un modelo de pendiente aleatoria para la variable wealth pero centrada.
mean(mydata$wealth)
## [1] 3.0082
mydata$wealthc<-mydata$wealth-mean(mydata$wealth)
mydata$meduc2 <- as.numeric(mydata$meduc2)
mydata$meduc3 <- as.numeric(mydata$meduc3)
fit <- glmer(antemed ~ magec + magecsq + meduc2 + meduc3 + wealthc + (1 + wealthc | comm), data = mydata, family = binomial("logit"), glmerControl(calc.derivs = FALSE))
summary(fit)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: antemed ~ magec + magecsq + meduc2 + meduc3 + wealthc + (1 +
## wealthc | comm)
## Data: mydata
## Control: glmerControl(calc.derivs = FALSE)
##
## AIC BIC logLik deviance df.resid
## 5987.0 6046.3 -2984.5 5969.0 5357
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3332 -0.6531 0.2905 0.6410 4.3985
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## comm (Intercept) 0.84208 0.9176
## wealthc 0.01427 0.1195 -1.00
## Number of obs: 5366, groups: comm, 361
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.4547939 0.0823497 -5.523 3.34e-08 ***
## magec -0.0001231 0.0065341 -0.019 0.985
## magecsq -0.0010561 0.0006816 -1.549 0.121
## meduc2 0.5471158 0.0846714 6.462 1.04e-10 ***
## meduc3 1.3092769 0.0968300 13.521 < 2e-16 ***
## wealthc 0.4050985 0.0301320 13.444 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) magec magcsq meduc2 meduc3
## magec -0.022
## magecsq -0.279 -0.489
## meduc2 -0.545 0.206 -0.047
## meduc3 -0.543 0.307 -0.069 0.553
## wealthc 0.008 -0.120 0.070 -0.159 -0.344
Comparamos el modelo fit con el modelo fita donde la pendiente de wealthc es fija
fita <- glmer(antemed ~ magec + magecsq + meduc2 + meduc3 + wealthc + (1 | comm), data = mydata, family = binomial("logit"), glmerControl(calc.derivs = FALSE))
-2*logLik(fita)
## 'log Lik.' 5979.303 (df=7)
-2*logLik(fit)
## 'log Lik.' 5968.988 (df=9)
-2*logLik(fita)-(-2*(logLik(fit)))
## 'log Lik.' 10.31497 (df=7)
La diferencia en los deviance es de 10.31 con grados de libertas = 7-9=2, que es estadísticamente significativo. Específicamente, el el valor p para \(\chi^2\) 10.31487 con 2 grados de libertad es 0.005756 que es menor que \(\alpha\) = .05.
pchisq(10.31497, df=2,lower.tail = FALSE)
## [1] 0.005756158
Por supuesto, obtenemos el mismo resultado si usamos la función anova para comparar ambos modelos.
anova(fita,fit)
## Data: mydata
## Models:
## fita: antemed ~ magec + magecsq + meduc2 + meduc3 + wealthc + (1 |
## fita: comm)
## fit: antemed ~ magec + magecsq + meduc2 + meduc3 + wealthc + (1 +
## fit: wealthc | comm)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fita 7 5993.3 6039.4 -2989.7 5979.3
## fit 9 5987.0 6046.3 -2984.5 5969.0 10.315 2 0.005756 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Means-as-outcomes con variables de Nivel 2
Añadir variables de nivel 2 en el contexto de un modelo logístico no difiere sustantivamente de un modelo normal. Las variables de nivel 2 se incluyen por sí solas o en conjunto con las variables de nivel 1. Hox (2010) recomienda agregar las variables de nivel 2 con las variables de nivel 1 que demostraron ser estadísticamente significativas.
#Modelo SIN urban
fita <- glmer(antemed ~ magec + (1 | comm), data = mydata, family = binomial("logit"),glmerControl(calc.derivs = FALSE))
summary(fita)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: antemed ~ magec + (1 | comm)
## Data: mydata
## Control: glmerControl(calc.derivs = FALSE)
##
## AIC BIC logLik deviance df.resid
## 6603.4 6623.2 -3298.7 6597.4 5363
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9757 -0.7431 0.3357 0.7190 3.2357
##
## Random effects:
## Groups Name Variance Std.Dev.
## comm (Intercept) 1.462 1.209
## Number of obs: 5366, groups: comm, 361
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.144604 0.071361 2.026 0.0427 *
## magec -0.032357 0.005163 -6.268 3.67e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## magec 0.009
#Modelo CON urban
fit <- glmer(antemed ~ magec + urban + (1 | comm), data = mydata, family = binomial("logit"),glmerControl(calc.derivs = FALSE))
summary(fit)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: antemed ~ magec + urban + (1 | comm)
## Data: mydata
## Control: glmerControl(calc.derivs = FALSE)
##
## AIC BIC logLik deviance df.resid
## 6492.1 6518.4 -3242.0 6484.1 5362
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1529 -0.7345 0.3229 0.7207 3.1805
##
## Random effects:
## Groups Name Variance Std.Dev.
## comm (Intercept) 0.9699 0.9848
## Number of obs: 5366, groups: comm, 361
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.346250 0.073817 -4.691 2.72e-06 ***
## magec -0.032481 0.005157 -6.299 3.00e-10 ***
## urban 1.494416 0.131532 11.362 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) magec
## magec 0.019
## urban -0.561 -0.018
Es importante considerar que en este modelo el coeficiente para magec corresponde a log odds para la edad de la madre cuando urban=0, es decir, cuando la comunidad es rural. El coeficiente para urban refleja la diferencia entre los log odds de rural y urbano para una madre de edad promedio.
La comparación entre ambos modelo muestra que incluir urban mejora el ajuste del modelo
anova(fita,fit)
## Data: mydata
## Models:
## fita: antemed ~ magec + (1 | comm)
## fit: antemed ~ magec + urban + (1 | comm)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fita 3 6603.4 6623.2 -3298.7 6597.4
## fit 4 6492.1 6518.4 -3242.0 6484.1 113.35 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Cross-level effects
La interacción entre urban y magec no es estadísticamente significativa.
#mydata$urban
#Modelo CON urban y cross-level effect
fit2 <- glmer(antemed ~ magec*wealthc + (1 | comm), data = mydata, family = binomial("logit"),glmerControl(calc.derivs = FALSE))
summary(fit2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: antemed ~ magec * wealthc + (1 | comm)
## Data: mydata
## Control: glmerControl(calc.derivs = FALSE)
##
## AIC BIC logLik deviance df.resid
## 6169.5 6202.4 -3079.7 6159.5 5361
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3570 -0.6803 0.3016 0.6689 3.1448
##
## Random effects:
## Groups Name Variance Std.Dev.
## comm (Intercept) 0.8004 0.8947
## Number of obs: 5366, groups: comm, 361
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.119839 0.057491 2.084 0.0371 *
## magec -0.029420 0.005326 -5.524 3.32e-08 ***
## wealthc 0.565705 0.026857 21.064 < 2e-16 ***
## magec:wealthc -0.003114 0.003848 -0.809 0.4183
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) magec welthc
## magec 0.009
## wealthc 0.001 -0.001
## magec:wlthc -0.012 0.111 -0.028