9.4 Modelo multinomial
En capítulos anteriores de este libro se desarrolló la manera de cómo hacer el análisis de variables aleatorias provenientes de encuestas complejas de tipo binario. Sin embargo, en encuestas complejas también existen 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 el modelo multinomial, el cual se desarrollará en esta sección.
El modelo de regresión logit multinomial es la extensión natural del modelo de regresión logística binomial simple para analizar variables con respuestas que tienen tres o más categorías distintas. Esta técnica es más apropiada para variables de encuesta con categorías de respuesta nominales.
Para el ajuste de un modelo multinomial se debe tener en cuenta que:
Su variable dependiente debe medirse en el nivel nominal.
Tiene una o más variables independientes que son continuas, ordinales o nominales (incluidas las variables dicotómicas).
Tener independencia de las observaciones y la variable dependiente debe tener categorías mutuamente excluyentes y exhaustivas
No debe haber multicolinealidad. La multicolinealidad ocurre cuando tiene dos o más variables independientes que están altamente correlacionadas entre sí.
Debe haber una relación lineal entre cualquier variable independiente continua y la transformación logit de la variable dependiente.
No debe haber valores atípicos, valores de apalancamiento elevados o puntos influyentes.
El modelo múltinomial está dado como:
\[ Pr\left(Y_{ik}\right)=Pr\left(y_{i}=k\mid\boldsymbol{x}_{i}:\boldsymbol{\beta}_{1},\dots\boldsymbol{\beta}_{m}\right)=\frac{\exp\left(\beta_{0k}+\boldsymbol{\beta}_{k}\boldsymbol{x}_{i}\right)}{\sum_{j=1}^{m}\exp\left(\beta_{0j}+\boldsymbol{\beta}_{j}\boldsymbol{x}_{i}\right)} \]
donde \(\boldsymbol{\beta}_k\) es el vector de coeficiente de \(\boldsymbol{X}\) para la k-ésima categoría.
Para la estimación de los parámetros del modelo logístico se utilizará la máxima verosimilitud como técnica de estimación. Para ello, implica maximizar la siguiente versión multinomial de la función de pseudoverosimilitud (Heeringa):
\[ 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}} \] donde,
\(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\).
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(x'\boldsymbol{\beta_{k}}\right)}{1+{ \sum_{k=2}^{K}exp\left(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 varianza-covarianza de los parámetros estimados, basada en la aplicación de Binder de la linealización de la serie 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:
diseno %>% group_by(Employment) %>%
summarise(Prop = survey_mean(vartype = c("se", "ci")))| Employment | Prop | Prop_se | Prop_low | Prop_upp |
|---|---|---|---|---|
| Unemployed | 0.0322524 | 0.0057006 | 0.0209647 | 0.0435401 |
| Inactive | 0.2734145 | 0.0106850 | 0.2522572 | 0.2945718 |
| Employed | 0.4141869 | 0.0138767 | 0.3867096 | 0.4416642 |
| NA | 0.2801462 | 0.0139885 | 0.2524476 | 0.3078449 |
Ahora bien, se filtran aquellas personas mayores a 15 años y se vuelve a estimar:
diseno %>% filter(Age >= 15)%>% group_by(Employment) %>%
summarise(Prop = survey_mean(vartype = c("se", "ci")))| Employment | Prop | Prop_se | Prop_low | Prop_upp |
|---|---|---|---|---|
| Unemployed | 0.0448040 | 0.0077955 | 0.0293683 | 0.0602398 |
| Inactive | 0.3798195 | 0.0150308 | 0.3500571 | 0.4095820 |
| Employed | 0.5753764 | 0.0130688 | 0.5494988 | 0.6012540 |
Con lo anterior, se observa un aumento en el porcentaje de personas empleadas en los hogares.
Ahora bien, se realiza la estimación usando el modelo multinomial con la función svy_vglm del paquete svyVGAM como se muestra a continuación:
library(svyVGAM)
diseno_15 <- diseno %>% filter(Age >= 15)
model_mul <- svy_vglm(
formula = Employment ~ Age + Sex + Region + Zone,
design = diseno_15,
crit = "coef",
family = multinomial(refLevel = "Unemployed"))
model_mul## Stratified 1 - level Cluster Sampling design (with replacement)
## With (238) clusters.
## Called via srvyr
## Sampling variables:
## - ids: PSU
## - strata: Stratum
## - weights: wk2
## 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), wk2 (dbl)
##
## Call:
## vglm(formula = formula, family = family, data = surveydata, weights = .survey.prob.weights,
## crit = "coef")
##
##
## Coefficients:
## (Intercept):1 (Intercept):2 Age:1 Age:2
## 2.29035224 2.09308478 0.02474836 0.02073894
## SexMale:1 SexMale:2 RegionSur:1 RegionSur:2
## -2.21949718 -0.55633376 -0.43624259 -0.27905123
## RegionCentro:1 RegionCentro:2 RegionOccidente:1 RegionOccidente:2
## 0.37132724 0.25576105 0.25356331 0.09279505
## RegionOriente:1 RegionOriente:2 ZoneUrban:1 ZoneUrban:2
## 0.61759957 0.47055482 -0.23461616 0.05674419
##
## Degrees of Freedom: 3758 Total; 3742 Residual
## Residual deviance: 2778.983
## Log-likelihood: -1389.492
##
## This is a multinomial logit model with 3 levelsAlgo a tener en cuenta es que la función broom::tidy(), que normalmente usamos para limpiar y estandarizar la salida del modelo, no puede ser empleada en este caso, sin embargo, en el link4 encuentra la función que utilizamos a continuación.
tidy.svyVGAM <- function(
x,
conf.int = FALSE,
conf.level = 0.95,
exponentiate = FALSE,
...
){
ret <- as_tibble(summary(x)$coeftable, rownames = "term")
colnames(ret) <- c("term", "estimate", "std.error", "statistic", "p.value")
coefs <- tibble::enframe(stats::coef(x), name = "term", value = "estimate")
ret <- left_join(coefs, ret, by = c("term", "estimate"))
if (conf.int){
ci <- broom:::broom_confint_terms(x, level = conf.level, ...)
ret <- dplyr::left_join(ret, ci, by = "term")
}
if (exponentiate){ret <- broom:::exponentiate(ret)}
ret %>%
tidyr::separate(term, into = c("term", "y.level"), sep = ":") %>%
arrange(y.level) %>%
relocate(y.level, .before = term)
}Utilizando el código anterior, obtenemos las estimaciones del modelo multinomial:
tab_model <- tidy.svyVGAM(model_mul,
exponentiate = FALSE,
conf.int = FALSE)
tab_modelA continuación, se grafica los intervalos de confianza para los coeficientes del modelo utilizando el formato definido para la CEPAL:
library(dotwhisker)
tab_model %>%
mutate(model = if_else(
y.level == 1,
"Inactive",
"Employed", ),
sig = gtools::stars.pval(p.value)
) %>%
dotwhisker::dwplot(
dodge_size = 0.3,
vline = geom_vline(xintercept = 1, colour = "grey60",
linetype = 2)
) +
guides(color = guide_legend(reverse = TRUE)) +
theme_bw() + theme(legend.position = "top")Una forma aternativa de ajustar los modelos de regresión multinomial es usando la función svy_vglm. Sin embargo, presenta limitaciones para hacer las predicciones con el modelo, por lo tanto, se puede usar como alternativa la función svymultinom. Cabe aclarar que esta función no está en el repositorio de R y hay que instalarla desde un repositorio de github. A continuación, se muestra cómo se realiza la instalación y posteriormente el uso de la función:
#devtools::install_github("BS1125/CMAverse")
library(CMAverse)
model_mul2 <- svymultinom(
formula = Employment ~ Age + Sex + Region + Zone,
weights = diseno_15$variables$wk2,
data = diseno_15$variables
)
summary(model_mul2)$summarydf| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| Inactive:(Intercept) | 2.2900932 | 0.5586694 | 4.0991924 | 0.0000432 |
| Inactive:Age | 0.0247522 | 0.0100125 | 2.4721273 | 0.0135199 |
| Inactive:SexMale | -2.2194689 | 0.3162459 | -7.0181738 | 0.0000000 |
| Inactive:RegionSur | -0.4361437 | 0.4258014 | -1.0242890 | 0.3058318 |
| Inactive:RegionCentro | 0.3715204 | 0.4910265 | 0.7566198 | 0.4493733 |
| Inactive:RegionOccidente | 0.2537092 | 0.4552794 | 0.5572605 | 0.5774164 |
| Inactive:RegionOriente | 0.6175485 | 0.5157976 | 1.1972692 | 0.2313540 |
| Inactive:ZoneUrban | -0.2346034 | 0.2906882 | -0.8070620 | 0.4197338 |
| Employed:(Intercept) | 2.0928541 | 0.5427128 | 3.8562827 | 0.0001190 |
| Employed:Age | 0.0207428 | 0.0096497 | 2.1495766 | 0.0317171 |
| Employed:SexMale | -0.5563127 | 0.3053398 | -1.8219460 | 0.0686234 |
| Employed:RegionSur | -0.2789895 | 0.4106178 | -0.6794386 | 0.4969444 |
| Employed:RegionCentro | 0.2559417 | 0.4678575 | 0.5470505 | 0.5844096 |
| Employed:RegionOccidente | 0.0929049 | 0.4438767 | 0.2093035 | 0.8342342 |
| Employed:RegionOriente | 0.4704651 | 0.5050940 | 0.9314407 | 0.3517463 |
| Employed:ZoneUrban | 0.0567466 | 0.2791795 | 0.2032619 | 0.8389525 |
El hacer uso de esta función permite obtener de manera simple la predicción de las probabilidades como se muestra a continuación:
library(tidyverse)
tab_pred <- predict(model_mul2, type = "probs") %>%
data.frame()
tab_pred %>% slice(1:15)| Unemployed | Inactive | Employed |
|---|---|---|
| 0.0387205 | 0.2236639 | 0.7376156 |
| 0.0150598 | 0.5948109 | 0.3901293 |
| 0.0310297 | 0.5550689 | 0.4139014 |
| 0.0907946 | 0.1854493 | 0.7237561 |
| 0.0134342 | 0.6005134 | 0.3860523 |
| 0.0317299 | 0.5537176 | 0.4145525 |
| 0.0466669 | 0.2157326 | 0.7376005 |
| 0.0172653 | 0.5878100 | 0.3949247 |
| 0.0790766 | 0.1920709 | 0.7288525 |
| 0.0294856 | 0.2349694 | 0.7355450 |
| 0.0095218 | 0.6169869 | 0.3734913 |
| 0.0621464 | 0.2031545 | 0.7346991 |
| 0.0687482 | 0.1985744 | 0.7326773 |
| 0.0839204 | 0.1892482 | 0.7268314 |
| 0.0730133 | 0.1958009 | 0.7311858 |
Las predicciones del modelo se realizan de la siguiente manera:
diseno_15$variables %<>%
mutate(predicion = predict(model_mul2))
diseno_15 %>% group_by(Employment) %>%
summarise(Prop = survey_mean(vartype = c("se", "ci")))| Employment | Prop | Prop_se | Prop_low | Prop_upp |
|---|---|---|---|---|
| Unemployed | 0.0448040 | 0.0077955 | 0.0293683 | 0.0602398 |
| Inactive | 0.3798195 | 0.0150308 | 0.3500571 | 0.4095820 |
| Employed | 0.5753764 | 0.0130688 | 0.5494988 | 0.6012540 |
De la anterior tabla se puede observar que, el porcentaje de personas en los hogares que están desempleados es del 4.48%, inactivos del 37.9% y empleados del 57%.