1 Ejemplo inicial: hipertensión pediátrica
Contexto. En esta sesión se desarrolla el ejemplo de
hipertensión pediátrica abordado en la lectura. Se plantea un modelo de
regresión poblacional que relaciona la presión sistólica
(SPB) \(Y\), el peso al
nacimiento \(X_1\) y la edad en días
\(X_2\).
En un segundo paso se revisará una observación con residual estandarizado extremo y alta influencia. Siguiendo el ejemplo didáctico, se ajustará un segundo modelo sin dicha observación y se compararán los resultados.
\[ SPB_i=\beta_0+\beta_1 Peso_i+\beta_2 Edad_i+\varepsilon_i, \]
donde \(\beta_1\) representa el cambio promedio esperado en la presión sistólica asociado con una unidad adicional de peso, manteniendo constante la edad, y \(\beta_2\) representa el cambio promedio asociado con una unidad adicional de edad, manteniendo constante el peso.
Importante. Las interpretaciones de los coeficientes son asociaciones ajustadas por las demás variables del modelo. El modelo de regresión, por sí solo, no establece causalidad.
2 Primer modelo de regresión
2.1 Captura y organización de los datos
Primeramente se captura la información completa y se forma un marco de datos:
datos_peso <- matrix(
c(
135,3,89, 120,4,90, 100,3,83, 105,2,77,
130,4,92, 125,5,98, 125,2,82, 105,3,85,
120,5,96, 90,4,95, 120,2,80, 95,3,79,
120,3,86, 150,4,97, 160,3,92, 125,3,88
),
ncol = 3,
byrow = TRUE
)
datos_peso2 <- data.frame(
Peso = datos_peso[, 1],
Edad = datos_peso[, 2],
SPB = datos_peso[, 3]
)
datos_peso22.2 Ajuste del modelo
#>
#> Call:
#> lm(formula = SPB ~ Peso + Edad, data = datos_peso2)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -4.0438 -1.3481 -0.2395 0.9688 6.6964
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 53.45019 4.53189 11.794 2.57e-08 ***
#> Peso 0.12558 0.03434 3.657 0.0029 **
#> Edad 5.88772 0.68021 8.656 9.34e-07 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 2.479 on 13 degrees of freedom
#> Multiple R-squared: 0.8809, Adjusted R-squared: 0.8626
#> F-statistic: 48.08 on 2 and 13 DF, p-value: 9.844e-07
El modelo explica aproximadamente 88.09% de la
variabilidad observada en SPB. El \(R^2\) ajustado es 86.26%.
La prueba global del modelo presenta un p-valor de 9.844e-07. Para los efectos parciales, el p-valor de Peso es 0.002896 y el de Edad es 9.342e-07.
Interpretación. Al nivel de significancia \(\alpha=0.05\), ambas variables presentan
evidencia de asociación con SPB después de ajustar por la
otra variable. Es importante distinguir \(R^2\) de \(R^2\) ajustado: el primero es la proporción
de variabilidad explicada por el modelo y el segundo penaliza por el
número de predictores.
3 Diagnóstico del primer modelo
El análisis de supuestos no debe basarse en una sola prueba. Se combinarán herramientas gráficas, residuales estandarizados y medidas de influencia.
3.1 Normalidad de los residuales
res <- rstandard(reg0)
qqnorm(
res,
main = "Gráfico Q-Q de residuales estandarizados"
)
qqline(res, col = "#6A0F35", lwd = 2)
#>
#> Shapiro-Wilk normality test
#>
#> data: res
#> W = 0.84847, p-value = 0.01294La prueba de Shapiro-Wilk produce un p-valor de 0.0129. Al ser menor que 0.05, existe evidencia contra el supuesto de normalidad de los residuales en este primer ajuste.
La prueba de Shapiro-Wilk es sensible al tamaño muestral. Su resultado debe leerse junto con el gráfico Q-Q y el análisis de observaciones potencialmente atípicas o influyentes.
3.2 Residuales, apalancamiento e influencia
Un residual estandarizado con valor absoluto mayor que 3 suele considerarse una señal de observación potencialmente atípica. Sin embargo, para decidir si una observación es problemática también conviene revisar su apalancamiento y la distancia de Cook.
diag0 <- data.frame(
Observacion = seq_len(nrow(datos_peso2)),
Residual_est = rstandard(reg0),
Residual_student = rstudent(reg0),
Leverage = hatvalues(reg0),
Cook = cooks.distance(reg0)
)
diag0$Clasificacion <- ifelse(
abs(diag0$Residual_est) > 3,
"Posible outlier",
"No outlier"
)
diag0La observación con el residual estandarizado de mayor magnitud es la 10, con residual estandarizado de aproximadamente 3.208 y distancia de Cook de 1.41.
Como referencia exploratoria, \(4/n=0.25\) puede utilizarse como umbral para señalar observaciones con influencia que ameritan inspección.
No debe eliminarse automáticamente una observación solo por ser extrema. Primero debe verificarse si existe un error de captura, una condición clínica especial o una explicación sustantiva. En esta sesión, siguiendo el ejemplo didáctico original y después de identificar su fuerte influencia, se continuará comparando el ajuste al excluir la observación 10.
3.3 Revisión gráfica de varianza constante
Para evaluar homocedasticidad, el gráfico de residuales frente a valores ajustados es más directo que analizar cada predictor por separado. Aun así, se conservan ambos gráficos por su utilidad didáctica.
plot_ly(
data = datos_peso3,
x = ~Peso,
y = ~Residuales,
color = ~Outlier,
colors = c("No outlier" = "#88898D", "Posible outlier" = "#6A0F35"),
type = "scatter",
mode = "markers"
) |>
layout(
xaxis = list(title = "Peso"),
yaxis = list(title = "Residual estandarizado"),
title = "Peso vs residuales estandarizados"
)plot_ly(
data = datos_peso3,
x = ~Edad,
y = ~Residuales,
color = ~Outlier,
colors = c("No outlier" = "#88898D", "Posible outlier" = "#6A0F35"),
type = "scatter",
mode = "markers"
) |>
layout(
xaxis = list(title = "Edad"),
yaxis = list(title = "Residual estandarizado"),
title = "Edad vs residuales estandarizados"
)plot_ly(
data = datos_peso3,
x = ~Predicho,
y = ~Residuales,
color = ~Outlier,
colors = c("No outlier" = "#88898D", "Posible outlier" = "#6A0F35"),
type = "scatter",
mode = "markers"
) |>
layout(
xaxis = list(title = "Valores ajustados"),
yaxis = list(title = "Residual estandarizado"),
title = "Valores ajustados vs residuales"
)Lectura del gráfico. Se busca una nube aproximadamente horizontal, centrada alrededor de cero y sin forma de embudo. Un patrón sistemático sugeriría problemas de especificación o heterocedasticidad.
4 Segundo modelo de regresión
4.1 Ajuste sin la observación 10
Siguiendo el planteamiento del ejemplo, se excluye la observación 10 y se ajusta nuevamente el modelo:
obs_excluir <- 10L
datos_peso_seg <- datos_peso2[-obs_excluir, ]
reg00 <- lm(SPB ~ Peso + Edad, data = datos_peso_seg)
summary(reg00)#>
#> Call:
#> lm(formula = SPB ~ Peso + Edad, data = datos_peso_seg)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.1850 -0.8061 0.2360 0.6790 1.9834
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 47.9377 2.3015 20.829 8.68e-11 ***
#> Peso 0.1832 0.0184 9.955 3.76e-07 ***
#> Edad 5.2825 0.3352 15.759 2.21e-09 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 1.177 on 12 degrees of freedom
#> Multiple R-squared: 0.9732, Adjusted R-squared: 0.9687
#> F-statistic: 217.5 on 2 and 12 DF, p-value: 3.741e-10El \(R^2\) aumenta de 88.09% en el modelo original a 97.32% en el modelo sin la observación 10. El \(R^2\) ajustado pasa de 86.26% a 96.87%.
El segundo modelo continúa siendo globalmente significativo, con p-valor 3.741e-10.
Un incremento de \(R^2\) después de eliminar una observación no constituye, por sí mismo, una justificación para eliminarla. La comparación es útil aquí porque la observación 10 ya había mostrado un residual extremo y alta influencia.
5 Selección de variables
Se comparan las estrategias backward, forward y
stepwise mediante el criterio AIC implementado por
step().
modelo_completo2 <- lm(SPB ~ Peso + Edad, data = datos_peso_seg)
modelo_reducido2 <- lm(SPB ~ 1, data = datos_peso_seg)
mod_backward <- step(
modelo_completo2,
scope = list(lower = modelo_reducido2, upper = modelo_completo2),
direction = "backward",
trace = 0
)
mod_forward <- step(
modelo_reducido2,
scope = list(lower = modelo_reducido2, upper = modelo_completo2),
direction = "forward",
trace = 0
)
mod_stepwise <- step(
modelo_reducido2,
scope = list(lower = modelo_reducido2, upper = modelo_completo2),
direction = "both",
trace = 0
)
seleccion_final <- data.frame(
Metodo = c("Backward", "Forward", "Stepwise"),
Modelo = c(
paste(deparse(formula(mod_backward)), collapse = ""),
paste(deparse(formula(mod_forward)), collapse = ""),
paste(deparse(formula(mod_stepwise)), collapse = "")
),
AIC = c(AIC(mod_backward), AIC(mod_forward), AIC(mod_stepwise))
)
seleccion_finalResultado de la selección. En este conjunto de datos, las tres estrategias conservan el modelo con Peso y Edad. Debido a que solo existen dos predictores y la muestra es pequeña, la selección automática aporta poco más que la comparación directa de modelos. En aplicaciones reales, la selección debe apoyarse también en conocimiento sustantivo y no únicamente en AIC.
6 Diagnóstico del segundo modelo
6.1 Normalidad de los residuales
res2 <- rstandard(reg00)
qqnorm(
res2,
main = "Gráfico Q-Q: segundo modelo"
)
qqline(res2, col = "#6A0F35", lwd = 2)
#>
#> Shapiro-Wilk normality test
#>
#> data: res2
#> W = 0.98095, p-value = 0.9756La prueba de Shapiro-Wilk produce un p-valor de 0.9756. En este segundo modelo no existe evidencia suficiente para rechazar la normalidad de los residuales al nivel \(\alpha=0.05\).
Es preferible decir que no se rechaza la normalidad en lugar de afirmar que se ha demostrado que los residuales son normales.
6.2 Residuales e influencia
obs_originales <- setdiff(seq_len(nrow(datos_peso2)), obs_excluir)
diag2 <- data.frame(
Observacion_original = obs_originales,
Residual_est = rstandard(reg00),
Residual_student = rstudent(reg00),
Leverage = hatvalues(reg00),
Cook = cooks.distance(reg00)
)
diag2$Clasificacion <- ifelse(
abs(diag2$Residual_est) > 3,
"Posible outlier",
"No outlier"
)
diag2El mayor residual estandarizado en valor absoluto es 2.13; por tanto, no existen observaciones con \(|r_i|>3\) en este segundo ajuste.
Usando \(4/n=0.267\) como referencia exploratoria para la distancia de Cook, las observaciones originales que ameritan una inspección adicional son: 12, 15.
Una observación influyente no es necesariamente errónea y no debe eliminarse de forma automática. La distancia de Cook señala cuánto puede cambiar el ajuste al retirar una observación; debe interpretarse junto con el contexto y otras medidas diagnósticas.
6.3 Varianza constante
plot_ly(
data = datos_peso4,
x = ~Peso,
y = ~Residuales,
color = ~Outlier,
colors = c("No outlier" = "#88898D", "Posible outlier" = "#6A0F35"),
type = "scatter",
mode = "markers"
) |>
layout(
xaxis = list(title = "Peso"),
yaxis = list(title = "Residual estandarizado"),
title = "Peso vs residuales: segundo modelo"
)plot_ly(
data = datos_peso4,
x = ~Edad,
y = ~Residuales,
color = ~Outlier,
colors = c("No outlier" = "#88898D", "Posible outlier" = "#6A0F35"),
type = "scatter",
mode = "markers"
) |>
layout(
xaxis = list(title = "Edad"),
yaxis = list(title = "Residual estandarizado"),
title = "Edad vs residuales: segundo modelo"
)plot_ly(
data = datos_peso4,
x = ~Predicho,
y = ~Residuales,
type = "scatter",
mode = "markers",
marker = list(color = "#6A0F35", size = 9)
) |>
layout(
xaxis = list(title = "Valores ajustados"),
yaxis = list(title = "Residual estandarizado"),
title = "Valores ajustados vs residuales: segundo modelo"
)Homocedasticidad. Si la dispersión de los residuales alrededor de cero es aproximadamente similar a lo largo del rango de valores ajustados y no aparece una forma de embudo clara, los gráficos son compatibles con el supuesto de varianza aproximadamente constante. Esto debe expresarse como evidencia gráfica, no como una demostración definitiva del supuesto.
7 Conclusiones
- El primer modelo con Peso y Edad presenta una alta proporción de variabilidad explicada y ambos predictores son estadísticamente significativos.
- La observación 10 presenta un residual estandarizado extremo y una distancia de Cook elevada, por lo que amerita revisión.
- Siguiendo el ejemplo didáctico y excluyendo dicha observación, el ajuste global aumenta de \(R^2=88.09\%\) a \(R^2=97.32\%\).
- Después de la exclusión, no existe evidencia suficiente para rechazar la normalidad de los residuales y no se observan residuales estandarizados con valor absoluto mayor que 3.
- La inspección de homocedasticidad debe basarse principalmente en residuales frente a valores ajustados y en la ausencia de patrones sistemáticos claros.
Sobre independencia. La independencia de los errores depende principalmente del diseño de recolección de los datos. No puede demostrarse únicamente a partir de los gráficos de residuales. Si las observaciones corresponden a individuos distintos y no existe agrupamiento, repetición de medidas o dependencia temporal/espacial, el supuesto puede considerarse razonable por diseño; en caso contrario debería modelarse explícitamente la estructura de dependencia.