EL CASO BINARIO

Regresión logística (matriz de confusión)

Departamento de Matemáticas y Estadística, Universidad del Norte (Barranquilla, Colombia)

hllinas@uninorte.edu.co

Biographical sketch

16/09/25

Abstract

La teoría mencionada puede revisarse en el volumen 8 de mis notas de clase que aparecen en el siguiente documento: 2.2. Regresión logística y en la referencia: LLinás (2006). En Rpubs:: toc se pueden ver otros documentos de posible interés.

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:

  • Rpubs::Mis paquetes.

  • Rpubs::Paquete lsm en R.

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.

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

3 Datasets

Para las aplicaciones, se 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

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

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

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

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

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

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}(\overline{y})\): 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}(\overline{y})\) 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}(\overline{y})\) 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 Modelo logístico

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

10.0.2 Supuestos

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

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

tiene rango completo \(Rg(C)=1+K\leq J\). 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á.

10.0.3 El modelo

Hypothesis 10.2 (Supuesto 4: modelo logístico) Para llegar a un modelo logístico se hace el supuesto adicional:

\[\begin{equation} \mbox{Logit}(p_j):= \ln\left(\frac{p_j}{1-p_j}\right) = \delta \;+\; \beta_1 \,x_{j1} \;+\;\cdots \;+\; \beta_K \,x_{jK} \tag{10.1} \end{equation}\]

Remark. Tenemos que:

1. \(\alpha =(\delta,\beta_1,\ldots,\beta_K)^T\) es el vector de parámetros en el modelo.

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

11 Riesgo

11.0.1 Definición

Definition 11.1 (Riesgo) En la práctica, la probabilidad \(p_j\) es conocida como riesgo.

Theorem 11.1 (Fórmula para el riesgo) Sea \(g_j:=\delta \;+\; \beta_1\,x_{1j} \;+\;\cdots \;+\; \beta_K \,x_{Kj}\). Entonces, la probabilidad

\[\begin{equation} p_j\; =\; P(Y_j=1|x_{j1}, \cdots, x_{jK}) \tag{11.1} \end{equation}\]

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

\[\begin{equation} p_j \; =\; \mbox{Logit}^{-1}(g_j) \;= \; \frac{e^{g_j}} {1 + e^{g_j}} \tag{11.2} \end{equation}\]

11.0.2 Ejemplo

Example 11.1 En la figura 11.1 aparece la gráfica de las dos funciones siguientes:

\[\begin{eqnarray*} p &=& \frac{e^{x}} {1 + e^{x}} \qquad \mbox{(Figura a, roja, con $\beta=1$)},\\ p &=&\frac{e^{-x}} {1 + e^{-x}} \qquad \mbox{(Figura b, verde, con $\beta=-1$)} \end{eqnarray*}\]

Observe que el signo de la pendiente (o sea, el valor numérico que multiplica a \(x\)) influye en el hecho si la gráfica es creciente o decreciente. La No.1 es creciente (por tener pendiente postiva) y la No.2, es decreciente (por tener pendiente negativa).

VI <- seq(-10, 10, 0.05)
Crec <- exp(VI)/(1+exp(VI))
Decr <- exp(-VI)/(1+exp(-VI))

ggplot() +
geom_point(mapping=aes(y=Crec, x = VI, color="Crec"), size=1.5) +
geom_point(mapping=aes(y=Decr, x = VI, color="Decr"),size=1.5) +  
labs(x="Variable explicativa X", y="Probabilidad de éxito p", fill= "") + 
#ylim(0, 100)+
facet_wrap(. ~ "Gráfica de dos logits con una sola variabla explicativa") +    
theme_bw(base_size = 12) +
#theme(legend.position = "none")+
  scale_color_manual(values = c("Crec" = "darkblue","Decr" = "red")) +
  labs(color = "Logits")+
scale_color_discrete(name = expression(paste("Pendiente", " ", beta, ":")), labels = c(expression(paste("(a)", " ", beta>0)), expression(paste("(b)", " ", beta<0)))) # Edit legend title and labels

Figure 11.1: Comparación de dos Logits

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

Theorem 12.1 (Log de la función de verosimilitud en el modelo logístico) Sea \(g_j:=\delta \;+\; \beta_1\,x_{1j} \;+\;\cdots \;+\; \beta_K \,x_{Kj}\). Reescribiendo \({\cal L}(p)\), dada en la ecuación (9.1), el logaritmo de la función de verosimilitud se puede escribir, en función de \(\alpha\), como:

\[\begin{eqnarray} {\cal L}({\alpha}) &=& \sum^{J}_{j=1}\left[z_j\ln\left(\frac{p_j}{1-p_j}\right)+ n_j\ln(1-p_j)\right]\\ &=& \sum^{J}_{j=1}z_j\,g_j \;-\; \sum^{J}_{j=1} n_j\ln \left[1+ e^{g_j}\right] \tag{12.1} \end{eqnarray}\]

En la primera expresión se observa que el paso de \(p_j\) hacia \(\ln\left(\frac{p_j} {1-p_j}\right)\) aparece de una manera natural.

13 Odds

Definition 13.1 (Odds) Un odds se define como la proporción entre las probabilidades de ocurrencia y no ocurrencia del evento que se relaciona con \(Y\) en la población \(j\). Es decir, es el cociente

