13 Verificación de modelos

¿Ajustó bien? ¿Predice bien? ¿Sobran variables?

Resumen

Ajustamos un modelo lineal, uno de Poisson y uno logístico y les aplicamos todas las pruebas de bondad de ajuste y de capacidad predictiva que se usan en epidemiología: R², AIC, ji cuadrada, sobredispersión, Hosmer-Lemeshow, curva ROC, sensibilidad, especificidad, MAE, RMSE y WIS. Usamos un spline para descubrir un efecto no lineal y efectos marginales para interpretar. Cada número viene con su lectura, redactada como se redactaría en un artículo.

library(tidyverse)
library(haven)
library(broom)
library(splines) # para los splines naturales
library(car) # VIF y pruebas de supuestos
library(AER) # prueba de sobredispersión
library(MASS) # binomial negativa
library(pROC) # curva ROC
library(ResourceSelection) # Hosmer-Lemeshow
library(marginaleffects) # efectos marginales
library(yardstick) # métricas de clasificación
library(performance) # batería de índices de ajuste
library(nnet) # regresión multinomial
library(scoringutils) # WIS
library(forecast) # ARIMA

# Evita que dplyr::select choque con MASS::select
select <- dplyr::select

El precio de usar muchos paquetes: los conflictos

Este capítulo carga quince paquetes, y varios exportan funciones con el mismo nombre. La que gana es la del paquete que se cargó al último.

El caso concreto que nos mordió aquí: forecast carga generics, que trae su propia función accuracy(). Como se carga después de yardstick, accuracy deja de ser la de yardstick y metric_set() truena con un mensaje que no ayuda nada.

La solución es escribir el paquete y dos puntos en las funciones que se prestan a confusión:

yardstick::accuracy(...) # así no hay duda de cuál estás usando

Para ver qué se está pisando en tu sesión:

conflicted::conflict_scout()

Usar muchos paquetes está bien y ahorra código, pero el precio son los conflictos de nombres. Si una función deja de funcionar de la nada, revisa primero si algo la está enmascarando.

Las dos preguntas que nunca hay que confundir

Verificar un modelo es en realidad dos trabajos distintos:

¿Ajusta bien? (bondad de ajuste). ¿El modelo es compatible con los datos que ya vio? ¿Se cumplen sus supuestos? Aquí viven la R², la ji cuadrada, la prueba de Hosmer-Lemeshow y los residuales.
¿Predice bien? (capacidad predictiva). ¿Le atina a datos que no vio? Aquí viven el MAE, el RMSE, el área bajo la curva ROC y el WIS.

Un modelo puede ajustar perfecto y predecir pésimo (le pasa a cualquier modelo con demasiadas variables), y puede predecir bien con supuestos violados. Son preguntas distintas y se contestan con herramientas distintas.

13.1 Los datos

Tres bases, cada una para un tipo de modelo:

# Ataques de asma: 120 pacientes (conteo -> Poisson)
asma <- read_csv("datasets/asthma.csv", show_col_types = FALSE) |>
  mutate(
    gender = factor(gender),
    res_inf = factor(res_inf, levels = c("no", "yes"))
  )

# Evento vascular cerebral: 226 pacientes (sí/no -> logística)
acv <- read_dta("datasets/stroke.dta") |>
  mutate(across(where(is.labelled), as_factor)) |>
  mutate(muerte = as.integer(status == "dead"))

# Dengue semanal en Puerto Rico (serie de tiempo)
dengue <- read_csv(
  "datasets/dengue_puerto_rico_semanal.csv",
  show_col_types = FALSE
)

Base	n	Desenlace	Modelo
`asma`	120	`attack`: ataques de asma al año (0 a 9)	Poisson
`acv`	226	`status`: vivo o muerto al egreso	Logística
`coronaria`	200	`pa_cat`: presión normal, elevada o hipertensión	Multinomial
`dengue`	1085	`casos`: casos semanales	ARIMA

Las dos primeras vienen de Data Analysis in Medicine and Health using R.

13.2 Modelo lineal: presión sistólica al ingreso

Empecemos con el más sencillo. ¿La presión sistólica al ingreso (sbp) depende de la edad, el sexo, la diabetes y el tipo de evento?

modelo_lineal <- lm(sbp ~ age + sex + dm + stroke_type, data = acv)

tidy(modelo_lineal, conf.int = TRUE) |>
  mutate(across(where(is.numeric), ~ round(.x, 3)))

# A tibble: 5 × 7
  term                   estimate std.error statistic p.value conf.low conf.high
  <chr>                     <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)             146.       11.5      12.7     0      123.      168.   
2 age                       0.287     0.167     1.72    0.086   -0.041     0.615
3 sexfemale                 5.12      4.63      1.11    0.27    -4.01     14.2  
4 dmyes                    -2.62      4.71     -0.556   0.579  -11.9       6.66 
5 stroke_typeHaemorrhag…   -3.11      4.89     -0.636   0.526  -12.7       6.53

13.2.1 `R²` y `R²` ajustada

glance(modelo_lineal) |>
  select(r.squared, adj.r.squared, statistic, p.value, AIC, BIC) |>
  mutate(across(everything(), ~ round(.x, 4)))

# A tibble: 1 × 6
  r.squared adj.r.squared statistic p.value   AIC   BIC
      <dbl>         <dbl>     <dbl>   <dbl> <dbl> <dbl>
1    0.0251        0.0074      1.42   0.228 2246. 2266.

Cómo se lee cada número

r.squared = 0.025. El modelo explica el 2.5% de la variabilidad de la presión sistólica. El otro 97.5% se debe a cosas que no medimos. Se lee así: “si te digo la edad, el sexo, la diabetes y el tipo de evento de un paciente, apenas mejoro un 2.5% mi capacidad de adivinar su presión respecto a decir siempre el promedio”.
adj.r.squared = 0.007. La R² ajustada castiga por cada variable que metes. Si al agregar una variable la R² sube pero la ajustada baja, esa variable no está aportando: sólo está usando grados de libertad. Aquí la caída de 0.025 a 0.007 es una señal de alarma.
statistic = 1.42 con p.value = 0.23. Es la prueba F global: “¿este modelo, en conjunto, sirve de algo comparado con no tener modelo?”. Con p = 0.23 la respuesta es no. Es la prueba más importante de la tabla y la que más se ignora.

13.2.2 Error de predicción: `MAE`, `MSE` y `RMSE`

res <- residuals(modelo_lineal)

