hllinas

1 Librerías

library(aplore3)     #Base de datos para los ejemplos
library(lsm)         #Base de datos para ejemplos y estimaciones del Log-verosimilitud
library(tidyverse)   #Incluye a dplyr y ggplot2

2 Introducción

Los métodos de regresión se han convertido en un componente integral de cualquier análisis de datos preocupado por describir la relación entre una variable de respuesta y una o variables más explicativas. Muy a menudo, la variable de resultado es discreta, tomando un valor de dos o más valores posibles. El modelo de regresión logística es el más modelo de regresión de mayor uso frecuente para el análisis de estos datos.

En el documento Rpbus :: Modelos lineales generalizados se explicó que estos modelos hacen parte de los modelos lineales generalizados y en Rpbus :: Regresión Logística binaria se explicó el caso binario. Allí se pueden ver los detalles correspondientes. En este documento se explicará el caso multinomial, en donde la variable de respuesta toma uno de tres valores posibles. Para conocer con profundidad estos modelos, también es importante estudiar los siguientes cuatro tipos de modelos:

  • Modelo de Bernoulli.

  • Modelo completo.

  • Modelo nulo.

  • Modelo saturado.

En el documento Rpbus :: Modelos completo, nulo y saturado se describieron sus propiedades, con los ejemplos correspondientes y en los documentos Rpbus :: Modelos logísticos (estimaciones), Rpbus :: Modelos logísticos (intervalos de confianza) y Rpbus :: Modelos logísticos (pruebas de comparación) se detalló todo lo relacionado con las estimaciones, intervalos de confianza para los parámetros logísticos y pruebas de comparación nde modelos. En este documento, se utilizarán las notaciones utilizados allá, así como los resultados encontrados en los ejemplos aplicados en esos documentos. A pesar de ello, se hará una breve descripción de los cuatro modelos mencionados arriba, los cuales serán base para la teoría que se explicará posteriormente.

3 Datasets

Para las aplicaciones, se utilizará la base datos hbsdemo (UCLA: Stat Consulting Group (2021)):

library(repmis)
source_data("https://github.com/hllinas/DatosPublicos/blob/main/hsbdemo.Rdata?raw=false")
Datos <- hsbdemo
attach(Datos)
names(Datos)
## [1] "hsbdemo"

El conjunto de datos contiene variables sobre 200 estudiantes. Los estudiantes que ingresan a la escuela secundaria hacen la elección de un programa de tres posibles: general, vocacional y académico. Su elección puede ser modelada usando algunas variables predictoras. A continuación, se describen las variables:

  1. Student y ID, el estudiante y su código de identificación.

  2. gender, el género del estudiante: female, male.

  3. ses, el estrato socioeconómico: low, middle, high.

  4. schtype, el tipo de escuela: private, public.

  5. prog, el tipo de programa elegido por el estudiante: 0=general, 1=vocacional y 2=académico, la cual será nuestra variable de respuesta.

  6. read, write, math, science y socst son variables continuas que representan los puntajes en lectura, esscritura, matemática, ciencia y sociales, respectivamente.

  7. Honors, estado de honores: enrolled, not enrolled.

  8. awards, número de premios recibidos: de 0 a 9.

  9. cid, puntaje no especificado: de 0 a 20.

  10. prog0, variable binaria: 1=si prog=general, 0= de otro modo.

  11. prog1, variable binaria: 1=si prog=vocation, 0= de otro modo.

  12. prog2, variable binaria: 1=si prog=academic, 0= de otro modo.

Las primeras 10 observaciones son:

Student id gender ses schtyp prog read write math science socst honors awards cid prog0 prog1 prog2
1 45 female low public vocation 34 35 41 29 26 not enrolled 0 1 0 1 0
2 108 male middle public general 34 33 41 36 36 not enrolled 0 1 1 0 0
3 15 male high public vocation 39 39 44 26 42 not enrolled 0 1 0 1 0
4 67 male low public vocation 37 37 42 33 32 not enrolled 0 1 0 1 0
5 153 male middle public vocation 39 31 40 39 51 not enrolled 0 1 0 1 0
6 51 female high public general 42 36 42 31 39 not enrolled 0 1 1 0 0
7 164 male middle public vocation 31 36 46 39 46 not enrolled 0 1 0 1 0
8 133 male middle public vocation 50 31 40 34 31 not enrolled 0 1 0 1 0
9 2 female middle public vocation 39 41 33 42 41 not enrolled 0 1 0 1 0
10 53 male middle public vocation 34 37 46 39 31 not enrolled 0 1 0 1 0

Algunas estadísticas descriptivas relacionadas con nuestra variable de respuesta de interés (prog) son las que se indican en el ejemplo siguiente.

Example 3.1 Tenemos:

  1. Tabulación cruzada entre prog y ses.
Frecuencia <- xtabs(~ prog + ses, data = Datos)
Frecuencia <- cbind(Frecuencia, apply(Frecuencia,1,sum))
Frecuencia <- rbind(Frecuencia, apply(Frecuencia,2,sum))
rownames(Frecuencia)[4] <- c('Total')
colnames(Frecuencia)[4] <- c('Total')
Frecuencia
1 2 3 4
prog_ses high low middle Total
academic 42 19 44 105
general 9 16 20 45
vocation 7 12 31 50
Total 58 47 95 200
  1. Medidas estadísticas de write dentro de cada nivel de prog:
Datos %>%   group_by(prog)  %>% 
            summarise(n = length(as.numeric(write)),
             Minimo = min(as.numeric(write)),
             Maximo = max(as.numeric(write)),
             Promedio = mean(as.numeric(write)),
             Desviacion = sd(as.numeric(write)))
1 2 3 4 5 6
Programa Total Mínimo Máximo Promedio Desviación
academic 105 33 67 56.2571 7.9433
general 45 31 67 51.3333 9.3978
vocation 50 31 67 46.76 9.3188

4 Paquetes R para modelado multinomial

