hllinas2023

1 Paquetes

1.0.1 Paquetes utilizados

En estas notas de clase utilizaremos los siguientes paquetes:

  • aplore3: paquete de Braglia (2016) y contiene bases de datos clásicas empleadas en el análisis de regresión logística, ideales para ilustrar conceptos teóricos con ejemplos reales.

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

  • lsm: paquete de autoría propia (Villalba JL, Llinás HJ y Fabregas OJ (2025)) que permite calcular de forma sencilla y precisa el log-verosímil del modelo saturado cuando la variable respuesta Y es dicotómica (0 o 1). Además, incluye conjuntos de datos de ejemplo útiles para ilustrar su aplicación.

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

1.0.2 Sobre el paquete lsm

  • El paquete lsm, desarrollado por Villalba JL, Llinás HJ y Fabregas OJ (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 dicotómica (es decir, toma valores de 0 o 1), la función lsm() 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 lsm no solo incluye funciones especializadas, sino también conjuntos de datos integrados que permiten realizar pruebas y ejemplos de manera práctica: chdage, icu, lowbwt, pros, survey y uis.

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

2 Introducción

Los métodos de regresión se han convertido en una herramienta fundamental en el análisis de datos, especialmente cuando se busca describir la relación entre una variable de respuesta y una o más variables explicativas. Con frecuencia, la variable de resultado es discreta, tomando uno de dos valores posibles. En estos casos, el modelo de regresión más utilizado es el modelo de regresión logística.

En el documento RPubs :: Modelos lineales generalizados se explicó que este tipo de modelos forma parte de la familia de los modelos lineales generalizados (GLM). Allí se abordan sus fundamentos y características principales. Para comprender en profundidad el funcionamiento y la lógica de la regresión logística, es importante estudiar los siguientes cuatro modelos clave:

  • Modelo de Bernoulli.

  • Modelo completo.

  • Modelo nulo.

  • Modelo saturado.

A continuación, se describirán sus propiedades y características principales, junto con ejemplos ilustrativos para facilitar su comprensión.

3 Datasets

Para las aplicaciones, se utilizaron bases de datos de las librerías aplore3 (creado por Braglia (2016)) y lsm (desarrollado por Villalba JL, Llinás HJ y Fabregas OJ (2025)).

library(aplore3) 
library(lsm)

Ambos paquetes incluyen, de manera no oficial, todos los conjuntos de datos utilizados en el texto de Hosmer, Lemeshow and Sturdivant (2013). En este link o en este otro se encuentran los nombres de los datasets con los respectivos detalles:

  • Descripción (description).

  • Uso (usage).

  • Formato (format).

  • Fuente (source) o Referencias (References).

Se resalta el hecho que lsm es un paquete que contiene:

  1. Otros datasets que son de autoría propia.

  2. La función lsm(), la cual nos permite estimar, entre otros, el logaritmo de la función de verosimilitud de los modelos completos, nulo, saturado y logístico. La teoría relacionada con los tres primeros modelos se van a explicar más adelante dentro de este documento. La que se relaciona con los modelos logísticos, en otros documentos.

4 Ejemplo introductorio

4.0.1 Descripción del datasets

Para los ejemplos, se utilizará chdage, el cual contiene datos recogidos con el fin de estudiar si la edad es un factor influyente en la presencia o no de enfermedades coronarias (CHD). Es un data frame con 100 observaciones. En aplore3 hay cuatro variables y en lsm, tres. A continuación, se describen cada una (entre paréntesis, el nombre de la variable en lsm):

  1. id (ID): código de identificación (1-100).

  2. chd (CHD): Presencia of CHD. En aplore3, los niveles son No, Si. En lsm son 0 (=No) y 1 (=Si).

  3. age (AGE): edad (en años) de los participantes.

  4. agegrp: edad agrupada en 8 niveles (20-39, 30-34, 35-39, 40-44, 45-49, 50-54, 55-59, 60-69).

4.0.2 Datasets desde lsm

Cargamos el data frame chdage de la librería lsm:

datos <- lsm::chdage
attach(datos)

head(datos, 10)

Las primeras 10 observaciones son:

ID AGE CHD
1 20 0
2 23 0
3 24 0
4 25 0
5 25 1
6 26 0
7 26 0
8 28 0
9 28 0
10 29 0

4.0.3 Datasets desde aplore3

Cargamos el data frame chdage de la librería aplore3:

df <- aplore3::chdage
attach(df)

head(df, 10)

Las primeras 10 observaciones son:

ID age agegrp chd
1 20 20-39 No
2 23 20-39 No
3 24 20-39 No
4 25 20-39 No
5 25 20-39 Yes
6 26 20-39 No
7 26 20-39 No
8 28 20-39 No
9 28 20-39 No
10 29 20-39 No

Observe que, en el datasets aplore3::chdage, la variable chd es binaria, pero no numérica. Al convertirla en integer, se codifica como 1=No y 2=Si. Por esta razón, esta variable se puede codificar como 0=No y 1=Si, de la siguiente manera:

df <- aplore3::chdage
df$chd_num <- as.integer(df$chd)-1

4.0.4 Ejemplo: enunciado

Example 4.1 Explore la relación entre las variables age, agegrp y chd siguiendo las instrucciones que se proponen abajo.

  1. Escriba un resumen de los datos.

  2. Realize dos diagramas: Uno de dispersión de age versus chd y otro que corresponda a uno de caja y bigotes para age dentro de cada valor de chd (incluya el diagrma de dispersión correspondiente). Responda las preguntas:

  1. ¿Tendencia para los individuos con evidencia de chd?
  2. ¿Describe claramente este diagrama la naturaleza dicotómica de chd?
  3. ¿Provee este diagrama una clara relación entre las dos variables?
  4. ¿Es pequeña la variabilidad en chd en todas las edades?
  5. ¿Qué es lo que dificulta describir la relación entre las dos variables?
  6. ¿Qué posible método podemos utilizar para remover alguna variación pero sin modificar la posible relación entre las dos variables?

  1. Construya una tabla de frecuencias paraagegrp.

  2. Realize una tabulación cruzada entre agegrp y chd. Por lo menos, para cada grupo etáreo, deben aparecer las cantidades de individuos con presencia y no presencia de enfermedades coronarias.

  3. Para cada grupo de edad, calcule la proporción de personas con enfermedades coronarias. Es decir, construir una tabla de frecuencias relativas para agegrp. ¿Es esta proporción creciente o decreciente a medida que la edad aumenta?

  4. Realize un diagrama de dispersión de los grupos etáreos (en el eje X) versus estas proporciones.

  1. ¿Provee este diagrama una clara relación entre las dos variables?
  2. ¿Cuál podría ser esta relación?

4.0.5 Ejemplo: solución

Inciso 1.

El resumen de los datos se puede encontrar así:

summary(chdage)
ID AGE CHD
Min. : 1.00 Min. :20.00 Min. :0.00
1st Qu.: 25.75 1st Qu.:34.75 1st Qu.:0.00
Median : 50.50 Median :44.00 Median :0.00
Mean : 50.50 Mean :44.38 Mean :0.43
3rd Qu.: 75.25 3rd Qu.:55.00 3rd Qu.:1.00
Max. :100.00 Max. :69.00 Max. :1.00

Inciso 2.

El diagrama de dispersión de age versus chd se puede visualizar en las Figuras 4.1 o 4.2.

datos <- lsm::chdage
attach(datos)
modelo <- lsm(CHD~AGE, data=datos)

plot(modelo, 
     type = "scatter",  
     color= "blue", 
     size =2, 
     title = "Diagrama de dispersión",
     xlab = "Age", ylab = "CHD")
Diagrama de dispersión para los datos chdage

Figure 4.1: Diagrama de dispersión para los datos chdage 

ggplot(chdage, aes(y = chd, x = age))+
  
geom_point(aes(color=chd), alpha = 1.9) +
  
labs(x="Edad", y="Presencia o no de CHD", fill= "")+ 
  
#guides(fill=FALSE)+                              #Etiquetas

facet_wrap(. ~ "Diagrama de dispersión") +    #Título
theme(legend.position = "none")  +
theme_bw(base_size = 12)
Diagrama de dispersión para los datos chdage

Figure 4.2: Diagrama de dispersión para los datos chdage

El diagrama de dispersión con el de cajas y bigotes incluido es:

ggplot(chdage, aes(y = chd, x = age))+
  
geom_boxplot(aes(fill= chd),color="black",size = .75)  + 
geom_jitter(alpha = .9, aes(color=chd)) +
  
labs(x="Edad", y="Presencia o no de CHD", fill= "")+ 
  
#guides(fill=FALSE)+                              #Etiquetas
scale_fill_manual(values = c("wheat", "yellow"))+ #Cambiar los colores

facet_wrap(. ~ "Diagramas de dispersión y de caja y bigotes ") +    #Título
theme() +
theme_bw(base_size = 12)  + 
theme(axis.text.x = element_text(angle = 0, hjust = 1, vjust = 1))

Inciso 3.

La tabla de frecuencias para agegrp es:

Tabla <- table(agegrp)
Tabla <- with(chdage, addmargins(Tabla))
Grupo etáreo Frecuencia
20-39 10
30-34 15
35-39 12
40-44 15
45-49 13
50-54 8
55-59 17
60-69 10
Sum 100

Inciso 4.

Esta es la tabla de frecuencias cruzada entre agegrp y chd:

with(chdage, addmargins(table(agegrp, chd)))
Grupo etáreo No Si Total
20-39 9 1 10
30-34 13 2 15
35-39 9 3 12
40-44 10 5 15
45-49 7 6 13
50-54 3 5 8
55-59 4 13 17
60-69 2 8 10
Total 57 43 100

Inciso 5.

La tabla de frecuencias relativas para agegrp es:

with(chdage, tapply(CHD, list(agegrp), mean))
Grupo etáreo Frecuencia relativa
20-39 0.1
30-34 0.1333
35-39 0.25
40-44 0.3333
45-49 0.4615
50-54 0.625
55-59 0.7647
60-69 0.8

Inciso 6.

El diagrama de dispersión de los grupos etáreos (en el eje X) versus las proporciones halladas en el punto anterior es:

Valor <- with(chdage, tapply(CHD, list(agegrp), mean))
Prueba <- as.data.frame(Valor)
Edad <- c("20-39", "30-34" ,"35-39", "40-44", "45-49", "50-54", "55-59", "60-69")
Tabla <- cbind(Edad, round(Prueba$Valor,4))

Tabla <- as.data.frame(Tabla)

ggplot(Tabla,aes(x = factor(Edad), y=as.numeric(V2),fill=factor(Edad))) +
geom_bar(width = 0.9,stat="identity", position = position_dodge()) +
geom_point(aes(fill=Edad), alpha = 1.9, size=4, color="black")+
labs(x="Edad", y="Propoción de personas con CHD", fill= "")+ 
ylim(0,1)+
facet_wrap(. ~ "Diagrama de dispersión") +    
theme_bw(base_size = 12) +
theme(legend.position = "none")

5 Modelo de Bernoulli

5.0.1 El modelo

La variable de interés \(Y\) es de Bernoulli. En símbolo, \(Y \sim {\cal B}(1,p)\), en donde la probabilidad \(p\) de que ocurra \(Y\) es:

\[p:=E(Y)=P(Y=1)\]

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

Se llega a un modelo estadístico de Bernoulli, en donde cada variable muestra \(Y_i\) también es de Bernoulli:

\[Y_i \;\sim\; {\cal B} (1,p_i), \quad i=1,\cdots,n\]

5.0.2 La función de verosimilitud

Fijando \(y=(y_1, \cdots ,y_n)^T\) obtenemos la función de verosimilitud en el parámetro \(p=(p_1, \cdots, p_n)^T\):

\[L(p)= \prod^n_{i=1}[p_i^{y_i}(1-p_i)^{1-y_i}]\]

Theorem 5.1 (Log-verosimilitud) En un modelo de Bernoulli, el logaritmo de la función de máxima verosimilitud será:

\[\begin{equation} {\cal L}(p):= \ln L(p)= \sum^n_{i=1}[y_i\ln p_i + (1-y_i)\ln (1-p_i)] \tag{5.1} \end{equation}\]

5.0.3 Propiedad de la función de verosimilitud

Remark. Dado que \(0 \le f(y, p) \le 1\), se deduce de la expresión (5.1) que la función de log-verosimilitud está acotada superiormente por cero:

\[-\infty \le {\cal L}(p) \le 0\]

Esta propiedad se puede visualizar en el siguiente gráfico, que muestra cómo la función de log-verosimilitud para un modelo Bernoulli depende de la proporción de éxitos observados. Cada curva representa un escenario con diferente número de éxitos (\(n_1\)) y fracasos (\(n_0\)). En todos los casos, el valor máximo de la función se alcanza cuando:

\[ \hat{p} = \frac{n_1}{n_1 + n_0} \]

Es decir, en la proporción observada de éxitos. Como es natural, este valor máximo nunca supera cero, lo que refleja la propiedad:

\[ \mathcal{L}(p) \le 0 \]

Esta visualización resulta útil para comprender por qué los modelos logísticos se construyen sobre la base de maximizar una función log-verosímil que, por definición, siempre es no positiva.

# Cargar paquetes necesarios
library(ggplot2)
library(dplyr)

# Crear diferentes escenarios con distintos n1 y n0
escenarios <- list(
  c(n1 = 2, n0 = 8),
  c(n1 = 5, n0 = 5),
  c(n1 = 7, n0 = 3),
  c(n1 = 9, n0 = 1)
)

# Vector de probabilidades
p <- seq(0.001, 0.999, length.out = 500)

# Generar la base de datos con todas las curvas
df_logL <- data.frame()

for (esc in escenarios) {
  n1 <- esc["n1"]
  n0 <- esc["n0"]
  logL <- n1 * log(p) + n0 * log(1 - p)
  temp <- data.frame(p = p,
                     logL = logL,
                     Escenario = paste0("n1 = ", n1, ", n0 = ", n0))
  df_logL <- rbind(df_logL, temp)
}

# Gráfico
ggplot(df_logL, aes(x = p, y = logL, color = Escenario)) +
  geom_line(size = 1.2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Función de log-verosimilitud en varios escenarios Bernoulli",
    x = expression(p),
    y = expression(log *L(p)),
    color = "Escenarios"
  ) +
  theme_minimal(base_size = 11)

6 Tipos de modelos Bernoulli

Dependiendo de la forma como se especifique el modelo y la información disponible, un modelo Bernoulli puede clasificarse en alguno de los siguientes casos:

  • Modelo completo: la probabilidad \(p\) de éxito se iguala exactamente al valor observado para cada unidad experimental.

  • Modelo nulo: se asume que todos las observaciones tienen la misma probabilidad \(p\), sin incluir predictores explicativos.

  • Modelo saturado: se permite un valor específico de \(p\) para cada combinación de predictores.

En las siguientes secciones se describirán con más detalle sus propiedades, diferencias y cómo se calculan sus respectivas funciones de log-verosimilitud.

7 Modelo completo

7.0.1 El modelo y sus estimaciones

Definition 7.1 El modelo completo es caracterizado por el supuesto de que todos \(p_i\), \(i=1, \cdots ,n\) son considerados como parámetros.

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

Theorem 7.1 (Estimaciones en el modelo completo) En el modelo completo, las ML-estimaciones de \(p_i\) son \(\hat{p}_i=Y_i\) con valores \(\hat{p}_i=y_i\), para todo \(i=1, \cdots ,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}(y) = 0, \quad L(y)= 1 \quad \mbox{y}\quad -2\,{\cal L}(y)=0\]