c(
  MAE = mean(abs(res)),
  MSE = mean(res^2),
  RMSE = sqrt(mean(res^2))
) |> round(2)

    MAE     MSE    RMSE 
  26.81 1149.09   33.90

Los tres errores, en español

Los tres miden lo mismo (qué tan lejos quedó el modelo) pero castigan distinto:

MAE (error absoluto medio) = 26.8 mmHg. El más fácil de explicar: en promedio, el modelo se equivoca por 26.8 mmHg. Está en las mismas unidades del desenlace, así que se puede poner tal cual en un resultado.
MSE (error cuadrático medio) = 1150. Eleva los errores al cuadrado, así que castiga mucho más los errores grandes. Su problema: está en mmHg², que no significa nada para nadie.
RMSE = 33.9 mmHg. Es la raíz del MSE, así que regresa a las unidades originales y sigue castigando los errores grandes.

La comparación entre MAE y RMSE te dice algo. Si RMSE es mucho mayor que MAE (aquí 33.9 contra 26.8), es porque hay unos pocos errores muy grandes jalando el promedio. Si fueran parecidos, los errores estarían repartidos parejo.

13.2.3 Supuestos: normalidad, homocedasticidad y colinealidad

# ¿Los residuales son normales?
shapiro.test(res)


    Shapiro-Wilk normality test

data:  res
W = 0.99046, p-value = 0.1436

# ¿La varianza es constante? (Breusch-Pagan)
ncvTest(modelo_lineal)

Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.0195577, Df = 1, p = 0.88878

# ¿Hay variables que se explican entre sí?
vif(modelo_lineal) |> round(3)

        age         sex          dm stroke_type 
      1.037       1.010       1.013       1.033

La lección más importante de este capítulo

Mira lo que acaba de pasar. Todos los supuestos se cumplen:

Prueba	Resultado	Lectura
Shapiro-Wilk	`p` = 0.14	No se rechaza la normalidad de los residuales
Breusch-Pagan	`p` = 0.89	La varianza es constante
`VIF`	todos ≈ 1	Sin colinealidad

Y sin embargo el modelo no sirve: R² de 0.025 y prueba F global con p = 0.23.

Que un modelo cumpla sus supuestos no quiere decir que sea útil. Los supuestos dicen si el modelo está bien construido; la R² y la F dicen si tiene algo que decir. Un modelo puede estar impecablemente construido sobre variables que no explican nada.

Cómo se reportaría esto en un artículo:

La presión sistólica al ingreso no se asoció de manera significativa con la edad, el sexo, el antecedente de diabetes ni el tipo de evento vascular (F(4,221) = 1.42; p = 0.23; R² ajustada = 0.007). Los supuestos de normalidad de residuales (Shapiro-Wilk p = 0.14) y homocedasticidad (Breusch-Pagan p = 0.89) se cumplieron.

Fíjate que se reporta el resultado nulo con todo y sus pruebas de supuestos. Un resultado nulo bien hecho es un resultado.

13.3 Modelo de Poisson: ataques de asma

Ahora un conteo: número de ataques de asma al año, según el puntaje de malestar psicológico GHQ-12, si hay infección respiratoria recurrente y el sexo.

modelo_poisson <- glm(attack ~ ghq12 + res_inf + gender,
  data = asma, family = poisson()
)

exp(cbind(IRR = coef(modelo_poisson), confint.default(modelo_poisson))) |>
  round(3)

              IRR 2.5 % 97.5 %
(Intercept) 0.730 0.509  1.045
ghq12       1.051 1.035  1.067
res_infyes  1.532 1.135  2.067
gendermale  0.959 0.754  1.219

Cómo se lee una razón de tasas de incidencia (IRR)

ghq12 = 1.051. Por cada punto más en la escala GHQ-12, el número esperado de ataques se multiplica por 1.051, es decir aumenta 5.1%. Como el intervalo (1.035–1.067) no incluye al 1, es estadísticamente significativo.
res_infyes = 1.532. Quienes tienen infección respiratoria recurrente tienen 53% más ataques que quienes no, ajustando por GHQ-12 y sexo.
gendermale = 0.959, intervalo 0.754–1.219: incluye al 1, así que no hay evidencia de diferencia por sexo.

En Poisson siempre se exponencia. Los coeficientes crudos están en escala logarítmica y no se interpretan directamente.

13.3.1 Bondad de ajuste: ¿se cumple el supuesto de Poisson?

La Poisson supone que la varianza es igual a la media. Hay que revisarlo siempre.

# Los datos crudos: ¿varianza = media?
c(media = mean(asma$attack), varianza = var(asma$attack)) |> round(3)

   media varianza 
   2.458    4.049

A primera vista hay problema: la varianza (4.05) es 1.65 veces la media (2.46). Pero eso es la dispersión cruda. Lo que le importa al modelo es la dispersión después de ajustar por las covariables:

# Tres formas de medir lo mismo
c(
  deviance_gl = modelo_poisson$deviance / modelo_poisson$df.residual,
  pearson_gl = sum(residuals(modelo_poisson, "pearson")^2) /
    modelo_poisson$df.residual
) |> round(3)

deviance_gl  pearson_gl 
      1.178       1.061

# Prueba formal de sobredispersión
dispersiontest(modelo_poisson)


    Overdispersion test

data:  modelo_poisson
z = 0.43048, p-value = 0.3334
alternative hypothesis: true dispersion is greater than 1
sample estimates:
dispersion 
  1.046941

Sobredispersión cruda ≠ sobredispersión condicional

Los datos crudos parecían sobredispersos (1.65) y el modelo no lo está (1.05, p = 0.33).

¿Por qué? Porque las covariables explicaron ese exceso de varianza. La gente con GHQ-12 alto tiene muchos ataques y la de GHQ-12 bajo tiene pocos; visto todo junto eso se ve como “mucha variabilidad”, pero una vez que el modelo sabe el GHQ-12 de cada quien, lo que sobra es justo lo que la Poisson espera.

Nunca decidas usar binomial negativa mirando la media y la varianza crudas. Se decide con los residuales del modelo ajustado.

13.3.2 La prueba ji cuadrada de bondad de ajuste

# ¿La devianza residual es compatible con una ji cuadrada?
pchisq(modelo_poisson$deviance, modelo_poisson$df.residual,
  lower.tail = FALSE
) |> round(4)

[1] 0.0922

Tip