Existe una amplia gama de paquetes R disponibles para el modelado multinomial, algunos de los cuales incluso permiten la incorporación de efectos aleatorios (este tema de los efectos aleatorios no se explicar[a en este documento).

4.0.1 Efectos fijos

Los siguientes paquetes se utilizan con frecuencia y se citan para situaciones en las que solo se utilizan efectos fijos:

  1. La función multinom en la librería nnet (Venables & Ripley 2002). Es el paquete de mayor uso y se basa en redes neuronales.

  2. La función vglm en la librería VGAM (Yee 2021). Se basa en modelos con vectores generalizados aditivos, con una amplia gama de distribuciones para la variable de respuesta.

  3. La función polr en la librería MASS (Venables & Ripley 2002). También es útil para Variables de respuesta categórica ordenadas.

  4. LA función glmnet en la librería glmnet (Friedman et al. 2010). Utiliza métodos de contracción donde las estimaciones de coeficientes se pueden reducir a cero durante el procesoel modelo adecuado.

4.0.2 Efectos fijos y aleatorios

Curiosamente, un número creciente de paquetes también se ajusta a modelos multinomiales tanto con efectos como efectos aleatorios (modelos multinomiales de efectos mixtos), por ejemplo:

  1. La función bayesx en la librería R2BayesX (Belitz et al. 2017, Umlauf et al. 2015), utiliza un enfoque bayesiano para la estimación de parámetros.

  2. La función npmlt en la librería mixcat (Papageorgiou y Hinde 2012).

  3. La función clmm en la librería *ordinal8 (Christensen 2018), para categorías ordenadas de la variable de respuesta.

5 Modelo de Bernoulli

La variable de interés \({Y}\) puede asumir tres niveles: \(0\), \(1\) o \(2\). Para cada \(r = 0, 1, 2\), sea

\[p_r: = P (Y = r)\]

la probabilidad de que \(Y\) tome el valor \(r\).

Haciendo \(n\) observaciones independientes de \({ Y}\), se obtiene la muestra \(Y=(Y_1, \ldots, Y_n)\) con los datos \(y_i\in\{0,1,2\}\), \(i=1, \cdots ,n\), donde \(y_i\) es un posible valor de \(Y_i\), las cuales son independientes entre sí.

Para construir la función de verosimilitud, debemos crear tres variables binarias con valores \(0\) y \(1\) e independientes, de la siguiente manera: \[U_{ri}=\left\{ \begin{array}{ll} 1, & \hbox{if $Y_i=r$;} \\ 0, & \hbox{de otra forma.} \end{array} \right. \]

donde \(r=0,1,2\) y \(i=1, \ldots, n\). Observe que \(U_{ri} \sim \mathcal{B}(1,p_{ri})\), siendo \(p_{ri}=P(Y_i=r)\).

En términos de las variables \(U_{ri}\), las variables muestrales serán \(Y_i = (U_{0i}, U_{1i}, U_{2i})\), con valores

\[y_i = (u_{0i}, u_{1i}, u_{2i}),\] siendo \(\sum\limits_{r=0}^2 u_{ri} =1\), para \(i\) fijo.

Example 5.1 En la Tabla 5.1 se ilustra hipotéticamente un conjunto de datos que contiene las las variables \(U_r\).

Table 5.1: Ilustración de un conjunto de datos agrupado en el modelo básico
1 2 3 4 5
\(Y\) \(X_1\) \(U_0\) \(U_1\) \(U_2\)
1 Bajo 0 1 0
0 Bajo 1 0 0
2 Bajo 0 0 1
0 Bajo 1 0 0
2 Bajo 0 0 1
1 Bajo 0 1 0
2 Bajo 0 0 1
1 Alto 0 1 0
2 Alto 0 0 1
0 Alto 1 0 0
1 Alto 0 1 0
1 Alto 0 1 0
General. \(Y\) es la variable de respuesta; \(X_1,\) una variable explicativa; \(U_r\) son las variables de respuestas dicotomizadas en el nivel \(r\).

Theorem 5.1 (Log-verosimilitud) Fijando \(y=(y_1, \cdots ,y_n)^T\), se llega a un modelo estadístico donde:

\[P(Y_i=y_i) = \prod^2_{r=0} p_{ri}^{u_{ri}}, \quad i=1, \ldots, n\]

Además, el logaritmo de la función de verosimilitud será:

\[\begin{equation} {\cal L}(p) = \sum^n_{i=1}[u_{0i} \ln p_{0i} + u_{1i} \ln p_{1i} + (1-u_{0i}-u_{1i}) \ln(1-p_{0i}-p_{1i})], \tag{5.1} \end{equation}\]

evaluada en el parámetro \(2n\)-dimensional

\[p=(p_{01}, p_{11}, \ldots, p_{0n}, p_{1n})^T\]

Remark. Hay varias situaciones que se pueden presentar en un modelo básico, como el que fue descrito anteriormente. Se dice que éste se puede identificar como alguno de los siguientes modelos: completo, nulo o saturado.

6 Modelo completo

Definition 6.1 El modelo completo es caracterizado por el supuesto de que todos \(p_{ri}\) (con \(r = 0, 1, 2\) y \(i = 1, \ldots, n\)) son considerados como parámetros.

Remark. El siguiente teorema describe las estimaciones en este modelo:

Theorem 6.1 (Estimaciones en el modelo completo) En el modelo completo, las ML-estimaciones de \(p_{ri}\) son \(\widehat{p}_{ri} = U_{ri}\) con valores \(\widehat{p}_{ri}=u_{ri}\) para \(r = 0, 1, 2\) y \(i = 1,\ldots, n\). Además, la estimación de la función de verosimilitud, de su logaritmo y de la llamada desviación tienen los siguientes valores:

\[{\cal L}(\widehat{p}) = 0, \quad L(\widehat{p})= 1 \quad \mbox{y}\quad -2\,{\cal L}(\widehat{p})=0\]

Example 6.1 Para los datos del archivo hsbdemo, en el modelo completo, se tiene que \({\cal L}(\widehat{p}) = -2\; {\cal L}(\widehat{p})= 0\) y \(L(\widehat{p})=1\). En R, se pueden verificar así:

Lp_i <-  ifelse(prog0 == 0 | prog1 == 0 | prog0+prog1 == 1, 0, 
                prog0*log(prog0)+prog1*log(prog1)+(1-prog0-prog1)*log(1-prog0-prog1)) 

LogCompleto <- sum(Lp_i)
L_Completo <- exp(LogCompleto)
DevLogCompleto <- -2*LogCompleto

7 Modelo nulo

Definition 7.1 El modelo nulo es caracterizado por el supuesto de que todos los \(p_{ri}\) (\(i=1, \cdots ,n\)) son considerados iguales; es decir, se tienen dos parámetros \(p_0\) y \(p_1\).

Theorem 7.1 (Log-verosimilitud en el modelo nulo) En este caso, (5.1) será:

\[\begin{eqnarray} {\cal L}(p)&=& n[\overline{u}_0\ln p_0 + \overline{u}_1\ln p_1+ (1- \overline{u}_0-\overline{u}_1)\ln (1-p_0-p_1)] \tag{7.1} \end{eqnarray}\]

El siguiente teorema describe las estimaciones en este modelo:

Theorem 7.2 (Estimaciones en el modelo nulo) En el modelo nulo, la ML-estimación de \(p_r\) es \(\hat{p}_r=\overline{U}_r\) con valor \(\hat{p}_r=\overline{u}_r\). Además,

\[{\cal L}(\overline{p})<0 \qquad \mbox{si y sólo si}\qquad 0 <\overline{u}_0 + \overline{u}_1 <1\]

8 Ejemplo (modelo nulo)

8.0.1 Enunciado

Para los datos del archivo hsbdemo, en el modelo nulo, encuentre:

  1. Las estimaciones de los parámetros del modelo nulo, utilizando las fórmulas planteadas en el documento.

  2. Esas estimaciones con la función multinom :: nnet.

  3. Esas estimaciones con la función vglm :: VGAM.

  4. El logaritmo de la función de versosimilitud \({\cal L}(\widehat{p})\) y el valor de la desviación residual \(-2{\cal L}(\widehat{p})\), utilizando las fórmulas planteadas en el documento.

  5. Las estimaciones pedidas en (d), con la función multinom :: nnet.

  6. Las estimaciones pedidas en (d), con la función vglm :: VGAM.

8.0.2 Solución parte (a)

Sea \(u_0:=\) prog0 y \(u_1:=\) prog1. Debido a que

\[\widehat{p}_0 \;=\; \overline{u_0}\;=\; 45/200 \;=\; 0.225, \qquad \widehat{p}_1 \;=\; \overline{u_1}\;=\; 50/200 \;=\; 0.25,\]

entonces, las estimaciones de los parámetros se pueden reunir en el siguiente vector:

\[\widehat{p} \;=\; (\widehat{p}_0, \widehat{p}_1)^T \;=\; (0.225, 0.25)^T\]

8.0.3 Solución parte (b)

La función summary::multinom() del paquete nnet entrega unas estimaciones que nos pueden ayudar a calcular lo solicitado, en especial, los valores de \(c\) y \(d\) que se muestran en el recuadro rojo de la figura 8.1). Es importante resaltar que se ha tomado como referencia a prog=academic:

## # weights:  6 (2 variable)
## initial  value 219.722458 
## final  value 204.096674 
## converged
library(nnet)
prog <- relevel(as.factor(prog), ref = "academic")

modelo <- multinom(prog ~ 1, data=Datos)
summary(modelo)
Modelo nulo. Fuente: Elaboración propia.

Figure 8.1: Modelo nulo. Fuente: Elaboración propia.

Con ayuda de los valores de \(c\) y \(d\), obtenemos los mismos resultados obtenidos en (a). En el documento Rpbus :: Regresión logística (estimaciones) se explican las fórmulas de abajo:

\[\begin{eqnarray*} \widehat{p}_0 &=& \frac{\exp\{c\}}{1 + \exp\{c\}+ \exp\{d\}} \;=\; \frac{\exp\{-0.847298\}}{1 + \exp\{-0.847298\}+ \exp\{-0.7419374\}} \;=\; 0.225\\ &&\\ \widehat{p}_1 &=& \frac{\exp\{d\}}{1 + \exp\{c\}+ \exp\{d\}\}} \;=\; \frac{\exp\{-0.7419374\}}{1 + \exp\{-0.847298\}+ \exp\{-0.7419374\}} \;=\; 0.25 \\ \end{eqnarray*}\]

En R:

c <- summary(modelo)$coefficients[1]
d <- summary(modelo)$coefficients[2]
p0 <- exp(c)/(1+ exp(c)+ exp(d))
p0
p1 <- exp(d)/(1+ exp(c)+ exp(d))
p1

8.0.4 Solución parte (c)

La función summary::vglm() del paquete VGAM entrega unas estimaciones que nos pueden ayudar a calcular lo solicitado, en especial, los valores de \(c\) y \(d\) que se muestran en el recuadro rojo de la figura 8.2). Es importante recalcar que se ha tomado como referencia a prog=academic (que es el nivel 1 en el datasets, por eso, refLevel=1 dentro de la familia multinomial()). Observe que son los mismos valores de \(c\) y \(d\) aplicados en la parte (b):