7.0.2 Ejemplo 1: modelo completo

Example 7.1 Para los datos del archivo chdage, en el modelo completo, se cumplen los valores que muetran en el teorema anterior.

En R (manual).

En R, se puede verificar así:

datos <- lsm::chdage

Y <- CHD
a <- ifelse(Y==0,log(1-Y),log(Y))
LogCompleto <- sum(a)
L_Completo <- exp(LogCompleto)
DevLogCompleto <- -2*LogCompleto

Con el paquete lsm

El paquete lsm también permite calcular directamente el valor del log-verosímil para el modelo completo, denotado como \(\mathcal{L}(y)\), que representa el caso en el que el modelo ajusta exactamente cada observación individual. En este modelo, la probabilidad asignada a cada \(y_i\) coincide con su valor observado. Podemos acceder a este valor de dos formas:

datos <- lsm::chdage
modelo <- lsm(CHD ~ datos$AGE, data = datos)
modelo                    # Opción 1
modelo$Log_Lik_Complete   # Opción 2

Con lsm: 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:
## lsm(formula = CHD ~ datos$AGE, data = datos)
## 
## Populations in Saturate Model: 43
## 
## Coefficients: 
##                  CoefB  Std.Error        ExpB
## (Intercept) -5.3094534 1.13365464 0.004944629
## datos$AGE    0.1109211 0.02405984 1.117306795
## 
## Log_Likelihood: 
##          Estimation
## Complete    0.00000
## Null      -68.33149
## Logit     -53.67655
## Saturate  -41.79938