Ésta es la prueba de bondad de ajuste por devianza. La hipótesis nula es que el modelo ajusta bien, así que aquí un p grande es buena noticia (al revés de lo acostumbrado).

Con p = 0.09 no se rechaza el ajuste, aunque está en la frontera. Se reporta así: “no se encontró evidencia de falta de ajuste (prueba de devianza, p = 0.09)”.

13.3.3 ¿Y si usáramos binomial negativa?

modelo_nb <- glm.nb(attack ~ ghq12 + res_inf + gender, data = asma)

c(
  AIC_Poisson = AIC(modelo_poisson),
  AIC_BinomialNegativa = AIC(modelo_nb)
) |> round(1)

         AIC_Poisson AIC_BinomialNegativa 
               417.7                419.7

Nota

El AIC de la Poisson (417.7) es menor que el de la binomial negativa (419.7), así que la Poisson gana. Tiene sentido: la binomial negativa agrega un parámetro extra para la sobredispersión, y aquí no hay sobredispersión que modelar, así que sólo paga el costo del parámetro sin ganancia.

Regla del AIC: entre más chico, mejor, y una diferencia menor a 2 se considera empate. Aquí la diferencia es exactamente 2, lo que refuerza quedarse con el modelo más simple.

13.3.4 Un spline: ¿el efecto del GHQ-12 es realmente lineal?

El modelo supone que cada punto de GHQ-12 pesa igual, tanto del 2 al 3 como del 30 al 31. Eso es una suposición fuerte y se puede probar. Un spline natural permite que la curva tome la forma que quiera:

modelo_spline <- glm(attack ~ ns(ghq12, df = 4) + res_inf + gender,
  data = asma, family = poisson()
)

anova(modelo_poisson, modelo_spline, test = "LRT")

Analysis of Deviance Table

Model 1: attack ~ ghq12 + res_inf + gender
Model 2: attack ~ ns(ghq12, df = 4) + res_inf + gender
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1       116     136.68                          
2       113     115.85  3   20.828 0.0001143 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

c(
  AIC_lineal = AIC(modelo_poisson),
  AIC_spline = AIC(modelo_spline)
) |> round(1)

AIC_lineal AIC_spline 
     417.7      402.9

Sí hay no linealidad, y es fuerte

La prueba de razón de verosimilitudes da p = 0.0001 y el AIC baja de 417.7 a 402.9 (casi 15 puntos: una diferencia enorme).

Conclusión: el efecto del GHQ-12 no es una línea recta. El modelo lineal estaba promediando formas distintas en un solo número.

Veamos la forma que tiene:

malla <- expand_grid(
  ghq12 = seq(0, 33, by = 0.5),
  res_inf = factor("no", levels = c("no", "yes")),
  gender = factor("female", levels = c("female", "male"))
)

pred_spl <- predictions(modelo_spline, newdata = malla)

ggplot(pred_spl, aes(x = ghq12, y = estimate)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
    fill = "#003f5c", alpha = 0.2
  ) +
  geom_line(color = "#003f5c", linewidth = 1) +
  geom_rug(data = asma, aes(x = ghq12, y = NULL), alpha = 0.3) +
  labs(
    x = "Puntaje GHQ-12 (malestar psicológico)",
    y = "Ataques de asma esperados al año",
    title = "El efecto del malestar psicológico no es lineal",
    caption = paste(
      "Mujeres sin infección respiratoria recurrente.",
      "Las marcas de abajo son los datos observados."
    )
  ) +
  theme_bw()

Efecto del GHQ-12 sobre el número esperado de ataques, con spline natural.

Y la pendiente local en distintos puntos de la escala:

puntos <- data.frame(
  ghq12 = c(3, 9, 15, 21, 27),
  res_inf = factor("no", levels = c("no", "yes")),
  gender = factor("female", levels = c("female", "male"))
)

slopes(modelo_spline, newdata = puntos, variables = "ghq12") |>
  as.data.frame() |>
  select(ghq12, estimate, conf.low, conf.high, p.value) |>
  mutate(across(-1, ~ round(.x, 4)))

  ghq12 estimate conf.low conf.high p.value
1     3  -0.0854  -0.1828    0.0119  0.0855
2     9   0.0724   0.0361    0.1086  0.0001
3    15   0.3219   0.1632    0.4805  0.0001
4    21  -0.0073  -0.2092    0.1945  0.9434
5    27   0.0024  -0.1430    0.1479  0.9741

Cómo se interpreta esto en un artículo

La curva y las pendientes locales cuentan una historia de umbral:

Por debajo de GHQ-12 ≈ 6 la curva es plana: el malestar psicológico bajo no se asocia a más ataques.
Entre 6 y 18 la curva sube con fuerza y las pendientes son significativas: ahí es donde el malestar sí importa.
Arriba de 18 vuelve a aplanarse: ya no empeora más.

Redactado como en un artículo:

Se modeló el puntaje GHQ-12 mediante un spline natural con 4 grados de libertad. El término no lineal resultó significativo frente al modelo lineal (prueba de razón de verosimilitudes, p = 0.0001; ΔAIC = 14.8), lo que indica una relación no lineal. La tasa de ataques permaneció estable en puntajes bajos, aumentó de forma marcada entre 6 y 18 puntos y se estabilizó posteriormente, lo que sugiere un efecto de umbral.

El modelo lineal habría reportado “5.1% más ataques por punto” y eso es engañoso: no aplica ni abajo ni arriba de la escala, sólo en medio.

El spline también cambió otra conclusión

Compara la IRR de infección recurrente en los dos modelos:

                  res_infyes 2.5 % 97.5 %
Modelo lineal          1.532 1.135  2.067
Modelo con spline      1.257 0.922  1.713

Pasó de 1.53 (significativa) a 1.26 (no significativa). Lo que ocurría es que el modelo lineal, al no capturar bien el efecto del GHQ-12, le estaba atribuyendo parte de ese efecto a la infección recurrente.

Especificar mal una variable de ajuste sesga los coeficientes de las demás. Es una forma de confusión residual, y es de las más comunes y menos revisadas.

13.3.5 Efectos marginales

Los coeficientes de un modelo de Poisson están en escala multiplicativa. Los efectos marginales los traducen a la escala del desenlace: “¿cuántos ataques más, en números?”.

avg_slopes(modelo_spline)


    Term      Contrast Estimate Std. Error     z Pr(>|z|)   S   2.5 % 97.5 %
 gender  male - female   0.1286     0.3117 0.413    0.680 0.6 -0.4823 0.7396
 ghq12   dY/dX           0.0152     0.0168 0.904    0.366 1.5 -0.0178 0.0482
 res_inf yes - no        0.5321     0.3492 1.524    0.128 3.0 -0.1524 1.2166