library(VGAM)

modelo <- vglm(prog ~ 1, multinomial(refLevel = 1), data=Datos)
summary(modelo)
Modelo nulo. Fuente: Elaboración propia.

Figure 8.2: Modelo nulo. Fuente: Elaboración propia.

8.0.5 Solución parte (d)

En el modelo nulo, la estimación del logaritmo de la función de verosimilitud es

\[{\cal L}(\widehat{p})\;=\; {\cal L}(0.225, 0.25) \;= \; -204.0967\]

y el valor de la desviación residual es:

\[-2{\cal L}(\widehat{p}) \;=\; 408.1933\]

En R, se puede verificar así:

n <- nrow(Datos)
n
u0_bar <- mean(prog0)
u0_bar
p0 <- u0_bar
u1_bar <- mean(prog1)
u1_bar
p1 <- u1_bar
u2_bar <- mean(prog2)
u2_bar
p2 <- u2_bar
LogNulo <- n*(u0_bar*log(p0)+u1_bar*log(p1)+(1-u0_bar-u1_bar)*log(1-p0-p1))
LogNulo
DevNulo <- -2*LogNulo

8.0.6 Solución parte (e)

La función multinom() del paquete nnet calcula directamente el valor de la desviación residual \(-2{\cal L}(\widehat{p})\) y, con ello, el valor de \({\cal L}(\widehat{p})\). En la salida de summary(), que se muestra en la figura 8.3, solo debe tenerse en cuenta los resultados que se indican en el recuadro rojo. Es importante recalcar que se ha tomado como referencia a prog=academic:

library(nnet)
prog <- relevel(as.factor(prog), ref = "academic")

modelo <- multinom(prog ~ 1, data=Datos)
summary(modelo)
Modelo nulo. Fuente: Elaboración propia.

Figure 8.3: Modelo nulo. Fuente: Elaboración propia.

8.0.7 Solución parte (f)

La función vglm() del paquete VGAM calcula directamente los valores de la desviación residual \(-2{\cal L}(\widehat{p})\) y de \({\cal L}(\widehat{p})\). En la salida de summary(), que se muestra en la figura 8.4, solo debe tenerse en cuenta los resultados que se indican en el recuadro rojo. Es importante recalcar que se ha tomado como referencia a prog=academic (que es el nivel 1 en el datasets, por eso, refLevel=1 dentro de la familia multinomial()):

library(VGAM)

modelo <- vglm(prog ~ 1, multinomial(refLevel = 1), data=Datos)
summary(modelo)
Modelo nulo. Fuente: Elaboración propia.

Figure 8.4: Modelo nulo. Fuente: Elaboración propia.

9 Modelo saturado

El modelo saturado está caracterizado por dos supuestos.

Hypothesis 9.1 (Supuesto 1 en el modelo saturado) Se supone que:

  1. Se tienen \(K\) variables explicativas \(X_1, \cdots, X_K\) (algunas pueden ser numéricas y otras categóricas) con valores \(x_{i1}, \cdots, x_{iK}\) para \(i=1, \cdots, n\) (fijadas u observadas por el estadístico, según sean variables determiní}sticas o aleatorias).

  2. Entre las \(n\) kuplas \((x_{i1}, \cdots, x_{iK})\) de los valores de la variable explicativa \(X\) haya \(J\) kuplas diferentes, definiendo las \(J\) poblaciones. Por tanto, \(J \le n\).

Remark. Para cada población \(j=1, \cdots ,J\) se denota:

  • El número de observaciones \(Y_{ij}\) (o de observaciones \(U_{rij}\) en la categoría \(r\)) en cada población \(j\) por \(n_j\), siendo \(n_1+\cdots +n_J=n\);

  • Para cada \(r=0,1,2\) fijo, la suma de las \(n_j\) observaciones \(U_{rij}\) en \(j\) por

\[Z_{rj}:=\sum\limits_{i=1}^{n_j}U_{rij} \quad \mbox{con valor}\quad z_{rj}=\sum\limits_{i=1}^{n_j} u_{rij},\quad \mbox{siendo}\quad \sum\limits^J_{j=1}z_{rj}= \sum\limits^n_{i=1} u_{ri}\]

En la Tabla 9.1 se ilustra hipotéticamente un conjunto de datos con \(J=2\) poblaciones.

Table 9.1: Ilustración de un conjunto de datos agrupado en \(J=2\) poblaciones
1 2 3 4 5 6 7 8 9 10 11 12 13 14
\(Y\) \(X_1\) \(X_2\) \(X_3\) \(X_4\) \(X_5\) \(U_0\) \(U_1\) \(U_2\) \(j\) \(n_j\) \(Z_1j\) \(Z_2j\) \(Z_3j\)
Población: Bajo, 80, Si, 170, Estrato 1
1 Bajo 80 Si 170 Estrato 1 0 1 0 \(j=1\) \(n_1=7\) \(Z_1=2\) \(Z_2=2\) \(Z_1=3\)
0 Bajo 80 Si 170 Estrato 1 1 0 0
2 Bajo 80 Si 170 Estrato 1 0 0 1
0 Bajo 80 Si 170 Estrato 1 1 0 0
2 Bajo 80 Si 170 Estrato 1 0 0 1
1 Bajo 80 Si 170 Estrato 1 0 1 0
2 Bajo 80 Si 170 Estrato 1 0 0 1
Población: Alto, 100, No, 180, Estrato 2
1 Alto 100 No 180 Estrato 2 0 1 0 \(j=2\) \(n_2=5\) \(Z_2=1\) \(Z_2=3\) \(Z_2=1\)
2 Alto 100 No 180 Estrato 2 0 0 1
0 Alto 100 No 180 Estrato 2 1 0 0
1 Alto 100 No 180 Estrato 2 0 1 0
1 Alto 100 No 180 Estrato 2 0 1 0
General. \(Y\) es la variable de respuesta; \(X_1, \cdots, X_5\) son las variables explicativas; \(U_r\) son las variables de respuestas dicotomizadas; \(j\) es la población; \(n_j\) es el tamaño de la población \(j\); \(Z_{rj}\) es el número de éxitos en la población \(j\), ubicado en el nivel \(r\).

Hypothesis 9.2 (Supuesto 2 en el modelo saturado) Para mayor simplicidad en la escritura, se abreviará la j-ésima población \((x_{j1}, \cdots ,x_{jK})\) por el símbolo \(\star\). Para cada \(r=0,1,2\) fijo, cada población \(j=1, \cdots ,J\) y cada observación \(i=1,\cdots,n\) en \(j\), se supone que:

  1. \((U_{rij}|\star)\) es de Bernoulli. Es decir,

\[(U_{rij}|\star) \sim {\cal B}(1,p_{rj})\]

  1. Las variables \((U_{rij}|\star)\) son independientes entre sí.

  2. La esperanza y la varianza son, respectivamente,

\[p_{rj}=P(U_{rij}=1|\star)=E(U_{rij}|\star), \qquad V(U_{rij}|\star)=p_{rj}(1-p_{rj})\]

A continuación, se oprimirá el símbolo \(\star\).

Remark. El supuesto 2 implica:

  1. Para cada \(r=0,1,2\) y cada poblaci'on \(j=1, \cdots ,J\), todos los \(p_{rij}\), \(i=1, \cdots ,n\) dentro de cada población \(j\) son iguales. Es decir, se tiene como parámetro el vector \(2J\)-dimensional:

\[p=(p_{01}, p_{11}, \ldots ,p_{0J},p_{1J})^T\]

  1. Para cada \(r=0,1,2\) y cada población \(j=1, \cdots ,J\):

    • La variable \(Z_{rj}\) es binomial. Es decir,

    \[Z_{rj}\sim{\cal B}(n_j,p_{rj})\]

    • Las variables \(Z_{rj}\) son independientes entre las poblaciones.

Theorem 9.1 (Log-verosimilitud en el modelo saturado) En el modelo saturado, el logaritmo de la función de máxima verosimilitud será

\[\begin{eqnarray} {\cal L}(p) &= & \sum^J_{j=1}\left[z_{0j}\ln p_{0j} \;+\; z_{1j}\ln p_{1j} \;+\; (n_j- z_{0j}-z_{1j})\ln (1-p_{0j}-p_{1j})\right] \tag{9.1} \end{eqnarray}\]

Theorem 9.2 (Estimaciones en el modelo saturado) En el modelo saturado, las ML-estimaciones de \(p_{rj}\) son \(\tilde{p}_{rj}=\frac{Z_{rj}}{n_j}\), con valores \(\tilde{p}_{rj}=\frac{z_{rj}}{n_j}\),\(j=1,\cdots ,J\). Además,

