1.1 Preparación del ambiente.

Cambiar el idioma a español:

Sys.setlocale("LC_ALL", "en_US.UTF-8")

## [1] "en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/C"

Cargar las librerías:

library(tidyverse)
library(MASS)
library(haven)
library(descr)

Carga de datos de Latinobarómetro, ola 2018:

datos <- read_sav("~/Dropbox/R/Latinobarometro_2018_Esp_Spss_v20190303.sav")

1.2 Variables ordinales.

Primer requisito para el uso de este tipo de modelos:

La variable explicada (\(Y\)) es de tipo categórica, medida en escala ordinal.

Características de una variable ordinal:

ordenación de las categorías de respuesta.
la ordenación permite identificar valores estructurados “mayores que”, “iguales que” y “menores que”.
la estructura de los valores pierde significado si no se respeta.
información particular del fenómeno debido a su ordenación.

Muchas de las variables utilizadas en las Ciencias Sociales son de tipo ordinal.

¿Se les ocurren ejemplos de variables ordinales?

1.3 Ejemplo de variable ordinal.

Ejemplos de este tipo de variables se encuentran en el libro de códigos de la base de datos de Latinobarómetro.

Distinción entre variable nominal y ordinal:

Libro de códigos de Latinobarómetro 2018.

¿Cuál de las variables está en escala nominal y cuál en escala ordinal?

1.4 Cálculo de probabilidades marginales y probabilidades acumuladas.

Los modelos logísticos trabajan con el cálculo de probabilidades.
Los modelos logit binominales o multinomiales calculan probabilidades marginales.
Los modelos logit ordinales calculan probabilidades acumuladas.

1.4.1 Probabilidades marginales.

Por ejemplo: A partir de la siguiente tabla se pueden calcular las probabilidades por categoría o marginales, así como probabilidades acumuladas de categorías.

descr::freq(datos $ P13STGBS.A, plot = FALSE)

## P13STGBS.A Grado de satisfacción con el funcionamiento de la democracia 
##       Frequency Percent Valid Percent
## 1          1376   6.811         7.093
## 2          3615  17.892        18.634
## 3          8879  43.947        45.768
## 4          5530  27.371        28.505
## NA's        804   3.979              
## Total     20204 100.000       100.000

Probabilidades por categoría:

\(P(Y=1) = \frac{n_1}{n} = 7.093\%\)

\(P(Y=2) = \frac{n_2}{n} = 18.634\%\)

\(P(Y=3) = \frac{n_3}{n} = 45.768\%\)

\(P(Y=4) = \frac{n_4}{n} = 1 - P(Y \leq 3) = 28.505\%\)

1.4.2 Probabilidades acumuladas.

Probabilidades acumuladas de categorías:

\(P(Y \leq 1) = \frac{n_1}{n} = 7.093\%\)

\(P(Y \leq 2) = \frac{n_1 + n_2}{n} = 7.093 + 18.634 = 25.727\%\)

\(P(Y \leq 3) = \frac{n_1 + n_2 + n_3}{n} = 7.093 + 18.634 + 45.768 = 71.495\%\)

\(P(Y \leq 4) = \frac{n_1 + n_2 + n_3 + n_4}{n} = 7.093 + 18.634 + 45.768 + 28.505 = 100\%\)

datos %>%
  ggplot(aes(P13STGBS.A)) +
  geom_bar()

2 Características del modelo logístico ordinal.