Type: response

Ojo con el efecto marginal promedio en modelos no lineales

El efecto marginal promedio del GHQ-12 sale no significativo (p ≈ 0.37) aunque el spline sí lo es (p = 0.0001). No es contradicción:

El efecto marginal promedio saca un promedio de la pendiente a lo largo de toda la curva. Como la curva sube fuerte en medio y es plana en los extremos, al promediarla los efectos se cancelan y queda un número chico y con mucha incertidumbre.

Cuando el efecto no es lineal, el efecto marginal promedio esconde justo lo interesante. Reporta la curva y las pendientes locales, no un solo número.

13.4 Modelo logístico: mortalidad hospitalaria

Ahora un desenlace binario: si el paciente falleció durante la hospitalización.

modelo_logistico <- glm(muerte ~ gcs + age + sex + dm + stroke_type,
  data = acv, family = binomial()
)

exp(cbind(OR = coef(modelo_logistico), confint.default(modelo_logistico))) |>
  round(3)

                           OR 2.5 % 97.5 %
(Intercept)             0.979 0.077 12.394
gcs                     0.720 0.645  0.803
age                     1.025 0.995  1.056
sexfemale               1.537 0.654  3.613
dmyes                   1.599 0.682  3.754
stroke_typeHaemorrhagic 3.522 1.506  8.238

Interpretación

gcs = 0.720. Por cada punto más en la escala de coma de Glasgow, los momios de morir se multiplican por 0.72, o sea bajan 28%. Es protector y su intervalo (0.645–0.803) no toca el 1.
stroke_typeHaemorrhagic = 3.522. Los eventos hemorrágicos tienen 3.5 veces los momios de muerte que los isquémicos, ajustando por lo demás. Intervalo 1.51–8.24: significativo pero muy ancho, porque son sólo 226 pacientes.
age, sex, dm: todos con intervalos que incluyen el 1. Sin evidencia de asociación.

13.4.1 Bondad de ajuste global

El paquete performance saca de un solo golpe casi todos los índices de ajuste, así no hay que recordar la fórmula de cada uno:

model_performance(modelo_logistico)

# Indices of model performance

AIC   |  AICc |   BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log
------------------------------------------------------------------------
171.4 | 171.7 | 191.9 |     0.410 | 0.332 |     1 |    0.353 |   -22.043

AIC   | Score_spherical |   PCP
-------------------------------
171.4 |           0.034 | 0.783

# Pseudo R2 de McFadden (no viene en la tabla de arriba)
1 - modelo_logistico$deviance / modelo_logistico$null.deviance

[1] 0.3646695

# ¿Qué variable aporta? Se agregan una por una
anova(modelo_logistico, test = "LRT")

Analysis of Deviance Table

Model: binomial, link: logit

Response: muerte

Terms added sequentially (first to last)

            Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                          225     250.83              
gcs          1   79.907       224     170.92 < 2.2e-16 ***
age          1    1.391       223     169.53  0.238311    
sex          1    0.783       222     168.75  0.376252    
dm           1    1.020       221     167.73  0.312478    
stroke_type  1    8.368       220     159.36  0.003818 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Qué es cada índice de model_performance

AIC / AICc / BIC: para comparar modelos. Entre más chico, mejor. El AICc corrige por muestras chicas y el BIC castiga más los parámetros.
Tjur's R2 (0.410): el pseudo-R² más intuitivo para logística. Es la diferencia entre la probabilidad promedio predicha a quienes murieron y a quienes no. Entre más separadas, mejor discrimina.
RMSE y Log_loss: error de las probabilidades predichas.
PCP (0.783): predicted correct probability, qué tan bien clasifica en promedio.

Nota

La pseudo R² de McFadden (0.365) no se lee como la R² de siempre. En regresión logística no existe una “proporción de varianza explicada”. Valores de 0.2 a 0.4 ya se consideran un ajuste excelente. Un 0.365 aquí es muy bueno.

La tabla de anova con test = "LRT" agrega las variables una por una y prueba si cada una aporta. Es la forma correcta de evaluar variables categóricas con varios niveles, donde mirar los p uno por uno no sirve.

13.4.2 Calibración: la prueba de Hosmer-Lemeshow

hoslem.test(modelo_logistico$y, fitted(modelo_logistico), g = 10)


    Hosmer and Lemeshow goodness of fit (GOF) test

data:  modelo_logistico$y, fitted(modelo_logistico)
X-squared = 5.5516, df = 8, p-value = 0.6973

Qué mide y cómo se lee

Hosmer-Lemeshow parte a los pacientes en 10 grupos según su riesgo predicho y compara, en cada grupo, cuántas muertes predijo el modelo contra cuántas hubo.

La hipótesis nula es que el modelo está bien calibrado, así que —igual que con la devianza— un p grande es buena noticia.

Aquí χ² = 5.55 con 8 grados de libertad y p = 0.70: no hay evidencia de mala calibración. Cuando el modelo dice “30% de riesgo”, muere alrededor del 30%.

Calibración y discriminación son cosas distintas. Un modelo puede ordenar perfectamente a los pacientes (buena discriminación) y aun así equivocarse en la magnitud del riesgo (mala calibración). Hay que reportar las dos.

13.4.3 Discriminación: la curva ROC

roc_obj <- roc(modelo_logistico$y, fitted(modelo_logistico), quiet = TRUE)

ggroc(roc_obj, color = "#003f5c", linewidth = 1) +
  geom_abline(intercept = 1, slope = 1, linetype = "dashed", color = "gray50") +
  annotate("text",
    x = 0.3, y = 0.2,
    label = paste0("AUC = ", round(auc(roc_obj), 3)), size = 5
  ) +
  labs(
    x = "Especificidad", y = "Sensibilidad",
    title = "¿Qué tan bien distingue el modelo quién muere?"
  ) +
  theme_bw()

ci.auc(roc_obj) |> round(4)

[1] 0.8435 0.8915 0.9395

El AUC y su intervalo

AUC = 0.892 (IC 95%: 0.844–0.940). Se interpreta así: si tomas al azar un paciente que murió y uno que sobrevivió, el modelo le asigna mayor riesgo al que murió el 89% de las veces.

