Regresión logística y multinomial

Objetivo del manual

• Expandir la regresión lineal a otros tipos de variable respuesta

• Introducir modelos que predicen variables categóricas

• Familiarizarse con el uso de estos modelos

Paquetes a utilizar en este manual:

# instalar/cargar paquetes

sketchy::load_packages(
  c("ggplot2", 
    "viridis", 
    "nnet",
    "caret"
   )
  )

Loading required package: caret

Loading required package: lattice

1 Regresión logística

La regresión logística es una técnica utilizada cuando la variable respuesta es binaria (por ejemplo, éxito/fracaso, sí/no). En lugar de modelar directamente la variable respuesta, la regresión logística modela la probabilidad de que ocurra un evento (es decir, P(Y=1)).

\(\hat{Y\ binaria} \sim{(enlace\ logit)}\ \beta_{o} + \beta_{1} * x_{1} + ... \beta_{n} * x_{n}\)

La función de enlace más comúnmente utilizada en la regresión logística es la función logit, que transforma las probabilidades en el rango [0,1] a valores en el rango (−∞,∞), permitiendo que la combinación lineal de predictores pueda tomar cualquier valor real.

1.1 Ejemplos de uso

En un estudio sobre comportamiento juvenil, los investigadores quieren predecir si los adolescentes se involucran en conductas de riesgo (por ejemplo, consumo de drogas) en función de factores como la influencia de pares, el rendimiento escolar, y el apoyo parental.

• Variable dependiente: Conducta riesgosa (1 = Sí, 0 = No).
• Variables independientes: Influencia de pares, rendimiento escolar, apoyo parental.
• Aplicación de regresión logística: Se usa para predecir la probabilidad de que un adolescente participe en conductas riesgosas.

Un psicólogo quiere predecir si un paciente tiene depresión (Sí/No) en función de factores como el nivel de estrés, la duración del sueño, el historial familiar de trastornos mentales, y la puntuación en un test de ansiedad.

• Variable dependiente: Depresión (1 = Sí, 0 = No).
• Variables independientes: Nivel de estrés, horas de sueño, historial familiar de trastornos, puntuación de ansiedad.
• Aplicación de regresión logística: Se usa para predecir si una persona tiene depresión a partir de los predictores.

1.2 Simular datos

Primero, vamos a simular un conjunto de datos donde la respuesta sea binaria. En este ejemplo, supongamos que tenemos dos variables predictoras: x1 y x2.

# Simulación de datos
set.seed(123)  # Para reproducibilidad
n <- 100  # Número de observaciones
x1 <- rnorm(n)  # Predictor 1: variable normal
x2 <- rnorm(n)  # Predictor 2: variable normal

# Coeficientes reales para la simulación
b0 <- -0.5  # Intercepto
b1 <- 3.4  # Coeficiente de x1
b2 <- -3  # Coeficiente de x2

# Calcular la probabilidad usando la función logística
logit_p <- b0 + b1 * x1 + b2 * x2
p <- 1 / (1 + exp(-logit_p))

# Generar la respuesta binaria
y <- rbinom(n, 1, p)

# Crear un data frame
datos <- data.frame(x1 = x1, x2 = x2, y = y)

# explorar datos
head(datos)

x1	x2	y
-0.56048	-0.71041	1
-0.23018	0.25688	0
1.55871	-0.24669	1
0.07051	-0.34754	1
0.12929	-0.95162	1
1.71506	-0.04503	1

Estos graficos nos muestran las relaciones entre x1, x2 y Y:

# graficar
ggplot(datos, aes(x = x1, y = y)) + 
  geom_point(color = "#3E4A89FF")

ggplot(datos, aes(x = x2, y = y)) + 
  geom_point(color =  "#1F9E89FF")

1.3 Ajustar el modelo

Para ajustar el modelo de regresión logística en R, utilizamos la función glm() con el argumento family = binomial. glm() es una función de R básico para ajustar modelos lineales generalizados.

# Ajustar el modelo de regresión logística
modelo_log <- glm(y ~ x1 + x2, data = datos, family = binomial)