Con lsm: resultado con la opción 2.

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

## [1] 0

El valor \(\mathcal{L}(y)\) del modelo completo sirve como valor de referencia máximo que puede tomar la log-verosimilitud, ya que corresponde al mejor ajuste posible. Todos los demás modelos (nulo, logístico, saturado) tendrán log-verosimilitudes menores o iguales a esta. Por eso, resulta fundamental para evaluar la calidad relativa del ajuste de un modelo específico frente al ajuste perfecto.

8 Modelo nulo

8.0.1 El modelo y sus estimaciones

Definition 8.1 El modelo nulo es caracterizado por el supuesto de que todos los \(p_i\), \(i=1, \cdots ,n\) son considerados iguales; es decir, se tiene un solo parámetro \(p=p_i, i=1, \cdots ,n\).

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

\[\begin{eqnarray} {\cal L}(p)&=& n[\overline{y}\ln p + (1- \overline{y})\ln (1-p)] \tag{8.1} \end{eqnarray}\]

El siguiente teorema describe las estimaciones en este modelo:

Theorem 8.2 (Estimaciones en el modelo nulo) En el modelo nulo, la ML-estimación de \(p\) es \(\hat{p}=\overline{Y}\) con valor \(\hat{p}=\overline{y}\). Además,

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

8.0.2 Visualización de la función de verosimilitud