\[\begin{equation} O_j\;=\;\frac{p_j}{1-p_j} \tag{13.1} \end{equation}\]

Remark (relación entre Riesgo y Odds). Se resalta el hecho de que los riesgos toman valores entre 0 y 1. En cambio, los odds toman valores entre 0 e infinito. Además, observe que:

\[\begin{equation} p_j\;=\;\frac{O_j}{1+O_j} \tag{13.2} \end{equation}\]

Example 13.1 En la Tabla 13.1 se presenta la relación entre riesgos y Odds. Se observa que cuando el riesgo aumenta, aumenta el odds. Esta propiedad se puede visualizar de manera más clara en la Figura 13.1.

Table 13.1: Equivalencia entre Riesgo y Odds.

1	2	3
Riesgo (\(p_j\))	No Riesgo (\(1-p_j\))	Odds (\(O_j\))
0.1	0.9	0.1/0.9 = 0.11
0.2	0.8	0.2/0.8 = 0.25
0.3	0.7	0.3/0.7 = 0.43
0.4	0.6	0.4/0.6 = 0.67
0.5	0.5	0.5/0.5 = 1.00
0.6	0.4	0.6/0.4 = 1.50
0.7	0.3	0.7/0.3 = 2.33
0.8	0.2	0.8/0.2 = 4.00
0.9	0.1	0.9/0.1 = 9.00

p <- seq(0.0001, 1, 0.00005)
Odds <- p/(1-p)
D <- data.frame(p,Odds)

ggplot(D, aes(y = Odds, x = p)) +
geom_point(aes(color=Odds), alpha = 1.9) +
labs(x="Riesgo p", y="Odds", fill= "") + 
ylim(0, 100)+
facet_wrap(. ~ "Riesgo versus Odds") +    
theme_bw(base_size = 12) +
theme(legend.position = "none")

Figure 13.1: Relación entre Riesgo y Odds

14 Riesgo relativo RR

Definition 14.1 (Riesgo relativo) El riesgo relativo se define como el cociente entre el riesgo en un grupo con un factor de exposición o de riesgo (población \(i\)) y el riesgo en un grupo de referencia, que no tiene el factor de exposición (población \(j\)). Es decir, es el siguiente cociente:

\[RR(i\; \mbox{vs} \; j) = \frac{\mbox{Incidencia acumulada en la población $i$}}{\mbox{Incidencia acumulada en la población $j$}}\]

O sea, el cociente entre las siguientes probabilidades de éxitos:

\[RR(i\; \mbox{vs} \; j) \;=\; \frac{p_i}{p_j}=\frac{P(\mbox{$Y=1$ | población $i$})}{P(\mbox{$Y=1$ | población $j$})}\]

Example 14.1 Una de las aplicaciones del cálculo del riesgo relativo es en el área de las ciencias de la salud, especificamente, en los estudios prospectivos (como, por ejemplo, el estudio de cohortes y el ensayo clínico). Para ello, de la población en estudio, se seleccionan dos muestras sin enfermedad, en donde una está expuesta al factor de riesgo (grupo \(i\)) y la otra, no (grupo \(j\)) y se hace un seguimiento del estudio en el tiempo. En la Figura 14.1 se muestra una estructura básica de un estudio de cohortes.

Figure 14.1: Estudio de cohorte. Fuente: Elaboración propia.

Supongamos que los datos encontrados en un estudio de cohorte son los que se muestran en la Tabla 14.1.

Table 14.1: Expuestos y no expuestos.

1	2	3	4
Factor	Expuesto (\(x=1\))	No expuesto (\(x=0\))	Total
Enfermo (\(y=1\))	\(a\)	\(c\)	\(a+c\)
Sano (\(y=0\))	\(b\)	\(d\)	\(b+d\)
Total	\(a+b\)	\(c+d\)	\(n\)

De cada muestra se calcula incidencia acumulada de expuestos y se halla su cociente (es decir, el riesgo relativo):