`AUC`	Lectura
0.5	El modelo no discrimina (equivale a una moneda)
0.7–0.8	Aceptable
0.8–0.9	Bueno
> 0.9	Excelente

Siempre reporta el intervalo de confianza. Un AUC de 0.89 con IC de 0.84–0.94 es sólido; el mismo 0.89 con IC de 0.60–0.98 no dice nada.

13.4.4 Sensibilidad, especificidad y compañía

Con yardstick no hace falta escribir ninguna fórmula a mano: se arma una tabla con la verdad, la probabilidad y la clase predicha, y el paquete calcula todo.

resultados <- tibble(
  verdad = factor(
    if_else(acv$muerte == 1, "muere", "vive"),
    levels = c("muere", "vive")
  ),
  prob = fitted(modelo_logistico)
) |>
  mutate(
    clase = factor(
      if_else(prob >= 0.5, "muere", "vive"),
      levels = c("muere", "vive")
    )
  )

conf_mat(resultados, truth = verdad, estimate = clase)

          Truth
Prediction muere vive
     muere    29   11
     vive     26  160

El error que le pasa a todo el mundo con yardstick

glm modela la probabilidad del segundo nivel del factor, mientras que yardstick toma el primero como el evento de interés.

Si no los alineas, la sensibilidad y la especificidad salen intercambiadas y el AUC sale por debajo de 0.5. Por eso arriba pusimos "muere" como primer nivel de forma explícita.

Si tu AUC sale menor a 0.5, casi nunca es que el modelo sea malísimo: es que invertiste las etiquetas.

metricas_clasificacion <- metric_set(
  yardstick::sensitivity, yardstick::specificity,
  yardstick::precision, yardstick::npv,
  yardstick::accuracy, yardstick::f_meas
)

metricas_clasificacion(resultados, truth = verdad, estimate = clase)

# A tibble: 6 × 3
  .metric     .estimator .estimate
  <chr>       <chr>          <dbl>
1 sensitivity binary         0.527
2 specificity binary         0.936
3 precision   binary         0.725
4 npv         binary         0.860
5 accuracy    binary         0.836
6 f_meas      binary         0.611

Aquí está la trampa más común en modelos de predicción clínica

La exactitud es 0.836: el modelo acierta en el 84% de los pacientes. Suena excelente.

Pero mira la sensibilidad: 0.527. De los 55 pacientes que murieron, el modelo sólo identificó a 29. Se le escapó casi la mitad de las muertes.

¿Cómo puede tener 84% de exactitud entonces? Porque el 76% de los pacientes sobrevivió, y el modelo es buenísimo detectando a los que sobreviven (especificidad 0.936). La exactitud está inflada por la clase mayoritaria.

En desenlaces poco frecuentes, la exactitud es una métrica engañosa. Si el 95% de tus pacientes sobrevive, un modelo que diga “todos viven” tiene 95% de exactitud y es completamente inútil. Reporta siempre sensibilidad y especificidad por separado.

Los nombres que vas a encontrar en la literatura:

En epidemiología	En ciencia de datos	Qué contesta
Sensibilidad	Recall	De los que murieron, ¿a cuántos detecté?
Especificidad	—	De los que vivieron, ¿a cuántos descarté bien?
Valor predictivo positivo	Precision	De los que marqué, ¿cuántos sí murieron?
Valor predictivo negativo	—	De los que descarté, ¿cuántos sí vivieron?

13.4.5 El punto de corte no es 0.5

El 0.5 que usamos es una convención, no un resultado estadístico. El índice de Youden busca el corte que maximiza sensibilidad + especificidad:

corte <- coords(roc_obj, "best", ret = c("threshold", "sensitivity", "specificity"))
corte |> mutate(across(everything(), ~ round(.x, 3)))

  threshold sensitivity specificity
1     0.195       0.855       0.813

resultados |>
  mutate(
    clase = factor(
      if_else(prob >= corte$threshold, "muere", "vive"),
      levels = c("muere", "vive")
    )
  ) |>
  metricas_clasificacion(truth = verdad, estimate = clase)

# A tibble: 6 × 3
  .metric     .estimator .estimate
  <chr>       <chr>          <dbl>
1 sensitivity binary         0.855
2 specificity binary         0.813
3 precision   binary         0.595
4 npv         binary         0.946
5 accuracy    binary         0.823
6 f_meas      binary         0.701

Compara las dos tablas

Corte	Sensibilidad	Especificidad	Exactitud
0.50	0.527	0.936	0.836
0.195 (Youden)	0.855	0.766	0.788

Bajando el corte, la sensibilidad sube de 53% a 86% a costa de bajar la especificidad y la exactitud global.

¿Cuál está bien? Depende de qué cuesta cada error, y ésa es una decisión clínica y de política, no estadística:

Si el modelo sirve para decidir a quién vigilar más de cerca, no detectar una muerte es gravísimo y una falsa alarma es barata → corte bajo.
Si sirviera para decidir a quién NO tratar, sería al revés.

Es la misma idea del capítulo de modelos: el punto de corte lo pones tú, y hay que justificarlo en la sección de métodos.

13.4.6 El puntaje de Brier

brier_class(resultados, truth = verdad, prob)

# A tibble: 1 × 3
  .metric     .estimator .estimate
  <chr>       <chr>          <dbl>
1 brier_class binary         0.110

Nota

El puntaje de Brier (0.110) es el error cuadrático medio de las probabilidades predichas. Va de 0 (perfecto) a 0.25 (inútil, cuando predices 0.5 siempre).

Su virtud es que combina calibración y discriminación en un solo número. Su defecto es que, justo por eso, no te dice cuál de las dos falló.

13.4.7 Efectos marginales del modelo logístico

avg_slopes(modelo_logistico)


        Term                        Contrast Estimate Std. Error      z
 age         dY/dX                            0.00262    0.00162  1.621
 dm          yes - no                         0.04963    0.04516  1.099
 gcs         dY/dX                           -0.03519    0.00420 -8.371
 sex         female - male                    0.04596    0.04639  0.991
 stroke_type Haemorrhagic - Ischaemic Stroke  0.15981    0.06157  2.595
 Pr(>|z|)    S     2.5 %   97.5 %
  0.10509  3.3 -0.000548  0.00579
  0.27175  1.9 -0.038877  0.13814
  < 0.001 54.0 -0.043425 -0.02695
  0.32177  1.6 -0.044957  0.13688
  0.00945  6.7  0.039127  0.28049

Type: response

De momios a puntos porcentuales