Remark. Algunas propiedades de la gráfica de \({\cal L}(p)\) son:

  1. Tiene un único mínimo en \(p=\frac{1}{2}\).

  2. Teniendo en cuenta que el dominio de \({\cal L}(p)\) es el intervalo abierto \((0,1)\), la gráfica de \({\cal L}(p)\):

    1. Es decreciente en \((0,\frac{1}{2})\) y creciente en \((\frac{1}{2},1)\).

    2. No tiene puntos de inflexión en \((0,1)\).

    3. Es cóncava hacia arriba.

    4. El punto mínimo es \(\left(\frac{1}{2},-n\log 2\right)\).

El siguiente teorema indica la forma de la gráfica de \({\cal L}(p)\), para \(0<p<1\).

Theorem 8.3 (Gráfica en el modelo nulo) La gráfica de la función \(\mathcal{L}(p)\) en el modelo nulo se presenta a continuación. Esta función representa la log-verosimilitud cuando se asume que todos los individuos tienen la misma probabilidad \(p\) de éxito, sin considerar variables explicativas.

n <- 100
p <- seq(0.0001, 1, 0.00005)
a <- n*(p*log(p)+(1-p)*log(1-p))
D <- data.frame(p,a)

ggplot(D, aes(y = a, x = p)) +
geom_point(aes(color=a), alpha = 1.9) +
labs(x="Probabilidad p", y="Log L(p)", fill= "") + 
ylim(-n*log(2), 0)+
facet_wrap(. ~ "Gráfica del Log L(p) en el modelo nulo con n=100") +    
theme_bw(base_size = 12) +
theme(legend.position = "none")