Los modelos logístico ordinales trabajan con las probabilidades acumuladas de \(Y\) cuando tiene \(k = n\) categorías.
Estos modelos buscan identificar cuál es la probabilidad acumulada de estar en diferentes combinaciones de categorías para la variable explicada (\(Y\)).
- \(F_{Y_k} = P(Y_k \leq i)\)
Estas probabilidades están condicionadas por la influencia de \(X\).
- \(F_k(X) = P(Y_k \leq i | X)\)
Las probabilidades son calculadas a partir de la razón de momios (odds).
- \(\frac{F_k(X)}{1-F_k(X)}\)
Transformación de los odds a unidades logarítmicas para el cálculo de probabilidades.
- \(L_k(X) = logit(F_k(X) = log\left(\frac{F_k(X)}{1-F_k(X)}\right)\)

2.1 Ecuación del modelo de regresión logística ordinal.

Ecuación:

\[L_k(X) = log\left(\frac{F_k(X)}{1-F_k(X)}\right) = \alpha_k - \beta_1X1 - ...-\beta_kX_k\]

Los signos negativos entre coeficientes son importantes en este tipo de modelos.
Se construirán \(k-1\) ecuaciones en total, por tanto hay \(\alpha_k\) interceptos.
Para todas las \(k-1\) ecuaciones, todos los coeficientes adoptan los mismos valores.

Razón del signo negativo en la ecuación:

Una vez que se inserta el signo negativo, una pendiente positiva (\(\beta\)) implica que mayores valores de \(X\) está asociado con un incremento en los odds de una categoría de respuesta mayor de \(Y\), en vez de una categoría de respuesta menos de \(Y\).
Esto facilita la interpretación de \(\beta\).

2.2 Preparación de las variables.

Es importante tener en cuenta que la variable explicada (\(Y\)) es de corte categórico y en escala de medida ordinal.
Este tipo de variable son reconocidas como del tipo factor en R.

Variables de interés:

Y = Satisfacción con la democracia, P13STGBS.A, k = 4 (1.- Muy satisfecho, 2.- Más bien satisfecho, 3.- No muy satisfecho, 4.- Nada satisfecho).
X1 = Satisfacción con el funcionamiento de la economía en el país, P13STGBS.B, k = 4 (1.- Muy satisfecho, 2.- Más bien satisfecho, 3.- No muy satisfecho, 4.- Nada satisfecho).

Revisión de la estructura de cada variable:

str(datos[, c(40, 41)])

## tibble [20,204 × 2] (S3: tbl_df/tbl/data.frame)
##  $ P13STGBS.A: dbl+lbl [1:20204] 3, 4, 3, 3, 3, 3, 3, 4, 3, 4, 4, 2, 4, 1, 3, 3, 2, 2...
##    ..@ label        : chr "P13STGBS.A Grado de satisfacción con el funcionamiento de la democracia"
##    ..@ format.spss  : chr "F2.0"
##    ..@ display_width: int 2
##    ..@ labels       : Named num [1:8] -4 -3 -2 -1 1 2 3 4
##    .. ..- attr(*, "names")= chr [1:8] "No preguntada" "No aplicable" "No responde" "No sabe" ...
##  $ P13STGBS.B: dbl+lbl [1:20204]  4,  4,  2,  3,  3,  3,  4,  4,  3,  4,  4,  3,  4, ...
##    ..@ label        : chr "P13STGBS.B Satisfacción con el funcionamiento de la economía en (país)"
##    ..@ format.spss  : chr "F2.0"
##    ..@ display_width: int 2
##    ..@ labels       : Named num [1:9] -5 -4 -3 -2 -1 1 2 3 4
##    .. ..- attr(*, "names")= chr [1:9] "Missing; no disponible" "No preguntado" "No aplicable" "No contesta" ...
##  - attr(*, "label")= chr "filelabel"

Se observa que las variables son del tipo “Named num”, pero se requiere transformarlas a tipo “factor” para lo que se usará el comando factor(). Aquí también se cambiaron los nombres de las variables para facilitar su manejo.

datos $ y <- factor (datos $ P13STGBS.A)
datos $ x1 <- factor (datos $ P13STGBS.B)
str(datos $ y)

##  Factor w/ 4 levels "1","2","3","4": 3 4 3 3 3 3 3 4 3 4 ...

str(datos $ x1)

##  Factor w/ 4 levels "1","2","3","4": 4 4 2 3 3 3 4 4 3 4 ...

También se requiere indicarle a R que la variable \(Y\) es de tipo ordinal, o que sus categorías están ordenadas. Para ello se usa el comando as.ordered().

datos $ y <- as.ordered(datos $ y)
str(datos $ y)

##  Ord.factor w/ 4 levels "1"<"2"<"3"<"4": 3 4 3 3 3 3 3 4 3 4 ...

Tras la ordenación, se identifica que el orden reconocido es “1”<“2”<“3”<“4”, pero como analistas sabemos que esto se debe invertir de manera que quede “4”<“3”<“2”<“1”. Para ello se usa el argumento rev dentro del comando levels().

datos $ y <- factor(datos $ y, levels = rev(levels(datos $ y)))
datos $ x1 <- factor(datos $ x1, levels = rev(levels(datos $ x1)))
str(datos $ y)

##  Ord.factor w/ 4 levels "4"<"3"<"2"<"1": 2 1 2 2 2 2 2 1 2 1 ...

str(datos $ x1)

##  Factor w/ 4 levels "4","3","2","1": 1 1 3 2 2 2 1 1 2 1 ...

2.3 Estadísticos descriptivos uni y bivariados.

Ahora es importante revisar la distribución de las variables.

Primero se solicita una distribución de frecuencias para cada variable por separado.

summary(datos[, c(396, 397)])

##     y           x1      
##  4   :5530   4   :7283  
##  3   :8879   3   :9005  
##  2   :3615   2   :2561  
##  1   :1376   1   : 766  
##  NA's: 804   NA's: 589

Posteriormente se revisa la asociación entre ambas variables mediante una tabla de contigencia.

tabla1 <- table(datos $ y, datos $ x1)
tabla1

##    
##        4    3    2    1
##   4 4442  833  139   42
##   3 1976 6136  564  122
##   2  495 1382 1538  134
##   1  174  466  262  443

Realización de un test de independencia de chi-cuadrada para evaluar la asociación entre las variables.

chisq.test (tabla1)

## 
##  Pearson's Chi-squared test
## 
## data:  tabla1
## X-squared = 12905, df = 9, p-value < 2.2e-16

A partir de los resultados y el “p-value” se puede rechazar la H0 (las variables son independientes), y se acepta que existe asociación estadística entre ambas variables.

3 Diseño del modelo logístico ordinal.

La realización de un modelo de regresión logístico ordinal se impelmenta con el comando polr()(proportional ordinal logistic regression), que es parte de la libraría MASS().
La estructura del comando en R se parece a la de los modelos de regresión lineal (OLS).

modelo1 <- polr(y ~ x1, data = datos, Hess = TRUE, na.action = na.omit, method = "logistic")
summary(modelo1)

## Call:
## polr(formula = y ~ x1, data = datos, na.action = na.omit, Hess = TRUE, 
##     method = "logistic")
## 
## Coefficients:
##     Value Std. Error t value
## x13 2.191    0.03610   60.68
## x12 3.790    0.05146   73.65
## x11 5.335    0.08875   60.11
## 
## Intercepts:
##     Value    Std. Error t value 
## 4|3   0.4170   0.0248    16.7957
## 3|2   3.2473   0.0364    89.1130
## 2|1   5.1950   0.0478   108.7592
## 
## Residual Deviance: 37845.35 
## AIC: 37857.35 
## (1056 observations deleted due to missingness)

3.1 Cálculo de “p-value” e intervalos de confianza.

Una de las limitaciones del comando polr() es que no arroja los valores del “p-value”, para ello se deben calcular aparte.

tabla_c <- coef(summary(modelo1))
p <- pnorm(abs(tabla_c[, "t value"]), lower.tail = FALSE) * 2
tabla_c <- cbind(tabla_c, "p value" = p)
tabla_c

##         Value Std. Error   t value      p value
## x13 2.1907172 0.03610285  60.67989 0.000000e+00
## x12 3.7902175 0.05146493  73.64661 0.000000e+00
## x11 5.3348968 0.08875338  60.10922 0.000000e+00
## 4|3 0.4169645 0.02482561  16.79574 2.622166e-63
## 3|2 3.2473122 0.03644038  89.11302 0.000000e+00
## 2|1 5.1949565 0.04776567 108.75920 0.000000e+00

Así se crea una tabla con los resultados del modelo1 en el que se añadió una columna con los valores del “p-value”.
Asimismo también se pueden calcular los intevalos de confianza para los coeficientes de regresión \(\beta\), para lo que se usa el comando confint().

ci <- confint(modelo1)

## Waiting for profiling to be done...

ci

##        2.5 %   97.5 %
## x13 2.120509 2.262061
## x12 3.690218 3.891991
## x11 5.160619 5.508597

A partir de los resultados del “output” del modelo de regresión se observan, inicialmente, 2 características:
- Se arrojan 3 valores de interceptos (\(\alpha\)), con los que se integran 3 ecuaciones.
- La categoría 4 de la variable \(x1\) es la categoría de referencia.
- Los coeficientes de regresión (\(\beta\)) son estadísticamente significativos para las 4 categorías evaluadas.

3.2 Ecuaciones de regresión e interpretación de los coeficientes.

Es importante recordar que estos modelos crearán varias ecuaciones de regresión, según la cantidad de \(Y_{k-1}\) categorías.

Ecuaciones: \(logit(\hat{P}(Y \leq 4)) = 0.42 - (2.19*x1_{k=3}) - (3.79*x1_{k=2}) - (5.33*x1_{k=1})\) \(logit(\hat{P}(Y \leq 3)) = 3.25 - (2.19*x1_{k=3}) - (3.79*x1_{k=2}) - (5.33*x1_{k=1})\) \(logit(\hat{P}(Y \leq 2)) = 5.19 - (2.19*x1_{k=3}) - (3.79*x1_{k=2}) - (5.33*x1_{k=1})\)

Interpretación de los coeficientes de regresión (\(\beta\)) en unidades logarítmicas.

Con un nivel de confianza del 95%, se espera un aumento de 2.19 unidades de log odds en la satisfacción con la democracia (\(Y\)) para aquellas personas que están no muy satisfecha con el funcionamiento de la economía (\(x1_{k=3}\)) en comparación con las personas que están nada satisfechas con el funcionamiento de la economía (\(x1_{k=4}\)), mantiendo el resto de las condiciones constantes.
Con un nivel de confianza del 95%, se espera un aumento de 3.79 unidades de log odds en la satisfacción con la democracia (\(Y\)) para aquellas personas que están más bien satisfechas con el funcionamiento de la economía (\(x1_{k=2}\)) en comparación con las personas que están nada satisfechas con el funcionamiento de la economía (\(x1_{k=4}\)), mantiendo el resto de las condiciones constantes.
Con un nivel de confianza del 95%, se espera un aumento de 5.33 unidades de log odds en la satisfacción con la democracia (\(Y\)) para aquellas personas que están muy satisfechas con el funcionamiento de la economía (\(x1_{k=1}\)) en comparación con las personas que están nada satisfechas con el funcionamiento de la economía (\(x1_{k=4}\)), mantiendo el resto de las condiciones constantes.

Las unidades logarítmicas se pueden transformar en unidades probabilísticas o porcentuales, lo que facilita su interpretación.

coef <- coefficients(modelo1)
exp(coef)

##        x13        x12        x11 
##   8.941624  44.266025 207.451341

Al solicitar los valores exponenciados de los coeficientes (con el comando exp()), estos se arrojan en razones de probabilidades (odds ratio) o momios de probabilidad. Y, a su vez, estos pueden ser transformados a porcentajes posteriormente.

(exp(coef) - 1) * 100

##        x13        x12        x11 
##   794.1624  4326.6025 20645.1341

A partir de estos valores se pueden interpretar los efectos de \(X1_k\) sobre \(Y\) de la siguiente manera:
- Las personas que se consideran no muy satisfechas con el funcionamiento de la economía (\(x1_{k=3}\)), comparados con las personas que están nada satisfechas con el funcionamiento de la economía (\(x1_{k=4}\)), están asociadas con un incremento de 8.9 veces con las categorías de mayor satisfacción con la democracia (\(Y\)), con un nivel de confianza del 95% y manteniendo constante el resto de las condiciones.
- Las personas que se consideran más bien satisfechas con el funcionamiento de la economía (\(x1_{k=2}\)), comparados con las personas que están nada satisfechas con el funcionamiento de la economía (\(x1_{k=4}\)), están asociadas con un incremento de 44.27 veces con las categorías de mayor satisfacción con la democracia (\(Y\)), con un nivel de confianza del 95% y manteniendo constante el resto de las condiciones.
- Las personas que se consideran muy satisfechas con el funcionamiento de la economía (\(x1_{k=1}\)), comparados con las personas que están nada satisfechas con el funcionamiento de la economía (\(x1_{k=4}\)), están asociadas con un incremento de 207.45 veces con las categorías de mayor satisfacción con la democracia (\(Y\)), con un nivel de confianza del 95% y manteniendo constante el resto de las condiciones.

3.3 Predicción del modelo.

Además de ayudar a medir el efecto de la variable explicativa sobre la variable explicada, en términos probabilísticos, estos modelos también permiten predecir las probabilidades de cada caso para ubicarse dentro de cada una de las categorías de \(Y\).

Para el cálculo de las probabilidades esperadas, se utiliza el comando predict().

pred <- predict(modelo1, datos[1:5, ], type = "prob")
round(pred, 2)

##      4    3    2    1
## 1 0.60 0.36 0.03 0.01
## 2 0.60 0.36 0.03 0.01
## 3 0.03 0.33 0.44 0.20
## 4 0.15 0.60 0.21 0.05
## 5 0.15 0.60 0.21 0.05

La tabla arroja un vistazo a las probabilidades esperadas para los primeros 1 casos en el data frame, así como las probabilidades esperadas para ubicarse en cada uno de las categorías de \(Y\).

3.4 Evaluación del modelo: matriz de confusión y tamaño de error del modelo.

Con el fin de evaluar la bondad de ajuste del modelo diseñado, comúnmente se revisan la matriz de confusión y el tamaño de error del modelo.

La matriz de confusión es una tabla cruzada en la que se combinan los casos registrados para la variable \(Y\) con respecto la ubicación esperada (o predicha) por el modelo para, así, ubicar a cada caso dentro de cada una de las categorías de \(\hat{Y}\).

pred1 <- predict(modelo1, datos)
tabla_conf <- table(pred1, datos $ y)
tabla_conf

##      
## pred1    4    3    2    1
##     4 4442 1976  495  174
##     3  833 6136 1382  466
##     2  139  564 1538  262
##     1   42  122  134  443

El supuesto de la matriz de confusión consiste en que si todos los casos fueran bien clasificados por el modelo, entonces estos se ubicarían en la diagonal descendente de la tabla (donde se intersectan las celdas 4-4, 3-3, 2-2, 1-1).
Los casos que se ubican fuera de dicha diagonal, son casos que el modelo clasificó mal. Estos son errores del modelo.
La medición del tamó de error del modelo se calcula a partir de conocer el porcentaje de casos que se ubicaron por fuera de la diagonal de la matriz de confusión.

1 - sum(diag(tabla_conf)) / sum(tabla_conf)

## [1] 0.344109

El tamaño de error del modelo es de 34% de los casos.

¿Este tamaño de error es aceptable o no?

Regresión logística ordinal.

0.1 Presentación.

1 Introducción.