Aquí está la razón por la que valen tanto la pena los efectos marginales. El OR del Glasgow era 0.72, que no le dice nada útil a un clínico. El efecto marginal promedio dice:

Cada punto adicional en la escala de Glasgow se asocia con una reducción absoluta de 3.5 puntos porcentuales en la probabilidad de morir (IC 95%: 2.7 a 4.3).

Eso sí se entiende y sí se puede usar para decidir. Es la traducción de la escala del modelo a la escala del mundo.

13.5 Modelo multinomial: categorías de presión arterial

Cuando el desenlace tiene más de dos categorías sin orden obligado, la logística se queda corta. La logística multinomial ajusta una ecuación por cada categoría contra una de referencia, tal como lo vimos en la forma de una regresión.

Usaremos la base coronary y clasificaremos la presión arterial en tres niveles:

coronaria <- read_dta("datasets/coronary.dta") |>
  mutate(across(where(is.labelled), as_factor)) |>
  mutate(
    pa_cat = case_when(
      sbp >= 140 | dbp >= 90 ~ "Hipertension",
      sbp >= 120 | dbp >= 80 ~ "Elevada",
      .default = "Normal"
    ),
    pa_cat = factor(
      pa_cat,
      levels = c("Normal", "Elevada", "Hipertension")
    )
  )

count(coronaria, pa_cat)

# A tibble: 3 × 2
  pa_cat           n
  <fct>        <int>
1 Normal          56
2 Elevada         70
3 Hipertension    74

Nota

El primer nivel del factor (Normal) será la categoría de referencia: todos los resultados se leen contra ella. Elegirla no es un detalle técnico, es una decisión de interpretación: se pone la categoría que sirve como comparación natural, casi siempre la sana o la más frecuente.

modelo_multi <- multinom(
  pa_cat ~ age + bmi + chol + gender + race,
  data = coronaria
)

tidy(modelo_multi, exponentiate = TRUE, conf.int = TRUE) |>
  filter(term != "(Intercept)") |>
  select(y.level, term, estimate, conf.low, conf.high, p.value) |>
  mutate(across(where(is.numeric), \(x) round(x, 3)))

# A tibble: 12 × 6
   y.level      term        estimate conf.low conf.high p.value
   <chr>        <chr>          <dbl>    <dbl>     <dbl>   <dbl>
 1 Elevada      age            1.04     0.942     1.14    0.46 
 2 Elevada      bmi            1.11     0.951     1.29    0.191
 3 Elevada      chol           1.31     0.923     1.86    0.131
 4 Elevada      genderman      0.634    0.304     1.32    0.224
 5 Elevada      racechinese    1.20     0.423     3.40    0.733
 6 Elevada      raceindian     1.07     0.188     6.07    0.941
 7 Hipertension age            1.08     0.967     1.20    0.176
 8 Hipertension bmi            0.863    0.74      1.01    0.06 
 9 Hipertension chol           1.81     1.24      2.64    0.002
10 Hipertension genderman      0.421    0.187     0.947   0.037
11 Hipertension racechinese    1.74     0.551     5.48    0.345
12 Hipertension raceindian     2.30     0.367    14.4     0.375

Cómo se lee una tabla multinomial

Hay dos bloques: uno para Elevada y otro para Hipertension, y los dos se comparan contra Normal. Cada OR contesta: “¿cuánto cambian los momios de estar en esta categoría, en vez de estar en Normal?”.

Lo más informativo es leer cada variable a lo largo de las dos categorías:

Variable	`Elevada` vs Normal	`Hipertension` vs Normal
`chol`	1.31	1.81
`age`	1.04	1.08
`genderman`	0.63	0.42

Fíjate en el gradiente: el colesterol pasa de 1.31 a 1.81 y la edad de 1.04 a 1.08. Es decir, el efecto es más fuerte para la categoría más severa.

Ese patrón creciente es lo que en epidemiología se llama una relación dosis-respuesta a lo largo del desenlace, y es uno de los argumentos clásicos a favor de que una asociación sea causal (criterio de Bradford Hill). Es justamente lo que un modelo multinomial deja ver y una logística binaria esconde: si hubiéramos juntado “Elevada” e “Hipertensión” en una sola categoría, habríamos reportado un solo número intermedio y perdido el gradiente.

13.5.1 Verificación: ¿clasifica bien?

Con yardstick las métricas multiclase salen igual de fácil que las binarias:

resultados_multi <- coronaria |>
  select(verdad = pa_cat) |>
  mutate(clase = predict(modelo_multi))

conf_mat(resultados_multi, truth = verdad, estimate = clase)

              Truth
Prediction     Normal Elevada Hipertension
  Normal           26      15            4
  Elevada          20      32           14
  Hipertension     10      23           56

metric_set(yardstick::accuracy, yardstick::kap)(
  resultados_multi,
  truth = verdad, estimate = clase
)

# A tibble: 2 × 3
  .metric  .estimator .estimate
  <chr>    <chr>          <dbl>
1 accuracy multiclass     0.57 
2 kap      multiclass     0.345

Exactitud y kappa

Exactitud = 0.570. El modelo acierta la categoría en 57% de los pacientes. Parece poco, pero el punto de comparación no es 100%: es lo que lograrías adivinando. Como la categoría más grande (Hipertension) es el 37% de la muestra, decir “todos hipertensos” daría 37%. El modelo mejora eso.
Kappa de Cohen = 0.345. Ésta es la métrica correcta cuando hay varias categorías, porque descuenta los aciertos que se explican por azar. Se lee así:

Kappa Acuerdo

< 0.20 Pobre

0.21–0.40 Aceptable

0.41–0.60 Moderado

0.61–0.80 Bueno

> 0.80 Muy bueno

Un 0.345 es aceptable: el modelo aporta información real, pero no alcanza para decidir tratamiento por sí solo.

Kappa	Acuerdo
< 0.20	Pobre
0.21–0.40	Aceptable
0.41–0.60	Moderado
0.61–0.80	Bueno
> 0.80	Muy bueno

También se puede calcular un AUC multiclase:

probabilidades <- predict(modelo_multi, type = "probs") |>
  as_tibble() |>
  rename_with(\(x) paste0(".pred_", x))

bind_cols(resultados_multi, probabilidades) |>
  roc_auc(truth = verdad, starts_with(".pred_"))

# A tibble: 1 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 roc_auc hand_till      0.719

Tip