Esta primera gráfica muestra el comportamiento de la función \(\mathcal{L}(p)\) en el modelo nulo cuando el tamaño de la muestra es \(n = 100\). Se observa que la función alcanza su valor máximo (menos negativo) cuando \(p = 0.5\), lo cual representa el caso de máxima entropía (incertidumbre máxima). A medida que p se aleja de 0.5, la log-verosimilitud decrece, confirmando que este valor central ofrece el mejor ajuste bajo el supuesto de homogeneidad en la probabilidad.

8.0.3 Comparación para distintos tamaños de muestra

A continuación se presenta la misma función \(\mathcal{L}(p)\) en el modelo nulo, pero evaluada para diferentes valores de \(n\). Esto permite observar cómo el tamaño muestral afecta la forma y el rango de la log-verosimilitud.

# Generar datos para todos los valores de n con orden correcto
n_vals <- c(10, 50, 100, 200)
p <- seq(0.0001, 0.9999, 0.0005)
D <- data.frame()

for (n in n_vals) {
  a <- n * (p * log(p) + (1 - p) * log(1 - p))
  temp <- data.frame(p = p, a = a, n = n)
  D <- rbind(D, temp)
}

# Convertir a factor con niveles ordenados
D$n <- factor(D$n, levels = n_vals)

# Graficar
ggplot(D, aes(x = p, y = a)) +
  geom_line(color = "steelblue", size = 1) +
  labs(
    x = "Probabilidad p",
    y = expression(log * L(p)),
    title = "Función log L en el modelo nulo para distintos valores de n"
  ) +
  facet_wrap(~ n, scales = "fixed") +
  theme_bw(base_size = 12)

8.0.4 Ejemplo 2: modelo nulo

Example 8.1 Para los datos del archivo chdage, en el modelo nulo:

  1. \(\widehat{p} \;=\; \overline{y}\;=\; 43/100 \;=\; 0.43\).

  2. \({\cal L}(\overline{y})\;=\; {\cal L}(0.43) \;= \; -68.33\).

  3. \(-2{\cal L}(\overline{y}) \;=\; 136.66\).

  4. El punto mínimo es \(\left(\frac{1}{2},-69.31\right)\).

En R (manual).

En R, se puede verificar así:

n <- length(Y)
YBarra <- mean(Y)
LogNulo <- n*(YBarra*log(YBarra)+(1-YBarra)*log(1-YBarra))
DevNulo <- -2*LogNulo
minimo <- -n*log(2)

Con el paquete lsm.

El paquete lsm 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 <- lsm::chdage
modelo <- lsm(CHD~datos$AGE, data=datos)
modelo               # Opción 1
modelo$Log_Lik_Null  # Opción 2

Con lsm: 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:
## lsm(formula = CHD ~ datos$AGE, data = datos)
## 
## Populations in Saturate Model: 43
## 
## Coefficients: 
##                  CoefB  Std.Error        ExpB
## (Intercept) -5.3094534 1.13365464 0.004944629
## datos$AGE    0.1109211 0.02405984 1.117306795
## 
## Log_Likelihood: 
##          Estimation
## Complete    0.00000
## Null      -68.33149
## Logit     -53.67655
## Saturate  -41.79938

Con lsm: resultado con la opción 2.

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

## [1] -68.33149

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.

9 Modelo saturado

9.0.1 Supuesto 1

El modelo saturado está caracterizado por dos supuestos.

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

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

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

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

  • El número de observaciones \(Y_{ij}\) en cada población \(j\) por \(n_j\), siendo \(n_1+\cdots +n_J=n\);

  • La suma de las \(n_j\) observaciones \(Y_{ij}\) en \(j\) por

\[Z_j:=\sum\limits_{i=1}^{n_j}Y_{ij} \quad \mbox{con valor}\quad z_j=\sum\limits_{i=1}^{n_j}y_{ij},\quad \mbox{siendo}\quad \sum\limits^J_{j=1}z_j \;= \; \sum\limits^n_{i=1}y_i\]

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

