Nota pedagógica

Este documento implementa computacionalmente, con detalle algebraico completo, los capítulos del texto de Kleinbaum & Klein. El énfasis central es en el álgebra matricial que subyace al vector de parámetros \(\hat{\boldsymbol{\beta}}\): desde su definición, pasando por la verosimilitud, la construcción de la matriz de información de Fisher, el algoritmo IRLS, hasta la varianza estimada y las pruebas de hipótesis. Cada resultado numérico de R tiene su correspondiente derivación matemática explícita.

1 La Regresión Logística

1.1 ¿Por qué necesitamos un modelo?

La pregunta fundamental en epidemiología es si una exposición \(E\) causa una enfermedad \(D\). Sin embargo, la asociación cruda \(E \to D\) rara vez es suficiente: la presencia de variables de confusión y posibles modificadores de efecto exigen un abordaje multivarido.

La regresión lineal ordinaria fallaría aquí porque modelar una probabilidad binaria \((0, 1)\) con una recta produciría predicciones fuera del rango \([0,1]\) y violaría los supuestos de homocedasticidad. La solución es aplicar una transformación que lleve la escala de las probabilidades al espacio de los reales.

Pregunta central: ¿Cuál es la probabilidad de que un individuo con perfil de exposición \(\mathbf{X} = (X_1, X_2, \ldots, X_k)\) desarrolle la enfermedad?

\[P(D = 1 \mid X_1, X_2, \ldots, X_k) = ?\]

1.2 La función logística — motivación y propiedades

1.2.1 Definición matemática

La función logística (o sigmoidea) es:

\[\boxed{f(z) = \frac{1}{1 + e^{-z}} = \frac{e^z}{1 + e^z}, \quad z \in \mathbb{R}}\]

El argumento \(z\) es el índice lineal de riesgo: una combinación lineal de los predictores \(z = \alpha + \beta_1 X_1 + \cdots + \beta_k X_k\).

Por qué esta función transforma lineales en probabilidades:

Cuando \(z \to -\infty\): \(e^{-z} \to +\infty\), entonces \(f(z) \to 0\)
Cuando \(z \to +\infty\): \(e^{-z} \to 0\), entonces \(f(z) \to 1\)
Siempre: \(0 < f(z) < 1\) para todo \(z\) finito
En \(z = 0\): \(f(0) = \frac{1}{1 + 1} = 0.5\) (punto de inflexión)
La función es monótonamente creciente: a mayor \(z\), mayor probabilidad

1.2.2 Derivada y curvatura

La derivada de \(f(z)\) tiene una forma elegante:

\[f'(z) = \frac{d}{dz} f(z) = f(z)\,[1 - f(z)]\]

Esta propiedad es fundamental para el algoritmo de estimación (IRLS), pues los pesos \(w_i = p_i(1 - p_i)\) provienen directamente de ella.

z   <- seq(-7, 7, length.out = 600)
fz  <- 1 / (1 + exp(-z))
dfz <- fz * (1 - fz)

df_plot <- tibble(z = z, fz = fz, dfz = dfz)

p1 <- ggplot(df_plot, aes(x = z, y = fz)) +
  geom_ribbon(aes(ymin = 0, ymax = fz), fill = COL_AZUL, alpha = 0.10) +
  geom_line(color = COL_AZUL, linewidth = 1.7) +
  geom_hline(yintercept = c(0, 0.5, 1), linetype = "dashed",
             color = COL_GRIS, linewidth = 0.4) +
  geom_vline(xintercept = 0, linetype = "dashed", color = COL_GRIS, linewidth = 0.4) +
  annotate("point", x = 0, y = 0.5, size = 4.5, color = COL_NARANJA) +
  annotate("text",  x = 1.4, y = 0.45, label = "f(0) = 0.5",
           size = 3.7, color = COL_NARANJA, fontface = "bold") +
  annotate("text", x = -5.5, y = 0.06, label = "f(z) → 0",
           color = COL_ROJO, size = 3.4, fontface = "italic") +
  annotate("text", x =  5.0, y = 0.94, label = "f(z) → 1",
           color = COL_VERDE, size = 3.4, fontface = "italic") +
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1)) +
  labs(title    = "Función Logística  f(z) = 1 / (1 + e⁻ᶻ)",
       subtitle = "Probabilidad como función del índice lineal de riesgo z",
       x = "z = α + β₁X₁ + ⋯ + βₖXₖ",
       y = "P(D = 1 | X) = f(z)",
       caption = "Fuente: Kleinbaum & Klein, Capítulo 1") +
  tema_doc

p2 <- ggplot(df_plot, aes(x = z, y = dfz)) +
  geom_ribbon(aes(ymin = 0, ymax = dfz), fill = COL_NARANJA, alpha = 0.12) +
  geom_line(color = COL_NARANJA, linewidth = 1.5) +
  geom_vline(xintercept = 0, linetype = "dashed", color = COL_GRIS, linewidth = 0.4) +
  annotate("point", x = 0, y = 0.25, size = 4.5, color = COL_ROJO) +
  annotate("text", x = 1.8, y = 0.23, label = "máx = 0.25\nen z = 0",
           size = 3.4, color = COL_ROJO) +
  scale_y_continuous(breaks = seq(0, 0.25, 0.05)) +
  labs(title    = "Derivada  f'(z) = f(z)·[1 − f(z)]",
       subtitle = "Pesos de la matriz W en el algoritmo IRLS",
       x = "z", y = "f'(z)") +
  tema_doc

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

Figura 1. Función logística y su derivada. La derivada máxima en z = 0 indica que el riesgo cambia más rápidamente alrededor del umbral p = 0.5.

1.2.3 La triada fundamental: probabilidad, odds y logit

El modelo logístico tiene tres representaciones equivalentes que el epidemiólogo debe dominar:

\[\underbrace{P(\mathbf{X})}_{\substack{\text{probabilidad}\\\text{(riesgo)}\\\in(0,1)}} \;\longleftrightarrow\; \underbrace{\frac{P(\mathbf{X})}{1 - P(\mathbf{X})}}_{\substack{\text{odds}\\\\\in(0,+\infty)}} \;\longleftrightarrow\; \underbrace{\ln\!\left(\frac{P(\mathbf{X})}{1-P(\mathbf{X})}\right)}_{\substack{\text{logit = log-odds}\\\\\in(-\infty,+\infty)}}\]

Lectura epidemiológica: Los odds expresan “cuántas veces más probable es el evento que su no ocurrencia”. Si \(P = 0.25\), el odds = \(0.25/0.75 = 1/3\) (por cada caso hay 3 no-casos). El logit simplemente toma el logaritmo natural del odds, convirtiendo la escala \((0, +\infty)\) en \((-\infty, +\infty)\), donde operan las regresiones lineales.

z_t  <- c(-5, -3, -2, -1, -0.5, 0, 0.5, 1, 2, 3, 5)
p_t  <- 1 / (1 + exp(-z_t))
o_t  <- p_t / (1 - p_t)
l_t  <- log(o_t)

tibble(
  `z (índice lineal)` = z_t,
  `P(X) = f(z)`       = round(p_t, 4),
  `Odds = P/(1−P)`    = round(o_t, 4),
  `Logit = ln(odds)`  = round(l_t, 4),
  `Interpretación`    = c(
    "Riesgo muy bajo (0.7%)",
    "Riesgo bajo (4.7%)",
    "Riesgo bajo (12%)",
    "Riesgo moderado-bajo (27%)",
    "Límite inferior neutro (38%)",
    "Punto de inflexión (50%)",
    "Límite superior neutro (62%)",
    "Riesgo moderado-alto (73%)",
    "Riesgo alto (88%)",
    "Riesgo muy alto (95%)",
    "Riesgo extremo (99.3%)"
  )
) |> tabla_doc(titulo = "Tabla 1.1 — Triada: z, probabilidad, odds y logit")

Tabla 1.1 — Triada: z, probabilidad, odds y logit
z (índice lineal)	P(X) = f(z)	Odds = P/(1−P)	Logit = ln(odds)	Interpretación
-5.0	0.0067	0.0067	-5.0	Riesgo muy bajo (0.7%)
-3.0	0.0474	0.0498	-3.0	Riesgo bajo (4.7%)
-2.0	0.1192	0.1353	-2.0	Riesgo bajo (12%)
-1.0	0.2689	0.3679	-1.0	Riesgo moderado-bajo (27%)
-0.5	0.3775	0.6065	-0.5	Límite inferior neutro (38%)
0.0	0.5000	1.0000	0.0	Punto de inflexión (50%)
0.5	0.6225	1.6487	0.5	Límite superior neutro (62%)
1.0	0.7311	2.7183	1.0	Riesgo moderado-alto (73%)
2.0	0.8808	7.3891	2.0	Riesgo alto (88%)
3.0	0.9526	20.0855	3.0	Riesgo muy alto (95%)
5.0	0.9933	148.4132	5.0	Riesgo extremo (99.3%)

1.3 El modelo logístico — forma canónica

1.3.1 Definición formal y notación vectorial

El modelo logístico múltiple para \(k\) covariables es:

\[\boxed{P(\mathbf{X}) = P(D = 1 \mid X_1, \ldots, X_k) = \frac{1}{1 + \exp\!\left[-\underbrace{\left(\alpha + \sum_{i=1}^{k} \beta_i X_i\right)}_{\text{índice lineal } z}\right]}}\]

Definiéndolo en notación vectorial:

\[P(\mathbf{X}) = \sigma(\boldsymbol{\beta}^{\top} \mathbf{x}) = \frac{1}{1 + e^{-\boldsymbol{\beta}^{\top}\mathbf{x}}}\]

donde:

\(\mathbf{x} = (1, X_1, X_2, \ldots, X_k)^{\top} \in \mathbb{R}^{k+1}\) — vector de covariables del individuo, con un 1 inicial para el intercepto
\(\boldsymbol{\beta} = (\alpha, \beta_1, \beta_2, \ldots, \beta_k)^{\top} \in \mathbb{R}^{k+1}\) — vector de parámetros a estimar
\(\sigma(\cdot)\) — función logística (sigmoidea)

En forma logit (lineal en los parámetros):

\[\text{logit}[P(\mathbf{X})] = \ln\!\left(\frac{P(\mathbf{X})}{1-P(\mathbf{X})}\right) = \alpha + \sum_{i=1}^{k} \beta_i X_i = \boldsymbol{\beta}^{\top} \mathbf{x}\]

Propiedad clave: El logit es lineal en los parámetros \(\boldsymbol{\beta}\). Aunque la probabilidad \(P(\mathbf{X})\) es una función no lineal (sigmoide) de \(\boldsymbol{\beta}\), el log-odds es exactamente \(\boldsymbol{\beta}^{\top}\mathbf{x}\). Esta linealidad es lo que hace posible la estimación por máxima verosimilitud con métodos numéricos eficientes.

1.3.2 Interpretación de los coeficientes \(\beta_i\)

\[\beta_i = \frac{\partial \,\text{logit}[P(\mathbf{X})]}{\partial X_i}\]

Cada \(\beta_i\) representa el cambio en el log-odds por cada unidad de aumento en \(X_i\), manteniendo las demás variables constantes (ceteris paribus). En la escala de odds:

\[e^{\beta_i} = \text{OR asociado a un incremento de 1 unidad en } X_i\]

1.3.3 Notación matricial del modelo para \(n\) individuos

Para un conjunto de datos con \(n\) individuos y \(k\) covariables, el modelo completo se escribe de forma compacta como:

\[\boldsymbol{\eta} = \mathbf{X}\boldsymbol{\beta}\]

donde:

Símbolo	Nombre	Dimensión	Descripción
\(\mathbf{X}\)	Matriz de diseño	\(n \times (k+1)\)	Columna 1 de unos + \(k\) covariables
\(\boldsymbol{\beta}\)	Vector de parámetros	\((k+1) \times 1\)	Intercepto \(\alpha\) + \(k\) coeficientes \(\beta_i\)
\(\boldsymbol{\eta}\)	Predictor lineal	\(n \times 1\)	Logit de cada individuo: \(\eta_i = \mathbf{x}_i^{\top}\boldsymbol{\beta}\)
\(\mathbf{p}\)	Vector de probabilidades	\(n \times 1\)	\(p_i = \sigma(\eta_i) = 1/(1+e^{-\eta_i})\)

La matriz de diseño tiene la siguiente estructura:

\[\mathbf{X} = \begin{pmatrix} 1 & X_{11} & X_{12} & \cdots & X_{1k} \\ 1 & X_{21} & X_{22} & \cdots & X_{2k} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & X_{n1} & X_{n2} & \cdots & X_{nk} \end{pmatrix}, \quad \boldsymbol{\beta} = \begin{pmatrix} \alpha \\ \beta_1 \\ \beta_2 \\ \vdots \\ \beta_k \end{pmatrix}\]

# ══════════════════════════════════════════════════════════════════
#  CONSTRUCCIÓN DE LA MATRIZ DE DISEÑO X
# ══════════════════════════════════════════════════════════════════
#
# model.matrix(~ ., data = df) construye automáticamente:
#   - Columna de 1s (intercepto)
#   - Las k covariables especificadas
#
# Internamente R usa esto para todo glm(): el usuario nunca lo ve
# directamente, pero es la base de TODA la estimación.

construir_X <- function(df) {
  model.matrix(~ ., data = df)
}

# ── Datos de demostración (Kleinbaum & Klein, Cap.1 CHD study) ────
demo_chd <- data.frame(
  CAT = c(1, 0, 1, 0, 1, 0, 1, 0),   # catecholamine high = 1
  AGE = c(45, 52, 61, 38, 55, 47, 63, 41),
  ECG = c(1, 0, 1, 0, 0, 1, 1, 0)    # ECG anomaly = 1
)

X_demo <- construir_X(demo_chd)

cat("══ Matriz de diseño X para el modelo CHD (n=8) ══\n\n")

#> ══ Matriz de diseño X para el modelo CHD (n=8) ══

cat("  Estructura: X = [1 | CAT | AGE | ECG]\n")

#>   Estructura: X = [1 | CAT | AGE | ECG]

cat("  Dimensión : n × (k+1) = 8 × 4\n\n")

#>   Dimensión : n × (k+1) = 8 × 4

print(X_demo)

#>   (Intercept) CAT AGE ECG
#> 1           1   1  45   1
#> 2           1   0  52   0
#> 3           1   1  61   1
#> 4           1   0  38   0
#> 5           1   1  55   0
#> 6           1   0  47   1
#> 7           1   1  63   1
#> 8           1   0  41   0
#> attr(,"assign")
#> [1] 0 1 2 3

cat("\n── Verificación de dimensiones ──\n")

#> 
#> ── Verificación de dimensiones ──

cat(sprintf("  nrow(X) = %d  (número de individuos)\n", nrow(X_demo)))

#>   nrow(X) = 8  (número de individuos)

cat(sprintf("  ncol(X) = %d  (1 intercepto + %d covariables)\n",
            ncol(X_demo), ncol(X_demo) - 1))

#>   ncol(X) = 4  (1 intercepto + 3 covariables)

¿Por qué el 1 en la primera columna?

El modelo es \(\text{logit}[P] = \alpha \cdot 1 + \beta_1 X_1 + \cdots + \beta_k X_k\). Al escribirlo como \(\boldsymbol{\beta}^{\top}\mathbf{x}\), necesitamos que \(\mathbf{x}\) tenga un componente constante de valor 1 que “multiplique” al intercepto \(\alpha\). Sin ese 1, no podríamos escribir el producto escalar de forma uniforme para todos los parámetros.

# ══════════════════════════════════════════════════════════════════
#  PREDICCIÓN: de β a probabilidades vía η = Xβ → p = σ(η)
# ══════════════════════════════════════════════════════════════════
#
# Pasos:
#   1.  η = X β          [producto matricial: n×(k+1) × (k+1)×1 = n×1]
#   2.  p = 1/(1+exp(-η)) [aplicar σ elemento a elemento]
#

# Parámetros de Kleinbaum & Klein, Capítulo 1
beta_chd <- c(alpha = -3.911, b_CAT = 0.652, b_AGE = 0.029, b_ECG = 0.342)

cat("── Vector de parámetros β̂ ──\n")

#> ── Vector de parámetros β̂ ──

cat("  (intercepto, CAT, AGE, ECG)\n")

#>   (intercepto, CAT, AGE, ECG)

print(beta_chd)

#>  alpha  b_CAT  b_AGE  b_ECG 
#> -3.911  0.652  0.029  0.342

# Paso 1: predictor lineal η = Xβ
eta_demo <- X_demo %*% beta_chd
cat("\n── Paso 1: η = Xβ (log-odds de cada individuo) ──\n")

#> 
#> ── Paso 1: η = Xβ (log-odds de cada individuo) ──

print(round(eta_demo, 4))

#>     [,1]
#> 1 -1.612
#> 2 -2.403
#> 3 -1.148
#> 4 -2.809
#> 5 -1.664
#> 6 -2.206
#> 7 -1.090
#> 8 -2.722

# Paso 2: probabilidades p = σ(η)
p_demo <- 1 / (1 + exp(-eta_demo))
cat("\n── Paso 2: p = σ(η) = 1/(1+e^{-η}) (probabilidades) ──\n")

#> 
#> ── Paso 2: p = σ(η) = 1/(1+e^{-η}) (probabilidades) ──

print(round(p_demo, 4))

#>     [,1]
#> 1 0.1663
#> 2 0.0829
#> 3 0.2409
#> 4 0.0568
#> 5 0.1592
#> 6 0.0992
#> 7 0.2516
#> 8 0.0617

# Tabla completa
tibble(
  Individuo = 1:8,
  CAT = demo_chd$CAT,
  AGE = demo_chd$AGE,
  ECG = demo_chd$ECG,
  `η = logit P` = round(as.vector(eta_demo), 4),
  `P̂(X) = σ(η)` = round(as.vector(p_demo), 4),
  `Odds = P/(1-P)` = round(as.vector(p_demo/(1-p_demo)), 4)
) |> tabla_doc(titulo = "Tabla 1.2 — Predictor lineal y probabilidades estimadas para 8 individuos (cálculo matricial)")

