hllinas2023

1 Paquetes utilizados

En estas notas de clase utilizaremos, de manera especial, los siguientes paquetes:

  • glsm: paquete de autoría propia (Llinás H, Villalba J, Borja J, Tilano J (2025)) que permite calcular de forma sencilla y precisa el log-verosímil del modelo saturado cuando la variable respuesta Y es multinomial (tiene más de dos niveles). Además, incluye conjuntos de datos de ejemplo útiles para ilustrar su aplicación.

  • tidyverse: una colección de paquetes útiles para la manipulación y visualización de datos (dplyr, ggplot2, tibble, entre otros), que permite mantener una sintaxis coherente y eficiente a lo largo del análisis.

library(glsm)      #Regresión logística (caso multinomial)
library(tidyverse) #Incluye a dplyr y ggplot2
library(plotly)
library(plyr)
library(jmv)       #Análisis descriptivos

2 Otros 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á en este documento). En este documento, utilizaremos algunos de ellos para algunos cálculos.

2.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.

2.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 ordinal (Christensen 2018), para categorías ordenadas de la variable de respuesta.

2.0.3 Sobre el paquete glsm

  • El paquete glsm, desarrollado por Llinás H, Villalba J, Borja J, Tilano J (2025), es una contribución propia que puede descargarse desde el repositorio oficial de CRAN. Su propósito principal es facilitar la estimación del log-verosimilitud del modelo saturado, herramienta clave en la evaluación comparativa de modelos logísticos.

  • Cuando la variable de respuesta \(Y\) es multinomial (tiene más de dos niveles), la función glsm() calcula automáticamente el valor del log-verosímil del modelo saturado. Este valor permite comparar otros modelos y evaluar su bondad de ajuste relativa.

  • El paquete glsm no solo incluye funciones especializadas, sino también un conjunto de dato integrado que permiten realizar pruebas y ejemplos de manera práctica: hsbdemo.

  • Más detalles sobre este paquete se describen en siguientes documentos:

3 Introducción

Los métodos de regresión se han convertido en una herramienta fundamental dentro del análisis de datos, ya que permiten describir la relación entre una variable de respuesta y una o más variables explicativas. Con frecuencia, la variable de respuesta es discreta, asumiendo dos o más categorías posibles.

En estos casos, el modelo de regresión logística es uno de los más utilizados, especialmente cuando el interés es analizar la probabilidad de ocurrencia de cada categoría de la variable dependiente.

En el documento Rpubs :: Modelos lineales generalizados se explicó que la regresión logística forma parte de los Modelos Lineales Generalizados (GLM), y en Rpubs :: Regresión Logística binaria se abordó el caso binario.

En este documento se estudiará el caso multinomial, en el cual la variable de respuesta puede asumir tres o más categorías. Para comprender en profundidad estos modelos, resulta esencial revisar cuatro casos particulares que sirven como base conceptual:

  • Modelo de Bernoulli.

  • Modelo completo.

  • Modelo nulo.

  • Modelo saturado.

A continuación, se presentarán sus propiedades junto con ejemplos ilustrativos que facilitarán la comprensión del modelo logístico multinomial.

4 Datasets

Para las aplicaciones, se utilizará la base datos hbsdemo de la librería glsm Llinás H, Villalba J, Borja J, Tilano J (2025).

datos <- glsm::hsbdemo
names(datos)
##  [1] "Student" "id"      "gender"  "ses"     "schtyp"  "prog"    "read"   
##  [8] "write"   "math"    "science" "socst"   "honors"  "awards"  "cid"    
## [15] "prog0"   "prog1"   "prog2"

El conjunto de datos contiene variables sobre 200 estudiantes. Los estudiantes que ingresan a la escuela secundaria hacen la elección de un programa (prog) de tres posibles: general (general), vocacional (vocation) y académico (academic). 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: general, vocaction y académic, la cual será nuestra variable de respuesta.

  6. read, write, math, science y socst son variables continuas que representan los puntajes en lectura, escritura, 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

5 Ejemplo 1: enunciado (exploratorio)