\[RR(i\; \mbox{vs} \; j) = \frac{\mbox{Incidencia acumulada en expuestos (población $i$)}}{\mbox{Incidencia acumulada en no expuestos (población $j$) }} = \frac{P_i(\mbox{Enfermo| Expuesto})}{P_j(\mbox{Enfermo | No expuesto })}=\frac{a/(a+b)}{c/(c+d)}\] ::: {.example #unnamed-chunk-58} Si RR=30, entonces, podemos interpretar de la siguiente manera: la probabilidad de que los expuestos adquieran la enfermedad es 30 veces la probabilidad de que los no expuestos la desarrollen. :::

15 Características del RR

1. Es adimensional y su valor se encuentra entre 0 e infinito.

2. Permite comparar la frecuencia de ocurrencia del evento entre los que tienen el factor de riesgo y los que no lo tienen. En este sentido, se puede considerar como medida de la magnitud o fuerza de la asociación. Algunas interpretaciones son las siguiente:

   • RR=1: No hay asociación entre la presencia del factor de riesgo y el evento.

   • RR>1: Existe asociación positiva. Es decir, la presencia del factor de riesgo se asocia a una mayor frecuencia de suceder el evento.

   • RR<1: Existe una asociación negativa. Es decir, no existe factor de riesgo, sino uno protector.

3. El riesgo relativo no puede utilizarse en los estudios de casos y controles o retrospectivos. Esto es así porque no es posible calcular las tasas de incidencia. En estos casos, se utilizará la razón de momios o razón odds (en inglés: odds ratio), concepto que se explicará a continuación.

16 Razón odds

En estudios de cohortes el RR se estima de forma directa ya que se conoce la incidencia de la enfermedad en expuestos y en no expuestos. Por el contrario, en los estudios de casos y controles no se puede calcular la incidencia, porque la población de estudio se selecciona a partir de individuos que ya han desarrollado la enfermedad. Por esta razón, en los estudio de casos y controles se calcula la razón de odds u odds ratio (OR). En la Figura 16.1 se muestra una estructura básica de un estudio de casos y controles.

Figure 16.1: Estudio de cohorte. Fuente: Notas metodológicas

Definition 16.1 (Razón odds) Una razón ODDS se define como el cociente entre dos odds:

\[OR(i\; \mbox{vs} \;j) \;=\; \frac{O_i}{O_j} \;=\; \frac{\frac{p_i}{1-p_i}}{\frac{p_j}{1-p_j}}\]

Theorem 16.1 (OR es exponencial de la pendiente) Siempre se cumple que \(OR(i\; \mbox{vs} \;j)\) es un número entre 0 e infinito. Además, en un modelo de regresión logística, se cumple que:

\[\begin{eqnarray} OR(i\; \mbox{vs} \;j) \;= \; e^{\beta_1(x_{i1}-x_{j1}) \; + \; \beta_2(x_{i2}-x_{j2}) \; + \; \cdots \; + \; \beta_K(x_{iK}-x_{jK})} \end{eqnarray}\]

Cuando \(x_{ik}-x_{jk}=1\) para todo \(k=1, \ldots, K\), entonces

\[OR:=OR(i\; \mbox{vs} \;j) \;= \;e^{\beta_1\;+\;\cdots \;+\;\beta_K}\]

Es decir, no depende de \(X_1, \ldots, X_K\) y muestra el cambio proporcional en la variable de respuesta cuando las variables independientes se incrementen en \(1\) unidad.

Example 16.1 Supongamos que los datos encontrados en un estudio de casos y controles son los que se muestran en la Tabla 16.1:

• Casos: se refiere aquellos individuos con la enfermedad al final del estudio.

• Controles: son los que no la padecen.

• Expuesto: los sujetos de estudio que tienen el factor de riesgo.

• No expuesto: los que no lo tienen.

Table 16.1: Expuestos y no expuestos.

1	2	3	4
Factor	Expuesto (\(x=1\))	No expuesto (\(x=0\))	Total
Caso (\(y=1\))	\(a\)	\(c\)	\(a+c\)
Control (\(y=0\))	\(b\)	\(d\)	\(b+d\)
Total	\(a+b\)	\(c+d\)	\(n\)

En este caso, en un estudio de casos y controles, la razón de momios es el cociente entre el odds de enfermedad en el grupo expuesto (o en el grupo tratado, población \(i\)) \(a/b\) y el odds de enfermedad en el grupo no expuesto (o no tratado, población \(j\)) \(c/d\):

\[OR(i\; \mbox{vs} \;j) \;=\; \frac{O_i}{O_j} \;=\; \frac{\frac{p_i}{1-p_i}}{\frac{p_j}{1-p_j}} = \frac{a/b}{c/d} = \frac{ad}{bc}\]

Example 16.2 Supongamos que en un caso particular, OR=3. Entonces, se puede interpretar así: la razón entre la presencia (casos) versus la no presencia de la enfermedad (controles) es 3 veces mayor en las personas expuestas al factor (población \(i\)) en comparación a las personas no expuestas (población \(j\)).

17 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 8 en LLinás (2006). En R, las funciones glm() y lsm() calculan estas estimaciones.

18 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\) o

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

4. En la práctica:

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

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

19 Ejemplo 3

19.0.1 Enunciado

Considere los datos del archivo chdage. Suponga que se quiere analizar un modelo de regresión logística, considerando a chd como variable dependiente y age como independiente.

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\), sin utilizar la función summary().

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

6. Utilizando las estimaciones halladas en el inciso anterior, escriba en el modelo correspondiente .

7. Haga la gráfica del riesgo de tener enfermedades coronarias versus la edad. ¿Es directa o indirecta esta relación?

8. Haga la gráfica del logit estimado versus la edad. ¿Qué tipo de relación hay entres estas dos variables?

9. Estime el logit estimado, para un sujeto con 50 años.

10. Estime la proporción estimada de personas con presencia de enfermedades coronarias a la edad de 50.

11. Halle los errores estándares estimados \(\widehat{S}_{\widehat{\beta}}\) y \(\widehat{S}_{\widehat{\delta}}\) de \(\widehat{\beta}\) y \(\widehat{\delta}\), respectivamente.

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

13. Encuentre el odds estimado para individuos con edad de 50 e interprételo.

14. Haga la gráfica del odds estimado versus la edad.

15. Halle la razón odds estimada cuando el incremento la edad se incrementa en 1 año. Inteprétela.

19.0.2 Solución

Consultar la solución del ejeplo 3 del documento: Rpubs::Regresión logística (estimaciones).

20 Usando matrices Esm y Elm

