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.

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) y Rpbus :: Modelos logísticos (intervalos de confianza) se detalló todo lo relacionado con las estimaciones e intervalos de confianza para los parámetros logísticos. 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.

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 Modelo básico

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

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

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

5.0.3 La función de verosimilitud

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

6.0.1 El modelo y sus estimaciones

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\]

6.0.2 Ejemplo 3: modelo completo

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

7 Modelo nulo

7.0.1 El modelo y sus estimaciones

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

9 Ejemplo 4: solución

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

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

9.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 9.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 el Inciso (b):

library(VGAM)

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

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

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

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

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

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

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

10 Modelo saturado

10.0.1 Supuesto 1

El modelo saturado está caracterizado por dos supuestos.

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

10.0.2 Ejemplo 5: tabla para el modelo saturado

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

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

10.0.3 Supuesto 2

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

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

10.0.5 Log L y estimaciones

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

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

10.0.6 Ejemplo 6: estimaciones en el modelo saturado

Example 10.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 10.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 10.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

11 Modelo logístico

11.0.1 Recordatorio: Rango y rango completo de una matriz

Definición

El rango de una matriz es el número máximo de columnas (o filas) linealmente independientes. En otras palabras, mide cuánta información nueva aportan las columnas de la matriz entre sí.

Una matriz tiene rango completo cuando su rango es igual al número de columnas (o al número de filas, si la matriz tiene más columnas que filas). Esto significa que no hay colinealidad exacta entre las columnas.

Ejemplo 1: Matriz con rango completo

\[ A = \begin{pmatrix} 1 & 2 \\ 3 & 4 \\ \end{pmatrix} \]

Las columnas de \(A\) no son múltiplos una de la otra. Por tanto, \(Rg(A) = 2\) y la matriz tiene rango completo. En este caso, como es una matriz cuadrada, también será invertible.

Ejemplo 2: Matriz con columnas linealmente dependientes

\[ B = \begin{pmatrix} 1 & 2 \\ 2 & 4 \\ \end{pmatrix} \]

Aquí, la segunda columna es el doble de la primera, es decir, hay colinealidad. Entonces \(Rg(B) = 1\), y no tiene rango completo.

11.0.2 Supuestos

Hypothesis 11.1 (Supuesto 3: matriz de diseño) Se hacen los supuestos 1 y 2 del modelo saturado (véase las hipótesis 10.1 y 10.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\). Esto significa que todas las columnas de la matriz \(C\) son linealmente independientes, es decir:

  • No hay redundancia entre las variables (por ejemplo, una covariable no puede ser combinación exacta de otras).

  • El modelo puede estimar todos los parámetros asociados a las covariables sin problemas de colinealidad.

  • El número de combinaciones observadas \(J\) debe ser mayor o igual al número de parámetros a estimar \((1 + K)\), para que haya información suficiente.

En la práctica:

Para poder ajustar el modelo logístico sin problemas, debe cumplirse que el número total de combinaciones observadas sea al menos igual al número de parámetros a estimar. Si no se cumple, el sistema de ecuaciones para los parámetros no tiene solución única y el modelo fallará.

11.0.3 El modelo

Hypothesis 11.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{11.1} \end{eqnarray}\]

Remark. Tenemos que:

  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\). Esto quiere decir que el parámetro \(\alpha\) se puede estimar de manera única a partir de los datos.

12 Riesgo

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

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

13 La función Log L en el modelo logístico

Theorem 13.1 (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{13.1} \end{eqnarray}\]

14 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) o el teorema 6 en Orozco, Llinás y Fonseca (2020). El paquete glsm (Llinás H, Villalba J, Borja J, Tilano J (2025)) calcula estas estimaciones.

15 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:

  1. \(n\to\infty\).

  2. \(J\to\infty\), caso que es únicamente es apropiado para datos no agrupados.

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

16 Ejemplo 5: 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 las funciones:

  • summary :: multinom :: nnet.

  • summary :: vglm :: VGAM.

  • summary :: glsm :: glsm.

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

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

  • multinom :: nnet.

  • glsm :: glsm.

  1. Halle las razones odds correspondientes. Utilize las funciones:

  • multinom :: nnet.

  • glsm :: glsm.

  1. Calcular las probabilidades predichas para cada uno de los niveles de la variable de respuesta. Utilize las funciones:

  • multinom :: nnet.

  • glsm :: glsm.

  1. Mantener fija la variable write=35 y examinar las probabilidades predichas para cada nivel de ses. Utilize las funciones:

  • multinom :: nnet.

  • glsm :: glsm.

  1. Mantener fija la variable write en su media y examinar las probabilidades predichas para cada nivel de ses. Utilize las funciones:

  • multinom :: nnet.

  • glsm :: glsm.

  1. Examinar los cambios en la probabilidad predicha asociados con cada nivel de ses, manteniendo write fija en su media. Utilize las funciones:

  • multinom :: nnet.

  • glsm :: glsm.

  1. Examinar las probabilidades predichas promediadas para diferentes valores de la variable predictora continua write dentro de cada nivel de ses. Utilize las funciones:

  • multinom :: nnet.

  • glsm :: glsm.

  1. 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. Utilize las funciones:

  • multinom :: nnet.

  • glsm :: glsm.

  1. Halle los errores estándares estimados de los estimadores de los parámetros logísticos, utilizando las funciones:

  • multinom :: nnet.

  • vglm :: VGAM.

  • glsm :: glsm.

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

  • multinom :: nnet.

  • vglm :: VGAM.

  • glsm :: glsm.

17 Ejemplo 5: Solución

17.0.1 Solución inciso (a)

Los vector de parámetros y 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.

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

17.0.3 Solución inciso (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 (11.1) o la (12.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}}} \]

17.0.4 Solución inciso (d)

Con la función summary :: multinom :: nnet

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 17.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 17.1: Estimaciones de los parámetros logísticos con nnet. Fuente: Elaboración propia.

Con la función summary :: vglm :: VGAM

En R, las estimaciones de los parámetros logísticos \(\delta_r\) y \(\beta_r\) también pueden obtenerse con vglm() del paquete VGAM. Antes de ajustar el modelo, es indispensable fijar el nivel de referencia de la respuesta; aquí usaremos prog = "academic". En la salida de summary(), deben interpretarse los coeficientes por categoría (en relación con la referencia) y sus errores estándar, z-values y \(p\)-valores. El enlace por defecto en multinomial() corresponde al baseline-category logit.

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 17.2: Estimaciones de los parámetrdos logísticos con nnet. Fuente: Elaboración propia.

Con la función summary :: glsm :: glsm

Como alternativa a VGAM y nnet, el paquete glsm implementa el modelo logístico multinomial (baseline-category logit) con una interfaz orientada a docencia y reporte. Al igual que antes, se debe fijar la categoría de referencia. La función glsm() entrega un resumen con coeficientes por categoría, errores estándar, estadísticos de Wald, \(p\)-valores, y medidas de verosimilitud útiles para comparar con los modelos nulo, completo y saturado.

datos <- glsm::hsbdemo

# Ajuste (cambia 'ref' si deseas otra categoría de referencia)
modelo <- glsm(prog ~ ses + write, data = datos, ref = "academic")

# Salida 1: Resumen tabular (coeficientes, SE, Wald, p, Exp(B)):
summary(modelo)

# Salida 2:  Coeficientes en escala logit:
coef(modelo)

Salida 1:

## 
## Call:
## glsm(formula = prog ~ ses + write, data = datos, ref = "academic")
## 
## Coefficients: 
##                        Coef(B) Std.Error    Exp(B)      Wald DF   P.value    
## (Intercept):general   1.689354  1.226940  5.415982  1.895809  1 0.1685481    
## (Intercept):vocation  4.235530  1.204650 69.098279 12.362147  1 0.0004381 ***
## seslow:general        1.162832  0.514220  3.198980  5.113707  1 0.0237376 *  
## seslow:vocation       0.982670  0.595530  2.671581  2.722757  1 0.0989270 .  
## sesmiddle:general     0.629541  0.465028  1.876749  1.832691  1 0.1758101    
## sesmiddle:vocation    1.274063  0.511065  3.575351  6.214826  1 0.0126685 *  
## write:general        -0.057928  0.021411  0.943718  7.320093  1 0.0068188 ** 
## write:vocation       -0.113603  0.022219  0.892613 26.141569  1 3.173e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Analysis of Deviance (Chi-squared): 
##                   Deviance  DF   P.value    
## Null vs Logit        48.23   6 1.063e-08 ***
## Logit vs Complete   359.96 392    0.8755    
## Logit vs Saturate   226.18 120 1.615e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Salida 2:

##  (Intercept):general (Intercept):vocation       seslow:general 
##           1.68935423           4.23552982           1.16283199 
##      seslow:vocation    sesmiddle:general   sesmiddle:vocation 
##           0.98267029           0.62954098           1.27406340 
##        write:general       write:vocation 
##          -0.05792841          -0.11360264

17.0.5 Solución inciso (e)

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 (11.1) o la (12.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*}\]

17.0.6 Solución inciso (f)

La gráfica correspondiente se muestra en la figura 17.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.

Con el paquete multinom :: nnet.

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"))

Con el paquete glsm :: glsm.

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

B <- as.matrix(coef(modelo))  

W <- seq(0, 100, 0.05)

g_0 <- B[1,1] + B[3,1] + B[7,1]*W
g_1 <- B[2,1] + B[4,1] + B[8,1]*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"))
Comparación de dos Logits

Figure 17.3: Comparación de dos Logits

17.0.7 Solución inciso (g)

En R, las estimaciones de los odds se pueden obtener con las funciones multinom :: nnet, vglm :: VGAM o glsm::glsm. Simplemente debe utilizar la función exp() a los coeficientes del modelo, como se resalta en el recuadro rojo de la figura 17.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().

Con el paquete multinom :: nnet.

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 17.4: Estimaciones de los odds. Fuente: Elaboración propia.

Con el paquete glsm :: glsm.

library(glsm)

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

exp(coef(modelo))     
##  (Intercept):general (Intercept):vocation       seslow:general 
##            5.4159821           69.0982789            3.1989800 
##      seslow:vocation    sesmiddle:general   sesmiddle:vocation 
##            2.6715806            1.8767489            3.5753512 
##        write:general       write:vocation 
##            0.9437175            0.8926126

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

17.0.8 Solución inciso (h)

Se pueden calcular las probabilidades predichas para cada uno de los niveles de la variable de respuesta.

Con el paquete multinom :: nnet.

Utilizando la función fitted():

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

head(fitted(modelo))
## # weights:  15 (8 variable)
## initial  value 219.722458 
## iter  10 value 179.983731
## final  value 179.981726 
## converged
Table 17.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

Con el paquete glsm :: glsm.

Utilizando la función modelo$Elm:

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

a <- modelo$Elm

kable(head(a[,c(1,2,7, 8, 9)]))
ses write p_academic_j p_general_j p_vocation_j
high 33 0.2917550 0.2336048 0.4746402
high 38 0.3966377 0.2377210 0.3656413
high 41 0.4631025 0.2332796 0.3036179
high 42 0.4852972 0.2307010 0.2840018
high 49 0.6324619 0.2004334 0.1671047
high 52 0.6876337 0.1831551 0.1292113

17.0.9 Solución inciso (i)

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

Con el paquete multinom :: nnet.

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

s<- summary(modelo)

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
## # weights:  15 (8 variable)
## initial  value 219.722458 
## iter  10 value 179.983731
## final  value 179.981726 
## converged
## [1] 0.3382509 0.5134769 0.1482721

Con el paquete glsm :: glsm.

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

B <- as.matrix(coef(modelo))  

W <- 35

g_0 <- B[1,1] + B[3,1] + B[7,1]*W
g_1 <- B[2,1] + B[4,1] + B[8,1]*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
##  (Intercept):general (Intercept):vocation  (Intercept):general 
##            0.3382488            0.5134731            0.1482781

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

Con el paquete multinom :: nnet.

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

Nuevo <- data.frame(ses = c("low", "middle", "high"), write = mean(write))
Tabla <- predict(modelo, newdata = Nuevo, "probs")
Tabla
## # weights:  15 (8 variable)
## initial  value 219.722458 
## iter  10 value 179.983731
## final  value 179.981726 
## converged
##    academic   general  vocation
## 1 0.4396813 0.3581915 0.2021272
## 2 0.4777451 0.2283359 0.2939190
## 3 0.7009046 0.1784928 0.1206026

Con el paquete glsm :: glsm.

#Error in UseMethod("predict"):  
#No applicable method for 'predict' applied to an object of class "glsm"

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

Nuevo <- data.frame(ses = c("low", "middle", "high"), write = mean(write))
Tabla <- predict(modelo, newdata = Nuevo, "probs")
Tabla

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

Con el paquete multinom :: nnet.

library(nnet)
prog <- relevel(as.factor(prog), ref = "academic")
modelo <- multinom(prog ~ ses + write, data=datos)
## # weights:  15 (8 variable)
## initial  value 219.722458 
## iter  10 value 179.983731
## final  value 179.981726 
## converged
s<- summary(modelo)

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
## [1] 0.3581915 0.2021272 0.4396813

Con el paquete glsm :: glsm.

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

B <- as.matrix(coef(modelo))  

W <- mean(write)

g_0 <- B[1,1] + B[3,1] + B[7,1]*W
g_1 <- B[2,1] + B[4,1] + B[8,1]*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
##  (Intercept):general (Intercept):vocation  (Intercept):general 
##            0.3581927            0.2021232            0.4396841

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

Con el paquete multinom :: nnet.

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

# 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))
## # weights:  15 (8 variable)
## initial  value 219.722458 
## iter  10 value 179.983731
## final  value 179.981726 
## converged
Table 17.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

Con el paquete glsm :: glsm.

#Error in UseMethod("predict") : 
#  No applicable method for 'predict' applied to an object of class "glsm"

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

# 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))

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

Con el paquete multinom :: nnet.

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

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")
## # weights:  15 (8 variable)
## initial  value 219.722458 
## iter  10 value 179.983731
## final  value 179.981726 
## converged

Con el paquete glsm :: glsm.

# Error: objeto 'Predichos' no encontrado

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

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")+
  ylim(0, 1)+
  theme_bw(base_size = 12) +
  labs(y="Probabilidad")

17.0.14 Solución inciso (n)

Con el paquete multinom :: nnet.

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 17.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 17.5: Estimaciones de los errores estándares con nnet. Fuente: Elaboración propia.

Con el paquete vglm :: VGAM.

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 17.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 17.6: Estimaciones de los errores estándares con VGAM. Fuente: Elaboración propia.

Con el paquete glsm :: glsm.

En R, las estimaciones de los errores estándares también se pueden obtener con la función glsm() del paquete glsm. Con la función summary, se pueden visualizar en el bloque de Coeffcients, en la clumna de Std.Error (opción 1). También se pueden obtener con modelo$Std.Error (opción 2).

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

summary(modelo)  # Opción 1
modelo$Std.Error # Opción 2

Opción 1.

## 
## Call:
## glsm(formula = prog ~ ses + write, data = datos, ref = "academic")
## 
## Coefficients: 
##                        Coef(B) Std.Error    Exp(B)      Wald DF   P.value    
## (Intercept):general   1.689354  1.226940  5.415982  1.895809  1 0.1685481    
## (Intercept):vocation  4.235530  1.204650 69.098279 12.362147  1 0.0004381 ***
## seslow:general        1.162832  0.514220  3.198980  5.113707  1 0.0237376 *  
## seslow:vocation       0.982670  0.595530  2.671581  2.722757  1 0.0989270 .  
## sesmiddle:general     0.629541  0.465028  1.876749  1.832691  1 0.1758101    
## sesmiddle:vocation    1.274063  0.511065  3.575351  6.214826  1 0.0126685 *  
## write:general        -0.057928  0.021411  0.943718  7.320093  1 0.0068188 ** 
## write:vocation       -0.113603  0.022219  0.892613 26.141569  1 3.173e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Analysis of Deviance (Chi-squared): 
##                   Deviance  DF   P.value    
## Null vs Logit        48.23   6 1.063e-08 ***
## Logit vs Complete   359.96 392    0.8755    
## Logit vs Saturate   226.18 120 1.615e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Opción 2.

