Sesión 4. EMED - MEB
Dr. Héctor de la Torre Gutiérrez
Fecha de modificación: 15-09-2023
1 Ejemplo inicial -Hipertensión pediátrica
En esta sesión se presentará y desarrollará el ejemplo respectivo a la hipertensión pediátrica (abordado en la lectura). Se plantea el modelo de regresión poblacional que relaciona la presión sistólica (por sus siglas en inglés, SPB) \(y\), el peso al nacimiento (\(x_{1}\)) y, la edad en días (\(x_{2}\)).
En un segundo paso, y dado que fue detectado un \(outlier\) en los datos, se ajustará un nuevo modelo de regresión sin considerar dicha observación, y se compararán los resultados.
2 Primer modelo de regresión
Primeramente, se captura la información completa, y se forma un marco de datos mediante las siguientes instrucciones:
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])A continuación, se ajusta el primer modelo de regresión:
reg0<-lm(SPB~Peso+Edad,data = datos_peso2)
summary(reg0)##
## 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-07Se observa del modelo anterior que ambas variables (Peso y Edad) resultaron ser significativas un nivel \(\alpha = 0.05\), el 86.26% de la varianza fue explicada por estas dos variables, y la regresión sí fue significativa ya que el p-value es muy pequeño (0.0000009844). Ahora, se analizará la varianza explicada por el modelo y por los residuales:
summary(aov(reg0))## Df Sum Sq Mean Sq F value Pr(>F)
## Peso 1 130.5 130.5 21.24 0.00049 ***
## Edad 1 460.5 460.5 74.92 9.34e-07 ***
## Residuals 13 79.9 6.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1De lo anterior se observa que, tal y como lo determinó la \(R^{2}\), el 88% de la varianza es explicada por el modelo.
El siguiente paso a seguir es probar los supuestos del modelo y el análisis de outliers.
3 Revisión de supuestos del primer modelo
3.1 Supuesto de Normalidad de los residuales
Se utilizará el gráfico Normal QQ y la prueba de Shapiro-Wilks para probar este supuesto. Dichos elementos se obtienen al ejecutar las siguientes instrucciones:
res<-rstandard(reg0) #Se obtienen los residuales estandarizados
qqnorm(res)
qqline(res)
shapiro.test(res) #Prueba de Shapiro-Wilks usando residuales estandarizados##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.84847, p-value = 0.01294Se observa de los elementos anteriores que NO se cumple el supuesto de Normalidad. Parece que existe un punto atípico.
3.2 Supuesto de Varianza Constante
Para probar este supuesto, se graficarán cada una de las variables independientes contra los valores de los residuales. Para lo anterior, nos auxiliaremos de los gráficos interactivos que se pueden elaborar en R:
library(plotly)
#Se genera una función la cual nos generará etiquetas sobre si un dato es outlier o no, según si es un valor de residual estandarizado mayor a 3
fun_outlier<-function(X){if(abs(X)>3){"Outlier"}else{"No outlier"}}
#Se aplica a los residuales obtenidos
outlierss<-sapply(res,FUN=fun_outlier)
#Se genera un nuevo marco de datos, agregándole una columna con los valores de los residuales estandarizados y las etiquetas de cada uno de estos
#Varianza constante
datos_peso3<-data.frame(datos_peso2,Residuales=res,Outlier=outlierss)
#El primer gráfico de Peso Vs residuales
plot_ly(data = datos_peso3, x = ~Peso, y = ~Residuales, color = ~Outlier)#El segundo gráfico de Peso Vs residuales
plot_ly(data = datos_peso3, x = ~Edad, y = ~Residuales, color = ~Outlier)## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levelsDe los dos gráficos anteriores se puede observar claramente que existe una observación tipo outlier, la cual se eliminará. Lo anterior se confirmar al correr la siguiente instrucción:
summary(res)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.80824 -0.58695 -0.10466 0.01442 0.42677 3.20842Eliminado dicha observación, se ajustarán nuevos modelos.
4 Segundo modelo de regresión
Ya eliminada la observación no. 10, se ajustará el siguiente modelo:
datos_peso_seg<-datos_peso2[-10,]
reg00<-lm(SPB~Peso+Edad,data=datos_peso_seg) Su análisis de varianza es:
summary(aov(reg00))## Df Sum Sq Mean Sq F value Pr(>F)
## Peso 1 258.8 258.8 186.7 1.12e-08 ***
## Edad 1 344.2 344.2 248.3 2.21e-09 ***
## Residuals 12 16.6 1.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El resumen de ajuste del nuevo modelo:
#¿Cuáles son las variables significativas?
summary(reg00) #Resumen del modelo ajustado##
## 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-10Se observa que, en comparación con el modelo que considera el outlier, las dos variables (Peso y Edad) siguen siendo significativas, la regresión es significativa, y el ajuste global del modelo aumentó sustancialmente, de 88% a 97%!
Se analizan los posibles tres modelos de selección de variables (modelos), backward, forward, step.
#############se define el modelo reducido###############
modelo_completo2<-lm(SPB~Peso+Edad,data = datos_peso_seg)
summary(modelo_completo2)##
## 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-10modelo_reducido2<-lm(SPB~1,data = datos_peso_seg)
summary(modelo_reducido2)##
## Call:
## lm(formula = SPB ~ 1, data = datos_peso_seg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.6 -5.1 0.4 4.4 10.4
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 87.600 1.718 51 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.653 on 14 degrees of freedom##################regresion backward###################
step(modelo_completo2,scope=list(lower=modelo_reducido2,upper=modelo_completo2),direction = "backward") # Se fija el argumento a backward## Start: AIC=7.55
## SPB ~ Peso + Edad
##
## Df Sum of Sq RSS AIC
## <none> 16.63 7.549
## - Peso 1 137.37 154.00 38.934
## - Edad 1 344.21 360.85 51.706##
## Call:
## lm(formula = SPB ~ Peso + Edad, data = datos_peso_seg)
##
## Coefficients:
## (Intercept) Peso Edad
## 47.9377 0.1832 5.2825##################Regresión forward###################
step(modelo_reducido2,scope=list(lower=modelo_reducido2,upper=modelo_completo2),direction = "forward") #Se cambia el argumento a forward## Start: AIC=57.82
## SPB ~ 1
##
## Df Sum of Sq RSS AIC
## + Edad 1 465.60 154.00 38.934
## + Peso 1 258.75 360.85 51.706
## <none> 619.60 57.815
##
## Step: AIC=38.93
## SPB ~ Edad
##
## Df Sum of Sq RSS AIC
## + Peso 1 137.37 16.632 7.549
## <none> 154.000 38.934
##
## Step: AIC=7.55
## SPB ~ Edad + Peso##
## Call:
## lm(formula = SPB ~ Edad + Peso, data = datos_peso_seg)
##
## Coefficients:
## (Intercept) Edad Peso
## 47.9377 5.2825 0.1832##################Regresión stepwise###################
step(modelo_reducido2,scope=list(lower=modelo_reducido2,upper=modelo_completo2),direction = "both") #Se cambia el argumento a “both”## Start: AIC=57.82
## SPB ~ 1
##
## Df Sum of Sq RSS AIC
## + Edad 1 465.60 154.00 38.934
## + Peso 1 258.75 360.85 51.706
## <none> 619.60 57.815
##
## Step: AIC=38.93
## SPB ~ Edad
##
## Df Sum of Sq RSS AIC
## + Peso 1 137.37 16.63 7.549
## <none> 154.00 38.934
## - Edad 1 465.60 619.60 57.815
##
## Step: AIC=7.55
## SPB ~ Edad + Peso
##
## Df Sum of Sq RSS AIC
## <none> 16.63 7.549
## - Peso 1 137.37 154.00 38.934
## - Edad 1 344.21 360.85 51.706##
## Call:
## lm(formula = SPB ~ Edad + Peso, data = datos_peso_seg)
##
## Coefficients:
## (Intercept) Edad Peso
## 47.9377 5.2825 0.1832De lo anterior podemos observar que el modelo que considera las dos variables es el que mejor ajuste ofrece.
Ahora, se procede a la validación de los supuestos del modelo y la detección de posibles outliers.
4.1 Supuesto de Normalidad
Se probará este supuesto con el gráfico QQ Normal y la prueba de Shapiro-Wilks:
res2<-rstandard(reg00) #Se obtienen los residuales estandarizados
qqnorm(res2)
qqline(res2)
shapiro.test(res2) #Prueba de Shapiro-Wilks usando residuales estandarizados##
## Shapiro-Wilk normality test
##
## data: res2
## W = 0.98095, p-value = 0.9756Se observa ahora que SÍ se cumple el supuesto de Normalidad de los residuales.
4.2 Supuesto de varianza constante
Se vuelve a generar gráficos interactivos para observar si varianza constante.
outlierss2<-sapply(res2,FUN=fun_outlier) #Se aplica la misma función fun_outlier, pero ahora a res2
datos_peso_seg## Peso Edad SPB
## 1 135 3 89
## 2 120 4 90
## 3 100 3 83
## 4 105 2 77
## 5 130 4 92
## 6 125 5 98
## 7 125 2 82
## 8 105 3 85
## 9 120 5 96
## 11 120 2 80
## 12 95 3 79
## 13 120 3 86
## 14 150 4 97
## 15 160 3 92
## 16 125 3 88#Varianza constante
datos_peso4<-data.frame(datos_peso_seg,Residuales=res2,Outlier=outlierss2)
library(plotly)
plot_ly(data = datos_peso4, x = ~Peso, y = ~Residuales, color = ~Outlier)## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levelsSe observa en el gráfico anterior que, en cuanto a la \(Peso\), sí existe una varianza constante. Ahora, respecto a la variable \(Edad\), se observa lo siguiente:
plot_ly(data = datos_peso4, x = ~Edad, y = ~Residuales, color = ~Outlier)## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levelsTambién se observa varianza constante en el gráfico respecto a \(Edad\).
En cuanto a la identificación de outliers, se buscará si existe alguno mayor a 3.0 o menor a -3.0
summary(res2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.13028 -0.75101 0.20818 -0.02241 0.67133 1.81346No se observó ningún residual estandarizado mayor a 3.0 o menor a -3.0, por lo que concluimos que no existen outliers.
En resumen, en los modelos iniciales se detectó un outlier y no se cumplió el supuesto de Normalidad. Removiendo dicho outlier, se mejoró el ajuste global del modelo y se cumplieron los supuestos del modelo.
El supuesto de independencia se omitió ya que, dado el contexto del problema estudiado, se entiende que las muestras obtenidas no guadan relación entre ellas.