El AUC multiclase (0.719) usa el método de Hand y Till: promedia el AUC de todas las parejas de categorías. Se interpreta igual que el binario, pero ojo: 0.719 con tres categorías es mejor de lo que suena, porque el azar aquí también da 0.5 y separar tres grupos es más difícil que separar dos.

13.5.2 ¿Sobran variables?

modelo_multi_reducido <- multinom(
  pa_cat ~ age + chol,
  data = coronaria,
  trace = FALSE
)

tibble(
  Modelo = c("Completo", "Sólo edad y colesterol"),
  AIC = c(AIC(modelo_multi), AIC(modelo_multi_reducido)),
  Parametros = c(modelo_multi$edf, modelo_multi_reducido$edf)
) |>
  mutate(AIC = round(AIC, 1))

# A tibble: 2 × 3
  Modelo                   AIC Parametros
  <chr>                  <dbl>      <dbl>
1 Completo                409.         14
2 Sólo edad y colesterol  410.          6

Nota

Aquí el AIC decide. Si el modelo reducido tiene un AIC menor o casi igual con menos parámetros, gana el reducido: estabas pagando grados de libertad por variables que no aportaban.

Redactado para un artículo:

Se ajustó un modelo de regresión logística multinomial con la categoría normal como referencia. El colesterol se asoció de forma creciente con las categorías de presión arterial (OR 1.31 para presión elevada y 1.81 para hipertensión, frente a presión normal), lo que sugiere una relación dosis-respuesta. La capacidad discriminativa del modelo fue moderada (AUC multiclase de Hand-Till 0.72; kappa de Cohen 0.35).

13.6 Series de tiempo: dengue

Para pronósticos las métricas son otras, porque hay que evaluar contra datos futuros que el modelo no vio.

H <- 52
entrena <- head(dengue, nrow(dengue) - H)
prueba <- tail(dengue, H)

serie <- ts(entrena$casos, frequency = 52)
fourier_dentro <- fourier(serie, K = 3)
fourier_fuera <- fourier(serie, K = 3, h = H)

modelo_arima <- auto.arima(serie, seasonal = FALSE, xreg = fourier_dentro)
pron_arima <- forecast(modelo_arima, xreg = fourier_fuera, level = c(50, 80, 90, 95))

# Modelo ingenuo: "va a pasar lo mismo que hace un año"
ingenuo <- tail(entrena$casos, H)

13.6.1 `MAE`, `MSE`, `RMSE` y `MAPE`

metricas_punto <- function(pred, obs) {
  c(
    MAE = mean(abs(pred - obs)),
    MSE = mean((pred - obs)^2),
    RMSE = sqrt(mean((pred - obs)^2)),
    MAPE = mean(abs((pred - obs) / obs)) * 100
  )
}

rbind(
  ARIMA = metricas_punto(as.numeric(pron_arima$mean), prueba$casos),
  Ingenuo = metricas_punto(ingenuo, prueba$casos)
) |> round(2)

           MAE      MSE   RMSE  MAPE
ARIMA    98.46 17480.86 132.22 64.00
Ingenuo 101.25 18811.75 137.16 63.04

Las métricas pueden contradecirse, y hay que saber cuál creer

Fíjate bien:

Por MAE (98.5 contra 101.3), MSE y RMSE, gana el ARIMA.
Por MAPE (64.0 contra 63.0), gana el modelo ingenuo.

¿Quién tiene razón? El problema es del MAPE: divide el error entre el valor observado, así que las semanas con pocos casos dominan el promedio. Errarle por 5 casos cuando hubo 10 cuenta como 50% de error; errarle por 50 cuando hubo 500 cuenta como 10%. En vigilancia epidemiológica, donde hay semanas con muy pocos casos, el MAPE es casi siempre mala idea.

Escoge la métrica antes de ver los resultados y justifícala. Si eliges la métrica después, siempre vas a encontrar una que favorezca a tu modelo.

13.6.2 El `WIS`: la métrica que sí toma en cuenta la incertidumbre

Todas las anteriores comparan un solo número predicho contra el observado, e ignoran por completo si el modelo dijo estar seguro o no. El WIS (Weighted Interval Score) evalúa el intervalo completo.

data.table::setDTthreads(1)

niveles <- c(0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.975)

q_arima <- cbind(
  pron_arima$lower[, 4], pron_arima$lower[, 3],
  pron_arima$lower[, 2], pron_arima$lower[, 1],
  pron_arima$mean,
  pron_arima$upper[, 1], pron_arima$upper[, 2],
  pron_arima$upper[, 3], pron_arima$upper[, 4]
)
sd_ingenuo <- sd(diff(entrena$casos, lag = 52), na.rm = TRUE)
q_ingenuo <- sapply(niveles, function(q) ingenuo + qnorm(q) * sd_ingenuo)

armar <- function(q, nombre) {
  tibble(
    modelo = nombre,
    semana = rep(prueba$semana, times = length(niveles)),
    observado = rep(prueba$casos, times = length(niveles)),
    nivel = rep(niveles, each = H),
    prediccion = as.vector(q)
  )
}

evaluacion <- bind_rows(armar(q_arima, "ARIMA"), armar(q_ingenuo, "Ingenuo")) |>
  as_forecast_quantile(
    forecast_unit = c("modelo", "semana"),
    observed = "observado", predicted = "prediccion",
    quantile_level = "nivel"
  ) |>
  score() |>
  summarise_scores(by = "modelo")

evaluacion |>
  select(
    modelo, wis, dispersion, underprediction, overprediction,
    interval_coverage_50, interval_coverage_90
  ) |>
  mutate(across(-1, ~ round(.x, 3)))

    modelo    wis dispersion underprediction overprediction
    <char>  <num>      <num>           <num>          <num>
1:   ARIMA 54.068     11.902          42.154          0.012
2: Ingenuo 64.826     11.423          53.396          0.006
   interval_coverage_50 interval_coverage_90
                  <num>                <num>
1:                0.269                0.673
2:                0.346                0.673

Cómo se lee el WIS y por qué es la mejor métrica para vigilancia

El WIS cobra tres cosas a la vez: qué tan ancho es el intervalo (para que no hagas trampa diciendo “entre 0 y un millón”), qué tanto se quedó corto y qué tanto se pasó. Entre más chico, mejor.