\[\begin{eqnarray} {\cal L}(\widetilde{p}) &=& \sum^J_{j=1} n_j[\tilde{p}_{0j}\ln \tilde{p}_{0j} \;+\; \tilde{p}_{1j}\ln \tilde{p}_{1j} \; + \;(1-\tilde{p}_{0j}-\tilde{p}_{1j})\ln(1-\tilde{p}_{0j}-\tilde{p}_{1j})] \end{eqnarray}\]

También se cumple que

\[{\cal L}(\widetilde{p})<0\quad \mbox{para}\quad 0< \tilde{p}_j <1\]

Example 9.1 Para los datos del archivo chdage, en el modelo saturado,hay \(J=43\) poblaciones y se cumple que \({\cal L}(\tilde{p})=-41.7991\), como se indica en la última fila de la Tabla 9.2:

Datos %>%
  group_by(gender,ses) %>%
  summarise(nj = n(),
            z0j = sum(prog0),
            z1j = sum(prog1),
            z2j = sum(prog2)) %>%
  mutate(p0j = round(z0j/nj,4),
         p1j = round(z1j/nj,4),
         p2j = round(z2j/nj,4),
         Lp_ref2 = ifelse(z0j==0 | z0j== nj| z1j==0 | z1j==nj |z0j+z1j==0 | z0j+z1j==nj, 0,
                     z0j*log(p0j)+z1j*log(p1j)+(nj-z0j-z1j)*log(1-p0j-p1j)),
        Lp_ref2 = round(Lp_ref2, 4)) -> saturado
        
L_saturado <- sum(saturado$Lp_ref2)
Table 9.2: Estimación en el modelo saturado: \({\cal L}(\tilde{p})= \sum\limits_{j=1}^J {\cal L}_j(\tilde{p}) =-41.7991\)
gender ses \(n_j\) \(z_{0j}\) \(z_{1j}\) \(z_{2j}\) \(p_{0j}\) \(p_{1j}\) \(p_{2j}\) \({\cal L}_j(\tilde{p})\)
female high 29 5 3 21 0.1724 0.1034 0.7241 -22.3736
female low 32 9 8 15 0.2812 0.2500 0.4688 -33.8722
female middle 48 10 16 22 0.2083 0.3333 0.4583 -50.4274
male high 29 4 4 21 0.1379 0.1379 0.7241 -22.6263
: : :
female middle 48 10 16 22 0.2083 0.3333 0.4583 -50.4274
male high 29 4 4 21 0.1379 0.1379 0.7241 -22.6263
male low 15 7 4 4 0.4667 0.2667 0.2667 -15.9090
male middle 47 10 15 22 0.2128 0.3191 0.4681 -49.3074

10 Modelo logístico

Hypothesis 10.1 (Supuesto 3: matriz de diseño) Se hacen los supuestos 1 y 2 del modelo saturado (véase las hipótesis 9.1 y 9.2), donde adicionalmente se supone que la matriz de diseño

\[C=\left(\begin{array}{cccc} 1 & x_{11} &\cdots &x_{1K}\\ 1 & x_{21} &\cdots &x_{2K}\\ \vdots &\vdots & &\vdots\\ 1 &x_{J1} &\cdots &x_{JK}\\ \end{array}\right)\]

tiene rango completo \(Rg(C)=1+K\leq J\).

Hypothesis 10.2 (Supuesto 4: modelo logístico) Para llegar a un modelo logístico se toma como referencia una de las categorías de la variable dependiente \(Y\), digamos 2, y se hace el supuesto adicional:

\[\begin{eqnarray} \mbox{Logit}(p_{0j}) &:=& \ln\left(\frac{p_{0j}}{p_{2j}}\right)= \delta_0 + \beta_{01} x_{j1} + \cdots + \beta_{0K} x_{jK} \\ \mbox{Logit}(p_{1j}) &:=& \ln\left(\frac{p_{1j}}{p_{2j}}\right)= \delta_1 + \beta_{11} x_{j1} + \cdots + \beta_{1K} x_{jK} \tag{10.1} \end{eqnarray}\]

Remark. Tenemos:

  1. El vector de los \(2(1+K)\) parámetros en el modelo es:

\[\alpha = (\beta_0,\beta_1)^T= (\delta_0,\beta_{01},\cdots,\beta_{0K}, \delta_1,\beta_{11},\cdots,\beta_{1K})^T\]

  1. Nótese que el supuesto sobre \(Rg(C)=1+K\), hace identificable al parámetro \(\alpha\).

11 Riesgo

Definition 11.1 (Riesgo) En la práctica, la probabilidad \(p_{rj}\) es conocida como riesgo.

Theorem 11.1 (Fórmula para el riesgo) Sea \(g_{rj}:=\delta_r \;+\; \beta_{r1}\,x_{j1} \;+\;\cdots \;+\; \beta_{rK} \,x_{jK}\). Entonces, la probabilidad

\[p_{rj}=P(Y_j=r|x_{j1}, \cdots, x_{jK})\]

de obtener el valor \(r\) en la población \(j=1, \ldots, J\), dado los valores \(x_{j1}, \cdots, x_{jK}\), viene dada por:

\[\begin{equation} p_{rj} \;= \; \frac{e^{g_{rj}}}{1 + e^{g_{0j}} + e^{g_{1j}}} \tag{11.1} \end{equation}\]

Theorem 11.2 (Log de la función de verosimilitud en el modelo logístico) Sea \(g_{rj}:=\delta_r \;+\; \beta_{r1}\,x_{j1} \;+\;\cdots \;+\; \beta_{rK} \,x_{jK}\). Entonces, el logaritmo de la función de verosimilitud se puede escribir, en función de \(\alpha\), como:

\[\begin{eqnarray} {\cal L}({\alpha}) &=& \sum\limits_{j=1}^J\left[z_{0j}g_{0j} + (n_j-z_{0j}-z_{1j})g_{1j} - n_j\ln\left(1 + e^{g_{0j}} + e^{g_{1j}} \right)\right] \tag{11.2} \end{eqnarray}\]

12 Método de estimación

El método que se propone para calcular las ML-estimaciones en un modelo logístico es el método iterativo de Newton-Raphson. Generalmente, el método requiere:

  1. Una estimación inicial para el valor que maximiza la función.

  2. La función es aproximada en una vecindad de aquella estimación por un polinomio de segundo grado.

  3. Entonces,la siguiente estimación se calcula como el máximo de dicho polinomio.

  4. Luego, se repite el proceso, usando esta estimación como la estimación inicial.

  5. De esta manera, el método genera una sucesión de estimaciones. Estas estimaciones convergen a la localización del máximo cuando la función es adecuada y/o la estimación inicial es buena.

Para más detalles, ver el teorema 5.1 en LLinás, Arteta y Tilano (2016).

13 Casos agrupado y no agrupado

  1. Cuando se trabaja con el modelo saturado, se tiene el caso de utilizar datos agrupados.

  2. Cuando se tiene el caso especial \(n_j=1\), para todo \(j\) (lo que implica que \(J=n\)) se habla de datos no agrupados.

  3. La distinción entre datos agrupados y no agrupados es importante por dos razones:

    1. Algunos métodos de análisis apropiados a datos agrupados no son aplicables a datos no agrupados.

    2. Las aproximaciones asintóticas pueden estar basados en uno de estos dos casos distintos: (i) \(n\to\infty\) o (ii) \(J\to\infty\), caso que es únicamente es apropiado para datos no agrupados.

  4. En la práctica:

    1. Cuando se tienen datos agrupados es importante tener en cuenta que \(J\) debe ser fijo. Por esta razón, debe tomarse como base el modelo saturado. Es decir, se empieza el análisis usando los vectores \(Z_{rj}\), \(j=1,\cdots,J\), \(r=0, 1,\ldots, R-1\).

    2. Si \(J\to\infty\) (por ejemplo, si \(J=n\)), entonces, en el modelo saturado no se puede considerar a \(J\) como fijo. Obsérvese que esta situación se presenta cuando se tienen datos no agrupados. En este caso, no se puede tomar como base el modelo saturado. Ahora se empezaría el análisis utilizando, de una vez, las observaciones \(Y_i\), \(i=1,\cdots, n\).

14 Ejemplo 1: Enunciado

