4  Regresión logística multinivel

En este capítulo se introduce la regresión logística multinivel, que extiende las ideas de la regresión lineal multinivel a contextos de variables de respuesta binarias y modelos lineales generalizados. El objetivo principal es estimar probabilidades de apoyo político a nivel estatal en EE.UU., utilizando encuestas nacionales con correcciones por no respuesta y combinando información demográfica con estructura geográfica.

La regresión logística multinivel se formula agrupando coeficientes en lotes y asignando distribuciones de probabilidad a cada lote. Alternativamente, puede entenderse como la adición de términos de error que representan distintas fuentes de variación en los datos.

4.1 Estimación de opiniones estatales a partir de encuestas nacionales

• Se busca estimar opiniones a nivel estatal usando encuestas nacionales previas a elecciones.

• Se consideran factores demográficos: sexo, etnicidad, edad y educación.

• Procedimiento:

  1. Ajustar un modelo de regresión para la respuesta individual \(y\), condicionada en demografía y estado.
  2. Usar datos del Censo de EE.UU. (\(N_l\)) para ponderar las estimaciones y obtener la media poblacional por estado:

\[ \theta_{j}=\frac{\sum_{l \in j} N_{l}\theta_{l}}{\sum_{l \in j} N_{l}}, \]

donde la suma recorre las categorías demográficas del estado \(j\). Este proceso se denomina postestratificación.

4.1.1 Preparación de datos y variables (carga, filtrado y renombrado)

Objetivo del bloque: cargar datos de encuestas, definir variables clave y dejar listo el subconjunto para el modelado.

# Carga de indicadores de región (objeto 'state' del paquete 'arm')
data(state)  # contiene información auxiliar de estados
state.abbr <- c(state.abb[1:8], "DC", state.abb[9:50])
dc <- 9
not.dc <- c(1:8, 10:51)

# Vector de regiones por estado (codificación usada en el ejemplo)
region <- c(3,4,4,3,4,4,1,1,5,3,3,4,4,2,2,2,2,3,3,1,1,1,2,2,3,2,4,2,4,1,1,4,1,3,2,2,3,4,1,1,3,2,3,3,4,1,3,4,1,2,4)

# Cargar encuestas (CBS 1988)
polls <- read_dta("../ARM_Data/election88/polls.dta")

# Subconjunto de la última encuesta (survey==8 en este dataset)
polls.subset <- polls %>%
  filter(survey == 8) %>%
  rename(y = bush)  # 1 si apoya a Bush, 0 si apoya a Dukakis

# Comprobación rápida
glimpse(polls.subset)

Rows: 2,193
Columns: 10
$ org    <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ year   <dbl> 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, …
$ survey <dbl+lbl> 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,…
$ y      <dbl> NA, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, NA, 1, 1, 0, 0, 0…
$ state  <dbl> 7, 39, 31, 7, 33, 33, 39, 20, 33, 40, 31, 33, 33, 31, 22, 7, 31…
$ edu    <dbl> 3, 4, 2, 3, 2, 4, 2, 2, 4, 1, 4, 3, 4, 4, 4, 2, 2, 1, 2, 3, 3, …
$ age    <dbl> 1, 2, 4, 1, 2, 4, 2, 4, 3, 3, 2, 4, 2, 2, 2, 1, 4, 4, 1, 2, 4, …
$ female <dbl> 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, …
$ black  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, …
$ weight <dbl> 923, 558, 448, 923, 403, 317, 1532, 896, 884, 4091, 279, 1212, …

4.1.2 Modelo simple con variación demográfica y geográfica

El modelo básico es:

\[ \begin{aligned} \Pr(y_{i}=1) &= \operatorname{logit}^{-1}\!\left(\alpha_{j[i]} + \beta^{\text{female}} \cdot female_{i} + \beta^{\text{black}} \cdot black_{i}\right), \\ \alpha_{j} &\sim \mathrm{N}(\mu_{\alpha}, \sigma_{\text{state}}^{2}), \quad j=1,\ldots,51. \end{aligned} \]

Donde:

• \(y_i\): apoyo al candidato republicano (1) o demócrata (0).
• \(\alpha_j\): intercepto específico del estado \(j\).
• \(\beta^{\text{female}}\), \(\beta^{\text{black}}\): efectos de ser mujer o afroamericano.

Objetivo del bloque: ajustar el modelo logístico multinivel con intercepto aleatorio por estado e interpretar coeficientes principales.

# Modelo logístico multinivel (intercepto aleatorio por estado)
M1 <- glmer(
  y ~ black + female + (1 | state),
  data   = polls.subset,
  family = binomial(link = "logit")
)