En el objeto generado por lsm(), las matrices Esm y Elm contienen resúmenes por población (\(j\)) en dos modelos distintos:

• Esm (Estimaciones en el modelo saturado) incluye, para cada población \(j\), las variables explicativas, el tamaño de la población (\(n_j\)), el número de éxitos (\(z_j\)), la probabilidad estimada \(\hat{p}_j\) sin restricciones logísticas, y la log-verosimilitud \(L_j\) correspondiente al modelo saturado.

• Elm (Estimaciones en el modelo logístico) presenta, para cada población \(j\), las mismas variables explicativas, \(n_j\) y \(z_j\), pero con las probabilidades estimadas \(\hat{p}_j\) del modelo logístico, el logit correspondiente (\(\text{Logit}(\hat{p}_j)\)), la log-verosimilitud \(L_j\) y la varianza del logit (Var.logit).

Estas matrices pueden unirse (por ejemplo, con cbind) para comparar directamente las estimaciones de ambos modelos en una misma tabla, lo que facilita el análisis y la extracción de métricas específicas para valores concretos de la variable explicativa.

library(knitr)
library(kableExtra)

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

# Matrices que el modelo entrega por población j
sm <- modelo$Esm
lm <- modelo$Elm

# Unimos 'Esm' con columnas relevantes de 'Elm' (Lp, gj, pj, Lj, Var.logit)
slm <- cbind(sm, lm[, 3:6])

# Cambiamos los nombres de columna
colnames(slm) <- c("Age", "nj", "zj", "pj_tilde", "L_saturado", 
                       "gj", "pj", "L_logistico", "Var_logit")

# Formatear una parte de la tabla
kable(slm[25:30, ]) %>%
kable_styling() %>%
kable_classic_2(full_width = TRUE) #%>%

	Age	nj	zj	pj_tilde	L_saturado	gj	pj	L_logistico	Var_logit
60	48	3	2	0.6666667	-1.909542	0.0147615	0.5036903	-2.0721425	0.0579086
63	49	3	1	0.3333333	-1.909542	0.1256826	0.5313794	-2.1482025	0.0607055
66	50	2	1	0.5000000	-1.386294	0.2366037	0.5588765	-1.4002572	0.0646601
68	51	1	0	0.0000000	0.000000	0.3475249	0.5860172	-0.8819309	0.0697725
69	52	2	1	0.5000000	-1.386294	0.4584460	0.6126455	-1.4383838	0.0760427
71	53	2	2	1.0000000	0.000000	0.5693672	0.6386171	-0.8969003	0.0834706

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

22 Comparación de modelos

Esta sección se refiere a resaltar que hay estadísticos para distintas pruebas de comparación de modelos:

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

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

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

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

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

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

23 Matriz de confusión

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

Table 23.1: Matrix de confusión.

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

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

1. La exactitud ( en inglés: accuracy) es la proporción de clasificaciones correctas. Es decir, mide el porcentaje global de predicciones correctas. Es útil como referencia general, pero puede ser engañosa en casos de clases desbalanceadas (por ejemplo, detectar enfermedades raras).

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

2. La sensitividad o tasa de verdaderos positivos (en inglés: sensivity, true positive rate, recall or hit rate) es la proporción de casos positivos observados que son correctamente clasificados. Evalúa qué tan bien el modelo detecta los casos positivos reales. Es clave en contextos donde perder un positivo es muy costoso (ejemplo: no diagnosticar un paciente enfermo).

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

3. La especificidad o tasa de verdaderos negativos (en inglés: specifity, selectivity or true negative rate) es la proporción de casos negativos observados que son correctamente clasificados. Indica la capacidad del modelo de identificar correctamente los casos negativos. Es importante cuando clasificar erróneamente un negativo como positivo tiene consecuencias graves (ejemplo: falsos positivos en pruebas judiciales).

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

4. La tasa de falsos positivos (1- especificidad) es la proporción de casos negativos observados que son incorrectamente clasificados. Cuantifica la proporción de negativos mal clasificados como positivos. Se usa para analizar el riesgo de alarmas falsas.

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

5. La precisión o valor predictivo positivo VPP (en inglés: precision, positive predictive value) es la proporción de casos positivos que son correctamente clasificados. Mide qué porcentaje de los casos predichos como positivos son realmente positivos. Es útil cuando se quiere confiar en que una predicción positiva es verdadera (ejemplo: spam filters).

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

6. El valor predictivo negativo VPN (en inglés: negative predictive value) es la proporción de casos negativos que son correctamente clasificados. Indica qué tan confiables son las predicciones de casos negativos. Es relevante cuando se necesita garantizar seguridad al descartar un positivo (ejemplo: pruebas médicas que descarten una enfermedad).

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

7. La prevalencia (en inglés, prevalence) es la proporción de casos positivos observados. Es decir, describe la proporción de positivos reales en la población. Sirve como contexto de referencia para interpretar las demás métricas.

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

8. La exactitud balanceada (en inglés: balanced accuracy) compensa el sesgo de la exactitud cuando hay clases desbalanceadas, promediando sensitividad y especificidad. Muy usada en problemas con datos desbalanceados.

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

9. La tasa de detección (en inglés: Detection rate) se define como la proporción de casos positivos observados y predichos. Representa la proporción de casos que son realmente positivos y además fueron correctamente clasificados. Da una idea de la cobertura real del modelo.