Considere los datos del archivo hsbdemo. Suponga que se quiere analizar un modelo de regresión logística, considerando a prog como variable dependiente y ses y write como independientes. Tome como referencia al valor prog=academic.

  1. Escriba, matemáticamente, el vector de parámetros logísticos y el de sus estimadores.

  2. Escriba, matemáticamente, La probabilidad estimada de que un individuo tenga enfermedades coronarias (chd\(=1\)), cuando tiene una edad determinada (digamos, age\(=x_j\)).

  3. Escriba, matemáticamente, el modelo logístico estimado.

  4. Obtenga las estimaciones \(\hat{\delta}\) y \(\hat{\beta}\) de los parámetros logísticos \(\delta\) y \(\beta\) utilizando la función summary :: multinom :: nnet.

  5. Obtenga las estimaciones \(\hat{\delta}\) y \(\hat{\beta}\) de los parámetros logísticos \(\delta\) y \(\beta\) utilizando la función summary :: vglm :: VGAM.

  6. Utilizando las estimaciones halladas en los incisos (d) o (e), escriba en el modelo correspondiente .

  7. Haga la gráfica del riesgo versus write, cuando ses=low. ¿Es directa o indirecta esta relación?

  8. Halle las razones odds correspondientes.

  9. Calcular las probabilidades predichas para cada uno de los niveles de la variable de respuesta.

  10. Mantener fija la variable write en su media y examinar las probabilidades predichas para cada nivel de ses.

  11. Examinar los cambios en la probabilidad predicha asociados con cada nivel de ses, manteniendo write fija en su media.

  12. Examinar las probabilidades predichas promediadas para diferentes valores de la variable predictora continua write dentro de cada nivel de ses.

  13. Usar las predicciones que se generó en el inciso anterior y grafique las probabilidades predichas contra la puntuación de write, para cada nivel de ses y para diferentes niveles de la variable de respuesta.

  14. Halle los errores estándares estimados de los estimadores de los parámetros logísticos, utilizando la función multinom :: nnet.

  15. Halle los errores estándares estimados de los estimadores de los parámetros logísticos, utilizando la función vglm :: VGAM.

  16. Calcule \({\cal L}(\hat{\alpha})\), la estimación del logaritmo de la función de máxima verosimilitud en el modelo logístico.

15 Ejemplo 1: Solución

15.0.1 Solución parte (a)

El vector de parámetros y el de sus estimadores son, respectivamente,

\[\alpha =(\delta_0,\beta_{01}, \beta_{02}, \delta_1,\beta_{11}, \beta_{12})^T, \qquad \hat{\alpha} =\left(\hat{\delta}_0, \hat{\beta}_{01}, \hat{\beta}_{02}, \hat{\delta}_1, \hat{\beta}_{11}, \hat{\beta}_{12}\right)^T\]

Aquí, \(T\) indica la transpuesta del vector.

15.0.2 Solución parte (b)

La probabilidad estimada de que un estudiante elija el programa \(r=0,1,2\), cuando tiene un estrato socioeconómico determinado (digamos, age\(=x_j\)), se puede escribir así:

\[\hat{p}_{rj} = \hat{P}(\mbox{prog}=r \,|\, \mbox{ses}=s_j, \, \mbox{write}=w_j )\]

15.0.3 Solución parte (c)

Sea

\[\hat{g}_{rj}:=\hat{\delta}_r + \hat{\beta}_{r1} s_j+ \hat{\beta}_{r2} w_j\]

donde \(s_j\) es un posible valor de la variable ses y \(w_j\) es un posible valor de la variable write. Sabiendo que \(\hat{p}_{rj}\) es como en el inciso anterior, el modelo estimado se puede escribir teniendo en cuenta la ecuación (10.1) o la (11.1):

\[\mbox{Logit}(\hat{p}_{rj}):= \ln\left(\frac{\hat{p}_{rj}}{\hat{p}_{2j}}\right) = \hat{g}_{rj}, \qquad \qquad \hat{p}_{rj} \;= \; \frac{e^{\hat{g}_{rj}}} {1 + e^{\hat{g}_{0j}}+ e^{\hat{g}_{1j}}} \]

15.0.4 Solución parte (d)

En R, las estimaciones de los parámetros logísticos \(\delta_r\) y \(\beta_r\) se pueden obtener con la función multinom() del paquete nnet. En la salida de summary(), que se muestra en la figura 15.1, solo debe tenerse en cuenta los resultados que se indican en el recuadro rojo. Primero, debemos elegir el nivel de referencia de nuestra variable de respuesta especificándolo con la función relevel. Es importante recalcar que se ha tomado como referencia a prog=academic:

library(nnet)
prog <- relevel(as.factor(prog), ref = "academic")

modelo <- multinom(prog ~ ses + write, data=Datos)
summary(modelo)
Estimaciones de los parámetros logísticos con nnet. Fuente: Elaboración propia.

Figure 15.1: Estimaciones de los parámetros logísticos con nnet. Fuente: Elaboración propia.

15.0.5 Solución parte (e)

En R, las estimaciones de los parámetros logísticos \(\delta_r\) y \(\beta_r\), también, se pueden obtener con la función vglm() del paquete VGAM. En la salida de summary(), que se muestra en la figura 15.2, solo debe tenerse en cuenta los resultados que se indican en el recuadro rojo. Como ya se mencionó anteriormenbte, primero, debemos elegir el nivel de referencia de nuestra variable de respuesta especificándolo con la función relevel. Es importante recalcar que se ha tomado como referencia a prog=academic (que es el nivel 1 en el datasets, por eso, refLevel=1 dentro de la familia multinomial()):

library(VGAM)

modelo <- vglm(prog ~ ses + write, multinomial(refLevel = 1), data=Datos)
summary(modelo)
Estimaciones de los parámetrdos logísticos con nnet. Fuente: Elaboración propia.

Figure 15.2: Estimaciones de los parámetrdos logísticos con nnet. Fuente: Elaboración propia.

15.0.6 Solución parte (f)

Sabemos que:

\[\hat{p}_{rj} = \hat{P}(\mbox{prog}=r \,|\, \mbox{ses}=s_j, \, \mbox{write}=w_j )\]

y

\[\hat{g}_{rj}\;:=\; \hat{\delta}_r + \hat{\beta}_{r1} s_j+ \hat{\beta}_{r2} w_j \]

donde \(s_j\) es un posible valor de la variable ses y \(w_j\) es un posible valor de la variable write. Además, el modelo estimado se puede escribir teniendo en cuenta la ecuación (10.1) o la (11.1):

\[\mbox{Logit}(\hat{p}_{rj}):= \ln\left(\frac{\hat{p}_{rj}}{\hat{p}_{2j}}\right) = \hat{g}_{rj}, \qquad \qquad \hat{p}_{rj} \;= \; \frac{e^{\hat{g}_{rj}}} {1 + e^{\hat{g}_{0j}} + e^{\hat{g}_{1j}}} \]

El modelo estimado se puede escribir utilizando una de las dos expresiones anteriores. Por ejemplo, para ses=low:

\[\begin{eqnarray*} \mbox{general} &:& \hat{g}_{0j} \;=\; \hat{\delta}_0 + \hat{\beta}_{01} s_j+ \hat{\beta}_{02} w_j \;=\; 1.6895 + 1.1628 - 0.0579 \,w_j \;=\; 2.8523 - 0.0579 \,w_j\\ \mbox{vocation} &:& \hat{g}_{1j}\;:=\; \hat{\delta}_1 + \hat{\beta}_{11} s_j+ \hat{\beta}_{12} \,w_j \;=\; 4.2356 + 0.9827 - 0.1136 \,w_j \;=\; 5.2183 - 0.1136 \,w_j \end{eqnarray*}\]

Por lo tanto, las ecuaciones correspondientes son:

\[\begin{eqnarray*} \mbox{general} &:& \hat{p}_{0j} \;= \; \frac{e^{\hat{g}_{0j}}} {1 + e^{\hat{g}_{0j}}+ e^{\hat{g}_{1j}}}\\ \mbox{vocation} &:& \hat{p}_{1j} \;= \; \frac{e^{\hat{g}_{1j}}} {1 + e^{\hat{g}_{0j}} + e^{\hat{g}_{1j}}} \end{eqnarray*}\]

15.0.7 Solución parte (g)

La gráfica correspondiente se muestra en la figura 15.3. Cuando ses=low, se observa que hay una relación indirecta entre write y prog, el riesgo de que un estudiante seleccione un programa:

library(nnet)
prog <- relevel(as.factor(prog), ref = "academic")

modelo <- multinom(prog ~ ses + write, data=Datos)
s<- summary(modelo)

W <- seq(0, 100, 0.05)

g_0 <- s$coefficients[1,1] + s$coefficients[1,2] + s$coefficients[1,4]*W
g_1 <- s$coefficients[2,1] + s$coefficients[2,2] + s$coefficients[2,4]*W