# Resumen del modelo
summary(modelo_log)

Call:
glm(formula = y ~ x1 + x2, family = binomial, data = datos)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -0.759      0.379   -2.00    0.045 *  
x1             2.885      0.699    4.13  3.7e-05 ***
x2            -3.515      0.771   -4.56  5.2e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 138.629  on 99  degrees of freedom
Residual deviance:  57.381  on 97  degrees of freedom
AIC: 63.38

Number of Fisher Scoring iterations: 7

Ahora podemos graficar los datos crudos junto a la curva de mejor ajuste. Para esto debemos estimar los valores predichos por el modelo primero:

# Crear un nuevo data frame con las predicciones para x1
nuevos_datos <-
  expand.grid(
    x1 = seq(min(datos$x1), max(datos$x1), length.out = 100),
    x2 = seq(min(datos$x2), max(datos$x2), length.out = 100)
  )

# anadir predicciones
nuevos_datos$pred_prob <-
  predict(object = modelo_log,
          newdata = nuevos_datos,
          type = "response")

# Crear el gráfico de puntos y la curva de mejor ajuste para x1
ggplot(datos, aes(x = x1, y = y)) +
  geom_point(alpha = 0.5, color = "#3E4A89FF") +  # Datos crudos
  geom_smooth(data = nuevos_datos, aes(y = pred_prob, x = x1), method = "glm", method.args = list(family = "binomial"), se = FALSE) +
  labs(y = "Probabilidad de y = 1") 

`geom_smooth()` using formula = 'y ~ x'

Warning in eval(family$initialize): non-integer #successes in a binomial glm!

# Crear el gráfico de puntos y la curva de mejor ajuste para x2
ggplot(datos, aes(x = x2, y = y)) +
  geom_point(alpha = 0.5, color = "#1F9E89FF") +  # Datos crudos
  geom_smooth(data = nuevos_datos, aes(y = pred_prob, x = x2), method = "glm", method.args = list(family = "binomial"), se = FALSE) +
  labs(y = "Probabilidad de y = 1") 

`geom_smooth()` using formula = 'y ~ x'

Warning in eval(family$initialize): non-integer #successes in a binomial glm!

1.4 Interpretación del modelo

En el modelo de regresión logística, los coeficientes estan dados como el logaritmo de los chances (log-odds) de que Y = 1:

# Obtener los log-chances
coefs <- coef(modelo_log)

coefs

(Intercept)          x1          x2 
    -0.7590      2.8846     -3.5149 

Los chances (“odds”, a veces traducido como probabilidades) se definen como la razón entre la probabilidad de la ocurrencia de un evento y la probabilidad de que ese evento no ocurra:

\[ \text{Odds} = \frac{\text{Probabilidad de éxito}}{\text{Probabilidad de fracaso}} = \frac{P(E)}{1 - P(E)} \]

El ‘log-odds’ es simplemente el logaritmo natural de ese cociente:

\[ \text{Log-Odds} = \log(\text{Odds}) = \log\left(\frac{P(E)}{1 - P(E)}\right) \]

Esto significa que por cada aumento de una unidad de x1, los ‘log-odds’ de que Y = 1 aumentan en 2.88.

Pueden interpretarse mas facilmente en términos de la razón de los chances (odds ratio). Para obtener las razones de chances, simplemente tomamos el exponente de los coeficientes:

# Obtener las razones de probabilidades
exp_coefs <- exp(coefs)

exp_coefs

(Intercept)          x1          x2 
   0.468134   17.896103    0.029752 

Esto quiere decir que el chance de Y=1 aumenta en 17.9 por cada aumento de una unidad de x1. Esto puede ser aun mas intuitivo si lo interpretamos como un porcentaje: un aumento de una unidad de x1 esta asociado con un aumento de aproximadamente 1690% en los chances de que Y=1.

Una forma intuitiva de evaluar el desempeño de modelos de clasificación (i.e. modelos que predicen categorías) es estimar el número (o proporción) de observaciones que fueron correctamente clasificadas. Este tipo de datos se pueden resumir en una matriz de confusión, que muestra cuántas observaciones fueron clasificadas correctamente y incorrectamente para cada combinación posible de categorías. Para esto usaremos la función confusionMatrix() del paquete caret:

# predecir valores
pred_vals <- predict(modelo_log, datos, type = "response")

# binarizar
pred_cat <- ifelse(pred_vals > 0.5, 1, 0)

# hacer la matriz de confusion
mat_conf <-
    confusionMatrix(factor(pred_cat), factor(datos$y))

# imprimir resultado
mat_conf$table

          Reference
Prediction  0  1
         0 43  6
         1  7 44

Esta función tambien estima otros descriptores del desempeño del modelo. Quizas el mas relevante es la precisión (accuracy), que es la proporción de observaciones que fueron clasificadas correctamente:

mat_conf$overall["Accuracy"]

Accuracy 
    0.87 

Mas adelante en el curso hablaremos de otros de estos descriptores. Una forma mas facil de interpretar la matriz de confusión es representarla gráficamente:

# convertir a data frame
conf_df <- as.data.frame(mat_conf$table)

# agregar totales por categoria
conf_df$total <-
    sapply(conf_df$Reference, function(x)
        sum(datos$y ==
                x))

# calcular proporciones
conf_df$proportion <- conf_df$Freq / conf_df$total

ggplot(conf_df, aes(x = Reference, y = Prediction, fill = proportion)) +
  geom_tile() + 
  coord_equal() + 
  scale_fill_distiller(palette = "Greens", direction = 1) + 
  geom_text(aes(label = round(proportion, 2)), color = "black", size = 3) + 
  labs(x = "Observado", y = "Predicho", fill = "Proporción") 

Como se mencionó antes la matriz de confusión nos permite observar cuántas veces el modelo predijo correctamente las categorías y cuántas veces falló. Los valores en la diagonal representan el número de aciertos (predicciones correctas), mientras que los valores fuera de la diagonal representan los errores de clasificación.

La regresión logística no hace suposiciones explícitas sobre la distribución de los residuos (errores), a diferencia de la regresión lineal, que asume que los residuos deben seguir una distribución normal. En la regresión logística, los residuos no son normales ni continuos, porque la variable dependiente es binaria (o multinomial en el caso de la regresión logística multinomial).

---

Estos modelos pueden ajustarse a estructuras mas complejas de forma similar a los modelos lineales. Por ejemplo en el siguiente modelo tenemos 2 variables predictoras y su interaccion:

# Simulación de datos
set.seed(123)  # Para reproducibilidad
n <- 100  # Número de observaciones
x1 <- rnorm(n)  # Predictor 1: variable normal
x2 <- rnorm(n)  # Predictor 2: variable normal

# Coeficientes reales para la simulación
b0 <- -0.5  # Intercepto
b1 <- 1.5  # Coeficiente de x1
b2 <- -2  # Coeficiente de x2
b3 <- 3 # interaccion
# Calcular la probabilidad usando la función logística
logit_p <- b0 + b1 * x1 + b2 * x2 + b3 * x1 * x2
p <- 1 / (1 + exp(-logit_p))

# Generar la respuesta binaria
y <- rbinom(n, 1, p)

# Crear un data frame
datos <- data.frame(x1 = x1, x2 = x2, y = y)

# explorar datos
head(datos)

x1	x2	y
-0.56048	-0.71041	0
-0.23018	0.25688	0
1.55871	-0.24669	0
0.07051	-0.34754	0
0.12929	-0.95162	1
1.71506	-0.04503	1

Estos graficos nos muestran las relaciones entre x1, x2 y Y:

# graficar
ggplot(datos, aes(x = x1, y = y)) + 
  geom_point(color = "#3E4A89FF")

ggplot(datos, aes(x = x2, y = y)) + 
  geom_point(color =  "#1F9E89FF")

… y ahora ajustamos el modelo de regresión logística conteniendo una interacción entre x1 y x2:

# Ajustar el modelo de regresión logística
modelo_log_int <- glm(y ~ x1 * x2, data = datos, family = binomial)

# Resumen del modelo
summary(modelo_log_int)