\[\mbox{Tasa de detección} \;=\; \frac{a}{n}\]

10. La prevalencia de detección (en inglés: Detection prevalence) se define como la proporción de casos positivos predichos. Indica la proporción de casos predichos como positivos (sean correctos o no). Permite comparar la tendencia del modelo a clasificar en positivo con la prevalencia real.

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

24 Indices de concordancias

24.0.1 Índice kappa de Cohen

Los índices de concordancia permiten medir el grado de acuerdo entre dos observadores o clasificadores. Sus valores oscilan entre 0 (total desacuerdo) y 1 (máximo acuerdo).. Un ejemplo es la exactitud (definido anteriormente), pero el más usado es el denominado índice kappa de Cohen (\(\kappa\)) y se define como:

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

donde \(P_E\) es la proporción de acuerdos esperados, dada por:

\[P_E \; = \; \frac{(\mbox{Total fila y=1)(Total columna y=1)} \;+\;\mbox{(Total fila y=0)(Total columna y=0)}}{\mbox{(Tamaño de la muestra)}^2}\] o, en términos de la tabla de contingencia:

\[P_E \; = \; \frac{(a+c)(a+b) + (b+d)(c+d)}{n^2}\]

24.0.2 ¿Para qué sirve?

• Sirve para responder a la pregunta: ¿hasta qué punto dos evaluadores (o dos métodos) clasifican de la misma forma a los mismos individuos?

• Es muy utilizado en medicina, psicología, educación y machine learning, donde se necesita medir la fiabilidad de una clasificación (por ejemplo, si dos médicos diagnostican lo mismo, si dos profesores califican igual, o si un modelo predice como lo haría un experto).

24.0.3 Escala de valoración para \(\kappa\)

El valor de \(\kappa\) se interpreta con la siguiente escala:

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

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

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

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

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

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

24.0.4 Ejemplo 4: \(\kappa\) (tabla de 2×2)

Enunciado

Supongamos dos evaluadores (A y B) que clasifican a 100 casos en Positivo o Negativo. La tabla de contingencia observada es:

• \(a = 50\): ambos dicen Positivo.

• \(b = 10\): A dice Positivo, B dice Negativo.

• \(c = 5\): A dice Negativo, B dice Positivo.

• \(d = 35\): ambos dicen Negativo.

• \(n = a+b+c+d = 100\)

Solución: Mostrar tabla y métricas básicas

library(knitr)
library(kableExtra)

# Datos (a, b, c, d)
a <- 50; b <- 10; c <- 5; d <- 35
n <- a + b + c + d

# Tabla de contingencia (A = filas, B = columnas)
tab <- matrix(c(a, b, c, d), nrow = 2, byrow = TRUE,
              dimnames = list(A = c("Positivo", "Negativo"),
                              B = c("Positivo", "Negativo")))

# Totales marginales
fila_pos <- a + b; fila_neg <- c + d      # Totales por filas (A)
col_pos  <- a + c; col_neg  <- b + d      # Totales por columnas (B)

# Exactitud observada (Po)
Po <- (a + d) / n

# Acuerdo esperado por azar (Pe)
Pe <- ((fila_pos * col_pos) + (fila_neg * col_neg)) / (n^2)

# Kappa de Cohen
kappa <- (Po - Pe) / (1 - Pe)

# Mostrar tabla y resultados
kable(addmargins(tab), align = "ccc",
      col.names = c("Positivo", "Negativo", "Total"),
      caption = "Tabla de contingencia A (filas) vs B (columnas)") %>%
  kable_styling(full_width = FALSE)

cat("Exactitud observada (Po):", round(Po, 4), "\n")
cat("Acuerdo esperado (Pe):", round(Pe, 4), "\n")
cat("Kappa de Cohen (κ):", round(kappa, 4), "\n")

Table 24.1: Table 24.2: Tabla de contingencia A (filas) vs B (columnas)

	Positivo	Negativo	Total
Positivo	50	10	60
Negativo	5	35	40
Sum	55	45	100

## Exactitud observada (Po): 0.85

## Acuerdo esperado (Pe): 0.51

## Kappa de Cohen (κ): 0.6939

Cálculos esperados (a mano)

\[\begin{eqnarray*} P_o &=& \frac{a+d}{n} \;=\; \frac{50+35}{100} \;=\; 0.85 \\ && \\ P_E &=& \frac{(a+c)(a+b) + (b+d)(c+d)}{n^2} \;=\; \frac{(55)(60) + (45)(40)}{100^2} \;=\; \frac{3300 + 1800}{10{,}000} \;=\; 0.51 \\ &&\\ \kappa &=& \frac{P_o - P_E}{1 - P_E} \;=\; \frac{0.85 - 0.51}{1 - 0.51} \;=\; \frac{0.34}{0.49} \;\approx\; 0.694 \end{eqnarray*}\]

Verificación con un paquete (opcional)

# Verificación con irr::kappa2 (necesita datos caso a caso)
# Generamos etiquetas sintéticas consistentes con a,b,c,d
library(irr)

A <- c(rep("Pos", a), rep("Pos", b), rep("Neg", c), rep("Neg", d))
B <- c(rep("Pos", a), rep("Neg", b), rep("Pos", c), rep("Neg", d))

