8.3 Modelo multinomial para variables con múltiples categorías
En las secciones anteriores se explicó cómo hacer un análisis de variables aleatorias provenientes de encuestas complejas de tipo binario usando modelos de regresión logística. Sin embargo, es muy común que en estas encuestas existan variables de tipo multinomial; es decir, que tengan más de dos niveles de respuesta. Para el análisis de este tipo de variables se utiliza un modelo multinomial, el cual se desarrollará en esta sección.
El modelo de regresión logístico multinomial es la extensión natural del modelo de regresión logístico binomial y se utiliza para analizar variables con respuestas que tienen tres o más categorías distintas. Esta técnica es más apropiada para variables de encuestas con categorías de respuesta nominales. Para el ajuste de un modelo multinomial se debe tener en cuenta que la variable dependiente debe medirse en el nivel nominal; se debe tener una o más variables independientes que pueden ser continuas o categóricas; la variable dependiente debe tener categorías mutuamente excluyentes y exhaustivas; no debe haber dos o más variables independientes que están altamente correlacionadas entre sí; ademas, debe haber una relación lineal entre cualquier variable independiente continua y la transformación logit de la variable dependiente. El modelo múltinomial está dado por la siguiente expresión:
\[ Pr\left(Y_{ik}\right)=Pr\left(y_{i}=k\mid\boldsymbol{X}_i\right)=\frac{\exp\left(\beta_{0k}+\boldsymbol{X}_{i}'\boldsymbol{\beta}_{k}\right)}{\sum_{j=1}^{m}\exp\left(\beta_{0j}+\boldsymbol{\boldsymbol{X}_{i}'\beta}_{j}\right)} \]
En donde \(Pr\left(y_{i}=k\mid\boldsymbol{X}_i\right)\) representa la probabilidad condicional de que la observación \(i\) pertenezca a la categoría \(k\) de la variable dependiente \(Y\), dado el conjunto de covariables \(\boldsymbol{X}_{i}\); \(\boldsymbol{\beta}_k\) es el vector de coeficientes de regresión para las variables \(\boldsymbol{X}\) en la \(k\)-ésima categoría. Para la estimación de los parámetros del modelo logístico es posible usar técnicas ajustadas de máxima pseudoverosimilitud, lo que implica maximizar la siguiente versión multinomial de la función de pseudoverosimilitud (Heeringa, West, y Berglund 2017):
\[ PL_{Mult}\left(\hat{\beta}\mid\boldsymbol{X}\right) = \prod_{i=1}^{n}\left\{ \prod_{i=1}^{k}\hat{\pi}_{k}\left(x_{i}\right)^{y_{i\left(k\right)}}\right\} ^{w_{i}} \]
En la anterior expresión, se tiene que
\(y_{i\left(k\right)}=1\) si \(y = k\) para la unidad muestreada \(i\), 0 en caso contrario. \(\hat{\pi}_{k}\left(x_{i}\right)\) es la probabilidad estimada de que \(y_{i}=k\mid x_{i}\). \(w_{i}\) es el peso de la encuesta para la unidad muestreada \(i\). David A. Binder (1983) afirma que la maximización implica la aplicación del algoritmo de Newton-Raphson para resolver el siguiente conjunto de ecuaciones de estimación \(\left(K-1\right)\times\left(p+1\right)\), suponiendo un diseño de muestra complejo con estratos indexados por \(h\) y conglomerados dentro de estratos indexados por \(\alpha\):
\[ S\left(\boldsymbol{\beta}\right)_{Mult} = \sum_{h}\sum_{\alpha}\sum_{i}\omega_{h\alpha i}\left(y_{h\alpha i}^{\left(k\right)}-\pi_{k}\left(\boldsymbol{\beta}\right)\right)x'_{h\alpha i}=0 \]
Donde,
\[ \pi_{k}\left(\boldsymbol{\beta}\right)=\frac{exp\left(\boldsymbol{x}'\boldsymbol{\beta_{k}}\right)}{1+{ \sum_{k=2}^{K}exp\left(\boldsymbol{x}'\boldsymbol{\beta_{k}}\right)}} \]
y
\[ \pi_{1\left(referencia\right)}\left(\boldsymbol{\beta}\right)=1-\sum_{k=2}^{K}\pi_{k}\left(\boldsymbol{\beta}\right) \]
Por último, la matriz de varianzas de los parámetros estimados basada en la aplicación de David A. Binder (1983) de la linealización de Taylor a las estimaciones derivadas utilizando la estimación de pseudomáxima verosimilitud es:
\[ \hat{Var}\left(\hat{\boldsymbol{\beta}}\right) = \left(J^{-1}\right)var\left[S\left(\hat{\boldsymbol{\beta}}\right)\right]\left(J^{-1}\right) \]
Para ejemplificar el uso de los modelos multinomiales en encuestas de hogares utilizando R se utilizan los siguientes códigos computacionales. Inicialemente, para efectos de comparación, estimemos el porcentaje de personas por desempleo:
| Employment | Prop | Prop_se | Prop_low | Prop_upp |
|---|---|---|---|---|
| Unemployed | 0.0427 | 0.0073 | 0.0283 | 0.0571 |
| Inactive | 0.3841 | 0.0153 | 0.3538 | 0.4143 |
| Employed | 0.5732 | 0.0140 | 0.5454 | 0.6010 |
Ahora bien, la estimación del modelo multinomial se realiza con la función svy_vglm del paquete svyVGAM como se muestra a continuación:
library(svyVGAM)
model_mul <- svy_vglm(
formula = Employment ~ Age + Sex + Region + Zone,
design = diseno,
crit = "coef",
family = multinomial(refLevel = "Unemployed")
)
summary(model_mul)## svy_vglm.survey.design(formula = Employment ~ Age + Sex + Region +
## Zone, design = diseno, crit = "coef", family = multinomial(refLevel = "Unemployed"))
## Stratified 1 - level Cluster Sampling design (with replacement)
## With (238) clusters.
## Called via srvyr
## Sampling variables:
## - ids: PSU
## - strata: Stratum
## - weights: qw
## Data variables: HHID (chr), Stratum (chr), NIh (int), nIh (dbl), dI (dbl),
## PersonID (chr), PSU (chr), Zone (chr), Sex (chr), Age (int), MaritalST (fct),
## Income (dbl), Expenditure (dbl), Employment (fct), Poverty (fct), dki (dbl),
## dk (dbl), wk (dbl), Region (fct), CatAge (ord), Age2 (I<dbl>), qw (dbl),
## Pobreza (dbl), Desempleo (dbl)
## Coef SE z p
## (Intercept):1 2.41129 0.78316 3.08 0.00208
## (Intercept):2 2.21128 0.63743 3.47 0.00052
## Age:1 0.02160 0.01066 2.03 0.04275
## Age:2 0.01823 0.00905 2.01 0.04406
## SexMale:1 -2.22316 0.32064 -6.93 4.1e-12
## SexMale:2 -0.59506 0.27908 -2.13 0.03299
## RegionSur:1 -0.43199 0.71209 -0.61 0.54409
## RegionSur:2 -0.27702 0.57377 -0.48 0.62923
## RegionCentro:1 0.36911 0.63111 0.58 0.55864
## RegionCentro:2 0.25684 0.54089 0.47 0.63490
## RegionOccidente:1 0.25319 0.63483 0.40 0.69002
## RegionOccidente:2 0.09176 0.51607 0.18 0.85888
## RegionOriente:1 0.61915 0.67153 0.92 0.35653
## RegionOriente:2 0.47226 0.61531 0.77 0.44277
## ZoneUrban:1 -0.22478 0.41830 -0.54 0.59102
## ZoneUrban:2 0.05425 0.37106 0.15 0.88376La anterior salida muestra los coeficientes estimados para cada variable independiente en el modelo, así como sus errores estándar y valores p asociados. Los valores de los interceptos para las clases Employed e Inactive del resultado multinomial indican el valor logarítmico de la proporción de probabilidades de pertenecer a la clase correspondiente en comparación con la clase de referencia que es Unemployed. Luego, se muestran los coeficientes asociados con la variable Age para las clases Employed e Inactive, que indican cómo el logaritmo de la proporción de probabilidades de pertenecer a una clase en comparación con la clase de referencia cambia por una unidad de cambio en la edad. A continuación, vienen los coeficientes asociados con la variable binaria ‘SexMale’ para las clases Employed e Inactive, que indican cómo el logaritmo de la proporción de probabilidades de pertenecer a una clase en comparación con la clase de referencia cambia cuando la variable ‘SexMale’ cambia de categoría. El resto de coeficientes en la salida se interpreta análogamente.
Las columnas SE, z y p proporcionan información sobre la significancia estadística de cada coeficiente. El valor p indica si el coeficiente es significativamente diferente de cero. En este caso, muchos de los coeficientes tienen valores de p bajos, lo que sugiere que son estadísticamente significativos.