Table 9.1: Ilustración de un conjunto de datos agrupado en \(J=4\) poblaciones
1 2 3 4 5 6 7 8 9
\(Y\) \(X_1\) \(X_2\) \(X_3\) \(X_4\) \(X_5\) \(j\) \(n_j\) \(Z_j\)
Población: Bajo, 80, Si, 170, Estrato 1
1 Bajo 80 Si 170 Estrato 1 \(j=1\) \(n_1=3\) \(Z_1=2\)
0 Bajo 80 Si 170 Estrato 1
1 Bajo 80 Si 170 Estrato 1
Población: Mediano, 100, Si, 150, Estrato 5
0 Mediano 100 Si 150 Estrato 5 \(j=2\) \(n_2=4\) \(Z_2=3\)
1 Mediano 100 Si 150 Estrato 5
1 Mediano 100 Si 150 Estrato 5
1 Mediano 100 Si 150 Estrato 5
Población: Mediano, 100, No, 180, Estrato 2
1 Mediano 100 No 180 Estrato 2 \(j=3\) \(n_3=2\) \(Z_3=1\)
0 Mediano 100 No 180 Estrato 2
Población: Alto, 100, No, 100, Estrato 4
0 Alto 100 No 100 Estrato 4 \(j=4\) \(n_4=3\) \(Z_4=2\)
1 Alto 100 No 100 Estrato 4
1 Alto 100 No 100 Estrato 4
General. \(Y\) es la variable de respuesta; \(X_1, \cdots, X_5\) son las variables explicativas; \(j\) es la población; \(n_j\) es el tamaño de la población \(j\); \(Z_j\) es el número de éxitos en la población \(j\).

9.0.2 Supuesto 2

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

  1. \((Y_{ij}|\star)\) es de Bernoulli. Es decir,

\[(Y_{ij}|\star) \sim {\cal B}(1,p_j)\]

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

  2. La esperanza y la varianza son, respectivamente,

\[p_j=P(Y_{ij}=1|\star)=E(Y_{ij}|\star), \qquad V(Y_{ij}|\star)=p_j(1-p_j)\]

9.0.3 Implicaciones del supuesto 2

Remark. A continuación, se oprimirá el símbolo \(\star\). El supuesto 2 implica:

  1. Todos los \(p_{ij}\), \(i=1, \cdots ,n\) dentro de cada población \(j\) son iguales. Es decir, se tiene como parámetro el vector \(p=(p_1, \cdots ,p_J)^T.\)

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

    • La variable \(Z_j\) es binomial. Es decir,

    \[Z_j\sim{\cal B}(n_j,p_j)\]

    • Las variables \(Z_j\) son independientes entre las poblaciones.

9.0.4 Log L y estimaciones

Theorem 9.1 (Log L 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(\sum_{i=1}^{n_j}[{y_{ij}}\ln p_j + (1- y_{ij})\ln (1-p_j)]\right)\nonumber\\ &=& \sum^J_{j=1}[{z_j}\ln p_j + (n_j- z_j)\ln (1-p_j)] \tag{9.1} \end{eqnarray}\]

Theorem 9.2 (Estimaciones en el modelo saturado) En el modelo saturado, las ML-estimaciones de \(p_j\) son \(\tilde{p}_j=\frac{Z_j}{n_j}\), con valores \(\tilde{p}_j=\frac{z_j}{n_j}\),\(j=1,\cdots ,J\). Además,

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

También se cumple que

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

9.0.5 Ejemplo 3: modelo saturado

Example 9.1 Para los datos del archivo chdage, en el modelo saturado,hay \(J=43\) poblaciones y se cumple que \({\cal L}(\tilde{p})=-41.7991\).

En R (manual).

En R, se puede verificar así (véase la última fila de la Tabla 9.2):

datos <- lsm::chdage

datos %>%
  group_by(AGE) %>%
  summarise(nj = n(),
            zj = sum(CHD)) %>%
  mutate(pj = zj/nj,
         pj = round(pj,3),
         Lp = ifelse(zj==0 | zj== nj, 0, zj*log(pj)+(nj-zj)*log(1-pj)),
         Lp = round(Lp, 4)
         ) -> saturado

Totales <- c("Total", sum(saturado$nj), sum(saturado$zj), "", sum(saturado$Lp))
Tabla <- rbind(saturado, Totales)
Table 9.2: Estimaciones en el modelo saturado.
AGE nj zj pj Lp
20 1 0 0 0
23 1 0 0 0
24 1 0 0 0
25 2 1 0.5 -1.3863
: : : : :
: : : : :
63 1 1 1 0
64 2 1 0.5 -1.3863
65 1 1 1 0
69 1 1 1 0
Total 100 43 -41.7991

Con el paquete lsm.

El paquete lsm 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 al ajustar perfectamente cada observación.