p_0 <- exp(g_0)/(1 + exp(g_0) + exp(g_1))
p_1 <- exp(g_1)/(1 + exp(g_0) + exp(g_1))
p_2 <- 1- (p_0 + p_1)

ggplot() +
geom_point(mapping=aes(y=p_0, x = W, color="p_0"),size=1.5 ) +
geom_point(mapping=aes(y=p_1, x = W, color="p_1"),size=1.5) +  
geom_point(mapping=aes(y=p_2, x = W, color="p_2"),size=1.5) +

labs(x="Write", y="Probabilidad de éxito", fill= "") + 

facet_wrap(. ~ "Gráfica dentro del grupo de ses=low") +    
theme_bw(base_size = 12) +
scale_color_discrete(name = "Tipo de programa", labels = c("(a) General", "(b) Vocation", "(c) Academic"))
Probabilidades de éxito cuando se mantiene constante a la variable ses

Figure 15.3: Probabilidades de éxito cuando se mantiene constante a la variable ses

15.0.8 Solución parte (h)

En R, las estimaciones de los odds se pueden obtener con las funciones multinom :: nnet o vglm :: VGAM. Simplemente debe utilizar la función exp() a los coeficientes del modelo, como se resalta en el recuadro rojo de la figura 15.4. Recuerde que, primero, debemos elegir el nivel de referencia de nuestra variable de respuesta (en nuestro caso, prog=academic y es el nivel 1 en el datasets). En nnet se especifica con la función relevel; y en VGAM, se coloca refLevel=1 dentro de la familia multinomial():

library(nnet)
prog <- relevel(as.factor(prog), ref = "academic")
modelo <- multinom(prog ~ ses + write, data=Datos)

library(VGAM)
modelo <- vglm(prog ~ ses + write, multinomial(refLevel = 1), data=Datos)

exp(coef(modelo))
Estimaciones de los odds. Fuente: Elaboración propia.

Figure 15.4: Estimaciones de los odds. Fuente: Elaboración propia.

Algunas interpretaciones (en el modelo del programa general relativo a académico):

  1. ses=low. Esta es la razón del riesgo relativo comparando ses low con ses high para la preferencia del programa general al académico, Para low relativo a high, el riesgo relativo de preferencia de general, en vez de académico, se esperaría que aumentara en un factor de 3.199, dado que las otras variables del modelo se mantienen fijas. En otras palabras, los estudiantes con ses=low son más probables que los hombres en seleccionar programa general sobre el académico.

  2. write. Esta es la razón del riesgo relativo para una unidad de incremento en los puntajes de write para la preferencia del programa general al académico, dado que las otras variables del modelo se mantienen fijas. Si un estudiante fuera a incrementar su puntaje en escritura en una unidad, el riesgo relativo para preferir el programa general (en vez del académico), se esperaría que decreciera en un factor de 0.9437, dado que las otras variables del modelo se mantienen fijas..

15.0.9 Solución parte (i)

Se pueden calcular las probabilidades predichas para cada uno de los niveles de la variable de respuesta, utilizando la función fitted():

head(fitted(modelo))
Table 15.1: Probabilidades predichas
1 2 3 4 5 6
Student ses write \(P_2\) (academic) \(P_0\) (general) \(P_1\) (vocation)
1 low 35 0.1483 0.3383 0.5135
2 middle 33 0.1202 0.1806 0.6992
3 high 39 0.4187 0.2368 0.3445
4 low 37 0.1727 0.3508 0.4765
5 middle 31 0.1001 0.1689 0.7309
6 high 36 0.3534 0.2378 0.4088

Observe que, por ejemplo, los valores de la primera fila se obtienen reemplazando ses=low y write=35, en las las fórmulas para las calcular las probabilidades estimadas, como se indica en la parte (c):

W <- 35
g_0 <- s$coefficients[1,1] + s$coefficients[1,2] + s$coefficients[1,4]*W
g_1 <- s$coefficients[2,1] + s$coefficients[2,2] + s$coefficients[2,4]*W
p_0 <- exp(g_0)/(1 + exp(g_0) + exp(g_1))
p_1 <- exp(g_1)/(1 + exp(g_0) + exp(g_1))
p_2 <- 1- (p_0 + p_1)
predichos <- c(p_0,p_1,p_2)
predichos

15.0.10 Solución parte (j)

Para examinar los cambios en la probabilidad predicha asociados con una de nuestras dos variables (en este caso, ses), se pueden crear pequeños conjuntos de datos que varíen la variable mientras se mantiene la otra constante (en este caso, sería write). Primero haremos esto manteniendo la variable write en su media y examinar las probabilidades predichas para cada nivel de ses.

Nuevo <- data.frame(ses = c("low", "middle", "high"), write = mean(write))
Tabla <- predict(modelo, newdata = Nuevo, "probs")
Tabla
##    academic   general  vocation
## 1 0.4396813 0.3581915 0.2021272
## 2 0.4777451 0.2283359 0.2939190
## 3 0.7009046 0.1784928 0.1206026

15.0.11 Solución parte (k)

Ahora sí se examinarán los cambios en la probabilidad predicha asociados con una de nuestras dos variables (en este caso, ses), manteniendo write fija en su media.

W <- mean(write)
g_0 <- s$coefficients[1,1] + s$coefficients[1,2] + s$coefficients[1,4]*W
g_1 <- s$coefficients[2,1] + s$coefficients[2,2] + s$coefficients[2,4]*W
p_0 <- exp(g_0)/(1 + exp(g_0) + exp(g_1))
p_1 <- exp(g_1)/(1 + exp(g_0) + exp(g_1))
p_2 <- 1- (p_0 + p_1)
predichos <- c(p_0,p_1,p_2)
predichos

15.0.12 Solución parte (l)

las probabilidades predichas promediadas para diferentes valores de la variable predictora continua write dentro de cada nivel de ses se pueden hallar siguiendo los pasos siguientes (que se ilustran en el código de abajo):

Paso 1: Se crea un data frame con diferentes valores de ses y write.

Paso 2: Se calculan y se almacenan las probabilidades predichas para cada nivel de ses.

Paso 3: Se calculan las probabilidades promedios dentro de cada nivel de ses.

# Paso 1:
  Nuevo <- data.frame(ses = rep(c("low", "middle", "high"), each = 51), write = rep(c(20:70),3))
  head(Nuevo)

# Paso 2: 
  Predichos <- cbind(Nuevo, predict(modelo, newdata = Nuevo, type = "probs", se = TRUE))
  Predichos

# Paso 3: 
  Predichos %>% 
  dplyr::group_by(ses) %>% 
  dplyr::summarise(P0_general=mean(general),
                   P1_vocation=mean(vocation),
                   P2_academic=mean(academic))
Table 15.2: Probabilidades predichas
1 2 3 4
ses \(P_0\) (general) \(P_1\) (vocation) \(P_2\) (academic)
high 0.1842 0.2898 0.526
low 0.313 0.3554 0.3316
middle 0.1875 0.4595 0.353

15.0.13 Solución parte (m)

Usando las predicciones que se generó en el inciso anterior a través del objeto Predichos, graficaremos las probabilidades predichas contra la puntuación de write, para cada nivel de ses y para diferentes niveles de la variable de respuesta. Utilizaremos la función reshape2::melt() con el fin de convertir la tabla de Predichos en una de tipo long y, así, utilizarla para ggplot2:

library(reshape2)
Tabla <- melt(Predichos, id.vars = c("ses", "write"), value.name = "probability")

ggplot(Tabla, aes(x = write, y = probability, colour = ses)) + 
  geom_line() + 
  facet_wrap(variable ~ ., scales = "free")

15.0.14 Solución parte (n)

En R, las estimaciones de los errores estándares se pueden obtener con la función multinom() del paquete nnet. En la salida de summary(), que se muestra en la figura 15.5, solo debe tenerse en cuenta los resultados que se indican en el recuadro rojo. Primero, debemos elegir el nivel de referencia de nuestra variable de respuesta especificándolo con la función relevel. Es importante recalcar que se ha tomado como referencia a prog=academic:

library(nnet)
prog <- relevel(as.factor(prog), ref = "academic")
modelo <- multinom(prog ~ ses + write, data=Datos)

summary(modelo)
Estimaciones de los errores estándares con nnet. Fuente: Elaboración propia.

Figure 15.5: Estimaciones de los errores estándares con nnet. Fuente: Elaboración propia.

15.0.15 Solución parte (o)

