Hasta ahora hemos desarrollado modelos multinivel para variables predichas continuas. Sin embargo, muchas variables de interés en ciencias sociales son categóricas nominales (aprobar vs reprobar un curso, votar por un candidato de derecha o izquierda) u ordinales (grado de acuerdo con una pregunta). Ahora Veremos el caso de variables de salida nominales binarias o dicotómicas donde sólo hay dos alternativas de respuesta que no tienen un orden inherente. La variable predictora puede ser continua o categórica. Para comenzar, veremos el caso de la regresión logística para ilustrar los principios que luego extenderemos al caso del análisis multinivel.
Regresión logística
Para esta parte del curso analizaremos los resultados de una encuesta nacional transeccional (BDHS2004). La variable dependiente es un indicador de atención prenatal (antemed) al menos una vez antes del nacimiento del hijo más reciente (dentro de los últimos 5 años). La variables analizadas y sus respectivos códigos se encuentran en la Tabla 1.
Tabla 1. Bangladesh Demographic and Health Survey (BDHS2004)
varnames
vardescr
comm
Community identifier
womid
Woman identifier
antemed
Received antenatal care at least once from a medically-trained provider, e.g. doctor, nurse or midwife (1 = yes, 0 = no)
bord
Birth order of child (ranges from 1 to 13)
mage
Mother’s age at the child’s birth (in years)
urban
Type of region of residence at survey (1 = urban, 0 = rural)
meduc
Mother’s level of education at survey (1 = none, 2 = primary, 3 = secondary or higher)
islam
Mother’s religion (1 = Islam, 0 = other)
wealth
Household wealth index in quintiles (1 = poorest to 5 = richest)
En total respondieron 5366 personas, de las cuales 2753 (51%) recibieron atención prenatal.
cbind(Freq = table(mydata$antemed), Perc = prop.table(table(mydata$antemed)),
Cum = cumsum(prop.table(table(mydata$antemed))))
## Freq Perc Cum
## 0 2613 0.4869549 0.4869549
## 1 2753 0.5130451 1.0000000
Media y varianza de variables binarias
El parámetro poblacional a estimar es \(\pi\) y lo estimaremos mediante la media de la variable dicotómica \(y\) que contiene sólo 1s y 0s.
\[\hat\pi=\overline{y}=Pr(y=1)=\frac{r}{N}\] donde r corresponde al número de individuos en la muestra que seleccionaron la opción “1”.1 La media y la desviación estándar se podrían obtener mediante summary(mydata$antemed)
La varianza muestral de y se puede obtener mediante la fórmula tradicional para obtener la varianza, que en el caso de variables dicotómicas (1 y 0) se puede simplificar a:
\[s^2=\hat\pi(1-\hat\pi)\]
Por lo tanto, la varianza de antemed se podría obtener como \(s^2=0.513(1-0.513)=0.248\) y la desviación estándar \(s=\sqrt{0.248}=0.4998\)
Se le llama distribución de Bernoulli a una distribución que tiene media \(\pi\) y varianza \(\pi(1-\pi)\).
En la ausencia de información adicional, el valor esperado para un individuo i es igual a la probabilidad de respuesta de la poblacióm, esto es \(\pi_i=\pi\). Si el 51% de la muestra dijo 1, entonces para cualquier individuo i la probabilidad de responder 1 es .51 más un error, esto es:
\[y_i=\pi+e_i\] Pero este es el caso si no hay ningún predictor. Si incluimos las variables \(X_1...X_p\) entonces \(\pi_i=\beta_0 + \beta_1X_1\), por lo que el valor esperado de \(y_i\) sería:
\[y_i=\pi+e_i=\beta_0 + \beta_1X_1+e_i\]
Visualizando la relación entre variables cuando una de ellas es binaria
plot(mydata$mage, mydata$antemed, xlab = "Edad de la madre al nacer", ylab =
"Atencion prenatal")
Una mejor manera de visualizar la relación entre ambas variables es agrupar la vd según la edad de la madres, de modo que cada punto represente la proporción de madres de x edad que recibieron atención prenatal.
propantemed <- tapply(mydata$antemed, mydata$mage, mean)
plot(names(propantemed), propantemed, xlim = c(10, 50), xlab = "Edad de la madre al nacer")
Evidentemente la probabilidad de recibir atención prenatal disminuye con la edad de la madre, pero esta relación no es lineal.
Consideremos si el lugar de residencia (urban) y el nivel de educación de la madre (meduc) como posibles predictores.
La proporción de madres que reciben atención prenatal viven preponderantemente en zonas urbanas (urban = 1) y tienen educación secundaria o superior (meduc = 3)
Regresión lineal para modelar atención prenatal
Hay varias dificultades asociadas al uso de la regresión lineal para modelar datos binarios, especialmente cuando los valores de \(\pi\) se aproximan a 1 o 0. Para ilustrar este punto, amalizaremos los datos anteriores usando la edad de la madre y su lugar de residencia. Previamente haremos algunas transformaciones (i.e., centrar la edad de la madre, agregar una variable para representar el efecto cuadrático de la edad centrada de la madres, y agregar dos variables dummy para representar el nivel de educación de la madre usando como categoría de referencia la falta de educación).
A cotinuación haremos un modelo para evaluar la relación entre magec y antemed pero teniendo la precaución de fijar los valores de los predictores del modelo en torno a su media. De esta manera, el coeficiente para magec se puede interpretar como la proporción de madres (de edad promedio) que recibieron atención prenatal que viven en una misma ubicación y que tienen el mismo nivel de educación.
mydata2 <- mydata2[order(mydata2$mage), ]
plot(mydata2$mage, mydata2$predantemed, xlim = c(10, 50), xlab = "Edad de la madre al nacimiento", ylab = "Pr(antemed)", type = "l")
Como se puede observar hay sólo dos valores posibles para cada error de predicción. Por ejemplo, para todas las madres cuya probabilidad de ser atendida es de .14 (que corresponde a los puntos más a la izquierda) el error de estimación es o -.31 o 1.91. Esto es así porque el valor observado para cada caso es o 1 o 0. El resultado es que el error de predicción no está normalmente distribuido.
Un segundo problema es que la varianza de \(\pi\) depende del valor de x. Particularmente, en los extremos de x la varianza de \(\pi\) va a ser menor.
Relacionado con lo anterior, la forma de la relación entre una variable continua y una variable binaria va a tener una forma de S. Sólo para los valores intermedios de \(\pi\) la relación con x va a ser lineal. Por ejemplo, la mayor intensidad del estudio (x) va a afectar la probabilidad de pasar el examen más para los alumnos que están en el centro de la distribución que los que están en los extremos.
Por último, un modelo de regresión lineal produce predicciones que están fuera del rango e \(\pi\), esto es, menores que 0 o mayores que 1, lo que claramente es inadecuado.
Por todas estas razones, es obligatorio hacer una transformación de los valores predichos de una regresión multiple que genere valores contenidos entre 0 y 1. Esta transformación se conoce como transformación logística.
El modelo logit/logístico
Una manera de representar la probabilidad de ocurrencia de un evento es mediante el cociente entre la probabilidad de que ocurra el evento \(\pi\) y la probabilidad de que no ocurra \(1-\pi\) que en inglés se llama odds.
\[odds=\frac{\pi}{1-\pi}\]
Este cociente tiene una distribución asimétrica positiva pero se puede transformar tomando su logaritmo natural3 Ants de continuar mira este video de Josh Starmer de Statquest
#Generar una distribución normal con media 0 y desviación estándar 1
bdist<-data.frame(rnorm(200, mean=0,sd=1))
#Etiquetar a primera columna como "z"
names(bdist)[1]<-"z"
#Transformar los valores z en odds
bdist$log<-exp(bdist$z)/(1+exp(bdist$z))
#Graficar los resultados
plot(bdist$z,bdist$log)
En esencia, vamos a interpretar los valores predichos de la regresión múltiple como un log(odds).
\[log\bigg{(}\frac{\pi}{1-\pi}\bigg{)}=\beta_0+\beta_1X_1\] ##Regresión logística con predictores continuos
Usando la edad de la madre (mage) como único predictor y como variable dependiente el recibir (o no) atención prenatal, podemos estimar un modelo de regresión logística usando la función glm y especificando que la distribución es binomial logística usando:
fit<-glm(antemed~mage, family = binomial(logit), data=mydata)
summary(fit)
##
## Call:
## glm(formula = antemed ~ mage, family = binomial(logit), data = mydata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.343 -1.208 1.045 1.134 1.501
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.785578 0.107847 7.284 3.24e-13 ***
## mage -0.031029 0.004415 -7.028 2.10e-12 ***
## ---
## 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: 7385.1 on 5364 degrees of freedom
## AIC: 7389.1
##
## Number of Fisher Scoring iterations: 4
El intercepto en el modelo anterior es 0.785578 y la pendiente -0.031029. Por lo tanto, el log odds de recibir atención prenatal para una madre de edad promedio (23.63) sería:
Para transformar el log odds resultante en la probabilidad de recibir atención prenatal hay que obtener el inverso de la función logit (o antilogit). La función logit del paquete VGAM convierte el logit en una probabilidad especificando inverse=TRUE
require(VGAM)
## Loading required package: VGAM
## Loading required package: stats4
## Loading required package: splines
logit(0.785578-0.031029*23.63436, inverse=TRUE)
## [1] 0.5130539
En consecuencia la probabilidad de que una madre de edad promedio reciba atención prenatal es de 51.31%. Podemos graficar los valores predicho para todos los niveles de mage.
#La opción `type = "response"` es importante para que la predicción
#ocurra en la escala de la VD
mydata$predprob <- predict(fit, type = "response")
plot(mydata$mage,mydata$predprob)
En el gráfico se observa que una madre de 45 años tiene una probabilidad mucho menor (aprox. 35%) de recibir atención prenatal que una de, por ejemplo, 15 años (58%).
Lo mismo pero usando la edad de la madre centrada
fit<-glm(antemed~magec, family = binomial(logit), data=mydata)
summary(fit)
##
## Call:
## glm(formula = antemed ~ magec, family = binomial(logit), data = mydata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.343 -1.208 1.045 1.134 1.501
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.052223 0.027440 1.903 0.057 .
## magec -0.031029 0.004415 -7.028 2.1e-12 ***
## ---
## 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: 7385.1 on 5364 degrees of freedom
## AIC: 7389.1
##
## Number of Fisher Scoring iterations: 4
Ahora el valor de \(X_1\) es 0, que es el promedio de magec. Sin embargo, la probabilidad resultante es exactamente la misma que en el modelo anterior salvo por el error de redondeo a partir del sexto decimal.
#La opción `type = "response"` es importante para que la predicción
#ocurra en la escala de la VD
mydata$predprob <- predict(fit, type = "response")
plot(mydata$mage,mydata$predprob)
Y una vez más pero incluyen magecsq
fit<-glm(antemed~magec+magecsq, family = binomial(logit), data=mydata)
summary(fit)
##
## Call:
## glm(formula = antemed ~ magec + magecsq, family = binomial(logit),
## data = mydata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.258 -1.238 1.099 1.115 1.973
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.1348006 0.0352034 3.829 0.000129 ***
## magec -0.0214217 0.0051078 -4.194 2.74e-05 ***
## magecsq -0.0021501 0.0005782 -3.718 0.000201 ***
## ---
## 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: 7370.9 on 5363 degrees of freedom
## AIC: 7376.9
##
## Number of Fisher Scoring iterations: 4
Se observa en el gráfico que dada la curvatura del modelo, la probabilidad para una madre de edad promedio (53.36%) mejora levemente en comparación con el modelo lineal (51.31%)
#La opción `type = "response"` es importante para que la predicción
#ocurra en la escala de la VD
mydata$predprob <- predict(fit, type = "response")
plot(mydata$mage,mydata$predprob)
Regresión logística con predictores categóricos
Tomando ahora el nivel de educación de la madre como variable predictora, podemos estimar el siguiente modelo:
La probabilidad de recibir atención prenatal en general es de 0.513, pero esta probabilidad varía según el nivel de educación de la madre. Las madres sin educación, con educación primaria y con educación secundaria o superior tienen una probabilidad de 0.318, 0.481, y 0.738 de recibir atención prenatal, respectivamente.
Por último, podemos comparar los odds de recibir atencion prenatal para los distintos niveles de educación. Por ejemplo, la razón de los odds (odds ratio) para las madres con educación primaria versus las madres sin educación sería:
(793/856)/(594/1272)
## [1] 1.98381
El 1.98 quiere decir que las madres con educación primaria tienen casi el doble de odds de recibir atención prenatal en comparación con las madres sin educación.
En resumen, tenemos que:
Educación
Probabilidad
Odds
Odds ratio
Sin educación
0.318
0.467
-
Primaria
0.481
0.926
1.984
Secundaria o superior
0.738
2.816
6.031
El modelo a estimar en este caso es: \[log\bigg{(}\frac{\pi}{1-\pi}\bigg{)}=\beta_0+\beta_1meduc1_i+\beta_2meduc3_i\] De nuevo, la sintaxis usa la función glm y se especifica a continuación:
fit <- glm(antemed ~ meduc2 + meduc3, data = mydata, family = binomial(logit))
summary(fit)
##
## Call:
## glm(formula = antemed ~ meduc2 + meduc3, family = binomial(logit),
## data = mydata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6367 -0.8754 0.7795 1.2100 1.5131
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.76147 0.04970 -15.323 <2e-16 ***
## meduc2 0.68502 0.06999 9.787 <2e-16 ***
## meduc3 1.79696 0.07255 24.768 <2e-16 ***
## ---
## 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: 6747.6 on 5363 degrees of freedom
## AIC: 6753.6
##
## Number of Fisher Scoring iterations: 4
Para obtener el efecto de educación podemos nuevamente usar la función logit del paquete VGAM. Como las variables meduc2 y meduc3 son variables dummy, hay que simplemente reemplazar por 1 y 0 para obtener los efectos que queremos determinar.
require(VGAM)
#Para madres sin educación (sólo intercepto)
logit(-0.76147+0.68502*0+1.79696*0, inverse=TRUE)
## [1] 0.3183272
#Para madres con educación primaria
logit(-0.76147+0.68502*1+1.79696*0, inverse=TRUE)
## [1] 0.4808968
#Para madres con educación secundaria o superior
logit(-0.76147+0.68502*0+1.79696*1, inverse=TRUE)
## [1] 0.7379789
Nótese que las probabilidades obtenidas calzan con las probabilidades de la tabla de resumen4 BAM!
#La opción `type = "response"` es importante para que la predicción
#ocurra en la escala de la VD
mydata$predprob <- predict(fit, type = "response")
plot(mydata$meduc,mydata$predprob)
Significancia estadística e intervalos de confianza
La significancia estadística de los coeficientes está incluida por defecto en el output de la función glm. Sin embargo, para obtener los intervalos de confianza para los coeficientes hay que usar la función confint() que también es parte del paquete VGAM.
En este caso, la hipótesis nula es que \(\beta=0\) que es lo mismo que decir que hay una probabilidad de un 50% de recibir atención prenatal. E ambos casos (meduc2 y meduc3) la probabilidad excede el 50%. Lo mismom se puede hacer (y de hecho me parece que es más común) usando los odds.
De nuevo, los valores reproducen lo que obtuvimos en la tabla de resumen pero ahora usando los odds en vez de las probabilidades. La interpretación es distinta. Un odd de 1 quiere decir que la chance de recibir el tratamiento es la misma que la de no recibirlo. En ambos caso, el intervalo de confianza de los odds exceden el valor de 1.
Por último, si el interés es evaluar el efecto de un set de predictores o la diferencia entre los coeficientes de dos predictores, podemos usar el test de Wald disponible en el paquete car5 ¡Precaución! La función logittambién es parte del paquete car. Para desambiguar su uso, de aquí en adelante hay que especificar si queremos usar la función logit del paquete VGAM VGAM::logit o car car::logit
require(car)
## Loading required package: car
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:VGAM':
##
## logit
#Evaluar un set de predictores
linearHypothesis(fit,c("meduc2","meduc3"),c(0,0))
El resultado del segundo test es interesante porque muestra que las madres con educación primaria y secundaria o superior difieren en cuanto a sus log odds de recibir atención prenatal.