##  (Intercept):general (Intercept):vocation       seslow:general 
##           1.22694022           1.20464972           0.51422012 
##      seslow:vocation    sesmiddle:general   sesmiddle:vocation 
##           0.59552969           0.46502835           0.51106550 
##        write:general       write:vocation 
##           0.02141082           0.02221890

17.0.15 Solución inciso (o)

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 (13.1):

\[ {\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, vglm::VGAM O glsm :: glsm. Con cualquier camino encontramos que \({\cal L}(\hat{\alpha})=-179.9817\).

Con el paquete multinom :: nnet.

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

\[{\cal L}(\hat{\alpha}) \;=\; \frac{\text{Residual Deviance}}{-2} \;=\; -179.9817\]

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 17.7: Estimaciones del Log-Likelihood con nnet. Fuente: Elaboración propia.

Con el paquete vglm :: VGAM.

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 17.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 17.8: Estimaciones del Log-Likelihood con VGAM. Fuente: Elaboración propia.

Con la función glsm: utilizar modelo$Elm

La matriz Elm contiene las sigientes variables.

library(glsm)
datos <- glsm::hsbdemo
modelo <- glsm(prog ~ ses + write, data = datos, ref = "academic")
a <- modelo$Elm
names(a)
##  [1] "ses"                "write"              "n"                 
##  [4] "z_academic_j"       "z_general_j"        "z_vocation_j"      
##  [7] "p_academic_j"       "p_general_j"        "p_vocation_j"      
## [10] "logit_p_general_j"  "logit_p_vocation_j" "v_academic_j"      
## [13] "v_general_j"        "v_vocation_j"

Con ayuda de algunas de estas variables, el objetivo es calcular el valor del log-verosímil \({\cal L}(\hat{\alpha})\) del modelo logístico ajustado. En el caso multinomial con 3 categorías, tomando como referencia la categoría academic, los dos logits corresponden a general y vocation. El log-verosímil se calcularía así:

\[ {\cal L}(\hat{\alpha}) \;=\; \sum_{j}\Big[z_{0j}\hat{g}_{0j} + z_{1j}\hat{g}_{1j} - n_j \log\!\left(1 + e^{\hat{g}_{0j}} + e^{\hat{g}_{1j}} \right)\Big]. \]

Por eso, de la matriz Elm se identifican las siguientes variables:

  • n ≡ \(n_j\)
  • z_general_j ≡ \(z_{0j}\)
  • z_vocation_j ≡ \(z_{1j}\)
  • logit_p_general_j ≡ \(\hat g_{0j}\)
  • logit_p_vocation_j ≡ \(\hat g_{1j}\)

O sea, de esa matriz, solo nos interesa las variables 1:6 y 10:11. La instrucción modelo$Elm[, c(1:6, 10,11)] devuelve una matriz con la información esencial para cada población \(j\):

library(glsm)
datos <- glsm::hsbdemo
modelo <- glsm(prog ~ ses + write, data = datos, ref = "academic")
modelo$Elm[, c(1:6, 10,11)]
ses write n z_academic_j z_general_j z_vocation_j logit_p_general_j logit_p_vocation_j
high 33 2 1 1 0 -0.2222833 0.4866427
high 38 1 1 0 0 -0.5119253 -0.0813704
high 41 2 2 0 0 -0.6857106 -0.4221784
high 42 1 1 0 0 -0.7436390 -0.5357810
high 49 3 1 0 2 -1.1491378 -1.3309995
high 52 6 5 0 1 -1.3229231 -1.6718074

De esta forma, el log-verosímil se calcula como:

# Calcular L desde los logits guardados en Elm
with(modelo$Elm, {
  g0 <- logit_p_general_j
  g1 <- logit_p_vocation_j
  z0 <- z_general_j
  z1 <- z_vocation_j
  nj <- n

  LL_desde_Elm_logit <- sum(z0 * g0 + z1 * g1 - nj * log(1 + exp(g0) + exp(g1)))
  LL_desde_Elm_logit
})
## [1] -179.9817

Con la función glsm: utilizar directamente modelo$Log_Lik_Logit

Con el código de abajo se extrae directamente el valor del log-verosímil mediante modelo$Log_Lik_Logit. Este valor corresponde a \({\cal L}(\hat{\alpha})\) calculado con la fórmula del log-verosímil del modelo logístico, y es el mismo que se obtiene previamente al sumar manualmente los términos individuales a partir de modelo$Elm.

#1. Con la función lsm:

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

modelo$Log_Lik_Logit
## [1] -179.9817

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

19 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 la Tabla 19.1 se presenta un resumen de las mencionadas pruebas. En las secciones siguientes se darán los detalles correspondientes. En especial, se observa que estas estadísticas tienen distribución asintótica chi-cuadrada.

Table 19.1: Pruebas de comparación de modelos con el modelo logístico.
1 2 3 4 5 6 7 8
Pruebas de Hipótesis No. de parámetros Vector de parámetros Logaritmo de la función de verosimilitud Estadístico de prueba (\(D\)) Distribución asintótica de \(D\). Grado de libertad (\(\nu\)) P-valor
Modelo logístico \(2(1+K)\) \(\hat{\alpha}\) \({\cal L}(\hat{\alpha})\)
\(H_0\): Logit vs \(H_1\): Completo \(2n\) \(\hat{p}=u\) \({\cal L}(\hat{p}) =0\) \(2[{\cal L}(\hat{p})-{\cal L}(\hat{\alpha})]= -2\,{\cal L}(\hat{\alpha})\) \(\chi^2\) \(2[n-(1+K)]\) \(P(\chi^2_{\nu} \geq D)\)
\(H_0\): Nulo vs \(H_1\): Logit \(2\) \(\hat{p}=\hat{\delta}_o=\overline{u}\) \({\cal L}(\hat{\delta}_o)\) \(2[{\cal L}(\hat{\alpha})-{\cal L}(\hat{\delta}_o)]\) \(\chi^2\) \(2K\) \(P(\chi^2_{\nu} \geq D)\)
\(H_0\): Logit vs \(H_1\): Saturado \(2J\) \(\tilde{p}\) \({\cal L}(\tilde{p})\) \(2[{\cal L}(\tilde{p})-{\cal L}(\hat{\alpha})]\) \(\chi^2\) \(2[J-(1+K)]\) \(P(\chi^2_{\nu} \geq D)\)
\(H_0\): Subm. vs \(H_1\): Logit \(2K\) \(\hat{\alpha}_s\) \({\cal L}(\hat{\alpha}_s)\) \(2[{\cal L}(\hat{\alpha})-{\cal L}(\hat{\alpha}_s)]\) \(\chi^2\) \(2\) \(P(\chi^2_{\nu} \geq D)\)
Fuente. Elaboración propia.
Notaciones. 1 Subm.: Submodelo logístico con una variable explicativa menos. 2 Otras notaciones en las secciones siguientes

19.0.1 Logístico vs Completo

Theorem 19.1 (Comparación de un modelo logístico con el modelo completo) Para la hipótesis

\(H_0\): el modelo logístico (con \(X_1, \ldots, X_K\)),

vs la alternativa

\(H_1\): el modelo completo (que no se basa en poblaciones)

la estadística de prueba es

\[D_C:=2\ln\left(\frac{L(\hat{p})}{L(\hat{\alpha})}\right)\;=\; 2[{\cal L}(\hat{p})-{\cal L}(\hat{\alpha})] \;= \; -2\,{\cal L}(\hat{\alpha})\]

y tiene distribución asintótica chi-cuadrada con \(v=2[n-(1+K)]\) grados de libertad cuando \(n\to\infty\).

Remark. Se espera que esta prueba no rechace \(H_0\) (p-valor alto), o sea, que los datos obtenidos no estén en contra del modelo logístico. Es decir, que al pasar del modelo completo al modelo logístico no se pierde información estadísticamente significativa.

19.0.2 Nulo vs Logístico

Theorem 19.2 (Comparación de un modelo logístico con el modelo nulo) Para la hipótesis

\(H_0\): el modelo nulo (sólo con el intercepto),

vs la alternativa

\(H_1\): el modelo logístico (con \(X_1,\ldots, X_K\))

la estadística de prueba es \[D_0 \;= \; 2\ln\left(\frac{L(\hat{\alpha})}{L(\hat{\delta}_o)}\right)\; =\; 2[{\cal L}(\hat{\alpha})-{\cal L}(\hat{\delta}_o)]\]

y tiene distribución asintótica chi-cuadrada con \(2K\) grados de libertad cuando \(J\to\infty\).

Aquí \(\hat{\delta}_o=\)logit\((\overline{U})\) es la estimación de \(\delta\) en el modelo nulo.

Remark. Tenga en cuenta que:

  1. La hipótesis es equivalente a la hipótesis \(H_0: \beta=0\).

  2. Esta prueba sólo es válida para el caso de datos no agrupados (\(J=n\)).

  3. En el trabajo práctico, se espera que esta prueba sí rechace \(H_0\) (p-valor bajo). Es decir, que las variables explicativas \(X_1, \ldots, X_K\) del modelo logístico, tiene una explicación más informativa que sólo el intercepto.

  4. En caso contrario (si la prueba no rechaza \(H_0\)), que no es muy común en problemas prácticos, se tendría que chequear otro modelo logístico con más o con otras variables.

19.0.3 Logístico vs Saturado

Theorem 19.3 (Comparación de un modelo logístico con el modelo saturado) La LR-estadística de prueba (según el método de cocientes de funciones de verosimilitud) para la hipótesis

\(H_0\): el modelo logístico (con \(X_1,\ldots, X_K\)),

vs la alternativa

\(H_1\): el modelo saturado correspondiente (con sus \(J\) poblaciones)

es equivalente a la llamada desviación que tiene el modelo logístico del modelo saturado

\[D_S:=2\ln\left(\frac{L(\tilde{p})}{L(\hat{\alpha})}\right)= 2[{\cal L}(\tilde{p})-{\cal L}(\hat{\alpha})]\]

la cual tiene distribución asintótica chi-cuadrada con \(v=2[J-(1+K)]\) grados de libertad cuando \(n\to\infty\) y \(J\) es fijo.

Remark. Tenga en cuenta lo siguiente:

  1. Aquí se requiere que \(J>1+K\). Para el caso en que \(J=1+K\), el análisis en un modelo logístico es el mismo que en el modelo saturado. En la sección 14 del documento Rpubs:: Regresión logística (estimaciones) se analizaron algunas relaciones entre los modelos logístico y saturado, en especial, estos dos casos.

  2. Esta prueba únicamente se cumple para datos agrupados porque \(J\) es fijo (lo que no sucede para el caso de datos no agrupados).

  3. Se espera que esta prueba no rechace \(H_0\) (p-valor alto), o sea, que los datos obtenidos no estén en contra del modelo logístico. Es decir, que al pasar del modelo saturado al modelo logístico no se pierde información estadísticamente significativa.

19.0.4 Submodelo vs Logístico

Theorem 19.4 (Comparación de un modelo logístico con algún submodelo) Para la hipótesis

\(H_0\): un submodelo logístico con \({X}_1,\cdots, {X}_{\tilde{K}}\),

vs la alternativa

\(H_1\): el modelo logístico con \({X}_1,\cdots,{X}_K\) con \(\tilde{K}<K\),

se cumplen los siguiente resultados:

Resultado 1:

La estadística de prueba es equivalente a la estadística

\[D_{L}:=2\ln\left(\frac{L(\hat{\alpha})}{L(\hat{\alpha}_s)}\right) = 2[{\cal L}(\hat{\alpha})-{\cal L}(\hat{\alpha}_s)]\]

Aquí

\[\hat{\alpha}=(\hat{\delta},\hat{\beta}_1,\cdots, \hat{\beta}_K)^T\]

es la ML-estimación en el modelo logístico de la alternativa \(H_1\) y

\[\hat{\alpha}_s=(\hat{\delta}_s,\hat{\beta}_{s1},\cdots, \hat{\beta}_{s\tilde{K}})^T\]

es la ML-estimación en el submodelo logístico de la hipótesis \(H_0\).

Resultado 2:

\(D_L\) tiene distribución asintótica chi-cuadrada con \(2[K-\tilde{K}]\) grados de libertad cuando \(J\to\infty\).

Remark. Para considerar:

  1. Esta prueba sólo es válida para datos no agrupados. Aunque, también, es posible realizarla teniendo en cuenta el modelo saturado. Pero, como en la prueba únicamente se considera el modelo logístico, no tiene mucho sentido comparar éste con un submodelo teniendo que pasar por el modelo saturado.

  2. Se espera que no se rechace la prueba (p-valor alto).

19.0.5 Submodelo con una X menos vs Logístico

Theorem 19.5 (Comparación de un modelo logístico con un submodelo que tiene una variable explicativa menos) Para la hipótesis

\(H_0\): el submodelo (con \({X}_1,\cdots,{X}_K\) sin un \({X}_k\)),

vs la alternativa

\(H_1\): el modelo logístico (con \({X}_1,\cdots,{X}_K\))

se puede tomar, alternativamente, una de las dos estadísticas de pruebas que se mencionan en los resultados siguientes:

Resultado 1:

Estadística equivalente a:

\[D_{L}:= 2\ln\left(\frac{L(\hat{\alpha})}{L(\hat{\alpha}_s)}\right)= 2[{\cal L}(\hat{\alpha})-{\cal L}(\hat{\alpha}_s)], \]

donde \(\hat{\alpha}_s\) es la estimación bajo \(H_0\). Se resalta que \(D_L\) tiene distribución asintótica chi-cuadrada con \(2\) grado de libertad cuando \(J\to\infty\).

Resultado 2:

Estadística de Wald:

\[W= \frac{\hat{\beta}^{2}_{rk}}{\hat{V}(\hat{\beta}_{rk})} = \left(\frac{\hat{\beta}_{rk}}{SE(\hat{\beta}_{rk})}\right)^2, \]

siendo \(\hat{\beta}_{rk}\) la estimación de \(\beta_{rk}\), para cada \(k=1,\cdots,K\) y \(r=0,1,2\) en el modelo (bajo \(H_1\)) con su varianza estimada \(\hat{V}(\hat{\beta}_{rk})\) y su error estándar

\[SE(\hat{\beta}_{rk})=\sqrt{\hat{V}(\hat{\beta}_{rk})}\]

También es importante anotar que \(W\) tiene distribución asintótica chi-cuadrada con \(2\) grado de libertad cuando \(J\to\infty\).

Remark. Para tener en cuenta:

  1. En este teorema se está considerando el caso de datos no agrupados.

  2. Con base en todas las pruebas parciales, para cada \(k=0,1,\cdots,K\) se eliminará la variable explicativa que menor aporte individual tenga en la explicación. Es decir, la variable que tenga p-valor parcial más alto. Así, se sigue eliminando variable tras variable hasta que se rechacen todas las pruebas parciales (todas las tengan p-valores bajos).

20 Ejemplo 6: 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. ¿Cuántos parámetros logíticos se deben estimar?

  2. Halle las estimaciones de los parámetros logísticos.

  3. Halle el valor estimado \({\cal L}(\hat{\alpha})\).

Luego, haga una prueba de comparación entre este modelo estimado y cada uno de los siguientes modelos que se indican abajo:

  1. Completo.

  2. Nulo.

  3. Saturado.

En cada uno de los incisos anteriores escriba las hipótesis correspondientes y calcule siempre: el número de parámetros que se deben estimar, la estimación del parámetro correspondiente y del logaritmo de la función de verosimilitud y el valor del estadístico de prueba. Además, diga cuál es la distribución de ese estadístico; halle el P-valor de la prueba y escriba, obviamente, la decisión.

  1. Resuma en una tabla algunos de los resultados encontrados en los incisos anteriores.

21 Ejemplo 6: Solución

21.0.1 Solución inciso (a)

Hay \(2\) variables explicativas (ses y write). Pero, ses es categórica con tres niveles (que son ses=low, ses=middle y ses=high), con la propiedad de que

\[\hat{\beta}_{ses=low} + \hat{\beta}_{ses=middle} + \hat{\beta}_{ses=high} =0\]

Es decir, para esta variable tenemos que considerar solo dos parámetros:

\[\hat{\beta}_{ses=low}, \qquad \hat{\beta}_{ses=middle}\]

ya que

\[\hat{\beta}_{ses=high} \;=\; -\hat{\beta}_{ses=low} - \hat{\beta}_{ses=middle} \]

Por esta razón, el número de parámetros logísticos que se deben estimar es \(2(1+K)=8\), a saber, los dos interceptos \(\delta_0\), \(\delta_1\) y las seis pendientes \(\beta_{0ses=low}\), \(\beta_{0ses=middle}\), \(\beta_{1ses=low}\), \(\beta_{1ses=middle}\), \(\beta_{0write}\), \(\beta_{1write}\). El vector de parámetros es (\(T\) indica la transpuesta del vector):

\[\alpha =(\delta_0,\beta_{0ses=low}, \beta_{0ses=middle},\beta_{0write}, \delta_1,\beta_{1ses=low}, \beta_{1ses=middle}, \beta_{1write})^T\]

21.0.2 Solución inciso (b)

En los incisos (d) o (e) del Ejemplo 5 se encuentran las estimaciones correspondientes de estos parámetros.

Con la función summary :: multinom :: nnet

Puede verse el recuadro rojo de la figura 21.1.

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 21.1: Estimaciones de los parámetros logísticos con nnet. Fuente: Elaboración propia.

Con la función summary :: glsm :: glsm

Como alternativa a VGAM y nnet, se puede implementar el paquete glsm (dos opciones posibles).

datos <- glsm::hsbdemo

# Ajuste (cambia 'ref' si deseas otra categoría de referencia)
modelo <- glsm(prog ~ ses + write, data = datos, ref = "academic")

# Salida 1: Resumen tabular (coeficientes, SE, Wald, p, Exp(B)):
summary(modelo)

# Salida 2:  Coeficientes en escala logit:
coef(modelo)

Opción 1.

## 
## Call:
## glsm(formula = prog ~ ses + write, data = datos, ref = "academic")
## 
## Coefficients: 
##                        Coef(B) Std.Error    Exp(B)      Wald DF   P.value    
## (Intercept):general   1.689354  1.226940  5.415982  1.895809  1 0.1685481    
## (Intercept):vocation  4.235530  1.204650 69.098279 12.362147  1 0.0004381 ***
## seslow:general        1.162832  0.514220  3.198980  5.113707  1 0.0237376 *  
## seslow:vocation       0.982670  0.595530  2.671581  2.722757  1 0.0989270 .  
## sesmiddle:general     0.629541  0.465028  1.876749  1.832691  1 0.1758101    
## sesmiddle:vocation    1.274063  0.511065  3.575351  6.214826  1 0.0126685 *  
## write:general        -0.057928  0.021411  0.943718  7.320093  1 0.0068188 ** 
## write:vocation       -0.113603  0.022219  0.892613 26.141569  1 3.173e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Analysis of Deviance (Chi-squared): 
##                   Deviance  DF   P.value    
## Null vs Logit        48.23   6 1.063e-08 ***
## Logit vs Complete   359.96 392    0.8755    
## Logit vs Saturate   226.18 120 1.615e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Opción 2.

##  (Intercept):general (Intercept):vocation       seslow:general 
##           1.68935423           4.23552982           1.16283199 
##      seslow:vocation    sesmiddle:general   sesmiddle:vocation 
##           0.98267029           0.62954098           1.27406340 
##        write:general       write:vocation 
##          -0.05792841          -0.11360264

21.0.3 Solución inciso (c)

En el inciso (o) del Ejemplo 5 se encuentra la estimación solicitada. Se obtiene que:

\[{\cal L}(\hat{\alpha}) \;=\; -179.9817\]

library(glsm)
datos <- glsm::hsbdemo
modelo <- glsm(prog ~ ses + write, data = datos, ref = "academic")
modelo$Log_Lik_Logit


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

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

21.0.4 Solución inciso (d)

Para mayor organización, primero, presentaremos los resultados globales y, luego, los resultados obtenidos con R.

Inciso (d): Resumen global

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

  2. El número de observaciones es \(n = \mbox{200}\). Entonces, el número de parámetros en el modelo completo que se deben estimar es \[2n = \mbox{400}\]

  3. Por el teorema 6.1, tenemos que

\[\hat{p}= u\]

  1. Por el teorema 6.1, se tiene que

\[{\cal L}(\hat{p})\;=\; {\cal L}(u)\;=\; \mbox{0}\]

  1. Por el teorema 19.1, un valor del estadístico de prueba es:

\[D_C \;= \; -2\,{\cal L}(\hat{\alpha}) \;=\; -2(\mbox{-179.9817}) \;=\; \mbox{359.9635}\]

  1. Para justificar la salida obtenida con glm y que se muestra en la figura 21.2, solo tiene que observar que:

\[{\cal L}(\hat{\alpha}) \;=\; \frac{D_C}{-2} \;=\; \frac{\mbox{359.9635}}{-2} \;=\; \mbox{-179.9817}\]

  1. Teniendo en cuenta que hay \(n=200\) observaciones, el estadístico \(D_C\) tiene distribución asintótica chi-cuadrada con los siguientes grados de libertad:

\[v\;=\; 2[n- (1+K)]\;=\; \mbox{400} - \mbox{8} \;=\; \mbox{392} \]

  1. El \(P\)-valor de la prueba es

\[P\mbox{-valor} \;=\;P(\chi^2_{\mbox{392}} > \mbox{359.9635})\;=\; \mbox{0.8755}\]

  1. No se rechaza \(H_0\) al nivel del \(5\%\). Es decir, el modelo logístico es válido.

Inciso (d): con glsm::modelo

El output mostrado corresponde a la estimación del modelo logístico mediante la función glsm(), incluyendo la información necesaria para contrastar el modelo logístico con el modelo completo. En la salida se observa la fórmula ajustada, el número de poblaciones en el modelo saturado, y los coeficientes estimados junto a sus errores estándar y odds ratios. En la sección de log-verosimilitud, se presentan los valores correspondientes al modelo completo, nulo, logístico y saturado. En particular, el valor de \({\cal L}(y)\) para el modelo completo es esencial para calcular el estadístico, sus grados de libertad y el valor-p, lo que permite evaluar si el modelo logístico explica adecuadamente los datos en comparación con el modelo completo.

library(glsm)
datos <- glsm::hsbdemo
modelo <- glsm(prog ~ ses + write, data = datos, ref = "academic")
modelo
## 
## Call:
## glsm(formula = prog ~ ses + write, data = datos, ref = "academic")
## 
## Populations in Saturated Model: 64
## 
## Coefficients: 
##                          Coef(B)  Std.Error     Exp(B)
## (Intercept):general   1.68935423 1.22694022  5.4159821
## (Intercept):vocation  4.23552982 1.20464972 69.0982789
## seslow:general        1.16283199 0.51422012  3.1989800
## seslow:vocation       0.98267029 0.59552969  2.6715806
## sesmiddle:general     0.62954098 0.46502835  1.8767489
## sesmiddle:vocation    1.27406340 0.51106550  3.5753512
## write:general        -0.05792841 0.02141082  0.9437175
## write:vocation       -0.11360264 0.02221890  0.8926126
## 
## Log Likelihood: 
##          Estimation
## Complete    0.00000
## Null     -204.09667
## Logit    -179.98173
## Saturate  -66.89303

Inciso (d): con glsm::modelo$...

En este bloque se extraen directamente, a partir del objeto modelo, los valores necesarios para contrastar el modelo logístico frente al modelo completo: número de parámetros, log-verosimilitud, estadístico de prueba, grados de libertad y valor-p. Estos resultados permiten evaluar si el modelo logístico presenta un ajuste significativamente diferente al modelo completo.

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

# Punto 2 (número de parámetros):
n <- nrow(chdage)
n

# Punto 4 (Log-verosimilitud):
Log_Lik_Complete <- modelo$Log_Lik_Complete  
Log_Lik_Complete

# Punto 5 (Estadístico de prueba):
Logit_vs_Complete <- modelo$Dev_Logit_vs_Complete 
Logit_vs_Complete

# Punto 7 (grados de libertad):
df_c <- modelo$Df_Logit_vs_Complete
df_c

# Punto 8 (P-valor):
pvalor_c <- modelo$P.v_Logit_vs_Complete
pvalor_c
## Número de parámetros (modelo completo): 100
## Log-verosimilitud (modelo completo): 0
## Estadístico de prueba (modelo completo): 359.9635
## Grados de libertad (modelo completo): 392
## P-valor (modelo completo): 0.8755302

Inciso (d): con glsm::summary

La prueba de comparación se puede llevar a cabo con la función glsm::summary(). Con ella obtenemos una matriz llamada Analysis of Deviance (Chi-squared) en donde se encuentran: el nombre de la prueba de comparación, el valor del estadístico de prueba (Deviance), los grados de libertad (DF) y el correspondiente P-valor de la prueba (P.value).

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

summary(modelo)
## 
## Call:
## glsm(formula = prog ~ ses + write, data = datos, ref = "academic")
## 
## Coefficients: 
##                        Coef(B) Std.Error    Exp(B)      Wald DF   P.value    
## (Intercept):general   1.689354  1.226940  5.415982  1.895809  1 0.1685481    
## (Intercept):vocation  4.235530  1.204650 69.098279 12.362147  1 0.0004381 ***
## seslow:general        1.162832  0.514220  3.198980  5.113707  1 0.0237376 *  
## seslow:vocation       0.982670  0.595530  2.671581  2.722757  1 0.0989270 .  
## sesmiddle:general     0.629541  0.465028  1.876749  1.832691  1 0.1758101    
## sesmiddle:vocation    1.274063  0.511065  3.575351  6.214826  1 0.0126685 *  
## write:general        -0.057928  0.021411  0.943718  7.320093  1 0.0068188 ** 
## write:vocation       -0.113603  0.022219  0.892613 26.141569  1 3.173e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Analysis of Deviance (Chi-squared): 
##                   Deviance  DF   P.value    
## Null vs Logit        48.23   6 1.063e-08 ***
## Logit vs Complete   359.96 392    0.8755    
## Logit vs Saturate   226.18 120 1.615e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Inciso (d): con glsm::summary$comparison_test

En este apartado se emplea la salida de summary(modelo) para acceder de forma directa a la tabla comparison_test, desde la cual se obtienen el estadístico de prueba, los grados de libertad y el valor-p correspondientes a la comparación entre el modelo logístico y el modelo completo. Esta forma de extracción facilita el cálculo y presentación de los resultados sin necesidad de acceder manualmente a cada elemento del objeto modelo.

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

LvsC <- summary(modelo)

# Punto 5 (Estadístico de prueba):
Logit_vs_Complete <- LvsC$`comparison test`[2,1]
Logit_vs_Complete

# Punto 7 (grados de libertad):
df_c <- LvsC$`comparison test`[2,2]
df_c

# Punto 8 (P-valor):
pvalor_c <- LvsC$`comparison test`[2,3]
pvalor_c
## Estadístico de prueba (modelo completo): 359.9635
## Grados de libertad (modelo completo): 392
## P-valor (modelo completo): 0.8755302

Inciso (d): con vglm::VGAM

Los valores obtenidos se pueden visualizar en el recuadro rojo en la figura 21.2.

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

K <- 3
n_alfa <- 2*(1+K)

# Punto 2 (número de parámetros):
c <- 2*nrow(datos)

# Punto 4 (Log-verosimilitud):
Lp_i <-  ifelse(prog0 == 0 | prog1 == 0 | prog0+prog1 == 1, 0, 
                prog0*log(prog0)+prog1*log(prog1)+(1-prog0-prog1)*log(1-prog0-prog1)) 

Log_Lik_Complete <- sum(Lp_i)

# Punto 5 (Estadístico de prueba):
Logit_vs_Complete <- deviance(modelo)  

# Punto 7 (grados de libertad)
df_c <- df.residual(modelo)

# Punto 8 (P-valor):
pvalor_c <- pchisq(deviance(modelo),df=df_c, lower.tail=FALSE)
Prueba de comparación del modelo logístico con el completo. Fuente: Elaboración propia.

Figure 21.2: Prueba de comparación del modelo logístico con el completo. Fuente: Elaboración propia.

21.0.5 Solución inciso (e)

Para mayor organización, primero, presentaremos los resultados globales y, luego, los resultados obtenidos con R.

Inciso (e): Resumen global

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

  2. El número de parámetros en el modelo nulo que se deben estimar es \[\mbox{2}\]

  3. Por la ecuación (7.1) , se tiene que

\[{\cal L}(\hat{p})\;=\; {\cal L}(\overline{u})\;=\; \mbox{-204.0967}\]

  1. Por el teorema 19.2, un valor del estadístico de prueba es: \[D_0 \;= \; 2[{\cal L}(\hat{\alpha})-{\cal L}(\hat{\delta})] \;=\;2[\mbox{-179.9817} - (\mbox{-204.0967})] \;=\; \mbox{48.2299}\]

  2. Teniendo en cuenta que el modelo nulo tiene \(n= \mbox{2}\) parámetros, el estadístico \(D_0\) tiene distribución asintótica chi-cuadrada con los siguientes grados de libertad:

\[v\;=\; 2(1+K)-2\;=\; \mbox{8} - \mbox{2} \;=\; \mbox{6} \]

  1. El \(P\)-valor de la prueba es

\[P\mbox{-valor} \;=\;P(\chi^2_{\mbox{6}} > \mbox{48.2299})\;=\; \mbox{0}\]

  1. No se rechaza \(H_0\) al nivel del \(5\%\). Es decir, el modelo logístico no es válido.

Inciso (e): con glsm::modelo

El output presentado contiene también la información necesaria para el contraste entre el modelo logístico y el modelo nulo. En la salida, el valor de la log-verosimilitud del modelo nulo \({\cal L}(\overline{y})\) se utiliza como referencia para calcular el estadístico de prueba que compara este modelo con el logístico. Junto con los grados de libertad y el valor-p obtenidos, esta comparación permite evaluar si la inclusión de la variable explicativa mejora significativamente el ajuste del modelo frente a la hipótesis de probabilidad constante.

library(glsm)
datos <- glsm::hsbdemo
modelo <- glsm(prog ~ ses + write, data = datos, ref = "academic")
modelo
## 
## Call:
## glsm(formula = prog ~ ses + write, data = datos, ref = "academic")
## 
## Populations in Saturated Model: 64
## 
## Coefficients: 
##                          Coef(B)  Std.Error     Exp(B)
## (Intercept):general   1.68935423 1.22694022  5.4159821
## (Intercept):vocation  4.23552982 1.20464972 69.0982789
## seslow:general        1.16283199 0.51422012  3.1989800
## seslow:vocation       0.98267029 0.59552969  2.6715806
## sesmiddle:general     0.62954098 0.46502835  1.8767489
## sesmiddle:vocation    1.27406340 0.51106550  3.5753512
## write:general        -0.05792841 0.02141082  0.9437175
## write:vocation       -0.11360264 0.02221890  0.8926126
## 
## Log Likelihood: 
##          Estimation
## Complete    0.00000
## Null     -204.09667
## Logit    -179.98173
## Saturate  -66.89303

Inciso (e): con glsm::modelo$...

En este bloque se extraen directamente, a partir del objeto modelo, los valores necesarios para contrastar el modelo logístico frente al modelo completo: número de parámetros, log-verosimilitud, estadístico de prueba, grados de libertad y valor-p. Estos resultados permiten evaluar si el modelo logístico presenta un ajuste significativamente diferente al modelo completo.

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

# Punto 2 (número de parámetros):
par_nulo <- 2
par_nulo

# Punto 3 (Log-verosimilitud):
Log_Lik_Nulo <- modelo$Log_Lik_Null  
Log_Lik_Nulo

# Punto 4 (Estadístico de prueba):
Logit_vs_Nulo <- modelo$Dev_Null_vs_Logit
Logit_vs_Nulo

# Punto 5 (grados de libertad):
df_0 <- modelo$Df_Null_vs_Logit
df_0

# Punto 6 (P-valor):
pvalor_0 <- modelo$P.v_Null_vs_Logit
pvalor_0
## Número de parámetros (modelo nulo): 2
## Log-verosimilitud (modelo nulo): -204.0967
## Estadístico de prueba (modelo nulo): 48.2299
## Grados de libertad (modelo nulo): 6
## P-valor (modelo nulo): 1.063001e-08
## NULL

Inciso (e): con glsm::summary

La prueba de comparación se puede llevar a cabo con la función glsm::summary(). Con ella obtenemos una matriz llamada Analysis of Deviance (Chi-squared) en donde se encuentran: el nombre de la prueba de comparación, el valor del estadístico de prueba (Deviance), los grados de libertad (DF) y el correspondiente P-valor de la prueba (P.value).

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

summary(modelo)
## 
## Call:
## glsm(formula = prog ~ ses + write, data = datos, ref = "academic")
## 
## Coefficients: 
##                        Coef(B) Std.Error    Exp(B)      Wald DF   P.value    
## (Intercept):general   1.689354  1.226940  5.415982  1.895809  1 0.1685481    
## (Intercept):vocation  4.235530  1.204650 69.098279 12.362147  1 0.0004381 ***
## seslow:general        1.162832  0.514220  3.198980  5.113707  1 0.0237376 *  
## seslow:vocation       0.982670  0.595530  2.671581  2.722757  1 0.0989270 .  
## sesmiddle:general     0.629541  0.465028  1.876749  1.832691  1 0.1758101    
## sesmiddle:vocation    1.274063  0.511065  3.575351  6.214826  1 0.0126685 *  
## write:general        -0.057928  0.021411  0.943718  7.320093  1 0.0068188 ** 
## write:vocation       -0.113603  0.022219  0.892613 26.141569  1 3.173e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Analysis of Deviance (Chi-squared): 
##                   Deviance  DF   P.value    
## Null vs Logit        48.23   6 1.063e-08 ***
## Logit vs Complete   359.96 392    0.8755    
## Logit vs Saturate   226.18 120 1.615e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Inciso (e): con glsm::summary$comparison_test

En este apartado se emplea la salida de summary(modelo) para acceder de forma directa a la tabla comparison_test, desde la cual se obtienen el estadístico de prueba, los grados de libertad y el valor-p correspondientes a la comparación entre el modelo logístico y el modelo nulo. Esta forma de extracción facilita el cálculo y presentación de los resultados sin necesidad de acceder manualmente a cada elemento del objeto modelo.

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

LvsC <- summary(modelo)

# Punto 5 (Estadístico de prueba):
Logit_vs_Nulo <- LvsC$`comparison test`[1,1]
Logit_vs_Nulo

# Punto 7 (grados de libertad):
df_0 <- LvsC$`comparison test`[1,2]
df_0

# Punto 8 (P-valor):
pvalor_0 <- LvsC$`comparison test`[1,3]
pvalor_0
## Estadístico de prueba (modelo nulo): 48.2299
## Grados de libertad (modelo nulo): 6
## P-valor (modelo nulo): 1.063001e-08

Inciso (e): con vglm::VGAM

Los valores obtenidos se pueden visualizar en el recuadro rojo en la figura 21.3):

library(VGAM)
modelo <- vglm(prog ~ ses + write, multinomial(refLevel = 1), data=datos)
Log_Lik_Logit <- deviance(modelo)/(-2)

modelo <- vglm(prog ~ 1, multinomial(refLevel = 1), data=datos)
summary(modelo)

# Punto 2 (número de parámetros):
n <- 2

# Punto 3 (Log-verosimilitud):
Log_Lik_Nulo <- deviance(modelo)/(-2)  

# Punto 4 (Estadístico de prueba):
Logit_vs_Nulo <- 2*(Log_Lik_Logit - Log_Lik_Nulo)

# Punto 5 (grados de libertad)
df_0 <- df.residual(modelo)

# Punto 6 (P-valor):
pvalor_0 <- pchisq(Logit_vs_Nulo, df=df_0, lower.tail=FALSE)
Prueba de comparación del modelo logístico con el nulo. Fuente: Elaboración propia.

Figure 21.3: Prueba de comparación del modelo logístico con el nulo. Fuente: Elaboración propia.

21.0.6 Solución inciso (f)

Para mayor organización, primero, presentaremos los resultados globales y, luego, los resultados obtenidos con R.

Inciso (f): Resumen global

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

  2. El número de poblaciones es \(J = \mbox{64}\). Entonces, el número de parámetros en el modelo saturado que se deben estimar es

\[2J = \mbox{128}\]

  1. Por el Teorema 10.2 , se tiene que

\[{\cal L}(\tilde{p})\;=\; \mbox{-66.8931}\]

  1. Por el teorema 19.3, un valor del estadístico de prueba es:

\[D_S \;= \; 2[{\cal L}(\tilde{p})-{\cal L}(\hat{\alpha})] \;=\;2[\mbox{-66.8931} - (\mbox{-179.9817})] \;=\; \mbox{226.1773}\]

  1. Teniendo en cuenta que el modelo nulo tiene \(n=\) 64 parámetros, el estadístico \(D_S\) tiene distribución asintótica chi-cuadrada con los siguientes grados de libertad:

\[v\;=\; 2[J- (1+K)]\;=\; \mbox{128} - \mbox{8} \;=\; \mbox{120}\]

  1. El \(P\)-valor de la prueba es

\[P\mbox{-valor} \;=\;P(\chi^2_{\mbox{120}} > \mbox{226.1773})\;=\; \mbox{0}\]

  1. Se rechaza \(H_0\) al nivel del \(5\%\). Es decir, el modelo logístico no es válido.

Inciso (f): con glsm::modelo

En este apartado se presenta el output generado por la función glsm() para el ajuste del modelo logístico, el cual resume información esencial para el contraste entre el modelo logístico y el modelo saturado. En la salida se incluye la fórmula del modelo ajustado, el número de poblaciones en el modelo saturado, los coeficientes estimados junto con sus errores estándar y odds ratios, así como los valores de la log-verosimilitud bajo diferentes configuraciones (modelo completo, nulo, logit y saturado). Estos resultados son la base para calcular los estadísticos de prueba, los grados de libertad y el valor-p que permiten evaluar la validez del modelo.

library(glsm)
datos <- glsm::hsbdemo
modelo <- glsm(prog ~ ses + write, data = datos, ref = "academic")
modelo
## 
## Call:
## glsm(formula = prog ~ ses + write, data = datos, ref = "academic")
## 
## Populations in Saturated Model: 64
## 
## Coefficients: 
##                          Coef(B)  Std.Error     Exp(B)
## (Intercept):general   1.68935423 1.22694022  5.4159821
## (Intercept):vocation  4.23552982 1.20464972 69.0982789
## seslow:general        1.16283199 0.51422012  3.1989800
## seslow:vocation       0.98267029 0.59552969  2.6715806
## sesmiddle:general     0.62954098 0.46502835  1.8767489
## sesmiddle:vocation    1.27406340 0.51106550  3.5753512
## write:general        -0.05792841 0.02141082  0.9437175
## write:vocation       -0.11360264 0.02221890  0.8926126
## 
## Log Likelihood: 
##          Estimation
## Complete    0.00000
## Null     -204.09667
## Logit    -179.98173
## Saturate  -66.89303

Inciso (f): con glsm::modelo$...

En este bloque se accede directamente a los elementos del objeto modelo para obtener la información necesaria en la comparación entre el modelo logístico y el modelo saturado. Se extraen el número de parámetros, la log-verosimilitud del modelo saturado, el estadístico de prueba, los grados de libertad y el valor-p, los cuales permiten evaluar si el modelo logístico presenta un ajuste estadísticamente equivalente al modelo saturado.

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

# Punto 2 (número de parámetros):
J <- modelo$Populations
2*J

# Punto 3 (Log-verosimilitud):
Log_Lik_Sat <- modelo$Log_Lik_Saturate  
Log_Lik_Sat

# Punto 4 (Estadístico de prueba):
Logit_vs_Sat <- modelo$Dev_Logit_vs_Saturate
Logit_vs_Sat

# Punto 5 (grados de libertad):
df_s <- modelo$Df_Logit_vs_Saturate
df_s

# Punto 6 (P-valor):
pvalor_s <- modelo$P.v_Logit_vs_Saturate
pvalor_s
## Número de parámetros (modelo saturado): 128
## Log-verosimilitud (modelo saturado): -66.89303
## Estadístico de prueba (modelo saturado): 226.1774
## Grados de libertad (modelo saturado): 120
## P-valor (modelo saturado): 1.614849e-08

Inciso (f): con glsm::summary

La prueba de comparación se puede llevar a cabo con la función glsm::summary(). Con ella obtenemos una matriz llamada Analysis of Deviance (Chi-squared) en donde se encuentran: el nombre de la prueba de comparación, el valor del estadístico de prueba (Deviance), los grados de libertad (DF) y el correspondiente P-valor de la prueba (P.value).

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

summary(modelo)
## 
## Call:
## glsm(formula = prog ~ ses + write, data = datos, ref = "academic")
## 
## Coefficients: 
##                        Coef(B) Std.Error    Exp(B)      Wald DF   P.value    
## (Intercept):general   1.689354  1.226940  5.415982  1.895809  1 0.1685481    
## (Intercept):vocation  4.235530  1.204650 69.098279 12.362147  1 0.0004381 ***
## seslow:general        1.162832  0.514220  3.198980  5.113707  1 0.0237376 *  
## seslow:vocation       0.982670  0.595530  2.671581  2.722757  1 0.0989270 .  
## sesmiddle:general     0.629541  0.465028  1.876749  1.832691  1 0.1758101    
## sesmiddle:vocation    1.274063  0.511065  3.575351  6.214826  1 0.0126685 *  
## write:general        -0.057928  0.021411  0.943718  7.320093  1 0.0068188 ** 
## write:vocation       -0.113603  0.022219  0.892613 26.141569  1 3.173e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Analysis of Deviance (Chi-squared): 
##                   Deviance  DF   P.value    
## Null vs Logit        48.23   6 1.063e-08 ***
## Logit vs Complete   359.96 392    0.8755    
## Logit vs Saturate   226.18 120 1.615e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Inciso (f): con glsm::summary$comparison_test

En este apartado se emplea la salida de summary(modelo) para acceder de forma directa a la tabla comparison_test, desde la cual se obtienen el estadístico de prueba, los grados de libertad y el valor-p correspondientes a la comparación entre el modelo logístico y el modelo saturado. Esta forma de extracción facilita el cálculo y presentación de los resultados sin necesidad de acceder manualmente a cada elemento del objeto modelo.

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

LvsC <- summary(modelo)

# Punto 5 (Estadístico de prueba):
Logit_vs_Sat <- LvsC$`comparison test`[3,1]
Logit_vs_Sat

# Punto 7 (grados de libertad):
df_s <- LvsC$`comparison test`[3,2]
df_s

# Punto 8 (P-valor):
pvalor_s <- LvsC$`comparison test`[3,3]
pvalor_s
## Estadístico de prueba (modelo saturado): 226.1774
## Grados de libertad (modelo saturado): 120
## P-valor (modelo saturado): 1.614849e-08

Inciso (f): con vglm::VGAM

Los valores obtenidos se pueden hallar así:

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

Log_Lik_Logit <- deviance(modelo)/(-2)

datos %>%
  group_by(ses, write) %>%
  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

# Punto 2 (número de poblaciones J y de parámetros):
J <- nrow(saturado)
N_parametros <- 2*J

# Punto 3 (Log-verosimilitud):
Log_Lik_Sat <- sum(saturado$Lp_ref2) 

# Punto 4 (Estadístico de prueba):
Logit_vs_Sat <- 2*(Log_Lik_Sat - Log_Lik_Logit)

# Punto 5 (grados de libertad)
df_s <- 2*J - n_alfa

# Punto 6 (P-valor):
pvalor_s <- pchisq(Logit_vs_Sat,df=df_s, lower.tail=FALSE)

21.0.7 Solución inciso (g)

En la Tabla 21.1 se presenta un resumen de los resultados encontrados en los ejercicios anteriores. Se ha tomado como base la tabla 19.1.

Table 21.1: Pruebas de comparación de modelos con el modelo logístico para los datos chdage.
1 2 3 4 5 6 7 8
Pruebas de Hipótesis Vector de parámetros No. de parámetros Logaritmo de la función de verosimilitud Estadístico de prueba (\(D\)) Distribución asintótica de \(D\). Grado de libertad (\(\nu\)) P-valor \(P(\chi^2_{\nu} \geq D)\)
Modelo logístico \(\hat{\alpha}\) \(2(1+K)=8\) \({\cal L}(\hat{\alpha})=-179.9817\)
\(H_0\): Logit vs \(H_1\): Completo \(\hat{p}=u\) \(2n=400\) \({\cal L}(\hat{p}) =0\) \(-2\,{\cal L}(\hat{\alpha})=359.9635\) \(\chi^2\) \(2[n-(1+K)]=392\) \(0.8755\)
\(H_0\): Nulo vs \(H_1\): Logit \(\hat{p}=\hat{\delta}\) \(2\) \({\cal L}(\hat{\delta})=-204.0967\) \(2[{\cal L}(\hat{\alpha})-{\cal L}(\hat{\delta})]=48.23\) \(\chi^2\) \(2(1+K)-2=6\) \(0\)
\(H_0\): Logit vs \(H_1\): Saturado \(\tilde{p}\) \(2J=128\) \({\cal L}(\tilde{p})=-66.8931\) \(2[{\cal L}(\tilde{p})-{\cal L}(\hat{\alpha})]=226.1773\) \(\chi^2\) \(2[J-(1+K)]=120\) \(0\)
\(H_0\): Subm. vs \(H_1\): Logit \(\hat{\alpha}_s\) \(2K=6\) \({\cal L}(\hat{\alpha}_s)\) \(2[{\cal L}(\hat{\alpha})-{\cal L}(\hat{\alpha}_s)]\) \(\chi^2\) \(2\)
Fuente. Elaboración propia.
Submodelo. 1 Subm.: Submodelo logístico con una variable explicativa menos; 2 En este ejemplo no hay submodelo.
Notaciones. a \(T:\) transpuesta de un vector

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) 10.1; (d) 10.2; (e) 12.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.