Call:
glm(formula = y ~ x1 * x2, family = binomial, data = datos)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -0.888      0.354   -2.51  0.01212 *  
x1             1.579      0.551    2.86  0.00418 ** 
x2            -2.092      0.575   -3.64  0.00028 ***
x1:x2          4.002      0.959    4.17    3e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 137.989  on 99  degrees of freedom
Residual deviance:  83.028  on 96  degrees of freedom
AIC: 91.03

Number of Fisher Scoring iterations: 7

Podemos nuevamente evaluar el desempeño del modelo con una matriz de confusión:

# predecir valores
pred_vals <- predict(object = modelo_log_int, newdata = datos, type = "response")

# binarizar
pred_cat <- ifelse(pred_vals > 0.5, 1, 0)

# hacer la matriz de confusion
mat_conf <-
    confusionMatrix(factor(pred_cat), factor(datos$y))

# imprimir resultado
mat_conf$table

          Reference
Prediction  0  1
         0 44 15
         1 10 31

mat_conf$overall["Accuracy"]

Accuracy 
    0.75 

# convertir a data frame
conf_df <- as.data.frame(mat_conf$table)

# agregar totales por categoria
conf_df$total <-
    sapply(conf_df$Reference, function(x)
        sum(datos$y ==
                x))

# calcular proporciones
conf_df$proportion <- conf_df$Freq / conf_df$total

ggplot(conf_df, aes(x = Reference, y = Prediction, fill = proportion)) +
  geom_tile() + 
  coord_equal() + 
  scale_fill_distiller(palette = "Greens", direction = 1) + 
  geom_text(aes(label = round(proportion, 2)), color = "black", size = 3) + 
  labs(x = "Observado", y = "Predicho", fill = "Proporción") 

1.5 Ejercicio 1

Los datos de sobrevivencia de pasajeros del Titanic son ampliamente utilizado para ilustrar modelos de predicción. Contiene información sobre los pasajeros del Titanic, incluyendo variables como:

• Sobrevivencia (Survived): indica si el pasajero sobrevivió o no
• Clase (Pclass): clase del pasajero (1, 2 o 3)
• Sexo (Sex): género del pasajero
• Edad (Age): edad del pasajero

Puedes cargarlo directamente en R usando el paquete titanic o manipularlo desde el conjunto de datos base. El siguiente código carga y da formato a los datos para que puedan ser utilizados en una regresión logística:

# cargar datos
data("Titanic")

# dar formato con una observacion por fila
datos_titanic <- as.data.frame(Titanic)
datos_titanic <- datos_titanic[datos_titanic$Freq > 0, ]

datos_tab_titanic <- do.call(rbind, lapply(1:nrow(datos_titanic), function(x) datos_titanic[rep(x, datos_titanic$Freq[x]),]))
datos_tab_titanic$Freq <- NULL

# explorar datos
head(datos_tab_titanic, 10)

	Class	Sex	Age	Survived
3	3rd	Male	Child	No
3.1	3rd	Male	Child	No
3.2	3rd	Male	Child	No
3.3	3rd	Male	Child	No
3.4	3rd	Male	Child	No
3.5	3rd	Male	Child	No
3.6	3rd	Male	Child	No
3.7	3rd	Male	Child	No
3.8	3rd	Male	Child	No
3.9	3rd	Male	Child	No

El siguiente codigo ajusta un modelo de regresión logística a estos datos con la variable “clase” como único predictor:

# Ajustar el modelo de regresión logística
modelo_log <- glm(Survived ~ Class, data = datos_tab_titanic, family = binomial)

# Resumen del modelo
summary(modelo_log)

Call:
glm(formula = Survived ~ Class, family = binomial, data = datos_tab_titanic)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)    0.509      0.115    4.44  8.8e-06 ***
Class2nd      -0.856      0.166   -5.16  2.5e-07 ***
Class3rd      -1.596      0.144  -11.11  < 2e-16 ***
ClassCrew     -1.664      0.139  -11.97  < 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: 2769.5  on 2200  degrees of freedom
Residual deviance: 2588.6  on 2197  degrees of freedom
AIC: 2597