# Resumen del modelo
summary(M1)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: y ~ black + female + (1 | state)
   Data: polls.subset

      AIC       BIC    logLik -2*log(L)  df.resid 
   2666.7    2689.1   -1329.3    2658.7      2011 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.7276 -1.0871  0.6673  0.8422  2.5271 

Random effects:
 Groups Name        Variance Std.Dev.
 state  (Intercept) 0.1692   0.4113  
Number of obs: 2015, groups:  state, 49

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.44523    0.10139   4.391 1.13e-05 ***
black       -1.74159    0.20954  -8.312  < 2e-16 ***
female      -0.09705    0.09511  -1.020    0.308    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
       (Intr) black 
black  -0.119       
female -0.551 -0.005

# Extraer coeficientes fijos y desviación entre estados
fixef(M1)

(Intercept)       black      female 
 0.44523389 -1.74159201 -0.09704648 

sqrt(unlist(VarCorr(M1)))

    state 
0.4113234 

Notas interpretativas (salida típica reportada):

• Intercepto promedio: \(0.44 \pm 0.1\)
• Efecto black: \(.1.74 \pm 0.2\)
• Efecto female: \(-0.10 \pm 0.1\)
• Desviación estándar entre estados: \(\sigma_{\text{state}} \approx 0.4\)
• No hay desviación estándar residual en logística (fijada por la familia binomial).

---

4.2 Modelo más completo con factores no anidados

Se incluyen todas las variables demográficas y sus interacciones relevantes, así como predictores a nivel estatal (región y voto previo republicano).

El modelo se expresa como:

\[ \begin{align*} \Pr(y_{i}=1)= & \operatorname{logit}^{-1}\!\Big(\beta^{0}+\beta^{\text{female}} \cdot female_{i}+\beta^{\text{black}} \cdot black_{i}+ \\ & +\beta^{\text{female.black}} \cdot female_{i}\cdot black_{i}+\alpha_{k[i]}^{\text{age}}+\alpha_{l[i]}^{\text{edu}}+\alpha_{k[i],l[i]}^{\text{age.edu}}+\alpha_{j[i]}^{\text{state}}\Big) \\ \alpha_{j}^{\text{state}} \sim & \mathrm{N}\!\left(\alpha_{m[j]}^{\text{region}}+\beta^{\text{v.prev}} \cdot v.prev_{j}, \sigma_{\text{state}}^{2}\right). \end{align*} \]

Los términos de variación adicional siguen distribuciones normales centradas en cero:

\[ \begin{aligned} \alpha_{k}^{\text{age}} &\sim N(0, \sigma_{\text{age}}^{2}), \\ \alpha_{l}^{\text{edu}} &\sim N(0, \sigma_{\text{edu}}^{2}), \\ \alpha_{k,l}^{\text{age.edu}} &\sim N(0, \sigma_{\text{age.edu}}^{2}), \\ \alpha_{m}^{\text{region}} &\sim N(0, \sigma_{\text{region}}^{2}). \end{aligned} \]

Objetivo del bloque: crear variables de interacción y predictores estatales, y ajustar el modelo ampliado.

# Cargar voto previo republicano y efectos de candidatos (para v.prev)
presvote <- read_dta("../ARM_Data/election88/presvote.dta")
v.prev   <- presvote$g76_84pr

candidate.effects <- read.table(
  "../ARM_Data/election88/candidate_effects.dat",
  header = TRUE
)

# Ajuste para D.C. y promedio de ciclos anteriores
v.prev[not.dc] <- v.prev[not.dc] +
  (candidate.effects$X76 + candidate.effects$X80 + candidate.effects$X84) / 3

# Número de categorías de educación
n.edu <- max(polls.subset$edu, na.rm = TRUE)

# Variables derivadas para el modelo completo
polls.subset <- polls.subset %>%
  mutate(
    age.edu     = n.edu * (age - 1) + edu,        # índice de interacción age x edu
    v.prev.full = v.prev[state],                  # predictor estatal expandido
    region.full = region[state]                   # región expandida a nivel individual
  )

# Modelo completo con múltiples efectos aleatorios
M2 <- glmer(
  y ~ black + female + black:female + v.prev.full +
    (1 | age) + (1 | edu) + (1 | age.edu) +
    (1 | state) + (1 | region.full),
  data   = polls.subset,
  family = binomial(link = "logit")
)

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.0357363 (tol = 0.002, component 1)

# Resumen del modelo completo
summary(M2)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: y ~ black + female + black:female + v.prev.full + (1 | age) +  
    (1 | edu) + (1 | age.edu) + (1 | state) + (1 | region.full)
   Data: polls.subset

      AIC       BIC    logLik -2*log(L)  df.resid 
   2649.5    2705.6   -1314.8    2629.5      2005 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.8649 -1.0432  0.6452  0.8401  2.8321 