kappa2(data.frame(A, B))  # Muestra κ y su IC

##  Cohen's Kappa for 2 Raters (Weights: unweighted)
## 
##  Subjects = 100 
##    Raters = 2 
##     Kappa = 0.694 
## 
##         z = 6.98 
##   p-value = 3.05e-12

Interpretación (opcional)

• Resultado: \(\kappa \approx 0.694\).

• Escala de interpretación: \(0.60 < \kappa \leq 0.80\). Entonces, es un valor sustancial.

• Conclusión práctica:

  • Hay un acuerdo sustancial entre los evaluadores A y B más allá del azar.

  • Esto sugiere que el procedimiento de evaluación es fiable: cuando A clasifica un caso, B tiende a coincidir en un grado alto.

Comentario adicional

• A diferencia de la exactitud \((0.85)\), \(\kappa\) ajusta por el acuerdo esperado al azar \((P_E = 0.51)\).

• Por eso, \(\kappa\) es una medida más justa del acuerdo real entre evaluadores o entre un modelo y un experto.

25 Ejemplo 5: enunciado

Considere los datos del archivo chdage. Suponga que se quiere analizar un modelo de regresión logística, considerando a chd como variable dependiente y age como independiente.

1. Calcule tanto las probabilidades predichas como los valores predichos de la variable de respuesta y reúnalos en el data frame chdage.

2. Teniendo en cuenta el data frame de (a), construya la matriz de confusión correspondiente.

3. Determine la exactitud del modelo.

4. A partir de la matriz de confusión, halle alguno de los estadísticos relacionados (exactitud, sensitividad, especificidad, etc.).

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

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

7. Interprete los resultados obtenidos.

26 Ejemplo 5: Solución

26.0.1 Solución inciso (a)

Abajo se muestra una parte del data frame solicitado:

Con el paquete glm

#1. Con la función glm:

df <- aplore3::chdage
modelo <- glm(chd~age, family=binomial(link = "logit"), data=df)

df$Prob_pred <- round(predict(modelo, newdata = df, "response"),4)
df$CHD_pred <- ifelse(df$Prob_pred > 0.5, "Yes", "No")

y_pred <- as.factor(df$CHD_pred)
y_obs <- as.factor(df$chd)
df

Con el paquete lsm

#2. Con la función lsm:

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

datos$Prob_pred <- round(predict(modelo, type = "response")[,1],4)
datos$CHD_pred <- ifelse(datos$Prob_pred > 0.5, "Yes", "No")

# Observado y predicho como factores DENTRO de 'datos'
datos$y_obs  <- factor(ifelse(datos$CHD == 1, "Yes", "No"), levels = c("No","Yes"))
datos$y_pred <- factor(datos$CHD_pred, levels = c("No","Yes"))

Table 26.1: Data frame con los predichos.

ID	AGE	CHD	Prob_pred	CHD_pred	y_obs	y_pred
1	20	0	0.0435	No	No	No
2	23	0	0.0596	No	No	No
3	24	0	0.0662	No	No	No
4	25	0	0.0733	No	No	No
5	25	1	0.0733	No	Yes	No
:	:	:	:	:
:	:	:	:	:
96	63	1	0.8427	Yes	Yes	Yes
97	64	0	0.8569	Yes	No	Yes
98	64	1	0.8569	Yes	Yes	Yes
99	65	1	0.8699	Yes	Yes	Yes
100	69	1	0.9125	Yes	Yes	Yes

26.0.2 Solución inciso (b)

La matriz de confusión correspondiente es:

y_obs <- datos$y_obs
y_pred<- datos$y_pred
Tabla <-  table(y_obs, y_pred)
Tabla

Table 26.2: Matriz de confusión.

	1	2	3
	y_pred=No	y_pred=Yes	Total
y_obs=No	45	12	57
y_obs=Yes	14	29	43
Total	59	41	100

26.0.3 Solución inciso (c)

La exactitud es del 74% Y ignifica que el modelo acierta 74 de cada 100 casos (y falla 26). Se puede calcular de dos maneras:

# Con el data frame de la parte (a):
E <- mean(y_pred == y_obs)*100

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

## [1] 74

26.0.4 Solución inciso (d)

En la tabla de abajo se presentan y calculan las métricas más usadas a partir de la matriz de confusión (filas = observados y_obs, columnas = predichos y_pred):

• TP (True Positives / Verdaderos Positivos): casos Yes correctamente predichos como Yes.

• TN (True Negatives / Verdaderos Negativos): casos No correctamente predichos como No.

• FP (False Positives / Falsos Positivos): casos No predichos como Yes.

• FN (False Negatives / Falsos Negativos): casos Yes predichos como No.

Las métricas que se calculan son:

• Accuracy (Exactitud): \((TP+TN)/\text{Total}\). Proporción de aciertos.

• Sensitivity (Recall, TPR): \(TP/(TP+FN)\). Cobertura de la clase positiva.

• Specificity (TNR): \(TN/(TN+FP)\). Cobertura de la clase negativa.

• Precision (PPV): \(TP/(TP+FP)\). Precisión sobre lo predicho como positivo.

• F1-score: media armónica entre Precision y Sensitivity.

• Baseline (accuracy nula): proporción de la clase mayoritaria observada (referencia de “predecir siempre la clase más frecuente”).