Number of Fisher Scoring iterations: 4

Podemos evaluar el desempeño del modelo con una matriz de confusión:

# predecir valores
pred_vals <- predict(object = modelo_log, newdata = datos_tab_titanic, type = "response")

# binarizar
pred_cat <- ifelse(pred_vals > 0.5, "Yes", "No")

# hacer la matriz de confusion
mat_conf <-
    confusionMatrix(factor(pred_cat), factor(datos_tab_titanic$Survived))

# imprimir resultado
mat_conf$table

          Reference
Prediction   No  Yes
       No  1368  508
       Yes  122  203

# convertir a data frame
conf_df <- as.data.frame(mat_conf$table)

# agregar totales por categoria
conf_df$total <-
    sapply(conf_df$Reference, function(x)
        sum(datos_tab_titanic$Survived ==
                x))

# calcular proporciones
conf_df$proportion <- conf_df$Freq / conf_df$total

# graficar
ggplot(conf_df, aes(x = Reference, y = Prediction, fill = proportion)) +
  geom_tile() + 
  coord_equal() + 
  scale_fill_distiller(palette = "Greens", direction = 1) + 
  geom_text(aes(label = round(proportion, 2)), color = "black", size = 3) + 
  labs(x = "Observado", y = "Predicho", fill = "Proporción") 

# precision
mat_conf$overall["Accuracy"]

Accuracy 
 0.71377 

• 1.1 Construya un modelo de regresión logística que prediga la sobrevivencia de los pasajeros del Titanic en función de su género. ¿Cuál es la precisión del modelo?

• 1.2 Añada una variable predictora mas a la vez (e.g. Survived ~ Sex; Survived ~ Sex + Class; Survived ~ Sex + Class + Age) y a cada nuevo modelo evalue su precisión.

• 1.3 ¿Qué modelo tiene la mejor precisión? ¿Qué variables parecen ser más importantes para predecir la sobrevivencia de los pasajeros del Titanic?

---

2 Regresión multinomial

La regresión multinomial es una técnica estadística utilizada cuando la variable dependiente tiene más de dos categorías y no hay un orden natural entre ellas. A diferencia de la regresión logística binaria, que se utiliza cuando la variable dependiente es dicotómica (es decir, tiene solo dos posibles resultados), la regresión multinomial maneja resultados con múltiples categorías nominales. Intenta modelar la probabilidad de pertenecer a cada categoría de una variable dependiente con múltiples clases. También se puede usar para predecir la categoría más probable para nuevas observaciones basadas en los valores de las variables independientes.

\(\hat{Y}\text{multinomial} \sim{(enlace\ logit)}\ \beta_{o} + \beta_{1} * x_{1} + ... \beta_{n} * x_{n}\)

La variable dependiente debe ser categórica y contener más de dos categorías. Además se asume que no existe un orden inherente entre las categorías de la variable dependiente (en cuyo caso, se utilizaría una regresión ordinal) y que las categorías son exhaustivas y mutuamente excluyentes; es decir, cada observación pertenece a una sola categoría.

2.1 Ejemplos de uso

Un psiquiatra quiere predecir a qué tipo de trastorno mental padece un paciente (por ejemplo, depresión, ansiedad o trastorno bipolar) en función de factores como el historial familiar, nivel de estrés, horas de sueño, y puntuaciones en diversos tests psicológicos.

• Variable dependiente: Tipo de trastorno mental (1 = Depresión, 2 = Ansiedad, 3 = Trastorno Bipolar).
• Variables independientes: Historial familiar, nivel de estrés, horas de sueño, puntuaciones en tests psicológicos.
• Aplicación de la regresión multinomial: Se utiliza para predecir la probabilidad de pertenecer a una de las tres categorías diagnósticas.

Una empresa de alimentos quiere predecir qué tipo de envase (plástico, vidrio o biodegradable) preferirán los consumidores para sus productos, en función de factores como el precio del producto, la percepción ambiental y la facilidad de uso.