Podemos obtener ambos valores de la siguiente manera:

datos <- lsm::chdage
modelo <- lsm(CHD~datos$AGE, data=datos)
modelo                                              # Opción 1
cbind(modelo$Populations, modelo$Log_Lik_Saturate)  # Opción 2

Con lsm: resultado con la opción 1.

Como se observa en el bloque de salida del modelo, los valores de \(J\) y \(\mathcal{L}(\tilde{p})\) se muestran como Populations in Saturate Model (parte superior) y Saturate (parte inferior), respectivamente:

## 
## Call:
## lsm(formula = CHD ~ datos$AGE, data = datos)
## 
## Populations in Saturate Model: 43
## 
## Coefficients: 
##                  CoefB  Std.Error        ExpB
## (Intercept) -5.3094534 1.13365464 0.004944629
## datos$AGE    0.1109211 0.02405984 1.117306795
## 
## Log_Likelihood: 
##          Estimation
## Complete    0.00000
## Null      -68.33149
## Logit     -53.67655
## Saturate  -41.79938

Con lsm: resultado con la opción 2.

Ambos valores pueden ser extraídos directamente usando modelo$Populations y modelo$Log_Lik_Saturate:

## Número de poblaciones (J): 43
## Log-verosimilitud del modelo saturado: -41.79938

El valor \(\mathcal{L}(\tilde{p})\) es útil para comparar con otros modelos como el nulo o el logístico, ya que representa el mejor ajuste posible a los datos. Todos los demás modelos tendrán log-verosimilitudes menores o iguales a esta

10 Ejercicios

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

  • Todos los datos mencionados aparecen en los paquetes 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.

  • Verifique cómo se obtienen las estimaciones correspondientes, los logaritmos de las funciones de máxima verosimilitud, las gráficas correspondientes, etc.

10.0.1 Ejercicios 1 a 3

  1. Demuestre los teoremas: (a) 7.1; (b) 8.2; (c) 9.2.

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

  3. Haga un listado de los paquetes de R que estimen el logaritmo de la función de máxima verosimilitud en el modelo saturado.

10.0.2 Ejercicio 4

  1. Los datos ICU corresponden a una muestra de 200 sujetos que hicieron parte de un estudio de supervivencia de pacientes que fueron remitidos a una unidad de cuidados intensivos (intensive care unit - ICU). La meta principal de este estudio fue desarrollar un modelo de regresión logística para predecir la probabilidad de supervivencia de estos pacientes en el hospital y estudiar los factores de riesgos asociados con el índice de mortalidad ICU. En estos datos tome a la variable AGE como independiente y STA como dependiente.

  1. ¿Qué características de STA nos pone a pensar que debamos considerar el modelo de regresión logística en vez del usual modelo de regresión lineal para describir la relación entre STA y AGE?

  2. Forme un diagrama de dispersión de STA contra AGE.

  3. Use los intervalos [15,24], [25,34], [35,44], [45,54], [55,64], [65,74], [75,84], [85,94] para AGE:

  • Calcule la media de STA de los sujetos dentro de cada intervalo.

  • Grafique estos valores de la media de STA contra el punto medio del intervalo de AGE usando el mismo conjunto de ejes que se utilizaron en la parte (b).

  1. Repita todos los análisis realizados en este documento.

10.0.3 Ejercicios 5 a 7

  1. Considere los datos ICU. Repita el ejercicio 4 utilizando la variable TYP (como variable dependiente) en vez de STA.

  2. Considere los datos ICU. Repita todos los análisis realizados en este documento, pero considerando ahora las variables AGE (como variable independiente) y STA (como variable dependiente).

  3. Considere los datos ICU. Haga el análisis correspondiente tomando a STA como variable dependiente y a AGE, SYS y HRA como independientes.

10.0.4 Ejercicios 8 a 9

  1. Los datos UIS se recogieron con el propósito de comparar dos programas de tratamiento A y B para reducir el abuso de la droga y prevenir sus riesgos. La descipción de los datos se puede ver también aquí. Detalle el análisis para estos datos, tomando a DFREE como variable dependiente y AGE, BECK y NDRUGTX como variables independientes.

  2. Los datos PROS corresponden a un estudio realizado pacientes con cáncer de próstata para determinar si las variables medidas en un examen básico pueden ser usadas para predecir si el tumor ha penetrado la cápsula prostática. Los datos fueron recogidos teniendo en cuenta 380 individuos, 153 de los cuales tuvieron un cáncer que penetró la cápsula prostática. En estos datos, una variable que se pensó que era particularmente predictiva para la penetración de cápsula es el nivel de antígeno prostático, PSA.

  • Repita los pasos del ejercicio 4 usando CAPSULE como variable dependiente y utilize para PSA, los intervalos: [0.0; 2.4], [2.5; 4.4], [4.5; 6.4], [6.5; 8.4], [8.5; 10.4], [10.5; 12.4], [12.5; 20.4], [20.5; 140].