• Accuracy 95% CI: intervalo de confianza del 95% para la exactitud (ensayo binomial).

# Descomposición clásica TP, TN, FP, FN
TN <- Tabla["No","No"];  FP <- Tabla["No","Yes"]
FN <- Tabla["Yes","No"]; TP <- Tabla["Yes","Yes"]

# Métricas
accuracy <- (TP + TN) / sum(Tabla)      # Exactitud
sens     <- TP / (TP + FN)              # Sensibilidad / Recall / TPR
spec     <- TN / (TN + FP)              # Especificidad / TNR
prec     <- TP / (TP + FP)              # Precisión / PPV
f1       <- 2 * prec * sens / (prec + sens)

# Baseline (clase mayoritaria observada)
baseline <- max(prop.table(rowSums(Tabla)))

# IC 95% para la exactitud (binomial)
ci_acc <- prop.test(sum(diag(Tabla)), sum(Tabla))$conf.int

# Tabla de resultados (sin interpretación)
res <- data.frame(
  Metric = c("Accuracy",
             "Baseline (majority class)",
             "Sensitivity (Recall, TPR)",
             "Specificity (TNR)",
             "Precision (PPV)",
             "F1-score",
             "Accuracy 95% CI [low, high]"),
  Value  = c(sprintf("%.4f (%.1f%%)", accuracy, 100*accuracy),
             sprintf("%.4f (%.1f%%)", baseline, 100*baseline),
             sprintf("%.4f (%.1f%%)", sens, 100*sens),
             sprintf("%.4f (%.1f%%)", spec, 100*spec),
             sprintf("%.4f (%.1f%%)", prec, 100*prec),
             sprintf("%.4f", f1),
             sprintf("[%.4f, %.4f]  ([%.1f%%, %.1f%%])",
                     ci_acc[1], ci_acc[2], 100*ci_acc[1], 100*ci_acc[2]))
)

kable(res, 
      align="lc",  
      col.names=c("Métrica", "Valor"), 
      caption = "Resultados (basado en una matriz de confusión)." 
      ) 

Table 26.3: Resultados (basado en una matriz de confusión).

Métrica	Valor
Accuracy	0.7400 (74.0%)
Baseline (majority class)	0.5700 (57.0%)
Sensitivity (Recall, TPR)	0.6744 (67.4%)
Specificity (TNR)	0.7895 (78.9%)
Precision (PPV)	0.7073 (70.7%)
F1-score	0.6905
Accuracy 95% CI [low, high]	[0.6410, 0.8203] ([64.1%, 82.0%])

26.0.5 Solución inciso (e)

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

library(caret)
confusionMatrix(data=y_pred, reference = y_obs,  positive="Yes")

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

26.0.6 Solución inciso (f)

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

library(gmodels)
CrossTable(y_pred, y_obs)

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

26.0.7 Solución inciso (g)

• Contexto de la matriz: Con 100 casos, el modelo acierta 74 (45 TN + 29 TP) y se equivoca 26 (12 FP + 14 FN). La clase positiva está definida como Yes.

• Exactitud = 0.74 (IC95%: 0.6427, 0.8226).
  El modelo clasifica correctamente aproximadamente 74% de los casos; el intervalo sugiere que, en población, la exactitud probable está entre 64% y 82%.

• NIR = 0.57; p-value (Acc > NIR) = 0.0003187. La exactitud es significativamente mayor que la tasa de la clase mayoritaria (57%). El modelo supera la línea base.

• Kappa de Cohen = 0.4666 (moderado).
  Indica acuerdo moderado más allá del azar entre predicción y realidad (ajusta la exactitud por el acuerdo esperado al azar).

• Sensibilidad (Recall, TPR) = 0.674.
  El modelo detecta 67.4% de los positivos reales (Yes). La tasa de falsos negativos es 1−0.674 ≈ 0.326.

• Especificidad (TNR) = 0.7895.
  El modelo reconoce 78.95% de los negativos reales (No). La tasa de falsos positivos es 1−0.7895 ≈ 0.2105.

• Valor predictivo positivo (PPV) = 0.7073.
  Cuando el modelo predice Yes, acierta aproximadamente 71% de las veces.

• Valor predictivo negativo (NPV) = 0.7632.
  Cuando el modelo predice No, acierta aproximadamente 76% de las veces.

• Balanced Accuracy = 0.7319.
  Promedio de sensibilidad y especificidad; útil con clases desbalanceadas. Coherente con la exactitud observada.

• Prevalence = 0.43 vs Detection Prevalence = 0.41.
  En los datos, 43% son positivos reales; el modelo predice positivos en 41% de los casos. Por lo tanto, ligera sub-predicción de la clase positiva.

• McNemar’s test p-value = 0.8445.
  No hay evidencia de asimetría sistemática entre errores tipo FP y FN (los errores no favorecen consistentemente a un lado).

En resumen, hay un desempeño superior al baseline, acuerdo moderado más allá del azar y buen balance entre sensibilidad y especificidad; los valores predictivos sugieren predicciones relativamente confiables para ambas clases.

27 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, estimaciones e intervalos de confianza para \(p_j\), ODDS, razones ODDS, intercepto, pendiente, etc.

27.0.1 Ejercicios 1 a 3

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

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

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