Example 5.1 Considerelos datos hsbdemo. Explore la relación entre prog y otras variables, siguiendo las instrucciones que se proponen abajo.

  1. Escriba un resumen de los datos.

  2. Construya una tabla de frecuencia para prog.

  3. Construya un diagrama de barras para prog.

  4. Construya una tabla de frecuencia cruzada (tablas de contingencia) entre ses y prog.

  5. Cargue el paquete jmv y utilice la función descriptives para conseguir datos descriptivos con algunas variables, digamos, prog, ses, math, science.

  6. Construya un diagrama de dispersión bivariado entre prog y awards, dentro de cada categoría de ses.

  7. Construya una tabla de frecuencias cruzadas de progr, honors y gender.

  8. Construya un diagrama de barras que muestre la distribución de honors dentro progr, para cada nivel de gender.

  9. Construya un diagrama circular que muestre el porcentaje de premios recibidos (awards) para cada tipo de programa (prog).

  10. Construya un diagrama dispersión bivariado entre los puntajes de lectura (read) y escritura (write) para cada tipo de programa (prog).

6 Ejemplo 1: solución

6.0.1 Solución inciso (a)

El resumen de los datos se puede encontrar así:

summary(datos)
Student id read write math science socst awards cid
Min. : 1.00 Min. : 1.00 Min. :28.00 Min. :31.00 Min. :33.00 Min. :26.00 Min. :26.00 Min. :0.00 Min. : 1.00
1st Qu.: 50.75 1st Qu.: 50.75 1st Qu.:44.00 1st Qu.:45.75 1st Qu.:45.00 1st Qu.:44.00 1st Qu.:46.00 1st Qu.:0.00 1st Qu.: 5.00
Median :100.50 Median :100.50 Median :50.00 Median :54.00 Median :52.00 Median :53.00 Median :52.00 Median :1.00 Median :10.50
Mean :100.50 Mean :100.50 Mean :52.23 Mean :52.77 Mean :52.65 Mean :51.85 Mean :52.41 Mean :1.67 Mean :10.43
3rd Qu.:150.25 3rd Qu.:150.25 3rd Qu.:60.00 3rd Qu.:60.00 3rd Qu.:59.00 3rd Qu.:58.00 3rd Qu.:61.00 3rd Qu.:2.00 3rd Qu.:15.00
Max. :200.00 Max. :200.00 Max. :76.00 Max. :67.00 Max. :75.00 Max. :74.00 Max. :71.00 Max. :7.00 Max. :20.00

6.0.2 Solución inciso (b)

La tabla de frecuencia es:

Tabla <- datos %>%
  dplyr::group_by(prog) %>%
  dplyr::summarise(Total = n()) %>%
  dplyr::mutate(Porcentaje = round(Total/sum(Total)*100, 2))
Tabla
prog Total Porcentaje
academic 105 52.5
general 45 22.5
vocation 50 25.0

6.0.3 Solución inciso (c)

El diagrama de barras correspondiente es:

ggplot(Tabla, aes(x = factor(prog), y=Total,fill=factor(prog)) ) + 
  
  geom_bar(width = 0.7,stat="identity", position = position_dodge())+
  
  geom_text(aes(label=paste0(Total," ", "", "\n(", Porcentaje, "%)")), vjust=-0.2, color="black",                     hjust=0.5, position = position_dodge(0.9),  angle=0, size=4.0) +
  
  ylim(c(0,150))+
  
  labs(x="Tipo de programa elegido", y= "Frecuencias \n (Porcentajes)")+ 
  labs(fill = "")+
   
  facet_wrap(~"Distribución del tipo de programa (n=159)") +
  
  theme_bw(base_size = 13)

6.0.4 Solución inciso (d)

Una tabla de frecuencias cruzada entre ses y prog es:

with(datos, addmargins(table(ses, prog)))
academic general vocation Sum
high 42 9 7 58
low 19 16 12 47
middle 44 20 31 95
Sum 105 45 50 200

6.0.5 Solución inciso (e)

Con jmv::descritptives podemos construir las tablas de frecuencias solicitadas:

jmv::descriptives(datos, vars = vars(prog, ses, math, science), freq = TRUE)
## 
##  DESCRIPTIVES
## 
##  Descriptives                                                  
##  ───────────────────────────────────────────────────────────── 
##                          prog    ses    math        science    
##  ───────────────────────────────────────────────────────────── 
##    N                      200    200         200         200   
##    Missing                  0      0           0           0   
##    Mean                                 52.64500    51.85000   
##    Median                               52.00000    53.00000   
##    Standard deviation                   9.368448    9.900891   
##    Minimum                              33.00000    26.00000   
##    Maximum                              75.00000    74.00000   
##  ─────────────────────────────────────────────────────────────

6.0.6 Solución inciso (f)

El diagrama solicitado es:

library(plyr)

resp <- datos
resp$prog <- revalue(resp$prog,c("general"="1","vocation"="2", "academic"="3"))

ggplot(resp, aes(y = prog, x=awards) ) +
  geom_point(aes(color = ses), size=2,  alpha=5) +
  
  labs(y="Tipo de programa elegido", x= "Número de premios recibido por cada estudiante")+ 
  labs(fill = "Estado socioeconómico") +
  
  facet_wrap(vars(ses)) +
  theme_bw(base_size = 13) +
  theme(strip.background = element_rect(fill="skyblue"))

6.0.7 Solución inciso (g)

La tabla de frecuencias solicitada es:

Tabla <- datos %>%
  dplyr::group_by(prog, honors, gender) %>%
  dplyr::summarise(Total = n()) %>%
  dplyr::mutate(Porcentaje = round(Total/sum(Total)*100, 2)) %>%
  dplyr::arrange(prog)
Tabla
Table 6.1: Tabla cruzada para prog, gender y honors.
prog honors gender Total Porcentaje
academic enrolled female 25 62.50
academic enrolled male 15 37.50
academic not enrolled female 33 50.77
academic not enrolled male 32 49.23
general enrolled female 5 71.43
general enrolled male 2 28.57
general not enrolled female 19 50.00
general not enrolled male 19 50.00
vocation enrolled female 5 83.33
vocation enrolled male 1 16.67
vocation not enrolled female 22 50.00
vocation not enrolled male 22 50.00

6.0.8 Solución inciso (h)

Un diagrama de barras que muestra la distribución de honors dentro progr, para cada nivel de gender es:

ggplot(Tabla, aes(x = factor(prog), y=Total,fill=factor(gender)) ) + 
  
  geom_bar(width = 0.9,stat="identity", position = position_dodge())+
  geom_text(aes(label=paste0(Total," ", "", "(", Porcentaje, "%)")), 
            color="black", 
            hjust=-0.15,
            position = position_dodge(0.9),  
            angle=90, 
            size=3.5)+
  
  ylim(c(0,50))+
  
  labs(x="Tipo de programa elegido", y= "Frecuencias \n (Porcentajes)")+ 
  labs(fill = "Genero") +
  
  facet_wrap(vars(honors)) +
  
  theme(strip.background = element_rect(fill="blanchedalmond"))

6.0.9 Solución inciso (i)

Un diagrama circular que muestra el porcentaje de premios recibidos (awards) para cada tipo de programa (prog) es:

library(dplyr)
library(plotly)

datos %>%
  mutate(prog = as.character(prog)) %>%
  dplyr::count(prog, name = "n") %>%
  plot_ly(labels = ~prog, values = ~n, type = "pie") %>%
  layout(title = "Distribución de estudiantes por tipo de programa")
library(dplyr)
library(plotly)

datos %>%
  group_by(prog) %>%
  summarise(awards = sum(awards, na.rm = TRUE), .groups = "drop") %>%
  plot_ly(labels = ~as.character(prog), values = ~awards, type = "pie") %>%
  layout(title = "Porcentaje del total de awards por tipo de programa")
library(dplyr)
library(plotly)

datos %>%
  group_by(prog) %>%
  summarise(prop_premiados = mean(awards > 0, na.rm = TRUE), .groups = "drop") %>%
  plot_ly(labels = ~as.character(prog), values = ~prop_premiados, type = "pie") %>%
  layout(title = "Proporción de estudiantes con ≥1 award por programa")