10.0.5 Ejercicio 10

  1. De todas las variables que aparecen en los datos PROS sólo considere a CAPSULE (como variable dependiente) y PSA (como variable independiente).

  1. ¿Qué características de la variable dependiente nos conduce a considerar la regresión logística como más apropiada que el modelo de regresión lineal para describir la relación entre las dos variables mencionadas anteriormente?

  2. Calcule:

  • \(\mathcal{L}(\widehat{p})\) en el modelo completo.

  • \(\mathcal{L}(\widehat{p})\) en el modelo nulo.

  • \(\mathcal{L}(\widetilde{p})\) en el modelo saturado.

  1. Halle las estimaciones de los siguientes parámetros e interprételos (justifique en forma clara y precisa todas sus afirmaciones):

  • \(P(CAPSULE=1 \, / \, PSA=11.2\) mg/ml\()\).

  • \(P(CAPSULE=0 \, / \, PSA=11.2\) mg/ml\()\).

10.0.6 Ejercicios 11 a 13

  1. Considere los datos PROS. Repita todos los análisis realizados en este documento tomando a CAPSULE como variable dependiente y VOL como variable independiente.

  2. Considere los datos PROS. Repita todos los análisis realizados en este documento tomando a CAPSULE como variable dependiente y AGE como variable independiente.

  3. Considere los datos PROS, tomando a CAPSULE como variable dependiente y AGE, PSA y VOL como variables independientes. Repita todos los análisis realizados en este documento

10.0.7 Ejercicio 14

  1. Los datos LOWBWT corresponden a un estudio realizado para identificar factores de riesgos asociados a nacimientos de bebés con bajo peso (peso menor que 2.500 gramos). Los datos fueron recogidos teniendo en cuenta 189 mujeres, 59 de las cuales tuvieron bebós con bajo peso y 130 de las cuales tuvieron bebés con peso normal. De todas las variables que aparecen sólo considere a LOW (como variable dependiente) y LWT (como variable independiente).

  1. ¿Qué características de la variable dependiente nos conduce a considerar la regresión logística como más apropiada que el modelo de regresión lineal para describir la relación entre las dos variables mencionadas anteriormente?

  2. Calcule:

  • \(\mathcal{L}(\widehat{p})\) en el modelo completo.

  • \(\mathcal{L}(\widehat{p})\) en el modelo nulo.

  • \(\mathcal{L}(\widetilde{p})\) en el modelo saturado.

  1. Halle las estimaciones para los siguientes parámetros e interprételos (justifique en forma clara y precisa todas sus afirmaciones):

  • \(P(LOW=1 \, / \, LWT=100.3\) libras\()\).

  • \(P(LOW=0 \, / \, LWT=100.3\) libras\()\).

10.0.8 Ejercicios 15 a 18

  1. Considere los datos LOWBWT, tomando a LOW como variable dependiente y AGE como variable independiente.

  2. Considere los datos LOWBWT, tomando a LOW como variable dependiente y LWT como variable independiente.

  3. Considere los datos LOWBWT, tomando a LOW como variable dependiente y BWT como variable independiente.

  4. Considere los datos LOWBWT, tomando a LOW como variable dependiente y AGE, LWT y BWT como variables independientes.

10.0.9 Ejercicios 19 a 21

  1. Considere los datos LOWBWTM11, tomando a LOW como variable dependiente y AGE como variable independiente.

  2. Considere los datos LOWBWTM11, tomando a LOW como variable dependiente y LWT como variable independiente.

  3. Considere los datos LOWBWTM11, tomando a LOW como variable dependiente y AGE y LWT como variables independientes.

10.0.10 Ejercicios 22

  1. En los datos CLSLOWBWT, entre otras variables, una que los físicos consideraron importante para el control del peso del bebé (variable dependiente LOW) fue el peso de la madre durante el último periodo menstrual (LWT).

  • Repita los pasos del ejercicio 4, pero para la parte (c) utilize los intervalos: [80,99], [100,109], [110,114], [115,119], [120,124], [125,129], [130,250].

  • La gráfica en la parte (c) no parece en forma de \(S\). La razón principal es que el rango de los valores graficados está aproximadamente entre 0.2 y 0.56.

  • Explique por qué un modelo para la probabilidad de LOW como una función de LWT pudiese ser el modelo de regresión logística.

10.0.11 Ejercicios 23 a 26

  1. Considere los datos CLSLOWBWT, tomando a LOW como variable dependiente y AGE como variable independiente.

  2. Considere los datos CLSLOWBWT, tomando a LOW como variable dependiente y LWT como variable independiente.

  3. Considere los datos CLSLOWBWT, tomando a LOW como variable dependiente y BWT como variable independiente.

  4. Considere los datos CLSLOWBWT, tomando a LOW como variable dependiente y AGE, LWT y BWT como variables independientes.

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.