ARIMA = 54.1 contra Ingenuo = 64.8: el ARIMA gana, y aquí sin ambigüedad (a diferencia del MAPE).
La descomposición dice por qué: casi todo el puntaje de ambos es subestimación (42.2 y 53.4), y casi nada es sobreestimación. Los dos se quedaron cortos, el ingenuo más.
La cobertura es lo más revelador: los intervalos al 90% sólo contuvieron el valor real el 67% de las veces. Deberían ser 90%.

Redactado como en un artículo:

El modelo ARIMA con términos de Fourier superó al modelo estacional ingenuo (WIS 54.1 frente a 64.8). La descomposición del WIS indicó que el error se debió predominantemente a subestimación. La cobertura empírica de los intervalos de predicción al 90% fue de 67%, por debajo del nivel nominal, lo que sugiere que los intervalos son demasiado estrechos y que la incertidumbre está subestimada.

13.7 Resumen: qué métrica para qué

Qué quiero saber	Herramienta	Modelo
¿El modelo explica algo?	R², prueba F global	Lineal
¿Sobran variables?	R² ajustada, AIC, BIC	Todos
¿Se cumplen los supuestos?	Shapiro-Wilk, Breusch-Pagan, residuales	Lineal
¿Hay variables redundantes?	VIF	Todos
¿Qué tan lejos quedé?	MAE, RMSE (MSE si comparas)	Lineal / series
¿El conteo está sobredisperso?	Devianza/gl, dispersiontest	Poisson
¿El riesgo predicho es correcto?	Hosmer-Lemeshow, Brier	Logística
¿Ordena bien a los pacientes?	AUC (curva ROC)	Logística
¿A cuántos enfermos detecto?	Sensibilidad (no exactitud)	Logística
¿Mi pronóstico es honesto?	WIS y cobertura	Series de tiempo

Los cinco errores más comunes al reportar modelos

Reportar R² sin la prueba F. Una R² de 0.30 con p = 0.4 no es nada.
Reportar exactitud en desenlaces poco frecuentes. Siempre da bien y siempre engaña. Usa sensibilidad y especificidad.
Decidir binomial negativa por la varianza cruda. Se decide con los residuales del modelo.
Interpretar el OR como riesgo. Si el desenlace es común, exagera (ver el capítulo de encuestas).
Suponer linealidad sin probarla. Un spline cuesta dos líneas de código y aquí cambió las conclusiones.

13.8 Ejercicios

13.8.1 Interpretación

El modelo lineal tuvo R² = 0.025 y todos sus supuestos se cumplieron. Escribe dos o tres renglones explicándole a alguien que no sabe estadística por qué “cumplir los supuestos” y “servir” no son lo mismo.
Un colega te dice: “mi modelo tiene 92% de exactitud para predecir muerte materna”. La mortalidad materna ocurre en menos del 0.1% de los embarazos. ¿Qué le preguntas antes de creerle?
Con la tabla de cortes (0.5 contra 0.195), ¿cuál usarías si el modelo fuera a servir para decidir a qué pacientes se les asigna cama de terapia intensiva? Justifica en términos de qué cuesta cada tipo de error.

13.8.2 Cálculo

Ajusta el modelo logístico sin gcs y compara el AUC con el que tiene. ¿Cuánto de la capacidad predictiva venía de esa sola variable?
Corre anova(modelo_logistico, test = "LRT") y decide, con ese criterio, qué variables quitarías. Ajusta el modelo reducido y compara los AIC.
Prueba splines con 3, 4 y 5 grados de libertad para el ghq12. ¿Cuál elige el AIC? ¿Cambia la forma de la curva?
Calcula el MAE y el RMSE del modelo de Poisson (compara attack observado contra fitted(modelo_poisson)). ¿Por qué el RMSE es mayor?

13.8.3 A mano

Con la matriz de confusión del corte 0.5, calcula con calculadora la sensibilidad, la especificidad y el valor predictivo positivo. Verifica con la función metricas.
Si un modelo tiene sensibilidad 0.90 y especificidad 0.90, y lo aplicas a una población donde el desenlace ocurre en el 1%, ¿cuál es el valor predictivo positivo? (Pista: usa el teorema de Bayes o arma una tabla con 10,000 personas.) ¿Te sorprende el resultado?

13.8.4 Con apoyo de una IA

Pásale a una IA la tabla de métricas de clasificación del corte 0.5 y pregúntale si el modelo es bueno. Si contesta que sí basándose en la exactitud, corrígela y pídele que reconsidere. Anota qué le faltó ver.
Pídele que te redacte la sección de métodos para el modelo de Poisson con spline. Verifica que mencione: los grados de libertad del spline, cómo se eligieron, la prueba contra el modelo lineal y el criterio de bondad de ajuste. Casi siempre se le olvida alguno.

13.1 Los datos

13.2 Modelo lineal: presión sistólica al ingreso

13.2.1 R² y R² ajustada

13.2.2 Error de predicción: MAE, MSE y RMSE

13.2.3 Supuestos: normalidad, homocedasticidad y colinealidad

13.3 Modelo de Poisson: ataques de asma

13.3.1 Bondad de ajuste: ¿se cumple el supuesto de Poisson?

13.3.2 La prueba ji cuadrada de bondad de ajuste

13.3.3 ¿Y si usáramos binomial negativa?

13.3.4 Un spline: ¿el efecto del GHQ-12 es realmente lineal?

13.3.5 Efectos marginales

13.4 Modelo logístico: mortalidad hospitalaria

13.4.1 Bondad de ajuste global

13.4.2 Calibración: la prueba de Hosmer-Lemeshow

13.4.3 Discriminación: la curva ROC

13.4.4 Sensibilidad, especificidad y compañía

13.4.5 El punto de corte no es 0.5

13.4.6 El puntaje de Brier

13.4.7 Efectos marginales del modelo logístico

13.5 Modelo multinomial: categorías de presión arterial

13.5.1 Verificación: ¿clasifica bien?

13.5.2 ¿Sobran variables?

13.6 Series de tiempo: dengue

13.6.1 MAE, MSE, RMSE y MAPE

13.6.2 El WIS: la métrica que sí toma en cuenta la incertidumbre

13.7 Resumen: qué métrica para qué

13.8 Ejercicios

13.8.1 Interpretación

13.8.2 Cálculo

13.8.3 A mano

13.8.4 Con apoyo de una IA

13.2.1 `R²` y `R²` ajustada

13.2.2 Error de predicción: `MAE`, `MSE` y `RMSE`

13.6.1 `MAE`, `MSE`, `RMSE` y `MAPE`

13.6.2 El `WIS`: la métrica que sí toma en cuenta la incertidumbre