En R, las estimaciones de los errores estándares se pueden obtener con la función multinom() del paquete nnet. En la salida de summary(), que se muestra en la figura 15.6, solo debe tenerse en cuenta los resultados que se indican en el recuadro rojo. Primero, debemos elegir el nivel de referencia de nuestra variable de respuesta especificándolo con la función relevel. Es importante recalcar que se ha tomado como referencia a prog=academic:

library(VGAM)
modelo <- vglm(prog ~ ses + write, multinomial(refLevel = 1), data=Datos)

summary(modelo)
Estimaciones de los errores estándares con VGAM. Fuente: Elaboración propia.

Figure 15.6: Estimaciones de los errores estándares con VGAM. Fuente: Elaboración propia.

15.0.16 Solución parte (p)

La estimación del logaritmo de la función de máxima verosimilitud en el modelo logístico se puede obtener de varias maneras. Una es reemplazando los valores correspondientes en la ecuación (11.2):

\[ {\cal L}(\hat{\alpha}) \;=\; \sum\limits_{j=1}^J\left[z_{0j}\hat{g}_{0j} + (n_j-z_{0j}-z_{1j})\hat{g}_{1j} - n_j\ln\left(1 + e^{\hat{g}_{0j}} + e^{\hat{g}_{1j}} \right)\right]\;=\; -179.9817\]

donde \[\hat{g}_{rj}:=\hat{\delta}_r \;+\; \hat{\beta}_{r1}\,x_{j1} \;+\;\cdots \;+\; \hat{\beta}_{rK} \,x_{jK}\]

Las otras maneras para calcular \({\cal L}(\hat{\alpha})\) son con las funciones multinom::nnet o vglm::VGAM. Con cualquier camino encontramos que \({\cal L}(\hat{\alpha})=-179.9817\).

Con la primera función se obtiene una salida (o al ejecutar summary) donde se obtiene, de alguna forma, ese valor. Ver recuadros rojos en la figura 15.7. Por ejemplo, en final value es el inverso aditivo y, en Residual Deviance se obtiene al dividir el valor de la salida por -2.

library(nnet)
prog <- relevel(as.factor(prog), ref = "academic")
modelo <- multinom(prog ~ ses + write, data=Datos)

summary(modelo)
Estimaciones del Log-Likelihood con nnet. Fuente: Elaboración propia.

Figure 15.7: Estimaciones del Log-Likelihood con nnet. Fuente: Elaboración propia.

Con la segunda función se obtiene una salida (o al ejecutar summary) donde se obtiene, de alguna forma, ese valor. Ver recuadros rojos en la figura 15.8. Por ejemplo, en Log-Likelihood es exactamente ese valor y, en Residual Deviance se obtiene al dividir el valor de la salida por -2.

library(VGAM)
modelo <- vglm(prog ~ ses + write, multinomial(refLevel = 1), data=Datos)

summary(modelo)
Estimaciones del Log-Likelihood con VGAM. Fuente: Elaboración propia.

Figure 15.8: Estimaciones del Log-Likelihood con VGAM. Fuente: Elaboración propia.

16 Intervalos de confianza

En el documento Rpbus :: Modelos logísticos (intervalos de confianza) se ha desarrollado esta teoría, acompañado de las aplicaciones correspondientes.

17 Comparación de modelos

En esta sección se presentan estadísticas para distintas pruebas de comparación de modelos:

  1. \(H_0\): Logístico versus \(H_1\): Completo.

  2. \(H_0\): Nulo versus \(H_1\): Logístico.

  3. \(H_0\): Logístico versus \(H_1\): Saturado.

  4. \(H_0\): Submodelo versus \(H_1\): Logístico.

  5. \(H_0\): Submodelo con una variable explicativa menos versus \(H_1\): Logístico.

En el documento Rpbus :: Modelos logísticos (pruebas de hipótesis) se ha desarrollado esta teoría, acompañada con las aplicaciones correspondientes.

18 Matriz de confusión

En el campo del Aprendizaje automático (Machine Learning), una matriz de confusión (también conocida como matriz de error), es una tabla de contingencia, con dos dimensiones (“real, observado” y “predicho”) que permite la visualización del rendimiento de un algoritmo, típicamente uno de aprendizaje supervisado. Por lo general, en aprendizaje no supervisado se denomina matriz de emparejamiento(en inglés: matching matrix). Cada fila de la matriz representa las ocurrencias en la clase observada, mientras que cada columna representa las ocurrencias en la clase predicha, o viceversa. El nombre proviene del hecho de que facilita determinar si el sistema está confundiendo dos clases (es decir, comúnmente etiquetado incorrectamente una como otra). En la Tabla 18.1 se muestra la estructura general de una matrix de confusión para culaquier clase \(r=0, 1, 2\):

Table 18.1: Matrix de confusión.
1 2 3 4 5 6
Condiciones Positiva predicha Negativa predicha Total Estadísticos Estadísticos
Positiva observada Verdadero positivo (\(a\)) Falso negativo (\(c\)) \(a+c\) Sensitividad = \(a/(a+c)\) Prevalencia \(=(a+c)/n\)
Negativa observada Falso positivo (\(b\)) Verdadero negativo (\(d\)) \(b+d\) Especificidad = \(d/(b+d)\) 1- Especificidad \(=b/(b+d)\)
Total \(a+b\) \(c+d\) Tamaño muestral (\(n=a+b+c+d)\)
Estadísticos Valor Predictivo positivo \(=a/(a+b)\) Valor predictivo negativo \(=d/(c+d)\) Exactitud = \((a+d)/n\)

Teniendo en cuenta la tabla anterior, se pueden definir los siguientes conceptos:

  1. La exactitud ( en inglés: accuracy) es la proporción de clasificaciones correctas.

\[Exactitud \; =\; \frac{a+d}{n}\]

  1. La sensitividad o tasa de verdaderos positivos (en inglés: sensivity, true positive rate, recall or hit rate) es la proporción de casos positivos observados que son correctamente clasificados:

\[ Sensitividad \;=\; \frac{a}{a+c} \]

  1. La especificidad o tasa de verdaderos negativos (en inglés: specifity, selectivity or true negative rate) es la proporción de casos negativos observados que son correctamente clasificados:

\[ Especificidad \;=\; \frac{d}{b+d} \]

  1. La tasa de falsos positivos (1- especificidad) es la proporción de casos negativos observados que son incorrectamente clasificados:

\[ 1- Especificidad \;=\; \frac{b}{b+d} \]

  1. La precisión o valor predictivo positivo VPP (en inglés: precision, positive predictive value) es la proporción de casos positivos que son correctamente clasificados:

\[ VPP \;=\; \frac{a}{a+b} \]

  1. El valor predictivo negativo VPN (en inglés: negative predictive value) es la proporción de casos negativos que son correctamente clasificados:

\[ VPN \;=\; \frac{d}{c+d} \]

  1. La prevalencia (en inglés, prevalence) es la proporción de casos positivos observados:

\[ Prevalencia \;=\; \frac{a+c}{n} \]

  1. La exactitud balanceada (en inglés: balanced accuracy) se define como:

\[\mbox{Exactitud balanceada} \;=\; \frac{\mbox{sensitividad} + \mbox{especificidad}}{2}\]

  1. La tasa de detección (en inglés: Detection rate) se define como la proporción de casos positivos observados y predichos:

\[\mbox{Tasa de detección} \;=\; \frac{a}{n}\] 8. La prevalencia de detección (en inglés: Detection prevalence) se define como la proporción de casos positivos predichos:

\[\mbox{Prevalencia de detección} \;=\; \frac{a+b}{n}\]

19 Indices de concordancias

Existen varios índices de concordancia propuestos y toman valores entre 0 (total desacuerdo) y 1 (máximo acuerdo). Un ejemplo es la exactitud (definido anteriormente), pero el más usado es el denominado índice kappa de Cohen (\(\kappa\)) y se define como:

\[\kappa \;=\; \frac{Exactitud - P_E}{1- P_E}\]

donde \(P_E\) es la proporción de acuerdos esperados. Está dado por:

\[\begin{eqnarray*} P_E &=& \frac{(\mbox{Suma de (total de fila $i$)$\cdot$(total de columna $i$)}}{\mbox{(Tamaño de la muestra)}^2} \\ && \\ &=& \frac{(total fila 1)(total columna 1) + $\cdots$ + (total última fila)(total última columna)}{n^2} \end{eqnarray*}\]

Se ha utilizado la siguiente escala de valoración del grado de acuerdo para \(\kappa\):

  1. \(\kappa \leq 0.00\): sin acuerdo.

  2. \(0.00 < \kappa \leq 0.20\): insignificante.

  3. \(0.20 <\kappa \leq 0.40\): discreto.

  4. \(0.40 < \kappa \leq 0.60\): moderado.

  5. \(0.60 < \kappa \leq 0.80\): sustancial.

  6. \(0.80 < \kappa \leq 1.00\): casi perfecto.

20 Ejemplo 2: enunciado

Considere los datos del archivo hsbdemo. Suponga que se quiere analizar un modelo de regresión logística, considerando a prog como variable dependiente y ses y write como independientes. Tome como referencia al valor prog=academic.

  1. Generar diagramas de barras o de caja y bigotes para explorar la relación entre prog y las otras dos variables.

  2. Crear el modelo logístico correspondiente.

  3. Calcule los valores predichos de la variable de respuesta.

  4. Halle la matriz de confusión.

  5. Determine la exactitud del modelo.

  6. Aplique la librería caret para hallar la matriz de confusión y los estadísticos relacionados (exactitud, sensitividad, especificidad, etc.).

  7. Aplique la librería gmodels para hallar la matriz de confusión y los estadísticos relacionados (exactitud, sensitividad, especificidad, etc.).

21 Ejemplo 2: Solución

21.0.1 Solución parte (a)

Abajo aparecen los diagramas solicitados:

p1 <- ggplot(Datos, aes(x=ses, fill=prog)) + geom_bar()  
p2 <- ggplot(Datos, aes(x=prog, y=write, fill=prog)) + geom_boxplot()  

library(gridExtra)
grid.arrange(p1, p2, ncol=2)

21.0.2 Solución parte (b)

El modelo logístico correspondiente se puede generar de la siguiente manera (en el ejemplo 1 se muestran detalles al respecto):

Datos <- hsbdemo

library(nnet)
prog <- relevel(as.factor(prog), ref = "academic")
modelo <- multinom(prog ~ ses + write, data=Datos)

21.0.3 Solución parte (c)

Pazra hallar los valores predichos de la variable de respuesta, se puede aplicar el siguiente código de R:

prog_pred <- predict(modelo, Datos, type="class")

21.0.4 Solución parte (d)

La matriz de confusión correspondiente es:

Tabla <-  table(prog, prog_pred)
Table 21.1: Matriz de confusión.
1 2 3 4
academic_pred general_pred vocation_pred Total
academic 92 4 9 105
general 27 7 11 45
vocation 23 4 23 50
Total 142 15 43 200

21.0.5 Solución parte (e)

La exactitud es del 61% y se puede calcular de dos maneras:

Tabla <-  table(prog, prog_pred)

# Primera forma:
E <- mean(prog_pred == prog)*100

# Con la matriz de confusión de la parte (d):
E <- round((sum(diag(Tabla))/sum(Tabla))*100,2)

21.0.6 Solución parte (f)

Con la función confusionMatrix de la librería caret se pueden hallar la matriz de confusión y los estadísticos solicitados (véase la figura 21.1):

library(caret)
confusionMatrix(data=prog_pred, reference = prog)
Matriz de confusión con el paquete caret. Fuente: Elaboración propia.

Figure 21.1: Matriz de confusión con el paquete caret. Fuente: Elaboración propia.

21.0.7 Solución parte (g)

Con la función CrossTable de la librería gmodels se pueden hallar la matriz de confusión y también los estadísticos solicitados (véase la figura 21.2):

library(gmodels)
CrossTable(prog_pred, prog)
Matriz de confusión con el paquete gmodels. Fuente: Elaboración propia.

Figure 21.2: Matriz de confusión con el paquete gmodels. Fuente: Elaboración propia.

22 Ejercicios

Para la solución de los siguientes ejercicios, téngase en cuenta los siguientes comentarios:

  • Todos los datos mencionados aparecen en los links mencionados en este documento.

  • Siempre debe detallar el análisis del conjunto de datos (con las variables especificadas) basado en lo explicado en este documento.

  • Donde sea posible, escoja siempre el mejor submodelo.

22.0.1 Ejercicios 1 a 3

  1. Demuestre los teoremas: (a) 6.1; (b) 7.2; (c) 9.1; (d) 9.2; (e) 11.1

  2. Haga un listado de los paquetes de R que estimen el logaritmo de la función de máxima verosimilitud en los modelos completo y nulo.

  3. Haga un listado de los paquetes de R que, en el caso multinomial, estimen el logaritmo de la función de máxima verosimilitud en los modelos saturado y logístico.

22.0.2 Ejercicio 4 a 7

  1. Considere los datos hbsdemo, tomando a prog como variable dependiente y gender como variable independiente. Repita todos los análisis realizados en este documento.

  2. Considere los datos hbsdemo, tomando a prog como variable dependiente y schtyp como variable independiente. Repita todos los análisis realizados en este documento.

  3. Considere los datos hbsdemo, tomando a prog como variable dependiente y honors como variable independiente. Repita todos los análisis realizados en este documento.

  4. Considere los datos hbsdemo, tomando a prog como variable dependiente y awards como variable independiente. Repita todos los análisis realizados en este documento.

22.0.3 Ejercicios 8 a 10

  1. Considere los datos hbsdemo, tomando a prog como variable dependiente y ses y schtyp como variables independientes. Repita todos los análisis realizados en este documento.

  2. Considere los datos hbsdemo, tomando a prog como variable dependiente y ses y honors como variables independientes. Repita todos los análisis realizados en este documento.

  3. Considere los datos hbsdemo, tomando a prog como variable dependiente y ses y awards como variables independientes. Repita todos los análisis realizados en este documento.

22.0.4 Ejercicio 11 a 13

  1. Considere los datos hbsdemo, tomando a prog como variable dependiente y gender y schtyp como variables independientes. Repita todos los análisis realizados en este documento.

  2. Considere los datos hbsdemo, tomando a prog como variable dependiente y gender y honors como variables independientes. Repita todos los análisis realizados en este documento.

  3. Considere los datos hbsdemo, tomando a prog como variable dependiente y gender y awards como variables independientes. Repita todos los análisis realizados en este documento.

22.0.5 Ejercicio 14 a 16

  1. Considere los datos hbsdemo, tomando a prog como variable dependiente y schtyp y honors como variables independientes. Repita todos los análisis realizados en este documento.

  2. Considere los datos hbsdemo, tomando a prog como variable dependiente y schtyp y awards como variables independientes. Repita todos los análisis realizados en este documento.

  3. Considere los datos hbsdemo, tomando a prog como variable dependiente y honors y awards como variables independientes. Repita todos los análisis realizados en este documento.

22.0.6 Ejercicios 17 a 19

  1. Considere los datos hbsdemo, tomando a prog como variable dependiente y gender, ses y schtyp como variables independientes. Repita todos los análisis realizados en este documento.

  2. Considere los datos hbsdemo, tomando a prog como variable dependiente y gender, ses y honors como variables independientes. Repita todos los análisis realizados en este documento.

  3. Considere los datos hbsdemo, tomando a prog como variable dependiente y gender, ses y awards como variables independientes. Repita todos los análisis realizados en este documento.

22.0.7 Ejercicios 20 a 21

  1. Considere los datos hbsdemo, tomando a prog como variable dependiente y gender, schtyp y honors como variables independientes. Repita todos los análisis realizados en este documento.

  2. Considere los datos hbsdemo, tomando a prog como variable dependiente y gender, schtyp y awards como variables independientes. Repita todos los análisis realizados en este documento.

22.0.8 Ejercicios 22 a 24

  1. Considere los datos hbsdemo, tomando a prog como variable dependiente y ses, schtyp y honors como variables independientes. Repita todos los análisis realizados en este documento.

  2. Considere los datos hbsdemo, tomando a prog como variable dependiente y ses, schtyp y awards como variables independientes. Repita todos los análisis realizados en este documento.

  3. Considere los datos hbsdemo, tomando a prog como variable dependiente y schtyp, honors y awards como variables independientes. Repita todos los análisis realizados en este documento.

22.0.9 Ejercicios 25 a 27

  1. Considere los datos hbsdemo, tomando a prog como variable dependiente y gender, ses, schtyp y honors como variables independientes. Repita todos los análisis realizados en este documento.

  2. Considere los datos hbsdemo, tomando a prog como variable dependiente y gender, ses, schtyp y awards como variables independientes. Repita todos los análisis realizados en este documento.

  3. Considere los datos hbsdemo, tomando a prog como variable dependiente y gender, ses, schtyp, honors y awards como variables independientes. Repita todos los análisis realizados en este documento.

Bibliografía

Consultar el documento RPubs :: Regresión logística (bibliografía).

 

  
If you found any ERRORS or have SUGGESTIONS, please report them to my email. Thanks.