27.0.2 Ejercicio 4

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

2. En estos datos tome a la variable AGE como independiente y STA como dependiente.

• Escriba la ecuación general para el modelo de regresión logística de STA contra AGE y para el riesgo estimado por este modelo.

• ¿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?

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

3. Partiendo del modelo original:

• Escriba una expresión para la función de verosimilitud y del logaritmo de esta función para el modelo de regresión logístico de (a) usando los 200 datos no agrupados.

• Obtenga una expresión para las dos ecuaciones de verosimilitud.

4. Partiendo del modelo original:

• Obtenga las estimaciones de los parámetros del modelo de regresión logístico de (a).

• Usando estas estimaciones, escriba las correspondientes ecuaciones para los valores ajustados.

• Grafique la ecuación para los valores ajustados utilizando los mismos ejes como en (b) y (c).

5. Resuma (describa en palabras) los resultados presentados en la gráfica obtenida en (a), (b) y (d).

6. Usando los resultados de la parte (d), verifique la significancia del coeficiente de AGE. ¿Qué supuestos se necesitan para realizar dicha prueba?

7. Usando los resultados de (d), halle un intervalo del 95% de confianza para la pendiente y la constante. Escriba una interpretación con respecto al intervalo encontrado para la pendiente.

8. Considere el modelo en (e).

• Obtenga la matriz de covarianzas estimada para el modelo en (e).

• Calcule el logit y la probabilidad logística estimada para una persona de 60 años.

• Calcule un intervalo del 95% de confianza para el logit.

• Calcule un intervalo del 95% de confianza para la probabilidad logística estimada. Interprete la probabilidad estimada y su intervalo de confianza.

1. Considere el modelo en (e).

• Obtenga el logit estimado y su error estándar para cada persona en el estudio ICU.

• Grafique el logit estimado y los límites del intervalo del 95 % de confianza versus AGE para cada persona.

• Explique (en palabras) similaridades y diferencias entre las apariencias de esta gráfica y una gráfica de una gráfica de un modelo de regresión ajustado y sus límites del intervalo del 95 % de confianza.

10. Realize las siguientes pruebas de comparación de modelos resumiendo en una tabla las pruebas realizadas, el valor y la distribución muestral del estadístico de prueba, los grados de libertad, el P-valor y su decisión.

• Nulo vs Logístico.

• Logístico vs Completo.

• Logístico vs Saturado

27.0.3 Ejercicios 5 a 7

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

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

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

27.0.4 Ejercicios 8 a 9

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

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

27.0.5 Ejercicio 10

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

1. Responda:

• ¿Cuál es la ecuación para el modelo de regresión logística?

• ¿Cuál es la ecuación para riesgo estimado por este modelo?

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

• \(\mathcal{L}(\widehat{\alpha})\) en el modelo logístico.

3. Construya intervalos de confianza del 95 % de confianza para los siguientes parámetros e interprételos (justifique en forma clara y precisa todas sus afirmaciones):

• La pendiente \(\beta\). ¿Es apropiado el modelo?

• El intercepto \(\delta\). ¿Pasa la curva de regresión logística por el origen?

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

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

• El odds cuando PSA=11.2.

• La razón odds OR. ¿Es PSA estadísticamente significativa en el modelo?

4. Realize las siguientes pruebas de comparación de modelos resumiendo en una tabla las pruebas realizadas, el valor y la distribución muestral del estadístico de prueba, los grados de libertad, el P-valor y su decisión.

• Nulo vs Logístico.

• Logístico vs Completo.

• Logístico vs Saturado

27.0.6 Ejercicios 11 a 13

11. Considere los datos PROS. Realize un análisis de regresión logística tomando a CAPSULE como variable dependiente y VOL como variable independiente.

12. Considere los datos PROS. Realize un análisis de regresión logística tomando a CAPSULE como variable dependiente y AGE como variable independiente.

13. Considere los datos PROS, tomando a CAPSULE como variable dependiente y AGE, PSA y VOL como variables independientes.

27.0.7 Ejercicio 14

14. 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. Responda:

• ¿Cuál es la ecuación para el modelo de regresión logística?

• ¿Cuál es la ecuación para riesgo estimado por este modelo?

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

• \(\mathcal{L}(\widehat{\alpha})\) en el modelo logístico.

3. Construya intervalos de confianza del 95 % de confianza para los siguientes parámetros e interprételos (justifique en forma clara y precisa todas sus afirmaciones):

• La pendiente \(\beta\). ¿Es apropiado el modelo?

• El intercepto \(\delta\). ¿Pasa la curva de regresión logística por el origen?

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

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

• El odds cuando LWT=100.3.

• La razón odds OR. ¿Es LWT estadísticamente significativa en el modelo?

4. Realize las siguientes pruebas de comparación de modelos resumiendo en una tabla las pruebas realizadas, el valor y la distribución muestral del estadístico de prueba, los grados de libertad, el P-valor y su decisión.

• Nulo vs Logístico.

• Logístico vs Completo.

• Logístico vs Saturado

27.0.8 Ejercicios 15 a 18

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

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

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

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

27.0.9 Ejercicios 19 a 21

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

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

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

27.0.10 Ejercicios 22

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

27.0.11 Ejercicios 23 a 26

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

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

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

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