Random effects:
 Groups      Name        Variance  Std.Dev. 
 state       (Intercept) 3.865e-02 0.1965901
 age.edu     (Intercept) 2.223e-02 0.1490936
 region.full (Intercept) 3.006e-02 0.1733877
 edu         (Intercept) 1.113e-02 0.1055169
 age         (Intercept) 8.635e-07 0.0009292
Number of obs: 2015, groups:  
state, 49; age.edu, 16; region.full, 5; edu, 4; age, 4

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -3.49589    1.03684  -3.372 0.000747 ***
black        -1.63511    0.32633  -5.011 5.42e-07 ***
female       -0.08939    0.09826  -0.910 0.362968    
v.prev.full   7.03156    1.85633   3.788 0.000152 ***
black:female -0.17747    0.42005  -0.423 0.672654    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) black  female v.prv.
black       -0.038                     
female      -0.052  0.182              
v.prev.full -0.991  0.022 -0.004       
black:femal  0.021 -0.762 -0.233 -0.008
optimizer (Nelder_Mead) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.0357363 (tol = 0.002, component 1)

# Varianzas de efectos aleatorios
sqrt(unlist(VarCorr(M2)))

       state      age.edu  region.full          edu          age 
0.1965901150 0.1490936479 0.1733876570 0.1055168764 0.0009292353 

Lectura rápida de la salida reportada en el texto:

• \(b_0 \approx -3.5 \ (\pm 1.0)\)
• black \(\approx -1.6\ (\pm 0.3)\)
• female \(\approx -0.1\ (\pm 0.1)\)
• v.prev.full \(\approx 7.0\ (\pm 1.7)\)
• Interacción black:female con alta incertidumbre.
• Desv. est. de efectos aleatorios \(\approx 0.0\) a \(0.2\) (según grupo).

4.2.1 Gráficos del modelo estimado (idea de linpred)

Para interpretar el modelo, se construye un predictor lineal combinado para las variables demográficas:

\[ \text{linpred}_{i} = \beta^{0}+\beta^{\text{female}} \cdot female_{i}+\beta^{\text{black}} \cdot black_{i}+ \beta^{\text{female.black}} \cdot female_{i}\cdot black_{i}+\alpha_{k[i]}^{\text{age}}+\alpha_{l[i]}^{\text{edu}}+\alpha_{k[i],l[i]}^{\text{age.edu}}. \]

Luego, la probabilidad de apoyo es:

\[ \Pr(y_{i}=1)=\operatorname{logit}^{-1}\!\left(\text{linpred}_{i}+\alpha_{j[i]}^{\text{state}}\right). \]

Objetivo del bloque: ilustrar cómo se podría construir el predictor lineal compuesto (cuando se cuente con simulaciones/coeficientes necesarios). Aquí se muestra un esquema conceptual.

# EJEMPLO CONCEPTUAL (requiere simulaciones/objetos del ajuste bayesiano si se usan)
# En este ejemplo didáctico no calculamos 'linpred' final porque
# demandaría extraer matrices completas de simulaciones.
# La idea es ilustrar el pipeline:

# pseudo-código (comentado):
# linpred <- with(polls.subset, {
#   b0  <- fixef(M2)["(Intercept)"]
#   bf  <- fixef(M2)["female"]
#   bb  <- fixef(M2)["black"]
#   bfb <- fixef(M2)["black:female"]
#
#   # Efectos aleatorios esperados (BLUPs)
#   ra  <- ranef(M2)
#   a_age     <- ra$age[as.character(age), "(Intercept)"]
#   a_edu     <- ra$edu[as.character(edu), "(Intercept)"]
#   a_ageedu  <- ra$`age.edu`[as.character(age.edu), "(Intercept)"]
#   a_state   <- ra$state[as.character(state), "(Intercept)"]
#
#   b0 + bf*female + bb*black + bfb*female*black + a_age + a_edu + a_ageedu + a_state
# })
#
# # Probabilidad predicha:
# p_hat <- plogis(linpred)

4.3 Comentario final: postestratificación y evaluación

• La postestratificación combina probabilidades por celda demográfica con frecuencias censales (\(N_l\)) para obtener estimaciones por estado:

  \[ y_{\text{state } j}^{\mathrm{pred}}=\frac{\sum_{l \in j} N_{l}\theta_{l}}{\sum_{l \in j} N_{l}}. \]

• En la comparación contra resultados electorales reales, el enfoque multinivel típicamente reduce el error absoluto medio frente a complete pooling y no pooling, logrando un equilibrio entre información global y heterogeneidad estatal.