6.0.10 Solución inciso (j)

Un diagrama dispersión bivariado entre los puntajes de lectura (read) y escritura (write) para cada tipo de programa (prog) es:

library(plotly)
plot_ly(data = datos, x = ~read, y = ~write, color = ~prog)  %>% 
     layout(title = 'Puntajes en lectura vs escritura, sectorizada por tipo de programa')

7 Modelo básico

7.0.1 El modelo

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.

7.0.2 Ejemplo 2: datos con una variable \(Y\) multinomial

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

Table 7.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\).

7.0.3 La función de verosimilitud

Theorem 7.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{7.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.

8 Modelo completo

8.0.1 El modelo y sus estimaciones

Definition 8.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 8.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\]

8.0.2 Ejemplo 3: modelo completo

Example 8.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 (manual).

En R, se puede verificar así:

datos <- glsm::hsbdemo

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

Con el paquete glsm

El paquete glsm también permite calcular directamente el valor del log-verosímil para el modelo completo, denotado como \(\mathcal{L}(y)\). Este modelo corresponde al caso en el que cada observación se ajusta de manera perfecta, es decir, la probabilidad asignada a cada \(y_i\) coincide exactamente con el valor observado.

Podemos acceder a este valor de dos formas:

datos <- glsm::hsbdemo
modelo <- glsm(prog ~ gender + ses, data = datos, ref = "academic")

modelo                    # Opción 1
modelo$Log_Lik_Complete   # Opción 2

Con glsm: resultado con la opción 1.

Como se observa, el valor de \(\mathcal{L}(y)\) se muestra como Complete en el bloque de salida del modelo:

## 
## Call:
## glsm(formula = prog ~ gender + ses, data = datos, ref = "academic")
## 
## Populations in Saturated Model: 6
## 
## Coefficients: 
##                         Coef(B) Std.Error    Exp(B)
## (Intercept):general  -1.6547758 0.4175354 0.1911349
## (Intercept):vocation -1.8469099 0.4478055 0.1577238
## gendermale:general    0.2216624 0.3716764 1.2481500
## gendermale:vocation   0.1098969 0.3599743 1.1161630
## seslow:general        1.4118732 0.5064763 4.1036349
## seslow:vocation       1.3534875 0.5548849 3.8709018
## sesmiddle:general     0.7550865 0.4561360 2.1277956
## sesmiddle:vocation    1.4430949 0.4709456 4.2337787
## 
## Log Likelihood: 
##          Estimation
## Complete     0.0000
## Null      -204.0967
## Logit     -195.5208
## Saturate  -194.5159

Con glsm: resultado con la opción 2.

Este valor también puede ser extraído directamente usando modelo$Log_Lik_Complete:

datos <- glsm::hsbdemo
modelo <- glsm(prog ~ gender + ses, data = datos, ref = "academic")
modelo$Log_Lik_Complete   # Opción 2
## [1] 0

El valor de \(\mathcal{L}(y)\) constituye un punto de referencia teórico, ya que representa el mejor ajuste posible: un modelo que reproduce exactamente las observaciones. Por esta razón, se utiliza como extremo superior para comparar la calidad de ajuste de otros modelos, como el modelo logístico o el modelo nulo.

9 Modelo nulo

9.0.1 El modelo y sus estimaciones

Definition 9.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 9.1 (Log-verosimilitud en el modelo nulo) En este caso, (7.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{9.1} \end{eqnarray}\]

El siguiente teorema describe las estimaciones en este modelo:

Theorem 9.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\]

10 Ejemplo 4: enunciado (modelo nulo)

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. Esas estimaciones con la función glsm :: glsm.

  5. 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.

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

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

  8. Las estimaciones pedidas en (e), con la función glsm :: glsm.

11 Ejemplo 4: solución

11.0.1 Solución inciso (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\]