• Variable dependiente: Tipo de envase preferido (1 = Plástico, 2 = Vidrio, 3 = Biodegradable).
• Variables independientes: Precio del producto, percepción ambiental, facilidad de uso.
• Aplicación de la regresión multinomial: Permite predecir qué tipo de envase es más probable que prefieran los consumidores, basado en sus actitudes y características.

2.2 Simular datos

Para ilustrar cómo funciona la regresión multinomial, vamos a simular un conjunto de datos con tres categorías posibles para la variable dependiente.

# Establecer una semilla para reproducibilidad
set.seed(123)

# Simular variables predictoras
n <- 500  # número de observaciones
x1 <- rnorm(n)
x2 <- rnorm(n)

# Probabilidades para cada categoría
p_y1 <- exp(1 + 2 * x1 - 5 * x2) / (1 + exp(1 + 2 * x1 - 5 * x2) + exp(-2 + x1 + 10 * x2))
p_y2 <- exp(-2 + x1 + 10 * x2) / (1 + exp(1 + 2 * x1 - 5 * x2) + exp(-2 + x1 + 10 * x2))
p_y3 <- 1 - p_y1 - p_y2  # Probabilidad de la tercera categoría

# Asignar categorías de acuerdo con las probabilidades
y <- sapply(1:n, function(i) sample(1:3, size = 1, prob = c(p_y1[i], p_y2[i], p_y3[i])))

# Crear un dataframe con las variables simuladas
datos_multin <- data.frame(y = factor(y, labels = c("a", "b", "c")), x1 = x1, x2 = x2)

# Visualizar las primeras filas de los datos simulados
head(datos_multin)

y	x1	x2
a	-0.56048	-0.60189
a	-0.23018	-0.99370
b	1.55871	1.02679
b	0.07051	0.75106
a	0.12929	-1.50917
a	1.71506	-0.09515

2.3 Explorar datos

Categorias de Y en el espacio X1 vs X2

ggplot(datos_multin, aes(x = x1, y = x2, color = y)) +
  geom_point() +
  scale_color_viridis_d(option = "G", begin = 0.4, end = 0.9) +
  labs(x = "x1", y = "x2", color = "Categoría")

2.4 Ajustar el modelo

El siguiente paso es ajustar un modelo de regresión multinomial a los datos simulados utilizando la función multinom() del paquete nnet.

# Ajustar el modelo de regresión multinomial
modelo_multinom <- multinom(y ~ x1 + x2, data = datos_multin)

# weights:  12 (6 variable)
initial  value 549.306144 
iter  10 value 124.300458
iter  20 value 123.855606
final  value 123.845494 
converged

# Resumen del modelo
summary(modelo_multinom)

Call:
multinom(formula = y ~ x1 + x2, data = datos_multin)

Coefficients:
  (Intercept)       x1     x2
b     -2.7457 -0.72491 14.052
c     -1.0776 -2.06886  4.611

Std. Errors:
  (Intercept)      x1      x2
b     0.49286 0.33465 1.67940
c     0.25239 0.31713 0.75636

Residual Deviance: 247.69 
AIC: 259.69 

Podemos agregar los valores de p para cada tamaño de efecto de esta forma:

summ_modelo_multinom <- summary(modelo_multinom)

coefs <- summ_modelo_multinom$coefficients

# anadir valores de p
p_x1 <- (1 - pnorm(abs(coefs[, 2]), 0, 1)) * 2
p_x2 <- (1 - pnorm(abs(coefs[, 3]), 0, 1)) * 2

coefs <- cbind(coefs, p_x1, p_x2)

coefs

  (Intercept)       x1     x2    p_x1       p_x2
b     -2.7457 -0.72491 14.052 0.46851 0.0000e+00
c     -1.0776 -2.06886  4.611 0.03856 4.0081e-06

 

 

Interpretación del modelo

 

Cuadro con los coeficientes (estimados):

  (Intercept)       x1     x2    p_x1       p_x2
b     -2.7457 -0.72491 14.052 0.46851 0.0000e+00
c     -1.0776 -2.06886  4.611 0.03856 4.0081e-06