Tabla 1.2 — Predictor lineal y probabilidades estimadas para 8 individuos (cálculo matricial)
Individuo	CAT	AGE	ECG	η = logit P	P̂(X) = σ(η	Odds = P/(1-P
1	1	45	1	-1.612	0.1663	0.1995
2	0	52	0	-2.403	0.0829	0.0904
3	1	61	1	-1.148	0.2409	0.3173
4	0	38	0	-2.809	0.0568	0.0603
5	1	55	0	-1.664	0.1592	0.1894
6	0	47	1	-2.206	0.0992	0.1101
7	1	63	1	-1.090	0.2516	0.3362
8	0	41	0	-2.722	0.0617	0.0657

1.4 Diseños de estudio y estimabilidad

tibble(
  `Diseño`              = c("Cohorte / Seguimiento",
                             "Casos y Controles",
                             "Transversal"),
  `Estima α`            = c("✓ Sí (tasa basal real)", "✗ No (α sesgado)", "✗ No"),
  `P̂(X) válida`         = c("✓ Sí", "✗ No", "✗ No"),
  `RR directo`          = c("✓ Sí", "✗ No", "✗ No"),
  `OR estimable`        = c("✓ Sí", "✓ Sí", "✓ Sí (prevalencia OR)"),
  `Vector β̂ completo`   = c("✓ Todos los β", "Solo β₁…βₖ", "Solo β₁…βₖ"),
  check.names = FALSE
) |>
  tabla_doc(titulo = "Tabla 1.3 — Estimabilidad del modelo logístico según diseño de estudio") |>
  column_spec(1, bold = TRUE)

Tabla 1.3 — Estimabilidad del modelo logístico según diseño de estudio
Diseño	Estima α	P̂(X) válida	\|RR directo	\|OR estimable	\|Vector β̂ complet	\|check.name
Cohorte / Seguimiento	✓ Sí (tasa basal real)	✓ Sí	✓ Sí	✓ Sí	✓ Todos los β	FALSE
Casos y Controles	✗ No (α sesgado)	✗ No	✗ No	✓ Sí	Solo β₁…βₖ	FALSE
Transversal	✗ No	✗ No	✗ No	✓ Sí (prevalencia OR)	Solo β₁…βₖ	FALSE

Nota crítica — casos y controles: En este diseño, el intercepto \(\alpha\) absorbe la fracción de muestreo de los controles y no es interpretable como tasa basal de enfermedad en la población. Solo los \(\beta_i\) (y por ende los OR = \(e^{\beta_i}\)) son estimables de forma válida. La razón es que el OR no depende de \(\alpha\) (el intercepto se cancela exactamente en la fórmula del OR, como se demuestra en el Capítulo 3).

1.5 Razón de Odds — Derivación algebraica completa

Para dos perfiles \(\mathbf{X}^{(1)}\) y \(\mathbf{X}^{(0)}\):

\[\text{OR} = \frac{P(\mathbf{X}^{(1)}) / [1 - P(\mathbf{X}^{(1)})]}{P(\mathbf{X}^{(0)}) / [1 - P(\mathbf{X}^{(0)})]}\]

Sustituyendo las expresiones logísticas. Para el numerador:

\[\frac{P(\mathbf{X}^{(1)})}{1 - P(\mathbf{X}^{(1)})} = \frac{\frac{1}{1+e^{-\boldsymbol{\beta}^{\top}\mathbf{x}^{(1)}}}}{\frac{e^{-\boldsymbol{\beta}^{\top}\mathbf{x}^{(1)}}}{1+e^{-\boldsymbol{\beta}^{\top}\mathbf{x}^{(1)}}}} = e^{\boldsymbol{\beta}^{\top}\mathbf{x}^{(1)}} = e^{\alpha + \sum_i \beta_i X_i^{(1)}}\]

Dividiendo numerador entre denominador:

\[\text{OR} = \frac{e^{\alpha + \sum_i \beta_i X_i^{(1)}}}{e^{\alpha + \sum_i \beta_i X_i^{(0)}}} = e^{\sum_i \beta_i (X_i^{(1)} - X_i^{(0)})}\]

\[\boxed{\text{OR} = \exp\!\left[\boldsymbol{\beta}^{\top}(\mathbf{x}^{(1)} - \mathbf{x}^{(0)})\right] = \exp\!\left[\sum_{i=1}^{k} \beta_i (X_i^{(1)} - X_i^{(0)})\right] = \prod_{i=1}^{k} e^{\beta_i (X_i^{(1)} - X_i^{(0)})}}\]

El intercepto \(\alpha\) se cancela exactamente: el OR depende solo de los \(\beta_i\) de las covariables, no del nivel basal de enfermedad. Este es el resultado algebraico que valida el uso del modelo logístico en estudios de casos y controles.

# ══════════════════════════════════════════════════════════════════
#  SUBRUTINA: Cálculo matricial del OR con trazabilidad completa
# ══════════════════════════════════════════════════════════════════
#
# Implementa OR = exp[β'(x1 - x0)] mostrando cada paso del álgebra
#
calc_OR_detallado <- function(x1, x0, beta_completo) {
  # Extraer solo los βi (sin intercepto)
  b <- beta_completo[-1]

  cat("═══════════════════════════════════════════════════\n")
  cat(" CÁLCULO DEL OR POR PRODUCTO ESCALAR MATRICIAL\n")
  cat("═══════════════════════════════════════════════════\n\n")

  cat("Perfil 1 (x1) — expuesto:     ", x1, "\n")
  cat("Perfil 0 (x0) — no expuesto:  ", x0, "\n\n")

  # Paso 1: diferencia de perfiles
  delta <- x1 - x0
  cat("Paso 1 — Δx = x1 - x0 (diferencia de perfiles):\n")
  nombres <- names(b)
  df_delta <- data.frame(Variable = nombres, `x1` = x1, `x0` = x0,
                          `Δx = x1-x0` = delta, check.names = FALSE)
  print(df_delta)

  # Paso 2: producto elemento a elemento β_i * Δx_i
  contrib <- b * delta
  cat("\nPaso 2 — Contribución de cada variable: β_i × Δx_i\n")
  df_contrib <- data.frame(Variable = nombres,
                             `β̂_i` = round(b, 4),
                             `Δx_i` = delta,
                             `β̂_i × Δx_i` = round(contrib, 4),
                             check.names = FALSE)
  print(df_contrib)

  # Paso 3: suma = log(OR)
  log_OR <- sum(contrib)
  cat(sprintf("\nPaso 3 — Suma = β'Δx = log(OR):  %.4f\n", log_OR))

  # Paso 4: exponenciar
  OR_val <- exp(log_OR)
  cat(sprintf("Paso 4 — OR = exp(%.4f)     = %.4f\n", log_OR, OR_val))
  cat("\nInterpretación: los individuos del perfil 1 tienen\n")
  cat(sprintf("  un OR de %.2f respecto al perfil 0.\n", OR_val))

  return(invisible(OR_val))
}

b_sin_alpha <- beta_chd[-1]
x1_cat <- c(CAT = 1, AGE = 40, ECG = 0)
x0_cat <- c(CAT = 0, AGE = 40, ECG = 0)

OR_cat <- calc_OR_detallado(x1_cat, x0_cat, beta_chd)

#> ═══════════════════════════════════════════════════
#>  CÁLCULO DEL OR POR PRODUCTO ESCALAR MATRICIAL
#> ═══════════════════════════════════════════════════
#> 
#> Perfil 1 (x1) — expuesto:      1 40 0 
#> Perfil 0 (x0) — no expuesto:   0 40 0 
#> 
#> Paso 1 — Δx = x1 - x0 (diferencia de perfiles):
#>     Variable x1 x0 Δx = x1-x0
#> CAT    b_CAT  1  0          1
#> AGE    b_AGE 40 40          0
#> ECG    b_ECG  0  0          0
#> 
#> Paso 2 — Contribución de cada variable: β_i × Δx_i
#>       Variable   β̂_i Δx_i β̂_i × Δx_i
#> b_CAT    b_CAT 0.652    1      0.652
#> b_AGE    b_AGE 0.029    0      0.000
#> b_ECG    b_ECG 0.342    0      0.000
#> 
#> Paso 3 — Suma = β'Δx = log(OR):  0.6520
#> Paso 4 — OR = exp(0.6520)     = 1.9194
#> 
#> Interpretación: los individuos del perfil 1 tienen
#>   un OR de 1.92 respecto al perfil 0.

2 El Modelo E, V, W — Confusión e Interacción

2.1 Marco conceptual

El modelo E, V, W es la formulación más general para estudios epidemiológicos analíticos:

\[\text{logit}[P(\mathbf{X})] = \alpha + \underbrace{\beta E}_{\substack{\text{exposición}\\\text{principal}}} + \underbrace{\sum_{j=1}^{p_1} \gamma_j V_j}_{\substack{\text{confusores}\\\text{(control)}}}\; + \underbrace{\sum_{l=1}^{p_2} \delta_l W_l \cdot E}_{\substack{\text{interacciones}\\\text{(modificación de efecto)}}}\]

Cuando hay interacción, el OR de \(E\) es:

\[\text{OR}_E = \exp\!\left(\beta + \sum_{l=1}^{p_2} \delta_l W_l\right)\]

El OR varía según el valor del modificador \(W_l\): esto es la modificación de efecto (heterogeneidad de efecto). El objetivo del epidemiólogo es:

Primero evaluar si hay interacción (prueba LR del término \(E \times W\))
Si hay interacción → reportar ORs estratificados por \(W\)
Si no hay interacción → evaluar confusión con el criterio del 10%

# ── Datos simulados con interacción E × tabaquismo ────────────────
set.seed(42)
n <- 800
datos_evw <- data.frame(
  E   = rbinom(n, 1, 0.50),
  AGE = rnorm(n, 50, 10),
  SMK = rbinom(n, 1, 0.40)
)
logit_p <- -2.0 + 0.8 * datos_evw$E + 0.03 * datos_evw$AGE +
            0.5 * datos_evw$SMK + 0.7 * datos_evw$E * datos_evw$SMK
datos_evw$D <- rbinom(n, 1, plogis(logit_p))

# Modelo E, V, W
mod_evw <- glm(D ~ E + AGE + SMK + E:SMK, data = datos_evw, family = binomial)

cat("── Coeficientes del modelo E, V, W ──\n")

#> ── Coeficientes del modelo E, V, W ──

tidy(mod_evw, conf.int = TRUE) |>
  mutate(
    OR         = round(exp(estimate), 4),
    `IC inf`   = round(exp(conf.low), 4),
    `IC sup`   = round(exp(conf.high), 4),
    estimate   = round(estimate, 4),
    std.error  = round(std.error, 4),
    p.value    = ifelse(p.value < 0.001, "<0.001", round(p.value, 4))
  ) |>
  rename(`β̂` = estimate, `SE(β̂)` = std.error, `p-valor` = p.value) |>
  tabla_doc(titulo = "Tabla 2.1 — Modelo E, V, W con interacción E × SMK")

Tabla 2.1 — Modelo E, V, W con interacción E × SMK
term	β	SE(β
(Intercept)	-2.0444	0.3959	-5.1640551	<0.001	-2.8313160	-1.2776477	0.1295	0.0589	0.2787
E	0.8028	0.1890	4.2467604	<0.001	0.4342270	1.1757456	2.2317	1.5438	3.2406
AGE	0.0311	0.0074	4.1848936	<0.001	0.0166820	0.0458604	1.0316	1.0168	1.0469
SMK	0.5834	0.2031	2.8725331	0.0041	0.1865652	0.9834333	1.7922	1.2051	2.6736
E:SMK	0.1544	0.3077	0.5019548	0.6157	-0.4459837	0.7609731	1.1670	0.6402	2.1404

# OR de E estratificado por SMK
coefs_evw <- coef(mod_evw)
b_E  <- coefs_evw["E"]
b_EW <- coefs_evw["E:SMK"]

cat("\n── OR de E según el valor del modificador SMK ──\n")

#> 
#> ── OR de E según el valor del modificador SMK ──

tibble(
  `SMK` = c(0, 1),
  `β_E + δ·SMK` = round(b_E + b_EW * c(0, 1), 4),
  `OR_E = exp(...)` = round(exp(b_E + b_EW * c(0, 1)), 4),
  `Interpretación` = c(
    "OR entre no fumadores: efecto de E solo",
    "OR entre fumadores: efecto de E potenciado por tabaquismo"
  )
) |> tabla_doc(titulo = "Tabla 2.2 — OR de E estratificado por tabaquismo (modificación de efecto)")

Tabla 2.2 — OR de E estratificado por tabaquismo (modificación de efecto)
SMK	β_E + δ·SMK	OR_E = exp(…)	Interpretación
0	0.8028	2.2317	OR entre no fumadores: efecto de E solo
1	0.9572	2.6043	OR entre fumadores: efecto de E potenciado por tabaquismo

3 Cálculo del Odds Ratio en la Práctica

3.1 Datos CHD — modelo de referencia

set.seed(100)
n_chd <- 609
datos_chd <- data.frame(
  CAT = rbinom(n_chd, 1, 0.30),
  AGE = round(rnorm(n_chd, 50, 10)),
  ECG = rbinom(n_chd, 1, 0.25)
)
logit_chd <- -3.911 + 0.652 * datos_chd$CAT +
              0.029 * datos_chd$AGE + 0.342 * datos_chd$ECG
datos_chd$CHD <- rbinom(n_chd, 1, plogis(logit_chd))

mod_chd <- glm(CHD ~ CAT + AGE + ECG, data = datos_chd, family = binomial)
cat("Modelo CHD ajustado con glm(). N =", nrow(datos_chd), "\n")

#> Modelo CHD ajustado con glm(). N = 609

cat("Eventos CHD =", sum(datos_chd$CHD), "\n")

#> Eventos CHD = 69

3.2 Tabla de ORs con incrementos clínicamente relevantes

# ══════════════════════════════════════════════════════════════════
#  SUBRUTINA: OR ajustados para incrementos personalizados
# ══════════════════════════════════════════════════════════════════
#
# Para variables continuas (como AGE) es más relevante el OR
# por 5 o 10 años que por 1 año. La fórmula es: OR(Δ) = exp(β·Δ)
#
tabla_OR_incremental <- function(mod, incrementos = NULL) {
  b      <- coef(mod)[-1]
  IC     <- tryCatch(confint(mod)[-1,], error = function(e) confint.default(mod)[-1,])
  pv     <- summary(mod)$coefficients[-1, 4]
  vars   <- names(b)

  if (is.null(incrementos)) incrementos <- setNames(rep(1, length(vars)), vars)

  lapply(vars, function(v) {
    d  <- if (!is.na(incrementos[v])) incrementos[v] else 1
    tibble(
      Variable    = v,
      `Δ (incremento)` = d,
      `β̂`         = round(b[v], 4),
      `β̂ × Δ`     = round(b[v] * d, 4),
      OR           = round(exp(b[v] * d), 4),
      `IC 95% inf` = round(exp(IC[v, 1] * d), 4),
      `IC 95% sup` = round(exp(IC[v, 2] * d), 4),
      `p-valor`    = ifelse(pv[v] < 0.001, "<0.001", as.character(round(pv[v], 4)))
    )
  }) |> bind_rows()
}

tabla_or <- tabla_OR_incremental(mod_chd,
                                  incrementos = c(CAT = 1, AGE = 5, ECG = 1))
tabla_doc(tabla_or,
          titulo = "Tabla 3.1 — OR ajustados: CAT y ECG por unidad; AGE por 5 años")

Tabla 3.1 — OR ajustados: CAT y ECG por unidad; AGE por 5 años
Variable	Δ (incremento)	β	β̂ ×
CAT	1	0.6764	0.6764	1.9668	1.1626	3.3066	0.0109
AGE	5	0.0479	0.2395	1.2706	1.1224	1.4426	<0.001
ECG	1	0.3929	0.3929	1.4813	0.8227	2.5871	0.1768

df_forest <- tibble(
  Variable = c("CAT (alto vs bajo)", "AGE (por 5 años)", "ECG (anormal vs normal)"),
  OR       = exp(coef(mod_chd)[-1] * c(1, 5, 1)),
  IC_L     = exp(confint.default(mod_chd)[-1, 1] * c(1, 5, 1)),
  IC_U     = exp(confint.default(mod_chd)[-1, 2] * c(1, 5, 1))
)

ggplot(df_forest, aes(x = OR, y = reorder(Variable, OR))) +
  geom_vline(xintercept = 1, linetype = "dashed", color = COL_GRIS, linewidth = 0.7) +
  geom_errorbarh(aes(xmin = IC_L, xmax = IC_U),
                 height = 0.20, color = COL_AZUL, linewidth = 0.9) +
  geom_point(size = 4.5, color = COL_AZUL, shape = 18) +
  geom_text(aes(label = sprintf("OR = %.2f\n[%.2f, %.2f]", OR, IC_L, IC_U)),
            hjust = -0.12, size = 3.2, color = COL_GRIS, fontface = "plain") +
  scale_x_log10(limits = c(0.5, 5), breaks = c(0.5, 1, 1.5, 2, 3)) +
  labs(title    = "Forest Plot — OR ajustados del modelo CHD",
       subtitle = "Escala logarítmica en el eje x",
       x = "Odds Ratio (escala log)",
       y = NULL,
       caption  = "IC 95% por perfil de verosimilitud. AGE expresado por 5 años.") +
  tema_doc +
  theme(panel.grid.major.y = element_line(color = "#e8f1f8"),
        panel.grid.major.x = element_line(color = "#e8f1f8"))

Figura 2. Forest plot de OR ajustados del modelo CHD. La línea punteada en OR = 1 representa la hipótesis nula.

4 CMáxima Verosimilitud y el Vector \(\hat{\boldsymbol{\beta}}\)

4.1 ¿Qué significa “estimar por Máxima Verosimilitud”?

La Máxima Verosimilitud (MV) responde a la pregunta: dado que observé los datos \(\mathcal{D} = \{(D_1, \mathbf{x}_1), \ldots, (D_n, \mathbf{x}_n)\}\), ¿qué valores de \(\boldsymbol{\beta}\) hacen que estos datos sean más probables de haber ocurrido?

Analogía epidemiológica: Si observamos que el 70% de los expuestos desarrolló la enfermedad y solo el 20% de los no expuestos, los parámetros que maximizan la verosimilitud son aquellos que predicen probabilidades cercanas al 70% y al 20% para esos perfiles. El MV “ajusta” \(\boldsymbol{\beta}\) hasta que las probabilidades predichas por el modelo sean lo más consistentes posible con los datos observados.

4.2 Función de verosimilitud en forma matricial

Dado que cada \(D_i \sim \text{Bernoulli}(p_i)\) con \(p_i = \sigma(\boldsymbol{\eta}_i) = \sigma(\mathbf{x}_i^{\top}\boldsymbol{\beta})\), la función de verosimilitud es:

\[L(\boldsymbol{\beta}) = \prod_{i=1}^{n} p_i^{D_i} (1 - p_i)^{1-D_i}\]

Tomando el logaritmo natural (log-verosimilitud):

\[\ell(\boldsymbol{\beta}) = \ln L(\boldsymbol{\beta}) = \sum_{i=1}^{n} \left[D_i \ln p_i + (1-D_i)\ln(1-p_i)\right]\]

Sustituyendo \(p_i = \sigma(\eta_i)\) y usando que \(\ln \sigma(\eta) = \eta - \ln(1+e^\eta)\) y \(\ln(1-\sigma(\eta)) = -\ln(1+e^\eta)\):

\[\ell(\boldsymbol{\beta}) = \sum_{i=1}^n \left[D_i \eta_i - \ln(1 + e^{\eta_i})\right]\]

Forma matricial compacta (con \(\boldsymbol{\eta} = \mathbf{X}\boldsymbol{\beta}\)):

\[\boxed{\ell(\boldsymbol{\beta}) = \mathbf{D}^{\top}\boldsymbol{\eta} - \mathbf{1}^{\top}\ln(\mathbf{1} + e^{\boldsymbol{\eta}}) = \mathbf{D}^{\top}\mathbf{X}\boldsymbol{\beta} - \mathbf{1}^{\top}\ln(\mathbf{1} + e^{\mathbf{X}\boldsymbol{\beta}})}\]

donde \(e^{\boldsymbol{\eta}}\) denota exponenciación elemento a elemento y \(\ln(\cdot)\) también es elemento a elemento.

4.3 Las ecuaciones de Score (gradiente)

Para encontrar el máximo de \(\ell(\boldsymbol{\beta})\), tomamos la derivada vectorial respecto a \(\boldsymbol{\beta}\) e igualamos a cero:

\[\frac{\partial \ell}{\partial \boldsymbol{\beta}} = \mathbf{X}^{\top}(\mathbf{D} - \mathbf{p}) = \mathbf{0}_{(k+1) \times 1}\]

Derivación detallada:

\[\frac{\partial \ell}{\partial \boldsymbol{\beta}} = \frac{\partial}{\partial \boldsymbol{\beta}} \left[\mathbf{D}^{\top}\mathbf{X}\boldsymbol{\beta} - \mathbf{1}^{\top}\ln(\mathbf{1} + e^{\mathbf{X}\boldsymbol{\beta}})\right]\]

Primer término: \(\frac{\partial}{\partial \boldsymbol{\beta}} \mathbf{D}^{\top}\mathbf{X}\boldsymbol{\beta} = \mathbf{X}^{\top}\mathbf{D}\)

Segundo término: \(\frac{\partial}{\partial \boldsymbol{\beta}} \mathbf{1}^{\top}\ln(\mathbf{1} + e^{\mathbf{X}\boldsymbol{\beta}}) = \mathbf{X}^{\top}\mathbf{p}\)

(porque \(\frac{\partial}{\partial \eta_i}\ln(1+e^{\eta_i}) = \frac{e^{\eta_i}}{1+e^{\eta_i}} = p_i\))

\[\therefore \quad \frac{\partial \ell}{\partial \boldsymbol{\beta}} = \mathbf{X}^{\top}\mathbf{D} - \mathbf{X}^{\top}\mathbf{p} = \mathbf{X}^{\top}(\mathbf{D} - \mathbf{p}) = \mathbf{0}\]

Interpretación geométrica de las ecuaciones de Score:

La condición \(\mathbf{X}^{\top}(\mathbf{D} - \mathbf{p}) = \mathbf{0}\) significa que el vector de residuales \((\mathbf{D} - \mathbf{p})\) debe ser ortogonal a cada columna de \(\mathbf{X}\). En otras palabras, la MV exige que los residuales “no tengan relación sistemática” con ninguno de los predictores. Son \(k+1\) ecuaciones simultáneas no lineales en \(k+1\) incógnitas \((\alpha, \beta_1, \ldots, \beta_k)\), sin solución analítica cerrada → requieren métodos iterativos.

4.4 La Hessiana y la Información de Fisher

La segunda derivada (curvatura) de \(\ell\) respecto a \(\boldsymbol{\beta}\):

\[\mathbf{H}(\boldsymbol{\beta}) = \frac{\partial^2 \ell}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^{\top}} = -\mathbf{X}^{\top}\mathbf{W}\mathbf{X}\]

donde \(\mathbf{W} = \text{diag}(w_1, w_2, \ldots, w_n)\) es la matriz de pesos diagonal con \(w_i = p_i(1-p_i) = f'(\eta_i)\).

Derivación:

\[\frac{\partial^2 \ell}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^{\top}} = \frac{\partial}{\partial \boldsymbol{\beta}^{\top}} \left[\mathbf{X}^{\top}(\mathbf{D} - \mathbf{p})\right] = -\mathbf{X}^{\top}\frac{\partial \mathbf{p}}{\partial \boldsymbol{\beta}^{\top}}\]

Como \(\frac{\partial p_i}{\partial \boldsymbol{\beta}} = p_i(1-p_i)\mathbf{x}_i = w_i \mathbf{x}_i\), entonces:

\[\frac{\partial \mathbf{p}}{\partial \boldsymbol{\beta}^{\top}} = \mathbf{W}\mathbf{X} \quad \Rightarrow \quad \mathbf{H} = -\mathbf{X}^{\top}\mathbf{W}\mathbf{X}\]

La información de Fisher esperada es \(\mathcal{I}(\boldsymbol{\beta}) = -\mathbb{E}[\mathbf{H}] = \mathbf{X}^{\top}\mathbf{W}\mathbf{X}\) (positiva semidefinida → el máximo es único bajo condiciones de identificabilidad).

4.5 El algoritmo IRLS — Mínimos Cuadrados Ponderados Iterativos

Dado que las ecuaciones de Score son no lineales, usamos el método de Newton-Raphson. Partiendo de un punto inicial \(\boldsymbol{\beta}^{(0)}\), la actualización Newton-Raphson es:

\[\boldsymbol{\beta}^{(t+1)} = \boldsymbol{\beta}^{(t)} - \left[\mathbf{H}^{(t)}\right]^{-1} \frac{\partial \ell}{\partial \boldsymbol{\beta}}\bigg|_{\boldsymbol{\beta}^{(t)}}\]

Sustituyendo \(\mathbf{H}^{(t)} = -\mathbf{X}^{\top}\mathbf{W}^{(t)}\mathbf{X}\) y \(\frac{\partial \ell}{\partial \boldsymbol{\beta}}\big|^{(t)} = \mathbf{X}^{\top}(\mathbf{D} - \mathbf{p}^{(t)})\):

\[\boldsymbol{\beta}^{(t+1)} = \boldsymbol{\beta}^{(t)} + \underbrace{\left(\mathbf{X}^{\top}\mathbf{W}^{(t)}\mathbf{X}\right)^{-1}}_{\text{inversa de }\mathcal{I}^{(t)}} \underbrace{\mathbf{X}^{\top}(\mathbf{D} - \mathbf{p}^{(t)})}_{\text{Score}^{(t)}}\]

¿Por qué se llama IRLS? En cada iteración se resuelve un sistema de mínimos cuadrados ponderados con respuesta de trabajo \(\mathbf{z}^{(t)}\):

\[\hat{\boldsymbol{\beta}}^{(t+1)} = \arg\min_{\boldsymbol{\beta}} \left(\mathbf{z}^{(t)} - \mathbf{X}\boldsymbol{\beta}\right)^{\top} \mathbf{W}^{(t)} \left(\mathbf{z}^{(t)} - \mathbf{X}\boldsymbol{\beta}\right)\]

donde la respuesta de trabajo es \(\mathbf{z}^{(t)} = \boldsymbol{\eta}^{(t)} + (\mathbf{W}^{(t)})^{-1}(\mathbf{D} - \mathbf{p}^{(t)})\).

La solución analítica de cada mínimos cuadrados ponderados es: \[\hat{\boldsymbol{\beta}}^{(t+1)} = \left(\mathbf{X}^{\top}\mathbf{W}^{(t)}\mathbf{X}\right)^{-1}\mathbf{X}^{\top}\mathbf{W}^{(t)}\mathbf{z}^{(t)}\]

# ══════════════════════════════════════════════════════════════════
#  IMPLEMENTACIÓN COMPLETA DEL ALGORITMO IRLS DESDE CERO
#
#  Pasos en cada iteración t:
#  (1) η^(t) = X β^(t)              predictor lineal
#  (2) p^(t) = σ(η^(t))             probabilidades actuales
#  (3) w_i^(t) = p_i^(t)(1-p_i^(t)) pesos diagonales de W
#  (4) X'W^(t)X                     matriz de información (k+1)×(k+1)
#  (5) Score = X'(D - p^(t))        gradiente de ℓ
#  (6) Δβ = (X'WX)^{-1} Score       paso Newton
#  (7) β^(t+1) = β^(t) + Δβ        actualización
#  (8) ℓ^(t) = D'η - 1'ln(1+exp(η)) log-verosimilitud
#  (9) ‖Δβ‖₂ < tol?                 criterio convergencia
# ══════════════════════════════════════════════════════════════════

IRLS_logistico <- function(X, y, tol = 1e-9, max_iter = 100,
                            verbose = TRUE) {
  n <- nrow(X); p <- ncol(X)

  # ── Inicialización: β⁰ = 0 ──────────────────────────────────────
  beta      <- rep(0, p)
  historial <- tibble(iter = integer(), log_lik = numeric(),
                      norm_delta = numeric())
  trayectoria <- matrix(NA_real_, nrow = max_iter + 1, ncol = p,
                        dimnames = list(NULL, colnames(X)))
  trayectoria[1,] <- beta

  if (verbose) {
    cat(sprintf("%-5s  %-16s  %-12s  %-s\n",
                "Iter", "log-lik ℓ(β)", "‖Δβ‖₂", "Estado"))
    cat(strrep("─", 58), "\n")
  }

  for (t in seq_len(max_iter)) {

    # (1) Predictor lineal: η = Xβ  [n × 1]
    eta   <- as.vector(X %*% beta)

    # (2) Probabilidades: p = σ(η), con clamp numérico
    p_hat <- pmin(pmax(1 / (1 + exp(-eta)), 1e-12), 1 - 1e-12)

    # (3) Pesos: w_i = p_i(1 - p_i) = f'(η_i)
    w <- p_hat * (1 - p_hat)

    # (4) Matriz de información: X'WX  [(k+1) × (k+1)]
    #     Forma eficiente: crossprod(X * w, X) evita crear W_nxn
    XtWX <- crossprod(X * w, X)

    # (5) Vector Score: X'(D - p)  [(k+1) × 1]
    score <- crossprod(X, y - p_hat)

    # (6) Paso Newton: (X'WX)^{-1} score
    delta_beta <- solve(XtWX, score)

    # (7) Actualización
    beta_nuevo <- beta + as.vector(delta_beta)

    # (8) Log-verosimilitud: ℓ(β) = Σ[D_i log p_i + (1-D_i) log(1-p_i)]
    log_lik <- sum(y * log(p_hat) + (1 - y) * log(1 - p_hat))

    # (9) Criterio de convergencia: ‖Δβ‖₂
    norm_d <- sqrt(sum(delta_beta^2))

    historial <- bind_rows(historial,
                           tibble(iter = t, log_lik = log_lik,
                                  norm_delta = norm_d))
    trayectoria[t + 1,] <- beta_nuevo

    if (verbose) {
      estado <- if (norm_d < tol) "✓ CONVERGIÓ" else ""
      cat(sprintf("  %-3d    %-16.6f  %-12.3e  %s\n",
                  t, log_lik, norm_d, estado))
    }

    beta <- beta_nuevo
    if (norm_d < tol) break
  }

  cat(sprintf("\n✓ Convergencia en iteración %d  (‖Δβ‖₂ = %.2e)\n",
              t, norm_d))

  return(list(beta_hat     = beta,
              iteraciones  = t,
              log_lik      = log_lik,
              XtWX_final   = XtWX,
              historial    = historial,
              trayectoria  = trayectoria[1:(t + 1),]))
}

# ── Aplicación al modelo CHD ──────────────────────────────────────
cat("═══ IRLS — Modelo CHD: CAT + AGE + ECG ═══\n\n")

#> ═══ IRLS — Modelo CHD: CAT + AGE + ECG ═══

X_chd <- model.matrix(CHD ~ CAT + AGE + ECG, data = datos_chd)
y_chd <- datos_chd$CHD

res_IRLS <- IRLS_logistico(X_chd, y_chd, verbose = TRUE)

#> Iter   log-lik ℓ(β)   ‖Δβ‖₂  Estado
#> ────────────────────────────────────────────────────────── 
#>   1      -422.126633       2.650e+00     
#>   2      -218.088738       1.583e+00     
#>   3      -205.217854       6.706e-01     
#>   4      -204.384091       7.070e-02     
#>   5      -204.376969       6.738e-04     
#>   6      -204.376968       6.132e-08     
#>   7      -204.376968       2.135e-15     ✓ CONVERGIÓ
#> 
#> ✓ Convergencia en iteración 7  (‖Δβ‖₂ = 2.14e-15)

df_h <- res_IRLS$historial

p_ll <- ggplot(df_h, aes(x = iter, y = log_lik)) +
  geom_line(color = COL_AZUL, linewidth = 1.5) +
  geom_point(color = COL_AZUL, size = 3, shape = 16) +
  labs(title    = "Log-verosimilitud  ℓ(β) por iteración",
       subtitle = "Sube monótonamente hasta el máximo global",
       x = "Iteración t", y = "ℓ(β̂⁽ᵗ⁾)") +
  tema_doc

p_nd <- ggplot(df_h, aes(x = iter, y = norm_delta)) +
  geom_line(color = COL_NARANJA, linewidth = 1.5) +
  geom_point(color = COL_NARANJA, size = 3, shape = 16) +
  scale_y_log10(labels = scales::scientific) +
  labs(title    = "Norma del paso  ‖Δβ‖₂  (escala log)",
       subtitle = "Convergencia cuadrática del método de Newton",
       x = "Iteración t", y = "‖Δβ‖₂  (log-escala)") +
  tema_doc

gridExtra::grid.arrange(p_ll, p_nd, ncol = 2)

Figura 4. Convergencia del algoritmo IRLS. Izquierda: la log-verosimilitud aumenta monótonamente hasta el máximo. Derecha: la norma del paso cae exponencialmente (convergencia cuadrática de Newton-Raphson).

# ── Verificación: IRLS debe reproducir exactamente a glm() ────────
cat("── Comparación IRLS vs glm() ─────────────────────────\n\n")

#> ── Comparación IRLS vs glm() ─────────────────────────

tibble(
  Parámetro  = colnames(X_chd),
  `β̂ (IRLS)` = round(res_IRLS$beta_hat, 7),
  `β̂ (glm)` = round(coef(mod_chd), 7),
  `|Δ|`      = round(abs(res_IRLS$beta_hat - coef(mod_chd)), 10),
  `Coincide` = ifelse(round(abs(res_IRLS$beta_hat - coef(mod_chd)), 6) < 1e-5,
                       "✓ Sí", "✗ No")
) |> tabla_doc(titulo = "Tabla 4.1 — Verificación: IRLS reproduce glm() hasta precisión de máquina")

Tabla 4.1 — Verificación: IRLS reproduce glm() hasta precisión de máquina
Parámetro	β̂ (IRLS	β̂ (gl
(Intercept)	-4.9089828	-4.9089828	✓ Sí
CAT	0.6764078	0.6764078	✓ Sí
AGE	0.0478915	0.0478915	✓ Sí
ECG	0.3928937	0.3928937	✓ Sí

cat("\nLog-verosimilitud IRLS: ", round(res_IRLS$log_lik, 6))

#> 
#> Log-verosimilitud IRLS:  -204.377

cat("\nLog-verosimilitud glm(): ", round(as.numeric(logLik(mod_chd)), 6), "\n")

#> 
#> Log-verosimilitud glm():  -204.377

4.6 La matriz de varianza-covarianza de \(\hat{\boldsymbol{\beta}}\)

Por la teoría asintótica del MV, el estimador \(\hat{\boldsymbol{\beta}}\) es asintóticamente normal:

\[\hat{\boldsymbol{\beta}} \xrightarrow{d} \mathcal{N}_{k+1}\!\left(\boldsymbol{\beta}_0,\; \mathcal{I}^{-1}(\boldsymbol{\beta}_0)\right)\]

La matriz de varianza-covarianza se estima evaluando la inversa de la información de Fisher en \(\hat{\boldsymbol{\beta}}\):

\[\widehat{\text{Var}}(\hat{\boldsymbol{\beta}}) = \left[\mathcal{I}(\hat{\boldsymbol{\beta}})\right]^{-1} = \left(\mathbf{X}^{\top}\hat{\mathbf{W}}\mathbf{X}\right)^{-1}\]

Los errores estándar de cada \(\hat{\beta}_i\) son:

\[SE(\hat{\beta}_i) = \sqrt{\left[\left(\mathbf{X}^{\top}\hat{\mathbf{W}}\mathbf{X}\right)^{-1}\right]_{ii}}\]

(raíz cuadrada del elemento diagonal \(i\)-ésimo de la matriz inversa)

# ══════════════════════════════════════════════════════════════════
#  CÁLCULO DE Var(β̂) = (X'ŴX)⁻¹ DESDE PRIMEROS PRINCIPIOS
# ══════════════════════════════════════════════════════════════════

cat("══ Matriz de Información de Fisher: I(β̂) = X'ŴX ══\n\n")

#> ══ Matriz de Información de Fisher: I(β̂) = X'ŴX ══

print(round(res_IRLS$XtWX_final, 4))

#>             (Intercept)       CAT         AGE      ECG
#> (Intercept)     58.7981   24.3566   3178.8890  16.4435
#> CAT             24.3566   24.3566   1295.2842   6.2651
#> AGE           3178.8890 1295.2842 178034.9779 885.1079
#> ECG             16.4435    6.2651    885.1079  16.4435

cat("\n── Matriz Var(β̂) = [X'ŴX]⁻¹ ──\n\n")

#> 
#> ── Matriz Var(β̂) = [X'ŴX]⁻¹ ──

V_hat <- solve(res_IRLS$XtWX_final)
print(round(V_hat, 8))

#>             (Intercept)         CAT         AGE         ECG
#> (Intercept)  0.52601545 -0.04361820 -0.00893287 -0.02856516
#> CAT         -0.04361820  0.07059467  0.00024861  0.00333900
#> AGE         -0.00893287  0.00024861  0.00016298  0.00006516
#> ECG         -0.02856516  0.00333900  0.00006516  0.08459999

SE_IRLS <- sqrt(diag(V_hat))
SE_glm  <- sqrt(diag(vcov(mod_chd)))

cat("\n── Errores estándar: raíces de los diagonales de Var(β̂) ──\n")

#> 
#> ── Errores estándar: raíces de los diagonales de Var(β̂) ──

tibble(
  Parámetro    = colnames(X_chd),
  `SE (IRLS)`  = round(SE_IRLS, 6),
  `SE (glm())`  = round(SE_glm, 6),
  `Diferencia`  = round(abs(SE_IRLS - SE_glm), 9)
) |> tabla_doc(titulo = "Tabla 4.2 — Errores estándar de β̂: IRLS vs vcov(glm)")

Tabla 4.2 — Errores estándar de β̂: IRLS vs vcov(glm)
Parámetro	SE (IRLS)	SE (glm())	Diferencia
(Intercept)	0.725269	0.725268	1.3e-06
CAT	0.265697	0.265696	3.0e-07
AGE	0.012767	0.012767	0.0e+00
ECG	0.290861	0.290861	3.0e-07

5 Inferencia: Pruebas y Errores Estándar

5.1 Prueba de Wald

Para cada \(\beta_i\), la prueba de Wald contrasta \(H_0: \beta_i = 0\) contra \(H_a: \beta_i \neq 0\):

\[Z_i = \frac{\hat{\beta}_i}{SE(\hat{\beta}_i)} = \frac{\hat{\beta}_i}{\sqrt{[\hat{\mathbf{V}}]_{ii}}} \xrightarrow{d} \mathcal{N}(0,1)\]

El p-valor bilateral es \(2\Phi(-|Z_i|)\) donde \(\Phi\) es la CDF de \(\mathcal{N}(0,1)\).

El intervalo de confianza de Wald al 95% es: \[\hat{\beta}_i \pm 1.96 \cdot SE(\hat{\beta}_i) \quad \Rightarrow \quad \text{IC}_{95\%}(\text{OR}_i) = \left[e^{\hat{\beta}_i - 1.96\,SE_i},\; e^{\hat{\beta}_i + 1.96\,SE_i}\right]\]

Limitación de Wald: En muestras pequeñas o con parámetros muy grandes (separación completa), la aproximación normal es pobre. En esos casos se prefieren los IC basados en el perfil de verosimilitud (ver más abajo).

# ══════════════════════════════════════════════════════════════════
#  SUBRUTINA: Prueba de Wald desde primeros principios
#  Usa directamente Var(β̂) = (X'ŴX)⁻¹ del IRLS
# ══════════════════════════════════════════════════════════════════

prueba_wald_completa <- function(beta_hat, V_hat, nombres = NULL, alfa = 0.05) {
  SE    <- sqrt(diag(V_hat))
  Z     <- beta_hat / SE
  pvals <- 2 * pnorm(-abs(Z))
  z_alfa <- qnorm(1 - alfa/2)   # cuantil crítico (1.96 para α=0.05)

  tibble(
    Parámetro     = if (!is.null(nombres)) nombres else paste0("β", seq_along(beta_hat)),
    `β̂`           = round(beta_hat, 4),
    `SE(β̂)`        = round(SE, 4),
    `Z = β̂/SE`    = round(Z, 3),
    `p-valor`      = ifelse(pvals < 0.001, "<0.001", as.character(round(pvals, 4))),
    `Sig.`         = case_when(
                       pvals < 0.001 ~ "***",
                       pvals < 0.01  ~ "**",
                       pvals < 0.05  ~ "*",
                       pvals < 0.10  ~ ".",
                       TRUE          ~ ""),
    OR             = round(exp(beta_hat), 4),
    `IC 95% inf`   = round(exp(beta_hat - z_alfa * SE), 4),
    `IC 95% sup`   = round(exp(beta_hat + z_alfa * SE), 4)
  )
}

tab_wald <- prueba_wald_completa(res_IRLS$beta_hat, V_hat,
                                  nombres = colnames(X_chd))
tabla_doc(tab_wald, titulo = "Tabla 5.1 — Prueba de Wald para todos los coeficientes (desde IRLS)")

Tabla 5.1 — Prueba de Wald para todos los coeficientes (desde IRLS)
Parámetro	β	SE(β
(Intercept)	-4.9090	0.7253	-6.768	<0.001	***	0.0074	0.0018	0.0306
CAT	0.6764	0.2657	2.546	0.0109		1.9668	1.1684	3.3107
AGE	0.0479	0.0128	3.751	<0.001	***	1.0491	1.0231	1.0756
ECG	0.3929	0.2909	1.351	0.1768		1.4813	0.8376	2.6195

5.2 Prueba de Razón de Verosimilitudes (LR Test)

La prueba LR compara dos modelos anidados: el completo (con \(q\) parámetros adicionales) vs. el reducido (sin esos parámetros):

\[LR = -2\left[\ell(\hat{\boldsymbol{\beta}}_{\text{reducido}}) - \ell(\hat{\boldsymbol{\beta}}_{\text{completo}})\right] \sim \chi^2_{q}\]

donde \(q = \dim(\boldsymbol{\beta}_{\text{completo}}) - \dim(\boldsymbol{\beta}_{\text{reducido}})\).

¿Por qué −2 × log? El factor −2 surge de la teoría asintótica: por el teorema de Wilks, la diferencia de log-verosimilitudes multiplicada por 2 converge a una distribución chi-cuadrado bajo \(H_0\). La prueba LR es generalmente más potente que la de Wald, especialmente cuando los parámetros son grandes.

# ══════════════════════════════════════════════════════════════════
#  SUBRUTINA: Prueba LR con diagnóstico completo
# ══════════════════════════════════════════════════════════════════

lr_test_completo <- function(mod_reducido, mod_completo, alfa = 0.05,
                              nombre_hipotesis = "H₀: β = 0") {
  ll_r  <- as.numeric(logLik(mod_reducido))
  ll_c  <- as.numeric(logLik(mod_completo))
  LR    <- -2 * (ll_r - ll_c)
  q     <- length(coef(mod_completo)) - length(coef(mod_reducido))
  pval  <- pchisq(LR, df = q, lower.tail = FALSE)
  chi_c <- qchisq(1 - alfa, df = q)   # valor crítico

  cat(sprintf("── Prueba de Razón de Verosimilitudes (%s) ──\n", nombre_hipotesis))
  cat(sprintf("  ℓ(modelo reducido):   %8.4f\n", ll_r))
  cat(sprintf("  ℓ(modelo completo):   %8.4f\n", ll_c))
  cat(sprintf("  Δℓ = ℓC - ℓR:         %8.4f\n", ll_c - ll_r))
  cat(sprintf("  LR = -2Δℓ:            %8.4f\n", LR))
  cat(sprintf("  Grados de libertad q: %8d\n",   q))
  cat(sprintf("  χ²(q=%.0f, α=%.2f):  %8.4f  (valor crítico)\n", q, alfa, chi_c))
  cat(sprintf("  p-valor:              %8.5f\n", pval))
  cat(sprintf("  LR > χ²_crit?         %s\n", ifelse(LR > chi_c, "SÍ → Rechazar H₀", "NO → No rechazar H₀")))

  return(invisible(list(LR = LR, df = q, pval = pval)))
}

# Test para CAT (¿aporta CAT al modelo más allá de AGE y ECG?)
mod_sin_cat <- glm(CHD ~ AGE + ECG, data = datos_chd, family = binomial)
lr_test_completo(mod_sin_cat, mod_chd, nombre_hipotesis = "H₀: β_CAT = 0")

#> ── Prueba de Razón de Verosimilitudes (H₀: β_CAT = 0) ──
#>   ℓ(modelo reducido):   -207.5303
#>   ℓ(modelo completo):   -204.3770
#>   Δℓ = ℓC - ℓR:           3.1533
#>   LR = -2Δℓ:              6.3066
#>   Grados de libertad q:        1
#>   χ²(q=1, α=0.05):    3.8415  (valor crítico)
#>   p-valor:               0.01203
#>   LR > χ²_crit?         SÍ → Rechazar H₀

5.3 Intervalos de Confianza — Perfil de Verosimilitud vs Wald

Los IC basados en el perfil de verosimilitud se construyen como el conjunto de valores de \(\beta_i\) compatibles con los datos:

\[\text{IC}_{95\%}(\beta_i) = \left\{\beta_i : -2\left[\ell(\beta_i^*) - \ell(\hat{\beta}_i)\right] \leq \chi^2_{1, 0.95}\right\}\]

No asumen simetría gaussiana → son más fiables que Wald con muestras pequeñas.

ic_comparativo <- function(mod) {
  IC_wald <- confint.default(mod)
  IC_perf <- tryCatch(confint(mod), error = function(e) IC_wald)
  b  <- coef(mod)

  tibble(
    Parámetro    = names(b),
    `β̂`          = round(b, 4),
    OR           = round(exp(b), 4),
    `OR Wald inf`  = round(exp(IC_wald[,1]), 4),
    `OR Wald sup`  = round(exp(IC_wald[,2]), 4),
    `OR Perfil inf`= round(exp(IC_perf[,1]), 4),
    `OR Perfil sup`= round(exp(IC_perf[,2]), 4)
  )
}

tabla_doc(ic_comparativo(mod_chd),
          titulo = "Tabla 5.2 — IC 95% Wald vs. Perfil de Verosimilitud (escala OR)")

Tabla 5.2 — IC 95% Wald vs. Perfil de Verosimilitud (escala OR)
Parámetro	β	O
(Intercept)	-4.9090	0.0074	0.0018	0.0306	0.0017	0.0295
CAT	0.6764	1.9668	1.1684	3.3107	1.1626	3.3066
AGE	0.0479	1.0491	1.0231	1.0756	1.0234	1.0760
ECG	0.3929	1.4813	0.8376	2.6195	0.8227	2.5871

6 Estrategia de Modelado: Confusión e Interacción

6.1 Estrategia jerárquica

La estrategia de modelado canónica opera en dos fases secuenciales:

FASE A — Evaluación de Interacción (primero siempre)

Ajustar modelo completo: \(E + V_j + E \cdot W_l\)
Prueba LR: ¿el término \(E \cdot W_l\) es significativo (\(p < 0.05\))?
- SÍ → retener interacción → reportar OR estratificados por \(W_l\)
- NO → pasar a Fase B

FASE B — Evaluación de Confusión (solo si no hay interacción)

Para cada \(V_j\): calcular OR de \(E\) con y sin \(V_j\)
¿\(|\text{OR}_{\text{con}} - \text{OR}_{\text{sin}}| / |\text{OR}_{\text{con}}| \geq 10\%\)?
- SÍ → \(V_j\) es confusor → retener en el modelo
- NO → \(V_j\) no confunde → puede excluirse

# ══════════════════════════════════════════════════════════════════
#  SUBRUTINA: Evaluación sistemática de confusión (criterio 10%)
# ══════════════════════════════════════════════════════════════════

evaluar_confusion <- function(datos, outcome, exposicion,
                               confusores) {
  f_completa <- as.formula(paste(outcome, "~",
                           paste(c(exposicion, confusores), collapse = " + ")))
  mod_c <- glm(f_completa, data = datos, family = binomial)
  OR_c  <- exp(coef(mod_c)[exposicion])

  lapply(confusores, function(v) {
    vars_r <- setdiff(confusores, v)
    f_r <- if (length(vars_r) == 0)
             as.formula(paste(outcome, "~", exposicion))
           else
             as.formula(paste(outcome, "~",
                         paste(c(exposicion, vars_r), collapse = " + ")))
    mod_r  <- glm(f_r, data = datos, family = binomial)
    OR_r   <- exp(coef(mod_r)[exposicion])
    cambio <- abs(OR_c - OR_r) / abs(OR_c) * 100

    tibble(
      `Confusor candidato` = v,
      `OR con v` = round(OR_c, 4),
      `OR sin v` = round(OR_r, 4),
      `|Cambio %|` = round(cambio, 2),
      Decisión = ifelse(cambio >= 10,
                        "✓ Confusor (≥10%)",
                        "✗ No confunde (<10%)")
    )
  }) |> bind_rows()
}

tab_conf <- evaluar_confusion(datos_chd, "CHD", "CAT", c("AGE", "ECG"))
tabla_doc(tab_conf, titulo = "Tabla 6.1 — Evaluación de confusión para CAT (criterio del 10%, Kleinbaum & Klein)")

Tabla 6.1 — Evaluación de confusión para CAT (criterio del 10%, Kleinbaum & Klein)
Confusor candidato	OR con v	OR sin v	\|Cambio %\|	Decisión
AGE	1.9668	1.8611	5.37	✗ No confunde (<10%)
ECG	1.9668	1.9423	1.25	✗ No confunde (<10%)

# ══════════════════════════════════════════════════════════════════
#  SUBRUTINA: Evaluación de interacción mediante prueba LR
# ══════════════════════════════════════════════════════════════════

evaluar_interaccion <- function(datos, outcome, exposicion,
                                 confusores, modificadores, alfa = 0.05) {
  f_base <- as.formula(paste(outcome, "~",
                        paste(c(exposicion, confusores), collapse = " + ")))
  mod_base <- glm(f_base, data = datos, family = binomial)

  lapply(modificadores, function(w) {
    f_inter <- as.formula(paste(outcome, "~",
                           paste(c(exposicion, confusores,
                                   paste0(exposicion, ":", w)), collapse = " + ")))
    mod_inter <- glm(f_inter, data = datos, family = binomial)
    LR <- as.numeric(-2 * (logLik(mod_base) - logLik(mod_inter)))
    pv <- pchisq(LR, df = 1, lower.tail = FALSE)

    tibble(
      `Término de interacción` = paste0(exposicion, " × ", w),
      `LR stat`  = round(LR, 4),
      df = 1,
      `p-valor` = ifelse(pv < 0.001, "<0.001", as.character(round(pv, 4))),
      Decisión  = ifelse(pv < alfa,
                          "✓ Interacción significativa → OR estratificados",
                          "✗ Sin interacción → evaluar confusión")
    )
  }) |> bind_rows()
}

tab_inter <- evaluar_interaccion(datos_chd, "CHD", "CAT", "AGE", "ECG")
tabla_doc(tab_inter, titulo = "Tabla 7.1 — Evaluación de interacción CAT × ECG mediante prueba LR")

Tabla 7.1 — Evaluación de interacción CAT × ECG mediante prueba LR
Término de interacción	LR stat	df	p-valor	Decisión
CAT × ECG	0.7224	1	0.3953	✗ Sin interacción → evaluar confusión

7 Bondad de Ajuste

7.1 Devianza, pseudo-R² y criterios de información

La devianza del modelo es \(D = -2\ell(\hat{\boldsymbol{\beta}})\). La devianza del modelo nulo (solo intercepto) es \(D_0 = -2\ell(\hat{\alpha}_0)\). La reducción \(D_0 - D\) tiene distribución \(\chi^2\) bajo \(H_0\) de que todos los \(\beta_i = 0\).

Los pseudo-\(R^2\) no son idénticos al \(R^2\) de la regresión lineal, pero miden la mejora relativa del ajuste:

Pseudo-R²	Fórmula	Rango típico
McFadden	\(1 - D/D_0 = 1 - \ell_M/\ell_0\)	[0, 1]; valores ≥ 0.2 = buen ajuste
Cox-Snell	\(1 - \exp[-(D_0 - D)/n]\)	Máximo < 1
Nagelkerke	Cox-Snell / máximo teórico	[0, 1]

bondad_ajuste_completa <- function(mod) {
  D0   <- mod$null.deviance;   gl0 <- mod$df.null
  Dm   <- mod$deviance;        glm_  <- mod$df.residual
  Delta_D  <- D0 - Dm
  gl_diff  <- gl0 - glm_
  p_chi    <- pchisq(Delta_D, df = gl_diff, lower.tail = FALSE)

  # Pseudo R²
  R2_mcf   <- 1 - Dm / D0
  R2_cs    <- 1 - exp(-Delta_D / length(mod$y))
  R2_cs_max <- 1 - exp(-D0 / length(mod$y))
  R2_nag   <- R2_cs / R2_cs_max

  tibble(
    Métrica = c("Devianza nula D₀", "Devianza modelo D",
                "Reducción ΔD = D₀ − D", "gl del test",
                "Prueba Chi² (p-valor)", "Pseudo-R² McFadden",
                "Pseudo-R² Cox-Snell", "Pseudo-R² Nagelkerke",
                "AIC", "BIC"),
    Valor = c(round(D0,4), round(Dm,4), round(Delta_D,4), gl_diff,
              ifelse(p_chi < 0.001, "<0.001", round(p_chi,5)),
              round(R2_mcf,4), round(R2_cs,4), round(R2_nag,4),
              round(AIC(mod),2), round(BIC(mod),2))
  ) |> tabla_doc(titulo = "Tabla 9.1 — Métricas de bondad de ajuste del modelo CHD")
}

bondad_ajuste_completa(mod_chd)

Tabla 9.1 — Métricas de bondad de ajuste del modelo CHD
Métrica	Valor
Devianza nula D₀	430.3933
Devianza modelo D	408.7539
Reducción ΔD = D₀ − D	21.6393
gl del test	3
Prueba Chi² (p-valor)	<0.001
Pseudo-R² McFadden	0.0503
Pseudo-R² Cox-Snell	0.0349
Pseudo-R² Nagelkerke	0.0689
AIC	416.75
BIC	434.4

7.2 Test de Hosmer-Lemeshow

El test de Hosmer-Lemeshow evalúa la calibración del modelo: ¿las probabilidades predichas corresponden a las frecuencias observadas? Se divide la muestra en \(g = 10\) grupos por deciles de riesgo predicho y se compara observados vs. esperados con un \(\chi^2\):

\[H_L = \sum_{j=1}^{g} \left[\frac{(O_{1j} - E_{1j})^2}{E_{1j}} + \frac{(O_{0j} - E_{0j})^2}{E_{0j}}\right] \sim \chi^2_{g-2}\]

Un p-valor no significativo (>0.05) indica buen ajuste (no se rechaza \(H_0\) de buen ajuste).

hl_res <- hoslem.test(mod_chd$y, fitted(mod_chd), g = 10)
cat(sprintf("Test Hosmer-Lemeshow: χ²(8) = %.4f,  p = %.4f\n",
            hl_res$statistic, hl_res$p.value))

#> Test Hosmer-Lemeshow: χ²(8) = 10.2831,  p = 0.2457

cat("Conclusión:", ifelse(hl_res$p.value > 0.05,
    "No se rechaza H₀ — calibración adecuada (p > 0.05).",
    "Se rechaza H₀ — evidencia de mal ajuste (p ≤ 0.05)."), "\n")

#> Conclusión: No se rechaza H₀ — calibración adecuada (p > 0.05).

# Tabla de calibración
probs_hl <- fitted(mod_chd)
grupos   <- cut(probs_hl, breaks = quantile(probs_hl, seq(0, 1, 0.1)),
                include.lowest = TRUE, labels = 1:10)

tab_hl <- lapply(levels(grupos), function(g) {
  idx <- which(grupos == g)
  if (length(idx) == 0) return(NULL)
  tibble(
    Grupo         = as.integer(g),
    n             = length(idx),
    `Obs. casos`  = sum(mod_chd$y[idx]),
    `Esp. casos`  = round(sum(probs_hl[idx]), 2),
    `Obs. ctrls`  = length(idx) - sum(mod_chd$y[idx]),
    `Esp. ctrls`  = round(length(idx) - sum(probs_hl[idx]), 2),
    `P̄ media`     = round(mean(probs_hl[idx]), 4)
  )
}) |>
  Filter(Negate(is.null), x = _) |>
  bind_rows()

tabla_doc(tab_hl, titulo = "Tabla 9.2 — Calibración Hosmer-Lemeshow (g = 10 grupos por decil de riesgo)")

Tabla 9.2 — Calibración Hosmer-Lemeshow (g = 10 grupos por decil de riesgo)
Grupo	n	Obs. casos	Esp. casos	Obs. ctrls	Esp. ctrls	P̄ medi
1	66	1	2.73	65	63.27	0.0414
2	56	1	3.15	55	52.85	0.0563
3	68	8	4.54	60	63.46	0.0668
4	59	5	4.79	54	54.21	0.0812
5	64	9	6.06	55	57.94	0.0947
6	59	5	6.41	54	52.59	0.1087
7	57	5	7.09	52	49.91	0.1244
8	59	12	8.47	47	50.53	0.1435
9	60	9	10.51	51	49.49	0.1751
10	61	14	15.24	47	45.76	0.2499

8 Curvas ROC y Discriminación

8.1 Área Bajo la Curva (AUC)

La curva ROC grafica la Sensibilidad vs. (1 − Especificidad) para todos los umbrales de clasificación. El AUC tiene una interpretación probabilística directa:

\[\text{AUC} = P(\hat{p}_{\text{caso}} > \hat{p}_{\text{control}})\]

Es decir, la probabilidad de que un caso tomado al azar tenga mayor riesgo predicho que un control tomado al azar.

AUC	Discriminación
0.50	Sin discriminación (azar)
0.60–0.70	Pobre
0.70–0.80	Aceptable
0.80–0.90	Buena
> 0.90	Excelente

y_obs   <- datos_chd$CHD
p_pred  <- fitted(mod_chd)

roc_obj  <- roc(y_obs, p_pred, quiet = TRUE)
AUC_val  <- as.numeric(auc(roc_obj))
IC_auc   <- ci(roc_obj)
coord_op <- coords(roc_obj, "best",
                   ret = c("threshold","sensitivity","specificity","ppv","npv"))

plot(roc_obj,
     main  = sprintf("Curva ROC — Modelo CHD   AUC = %.3f  (IC 95%%: %.3f – %.3f)",
                     AUC_val, IC_auc[1], IC_auc[3]),
     col   = COL_AZUL, lwd = 2.5,
     print.auc = FALSE,
     auc.polygon = TRUE, auc.polygon.col = COL_PALIDO,
     print.thres = "best",
     print.thres.pattern = "Umbral = %.3f",
     print.thres.col = COL_NARANJA,
     identity.lty = 2, identity.col = "gray60",
     cex.main = 0.9, font.main = 2)

Figura 5. Curva ROC del modelo CHD. El área coloreada representa el AUC. El punto naranja indica el umbral óptimo según el índice de Youden (máxima Sensibilidad + Especificidad − 1).

# Tabla de métricas de desempeño
tibble(
  Métrica = c("AUC", "IC 95% inf", "IC 95% sup",
              "Umbral óptimo (Youden)",
              "Sensibilidad", "Especificidad",
              "VPP (Valor Predictivo Positivo)",
              "VPN (Valor Predictivo Negativo)"),
  Valor = c(round(AUC_val, 4),
            round(IC_auc[1], 4), round(IC_auc[3], 4),
            round(coord_op$threshold, 4),
            round(coord_op$sensitivity, 4),
            round(coord_op$specificity, 4),
            round(coord_op$ppv, 4),
            round(coord_op$npv, 4))
) |> tabla_doc(titulo = "Tabla 10.1 — Métricas de desempeño discriminatorio del modelo CHD")

Tabla 10.1 — Métricas de desempeño discriminatorio del modelo CHD
Métrica	Valor
AUC	0.6647
IC 95% inf	0.6010
IC 95% sup	0.7283
Umbral óptimo (Youden)	0.0899
Sensibilidad	0.7826
Especificidad	0.4667
VPP (Valor Predictivo Positivo)	0.1579
VPN (Valor Predictivo Negativo)	0.9438

9 Logística Condicional (Datos pareados)

9.1 ¿Por qué la logística estándar falla con pareamiento?

En estudios con apareamiento 1:M, cada par/estrato \(s\) tiene su propio intercepto \(\alpha_s\) (efecto del estrato). Si estimamos \(\hat{\boldsymbol{\beta}}\) con todos los \(\alpha_s\), el número de parámetros crece con \(n\) → inconsistencia de Neyman-Scott: los \(\hat{\beta}_i\) son sesgados.

La solución es la logística condicional, que elimina los \(\alpha_s\) como parámetros de molestia (nuisance parameters) condicionando la verosimilitud sobre el total de casos en cada estrato. La log-verosimilitud condicional es:

\[\ell_c(\boldsymbol{\beta}) = \sum_{s=1}^{S} \ln\!\left[\frac{\exp(\boldsymbol{\beta}^{\top}\mathbf{x}_{\text{caso},s})}{\sum_{j \in R_s} \exp(\boldsymbol{\beta}^{\top}\mathbf{x}_{js})}\right]\]

donde \(R_s\) es el conjunto de todos los individuos del par/grupo \(s\) (caso + controles).

set.seed(99)
n_pares  <- 300
datos_ap <- data.frame(
  par_id = rep(1:n_pares, each = 2),
  caso   = rep(c(1, 0), n_pares),
  E      = c(rbinom(n_pares, 1, 0.45), rbinom(n_pares, 1, 0.28)),
  AGE    = round(rnorm(n_pares * 2, 52, 10)),
  SMK    = rbinom(n_pares * 2, 1, 0.40)
)

mod_cond <- clogit(caso ~ E + AGE + SMK + strata(par_id), data = datos_ap)

tidy(mod_cond, exponentiate = TRUE, conf.int = TRUE) |>
  mutate(across(c(estimate, conf.low, conf.high, p.value), ~round(., 4))) |>
  rename(OR = estimate, `IC inf` = conf.low, `IC sup` = conf.high,
         `p-valor` = p.value) |>
  tabla_doc(titulo = "Tabla 11.1 — OR del modelo logístico condicional (apareado 1:1)")

Tabla 11.1 — OR del modelo logístico condicional (apareado 1:1)
term	OR	std.error	statistic	p-valor	IC inf	IC sup
E	0.8170	0.1743412	-1.1593498	0.2463	0.5805	1.1498
AGE	0.9955	0.0079702	-0.5626742	0.5737	0.9801	1.0112
SMK	0.9974	0.1651279	-0.0154969	0.9876	0.7217	1.3786

10 Logística Politómica (Desenlace Multinomial)

10.1 Estructura del modelo

Para desenlaces con \(J \geq 3\) categorías no ordenadas, se ajustan \(J-1\) ecuaciones logísticas simultáneas contra la categoría de referencia \(j_0\):

\[\ln\!\left(\frac{P(Y = j \mid \mathbf{X})}{P(Y = j_0 \mid \mathbf{X})}\right) = \alpha_j + \boldsymbol{\beta}_j^{\top}\mathbf{X}, \quad j \neq j_0\]

Cada ecuación produce su propio vector \(\boldsymbol{\beta}_j\) (efectos de las covariables diferentes para cada categoría vs. referencia). Las probabilidades se recuperan mediante:

\[P(Y = j \mid \mathbf{X}) = \frac{e^{\alpha_j + \boldsymbol{\beta}_j^{\top}\mathbf{X}}}{1 + \sum_{j' \neq j_0} e^{\alpha_{j'} + \boldsymbol{\beta}_{j'}^{\top}\mathbf{X}}}\]

set.seed(55)
n_pol <- 700
datos_pol <- data.frame(
  AGE  = rnorm(n_pol, 55, 12),
  SMK  = rbinom(n_pol, 1, 0.40),
  TIPO = sample(c("Control","Tipo_A","Tipo_B"), n_pol, replace = TRUE,
                prob = c(0.5, 0.3, 0.2))
)
datos_pol$TIPO <- relevel(factor(datos_pol$TIPO), ref = "Control")

mod_pol <- multinom(TIPO ~ AGE + SMK, data = datos_pol, trace = FALSE)

tidy(mod_pol, exponentiate = TRUE, conf.int = TRUE) |>
  mutate(across(c(estimate, conf.low, conf.high), ~round(., 3))) |>
  rename(`Categoría vs Control` = y.level, OR = estimate,
         `IC inf` = conf.low, `IC sup` = conf.high) |>
  tabla_doc(titulo = "Tabla 12.1 — OR del modelo politómico (referencia: Control)")

Tabla 12.1 — OR del modelo politómico (referencia: Control)
Categoría vs Control	term	OR	std.error	statistic	p.value	IC inf	IC sup
Tipo_A	(Intercept)	0.369	0.4248937	-2.3460670	0.0189727	0.160	0.849
Tipo_A	AGE	1.008	0.0074775	1.0856485	0.2776346	0.993	1.023
Tipo_A	SMK	0.985	0.1782506	-0.0875834	0.9302078	0.694	1.396
Tipo_B	(Intercept)	0.438	0.4649810	-1.7777166	0.0754504	0.176	1.088
Tipo_B	AGE	1.004	0.0082841	0.5212529	0.6021906	0.988	1.021
Tipo_B	SMK	0.507	0.2081621	-3.2621777	0.0011056	0.337	0.763

11 Logística Ordinal

11.1 Modelo de odds proporcionales

Para desenlaces ordinales con \(J\) niveles ordenados, el modelo de odds proporcionales (McCullagh, 1980) es:

\[\text{logit}[P(Y \leq j \mid \mathbf{X})] = \alpha_j - \boldsymbol{\beta}^{\top}\mathbf{X}, \quad j = 1, \ldots, J-1\]

Supuesto clave: \(\boldsymbol{\beta}\) es el mismo para todos los cortes \(j\). Solo los interceptos \(\alpha_j\) varían (representan los umbrales de corte en la escala latente). El OR de una covariable \(X_i\) es:

\[\text{OR}_i = e^{\beta_i}\]

y mide cuánto mayor es el odds de estar en una categoría más severa por unidad de \(X_i\), independientemente del punto de corte.

set.seed(77)
n_ord  <- 700
AGE_z  <- scale(rnorm(n_ord, 50, 10))[,1]
SMK_s  <- rbinom(n_ord, 1, 0.40)
DOSE_z <- scale(rnorm(n_ord, 0, 1))[,1]
indice <- 0.4 * AGE_z + 0.5 * SMK_s + 0.8 * DOSE_z + rnorm(n_ord, 0, 1.2)
SEV    <- cut(indice, breaks = quantile(indice, c(0,1/3,2/3,1)),
              labels = c("Leve","Moderada","Severa"), include.lowest = TRUE)
datos_ord <- data.frame(AGE_z = as.numeric(AGE_z), SMK = SMK_s,
                         DOSE_z = as.numeric(DOSE_z),
                         SEVERIDAD = ordered(SEV, levels = c("Leve","Moderada","Severa")))

cat("Distribución del desenlace ordinal:\n")

#> Distribución del desenlace ordinal:

print(table(datos_ord$SEVERIDAD))

#> 
#>     Leve Moderada   Severa 
#>      234      233      233

mod_ord <- polr(SEVERIDAD ~ AGE_z + SMK + DOSE_z,
                data = datos_ord, Hess = TRUE, method = "logistic")

tidy(mod_ord, exponentiate = TRUE, conf.int = TRUE) |>
  mutate(across(where(is.numeric), ~round(., 4))) |>
  rename(OR = estimate) |>
  tabla_doc(titulo = "Tabla 13.1 — Modelo de odds proporcionales (logística ordinal)")

Tabla 13.1 — Modelo de odds proporcionales (logística ordinal)
term	OR	std.error	statistic	conf.low	conf.high	coef.type
AGE_z	1.7086	0.0793	6.7526	1.4650	1.9999	coefficient
SMK	2.4329	0.1595	5.5728	1.7828	3.3332	coefficient
DOSE_z	3.3124	0.0918	13.0436	2.7764	3.9802	coefficient
Leve\|Moderada	0.5463	0.1052	-5.7450	NA	NA	scale
Moderada\|Severa	3.5161	0.1143	10.9985	NA	NA	scale

12 Datos Correlacionados: GEE

12.1 Estimación GEE — ecuaciones y matriz de correlación

Las Ecuaciones de Estimación Generalizadas (GEE; Liang & Zeger 1986) producen estimaciones poblacionales marginales robustas cuando hay datos correlacionados (medidas repetidas, conglomerados).

La ecuación de estimación GEE generaliza el score del MV:

\[\sum_{i=1}^{n} \frac{\partial \boldsymbol{\mu}_i^{\top}}{\partial \boldsymbol{\beta}} \mathbf{V}_i^{-1}(\mathbf{Y}_i - \boldsymbol{\mu}_i) = \mathbf{0}\]

donde \(\mathbf{V}_i = \phi \mathbf{A}_i^{1/2} \mathbf{R}(\boldsymbol{\alpha}) \mathbf{A}_i^{1/2}\) con \(\mathbf{R}(\boldsymbol{\alpha})\) = matriz de correlación de trabajo.

Estructura de correlación	Descripción	Cuándo usar
Independiente	\(R = \mathbf{I}\)	Sin correlación esperada
Intercambiable	\(\text{Cor}(Y_{it}, Y_{it'}) = \alpha\)	Grupos/clusters
AR(1)	\(\text{Cor}(Y_{it}, Y_{it'}) = \alpha^{\\|t-t'\\|}\)	Series de tiempo
No estructurada	\(R\) sin restricciones	Máxima flexibilidad

set.seed(88)
n_s <- 120; n_v <- 4
datos_long <- data.frame(
  ID   = rep(1:n_s, each = n_v),
  TRAT = rep(rbinom(n_s, 1, 0.5), each = n_v),
  TIEM = rep(0:(n_v-1), times = n_s),
  AGE  = rep(round(rnorm(n_s, 50, 10)), each = n_v)
)
logit_l <- -1.5 + 0.9 * datos_long$TRAT + 0.03 * datos_long$AGE -
            0.25 * datos_long$TIEM
datos_long$Y <- rbinom(nrow(datos_long), 1, plogis(logit_l))

mod_gee <- geeglm(Y ~ TRAT + AGE + TIEM, data = datos_long,
                   id = ID, family = binomial, corstr = "exchangeable")

tidy(mod_gee, exponentiate = TRUE) |>
  mutate(across(c(estimate, std.error, statistic, p.value), ~round(., 4))) |>
  rename(OR = estimate) |>
  tabla_doc(titulo = "Tabla 14.1 — Modelo GEE con correlación intercambiable (efectos marginales poblacionales)")

Tabla 14.1 — Modelo GEE con correlación intercambiable (efectos marginales poblacionales)
term	OR	std.error	statistic	p.value
(Intercept)	0.5226	0.5154	1.5852	0.2080
TRAT	2.9556	0.1970	30.2681	0.0000
AGE	1.0115	0.0094	1.5047	0.2200
TIEM	0.8297	0.0794	5.5252	0.0187

cat("\nCorrelación de trabajo estimada α̂ =",
    round(mod_gee$geese$alpha, 4),
    "\n(dos medidas del mismo sujeto tienen correlación", round(mod_gee$geese$alpha, 4), ")\n")

#> 
#> Correlación de trabajo estimada α̂ = 0.0117 
#> (dos medidas del mismo sujeto tienen correlación 0.0117 )

13 Formulario Matricial — Compendio Técnico

Formulario Completo — Regresión Logística con Álgebra Matricial
Componente	Expresión matemática
Vector de parámetros β	β = (α, β₁, …, βₖ)ᵀ ∈ ℝᵏ⁺¹
Matriz de diseño X	X = [1 \| X₁ \| ⋯ \| Xₖ], dimensión n×(k+1)
Predictor lineal η	η = Xβ [n×1]; ηᵢ = βᵀxᵢ = α + Σβⱼxᵢⱼ
Función logística (probabilidades p)	p = σ(η) = 1/(1+e⁻ᵑ) [n×1]
Transformación logit	logit[P(X)] = ln[p/(1−p)] = βᵀx
Log-verosimilitud ℓ(β)	ℓ(β) = Dᵀη − 1ᵀln(1+eᵑ) [escalar]
Score (gradiente ∂ℓ/∂β)	∂ℓ/∂β = Xᵀ(D−p) = 0 [(k+1)×1]
Matriz de pesos W	W = diag(pᵢ(1−pᵢ)) [n×n]
Información de Fisher I(β)	I(β) = XᵀWX [(k+1)×(k+1)]
Varianza estimada Var(β̂)	\|Var(β̂) = I(β̂)⁻¹ = (XᵀŴX)⁻¹ [(k+1)×(k+1)]
Actualización IRLS (Newton)	β̂⁽ᵗ⁺¹⁾ = β̂⁽ᵗ⁾ + (XᵀWX)⁻¹ Xᵀ(D−p)
OR para comparación x₁ vs x₀	OR = exp[βᵀ(x₁−x₀)] = exp[Σβᵢ(x₁ᵢ−x₀ᵢ)]
OR para variable binaria Xᵢ	OR = exp(βᵢ) [si Xᵢ ∈ {0,1}]
OR para variable continua (por Δ unidades)	OR(Δ) = exp(βᵢ·Δ)
Prueba de Wald (por coef.)	Z = β̂ᵢ / SE(β̂ᵢ) ~ N(0,1); SE = √[Var(β̂)ᵢᵢ]
Prueba LR (test de modelos anidados)	LR = −2[ℓ(β̂ᵣ) − ℓ(β̂꜀)] ~ χ²_q (q = Δparams)
IC Wald 95% para OR	[exp(β̂ᵢ − 1.96·SEᵢ), exp(β̂ᵢ + 1.96·SEᵢ)]
IC Perfil de verosimilitud	{β : −2[ℓ(β*) − ℓ(β̂)] ≤ χ²₁,₀.₉₅}

Regresión Logística

19 de junio de 2026