11.0.2 Solución inciso (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 11.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 11.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

11.0.3 Solución inciso (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 11.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 11.2: Modelo nulo. Fuente: Elaboración propia.

11.0.4 Solución inciso (d)

El modelo nulo corresponde al caso en el que la variable de respuesta prog se modela únicamente con interceptos, es decir, sin predictores adicionales. En este escenario, cada categoría de prog tiene asociada una constante que refleja su proporción en la muestra. Al ejecutar código de abajo, se obtiene la estimación de los parámetros del modelo nulo:

#PENDIENTE. NO FUNCIONA. 
datos <- glsm::hsbdemo

modelo <- glsm(prog ~ 1, data = datos, ref = "academic")
modelo

11.0.5 Solución inciso (e)

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

11.0.6 Solución inciso (f)

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 11.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 11.3: Modelo nulo. Fuente: Elaboración propia.

11.0.7 Solución inciso (g)

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 11.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 11.4: Modelo nulo. Fuente: Elaboración propia.

11.0.8 Solución inciso (h)

El paquete glsm permite calcular directamente el valor del log-verosímil para diferentes modelos, incluyendo el modelo nulo, denotado como \({\cal L}(\overline{y})\), que representa el valor de la log-verosimilitud bajo el supuesto de que todos los casos tienen la misma probabilidad de éxito. Podemos acceder a este valor de dos formas:

datos <- glsm::hsbdemo
modelo <- glsm(prog ~ gender + ses, data = datos, ref = "academic")

modelo               # Opción 1
modelo$Log_Lik_Null  # Opción 2

Con glsm: resultado con la opción 1.

Como se observa, el valor de \(\mathcal{L}(\overline{y})\) se muestra como Null en el bloque de salida del modelo:

## 
## Call:
## glsm(formula = prog ~ gender + ses, data = datos, ref = "academic")
## 
## Populations in Saturated Model: 6
## 
## Coefficients: 
##                         Coef(B) Std.Error    Exp(B)
## (Intercept):general  -1.6547758 0.4175354 0.1911349
## (Intercept):vocation -1.8469099 0.4478055 0.1577238
## gendermale:general    0.2216624 0.3716764 1.2481500
## gendermale:vocation   0.1098969 0.3599743 1.1161630
## seslow:general        1.4118732 0.5064763 4.1036349
## seslow:vocation       1.3534875 0.5548849 3.8709018
## sesmiddle:general     0.7550865 0.4561360 2.1277956
## sesmiddle:vocation    1.4430949 0.4709456 4.2337787
## 
## Log Likelihood: 
##          Estimation
## Complete     0.0000
## Null      -204.0967
## Logit     -195.5208
## Saturate  -194.5159

Con glsm: resultado con la opción 2.

Aquí puede ser extraído directamente usando modelo$Log_Lik_Null:

## [1] -204.0967

Este valor es útil para comparar con la log-verosimilitud de otros modelos (como el modelo logístico o el modelo saturado) y así evaluar la calidad del ajuste. Cuanto mayor (menos negativo) sea el valor de la log-verosimilitud, mejor será el ajuste del modelo a los datos observados.

12 Modelo saturado

12.0.1 Supuesto 1

El modelo saturado está caracterizado por dos supuestos.

Hypothesis 12.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}\]

12.0.2 Ejemplo 5: tabla para el modelo saturado

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

Table 12.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\).

12.0.3 Supuesto 2

Hypothesis 12.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})\]

12.0.4 Implicaciones del supuesto 2

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.

12.0.5 Log L y estimaciones

Theorem 12.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{12.1} \end{eqnarray}\]

Theorem 12.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\]

12.0.6 Ejemplo 6: estimaciones en el modelo saturado

Example 12.1 Para los datos del archivo hsbdemo, en el modelo saturado, hay \(J=6\) poblaciones y se cumple que \({\cal L}(\tilde{p})=-194.5159\).

En R (manual).

En R, se puede verificar así (los valores obtenidos se pueden ver en la Tabla 12.2):

datos <- glsm::hsbdemo

datos %>%
  dplyr::group_by(gender,ses) %>%
  dplyr::summarise(nj = n(),
            z0j = sum(prog0),
            z1j = sum(prog1),
            z2j = sum(prog2)) %>%
  dplyr::mutate(p0j = round(z0j/nj,4),
         p1j = round(z1j/nj,4),
         p2j = round(z2j/nj,4),
         LogL_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)),
         LogL_ref2 = round(LogL_ref2, 4)) -> saturado
        