• El modelo encontró que \(\beta_1b\) es -0.72491 y que no es significativamente diferente de 0 (p = 0.46851). Podemos expresarlo como chances al calcular exp(coefs[1, 2]) que es 0.48. Es decir, la probabilidad de pertenecer a la categoría “b” aumenta en 0.48 comparado con la probabilidad de pertenecer categoría “a” con un aumento de una unidad de x1. En este caso el aumento es menor a uno lo que quiere decir que la probabilidad de ocurrencia de “b” es menor a la de “a” con un aumento de x1. Si mas bien calculamos 1 / exp(coefs[1, 2]) podemos decir que la probabilidad de pertenecer a la categoría “a” aumenta en 2.06 comparado con la probabilidad de pertenecer categoría “b” con un aumento de una unidad de x1.

• El modelo encontró que \(\beta_1c\) es -2.06886 y que es significativamente diferente de 0 (p = 0.03856). Al calcular los chances (exp(coefs[2, 2])) podemos decir que la probabilidad de pertenecer a la categoría “c” aumenta en 0.48 comparado con la probabilidad de pertenecer a la categoría “a”.

• También se encontró que el \(\beta_2b\) es 14.05159 y que es significativamente diferente de 0 (p = 0). Al calcular los chances (exp(coefs[1, 3])) podemos decir que la probabilidad de pertenecer a la categoría “b” aumenta en 1.26627^{6} comparado con la probabilidad de pertenecer a la categoría “a”.

• También se encontró que el \(\beta_2c\) es 4.61096 y que también es significativamente diferente de 0 (p = 4.0081^{-6}). Al calcular los chances (exp(coefs[2, 3])) podemos decir que la probabilidad de pertenecer a la categoría “c” aumenta en 100.58 comparado con la probabilidad de pertenecer a la categoría “a”.

2.5 Visualización de resultados

Una forma útil de visualizar los resultados del modelo es graficar las probabilidades predichas para cada categoría en función de los predictores. Esto permite observar cómo cambian las probabilidades a medida que cambian los valores de las variables independientes.

En el siguiente ejemplo, graficaremos las probabilidades predichas de pertenecer a cada categoría con respecto a una de las variables independientes (x1):

# Crear una secuencia de valores para x1
x1_nuevo <- seq(min(datos_multin$x1), max(datos_multin$x1), length.out = 100)

# Generar un nuevo conjunto de datos con x1 variando y x2 fijo
nuevo_data <- data.frame(x1 = x1_nuevo, x2 = 0)

# Obtener las probabilidades predichas para el nuevo conjunto de datos
prob_predichas <- predict(modelo_multinom, newdata = nuevo_data, type = "probs")

# definir colores
cols <- mako(3, begin = 0.4, end = 0.9)

# Graficar las probabilidades predichas
matplot(x1_nuevo, prob_predichas, type = "l", lty = 1, col = cols,
        xlab = "x1", ylab = "Probabilidad predicha",
        main = "Probabilidades predichas para cada categoría")
legend("topright", legend = levels(datos_multin$y), col = cols, lty = 1)

Este gráfico muestra cómo varían las probabilidades predichas de pertenecer a cada categoría conforme el valor de x1 cambia, con x2 fijo en 0. Las líneas de diferentes colores corresponden a las diferentes categorías de la variable dependiente.

2.5.1 Matriz de confusión

Al ser este en escencia un modelo de clasificación, podemos evaluar su desempeño con una matriz de confusión:

# predecir valores
pred_cat <- predict(modelo_multinom, datos_multin, type = "class")

# hacer la matriz de confusion
mat_conf <-
    confusionMatrix(factor(pred_cat), factor(datos_multin$y))

# imprimir resultado
mat_conf$table

          Reference
Prediction   a   b   c
         a 238   7  17
         b   7 186  13
         c   7   1  24

Esta función tambien estima otros descriptores del desempeño del modelo. Quizas el mas relevante es la precisión (accuracy), que es la proporción de observaciones que fueron clasificadas correctamente:

mat_conf$overall["Accuracy"]

Accuracy 
   0.896 

Mas adelante en el curso hablaremos de otros de estos descriptores. Una forma mas facil de interpretar la matriz de confusión es representarla gráficamente:

# convertir a data frame
conf_df <- as.data.frame(mat_conf$table)

# agregar totales por categoria
conf_df$total <-
    sapply(conf_df$Reference, function(x)
        sum(datos_multin$y ==
                x))

# calcular proporciones
conf_df$proportion <- conf_df$Freq / conf_df$total

ggplot(conf_df, aes(x = Reference, y = Prediction, fill = proportion)) +
  geom_tile() + 
  coord_equal() + 
  scale_fill_distiller(palette = "Greens", direction = 1) + 
  geom_text(aes(label = round(proportion, 2)), color = "black", size = 3) + 
  labs(x = "Observado", y = "Predicho", fill = "Proporción") 

2.6 Ejercicio 2

Para este ejercicio puedo utilizar el código del ejercicio 1 para la evaluación de los modelos.

2.1 Utilice los datos iris para ajustar una regresión multinomial utilizando como respuesta el nombre de la especie (“Species”) y el largo del pétalo como predictor. ¿Cuál es la precisión del modelo?

2.2 Añada un segundo predictor (e.g. largo del sépalo) y evalúe la precisión del modelo.

2.3 Ahora use todas las variables como predictoras y evalúe la precisión del modelo.

2.4 ¿Cual es el modelo con el mejor desempeño? ¿Qué variables parecen ser más importantes para predecir la especie de la flor?

---

Información de la sesión

R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Costa_Rica
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] caret_6.0-94      lattice_0.20-45   nnet_7.3-17       viridis_0.6.5    
[5] viridisLite_0.4.2 ggplot2_3.5.1     knitr_1.48       

loaded via a namespace (and not attached):
 [1] gtable_0.3.5         xfun_0.47            htmlwidgets_1.6.4   
 [4] recipes_1.0.9        remotes_2.5.0        vctrs_0.6.5         
 [7] tools_4.3.2          generics_0.1.3       stats4_4.3.2        
[10] parallel_4.3.2       proxy_0.4-27         tibble_3.2.1        
[13] fansi_1.0.6          ModelMetrics_1.2.2.2 pkgconfig_2.0.3     
[16] Matrix_1.6-5         data.table_1.14.10   RColorBrewer_1.1-3  
[19] lifecycle_1.0.4      farver_2.1.2         compiler_4.3.2      
[22] stringr_1.5.1        munsell_0.5.1        codetools_0.2-18    
[25] sketchy_1.0.3        htmltools_0.5.8.1    class_7.3-20        
[28] yaml_2.3.10          prodlim_2023.08.28   pillar_1.9.0        
[31] crayon_1.5.3         MASS_7.3-55          gower_1.0.1         
[34] iterators_1.0.14     rpart_4.1.16         foreach_1.5.2       
[37] nlme_3.1-155         parallelly_1.38.0    lava_1.7.3          
[40] tidyselect_1.2.1     packrat_0.9.2        digest_0.6.37       
[43] stringi_1.8.4        future_1.34.0        reshape2_1.4.4      
[46] purrr_1.0.2          dplyr_1.1.4          listenv_0.9.1       
[49] labeling_0.4.3       splines_4.3.2        fastmap_1.2.0       
[52] grid_4.3.2           colorspace_2.1-1     cli_3.6.3           
[55] magrittr_2.0.3       survival_3.2-13      utf8_1.2.4          
[58] e1071_1.7-16         future.apply_1.11.2  withr_3.0.1         
[61] scales_1.3.0         timechange_0.2.0     lubridate_1.9.3     
[64] rmarkdown_2.28       globals_0.16.3       timeDate_4032.109   
[67] gridExtra_2.3        evaluate_1.0.0       hardhat_1.3.0       
[70] mgcv_1.8-39          rlang_1.1.4          Rcpp_1.0.13         
[73] glue_1.8.0           xaringanExtra_0.8.0  pROC_1.18.5         
[76] ipred_0.9-14         rstudioapi_0.16.0    jsonlite_1.8.9      
[79] R6_2.5.1             plyr_1.8.9