LogL_saturado <- sum(saturado$LogL_ref2)
## [1] -194.5159
Table 12.2: Estimación en el modelo saturado: \({\cal L}(\tilde{p})= \sum\limits_{j=1}^J {\cal L}_j(\tilde{p}) =-194.5159\)
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
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

Con el paquete glsm.

El paquete glsm permite calcular directamente los valores de:

  • \(J\): número de poblaciones únicas (combinaciones de predictores) en el modelo saturado.
  • \(\mathcal{L}(\tilde{p})\): log-verosimilitud del modelo saturado, es decir, el valor máximo que puede alcanzar la log-verosimilitud cuando el modelo reproduce exactamente cada observación.

Podemos obtener ambos valores (o visualizar la tabla completa) podemos proceder de la siguiente manera:

datos  <- glsm::hsbdemo
modelo <- glsm(prog ~ gender + ses, data = datos, ref = "academic")

modelo                                             # Opción 1
cbind(modelo$Populations, modelo$Log_Lik_Saturate) # Opción 2
modelo$Esm                                         # Opción 3

Con glsm: resultado con la opción 1.

En la salida del modelo, los valores de \(J\) y \(\mathcal{L}(\tilde{p})\) aparecen como Populations in Saturate Model (parte superior) y Saturate (en el bloque de log-likelihood):

## 
## Call:
## glsm(formula = prog ~ gender + ses, data = datos, ref = "academic")
## 
## Populations in Saturated Model: 6
## 
## Coefficients: 
##                         Coef(B) Std.Error    Exp(B)
## (Intercept):general  -1.6547758 0.4175354 0.1911349
## (Intercept):vocation -1.8469099 0.4478055 0.1577238
## gendermale:general    0.2216624 0.3716764 1.2481500
## gendermale:vocation   0.1098969 0.3599743 1.1161630
## seslow:general        1.4118732 0.5064763 4.1036349
## seslow:vocation       1.3534875 0.5548849 3.8709018
## sesmiddle:general     0.7550865 0.4561360 2.1277956
## sesmiddle:vocation    1.4430949 0.4709456 4.2337787
## 
## Log Likelihood: 
##          Estimation
## Complete     0.0000
## Null      -204.0967
## Logit     -195.5208
## Saturate  -194.5159

Con glsm: resultado con la opción 2.

Ambos valores también pueden extraerse directamente desde el objeto del modelo (usando modelo$Populations y modelo$Log_Lik_Saturate):

## Número de poblaciones (J): 6
## Log-verosimilitud del modelo saturado (LogL): -194.5159

El valor \(\mathcal{L}(\tilde{p})\) sirve como referencia superior para evaluar la calidad del ajuste de otros modelos (por ejemplo, el nulo o el logístico): cualquier modelo ajustado tendrá una log-verosimilitud menor o igual a la del modelo saturado.

Con glsm: resultado con la opción 3.

gender ses n z_academic_j z_general_j z_vocation_j p_academic_j_tilde p_general_j_tilde p_vocation_j_tilde Lp
female high 29 21 5 3 0.7241379 0.1724138 0.1034483 -22.37358
female low 32 15 9 8 0.4687500 0.2812500 0.2500000 -33.87224
female middle 48 22 10 16 0.4583333 0.2083333 0.3333333 -50.42744
male high 29 21 4 4 0.7241379 0.1379310 0.1379310 -22.62625
male low 15 4 7 4 0.2666667 0.4666667 0.2666667 -15.90903
male middle 47 22 10 15 0.4680851 0.2127660 0.3191489 -49.30740

13 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.

  • Haga un análisis exploratorio diferente de con los datos mencionados.

13.0.1 Ejercicios 1 a 3

  1. Demuestre los teoremas: (a) 8.1; (b) 9.2; (c) 12.2.

  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 estimen el logaritmo de la función de máxima verosimilitud en el modelo saturado.

13.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.

13.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.

13.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.

13.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.

13.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.

13.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.

13.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.

13.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.