Modelos Estadisticos Aplicados I - Taller 4
Nelson Abad (nsabad@espol.edu.ec)
Dataset: FEV - Ejercicio 2
El dataset \(\texttt{fev.txt}\) incluye observaciones de 654 niños a los que se les mide el Volumen Espiratorio Forzado (FEV: forced expiratory volume). En este ejercicio estamos interesados en modelar la relación de la variable de respuesta FEV y las siguientes variables explicativas.
age: la edad del niño en años
height: la altura del niño en pulgadas
# Importamos los datos.
FEV_data = read.table("fev.txt", header = TRUE)
# Seleccionamos las variables a utilizar.
FEV_subset = FEV_data[,1:3]
# Adjuntamos las variables del dataframe.
attach(FEV_subset)Veamos un primer vistazo a los datos:
head(FEV_subset)summary(FEV_subset)## age FEV height
## Min. : 3.000 Min. :0.791 Min. :46.00
## 1st Qu.: 8.000 1st Qu.:1.981 1st Qu.:57.00
## Median :10.000 Median :2.547 Median :61.50
## Mean : 9.931 Mean :2.637 Mean :61.14
## 3rd Qu.:12.000 3rd Qu.:3.119 3rd Qu.:65.50
## Max. :19.000 Max. :5.793 Max. :74.00Revisemos los histogramas de cada variable para poder observar si tienen algun sesgo (especialmente la variable de respuesta).
par(mfrow = c(2, 2))
hist(FEV, freq = FALSE)
hist(age, freq = FALSE)
hist(height, freq = FALSE)
Se puede observar que la variable FEV tiene un sesgo a la derecha.
Pregunta 1. Realice una exploración gráfica de la relación entre FEV y la edad y la altura del niño. ¿Se ajustaría un modelo aditivo de primer orden o algún orden superior?
plot(FEV_subset)
# Graficamos las curvas para poder tener una mejor apreciación.
pairs(FEV_subset,
panel = function(x, y) {
points(x, y)
lines(lowess(x, y), col = "red")
})
La relación entre FEV tanto con las variables age y altura contienen tiene una cierta curvatura, se podría ajustar un modelo aditivo de orden superior.
Pregunta 2. Ajuste el siguiente modelo aditivo de primer orden, analice los residuos del modelo y detecte la curvatura y heterocedasticidad en los residuos:
\[ FEV_i = \beta_0 + \beta_1 age_i + \beta_2 height_i + \epsilon_i \]
# Modelo aditivo de primer orden.
modelo_orden1 = lm(FEV~height+age, data = FEV_subset)
summary(modelo_orden1)##
## Call:
## lm(formula = FEV ~ height + age, data = FEV_subset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.50533 -0.25657 -0.01184 0.24575 2.01914
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.610466 0.224271 -20.558 < 2e-16 ***
## height 0.109712 0.004716 23.263 < 2e-16 ***
## age 0.054281 0.009106 5.961 4.11e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4197 on 651 degrees of freedom
## Multiple R-squared: 0.7664, Adjusted R-squared: 0.7657
## F-statistic: 1068 on 2 and 651 DF, p-value: < 2.2e-16par(mfrow = c(2, 2))
plot(modelo_orden1)
Basándonos en el los gráficos del lado izquierdo, que nos ayudan a analizar los residuos, podemos observar que existe heterocedasticidad ademas de no linealidad.
Pregunta 3. Ajuste el siguiente modelo aditivo de segundo orden, analice los residuos del modelo:
\[ FEV_i = \beta_0 + \beta_1 age_i + \beta_2 height_i + \beta_3 (age_i - \bar{age})^2 + \beta_4 (height_i - \bar{height})^2 + \epsilon_i \]
# Modelo aditivo de segundo orden
modelo_orden2 = lm(FEV~height+age+
I((height-mean(height))^2)+
I((age-mean(age))^2), data = FEV_subset)par(mfrow = c(2, 2))
plot(modelo_orden2)
Analizando gráficamente podemos observar que el modelo se ajusta mucho mejor a los datos y existe una disminución en la heterocedasticidad.
Pregunta 4. Realice una prueba de hipótesis para los efectos cuadráticos de age y height mejoran el ajuste del modelo. Escriba el contraste de hipótesis, el estadístico de prueba, el valor p de la prueba y su conclusión.
\[ H_0: \beta_3 = \beta_4 = 0 \quad \quad \textrm{vs.} \quad \quad H_1: \beta_3 \neq 0 \quad \lor \quad \beta_4 \neq 0 \]
anova(modelo_orden1, modelo_orden2)Estadístico \(F_0 = 35.142\). Valor p: \(3.229 \cdot10^{-15}\). Rechazamos \(H_0\). Concluimos entonces que al menos uno de los efectos cuadráticos explica la variabilidad de FEV y mejoran el ajuste del modelo.
Pregunta 5. ¿Alguno de los efectos cuadráticos no es significativo, basado en los t-test marginales?. Elimine el efecto que no es significativo y ajuste el modelo reducido. Analice los residuos, comente qué supuesto parece no cumplirse.
# Resumen del segundo modelo.
summary(modelo_orden2)##
## Call:
## lm(formula = FEV ~ height + age + I((height - mean(height))^2) +
## I((age - mean(age))^2), data = FEV_subset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.65216 -0.22852 0.00417 0.24157 1.87632
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.9552802 0.2359057 -21.005 < 2e-16 ***
## height 0.1134475 0.0051600 21.986 < 2e-16 ***
## age 0.0547062 0.0109636 4.990 7.77e-07 ***
## I((height - mean(height))^2) 0.0031505 0.0004884 6.451 2.17e-10 ***
## I((age - mean(age))^2) 0.0011300 0.0017611 0.642 0.521
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3993 on 649 degrees of freedom
## Multiple R-squared: 0.7892, Adjusted R-squared: 0.7879
## F-statistic: 607.6 on 4 and 649 DF, p-value: < 2.2e-16Podemos observar que el efecto cuadrático de la variable age no explica significativamente la variable FEV. El valor p es de \(0.521\), por lo que no podemos rechazar \(H_0: \beta_3 = 0\).
# Modelo reducido.
modelo_orden2_red = lm(FEV~height+age+
I((height-mean(height))^2), data = FEV_subset)
summary(modelo_orden2_red)##
## Call:
## lm(formula = FEV ~ height + age + I((height - mean(height))^2),
## data = FEV_subset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.65645 -0.23224 0.00792 0.23676 1.88121
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.8944923 0.2159506 -22.665 < 2e-16 ***
## height 0.1118203 0.0044918 24.894 < 2e-16 ***
## age 0.0590025 0.0086776 6.799 2.38e-11 ***
## I((height - mean(height))^2) 0.0033316 0.0003984 8.363 3.73e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3991 on 650 degrees of freedom
## Multiple R-squared: 0.7891, Adjusted R-squared: 0.7881
## F-statistic: 810.7 on 3 and 650 DF, p-value: < 2.2e-16par(mfrow = c(2, 2))
plot(modelo_orden2_red)
Aun existe cierta heterocedasticidad, a pesar de que hubo una considerable disminución respecto a los modelos anteriores, entonces el supuesto que no se cumple seria el de varianza constante. Tambien basándonos en el gráfico Q-Q plot podemos ver que a los costados existe cierto desvió, por lo tanto se podría sospechar que los residuos no siguen una distribución normal.
Pregunta 6. Compare el modelo 2 y 3 en base al coeficiente de determinación ajustado. ¿Qué modelo prefiere?
data.frame(
"Modelo_orden2" = summary(modelo_orden2)$adj.r.squared,
"Modelo_orden2_reducido" = summary(modelo_orden2_red)$adj.r.squared
)Vemos que el modelo de orden 2 reducido tiene un mayor \(R^2_{a}\), por lo tanto explica mayor variabilidad de la variable de respuesta FEV.
Pregunta 7. Para corregir el problema de heterocedasticidad, estimemos el modelo 3 usando los mínimos cuadrados ponderados. Analice los residuos \(r^{(w)}\). ¿Existe violación grave de algún supuesto?
modelo_orden2_red_ws = lm(abs(residuals(modelo_orden2_red)) ~ fitted.values(modelo_orden2_red))
ws = 1/fitted.values(modelo_orden2_red_ws)^2
modelo_corregido = lm(FEV~height+age+
I((height-mean(height))^2),
data = FEV_subset,
weights = ws)
summary(modelo_corregido)##
## Call:
## lm(formula = FEV ~ height + age + I((height - mean(height))^2),
## data = FEV_subset, weights = ws)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -4.6818 -0.8155 0.0202 0.7848 5.1240
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.922580 0.231615 -21.253 < 2e-16 ***
## height 0.114427 0.004724 24.221 < 2e-16 ***
## age 0.047012 0.008627 5.449 7.18e-08 ***
## I((height - mean(height))^2) 0.002887 0.000336 8.590 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.289 on 650 degrees of freedom
## Multiple R-squared: 0.8132, Adjusted R-squared: 0.8123
## F-statistic: 943.2 on 3 and 650 DF, p-value: < 2.2e-16par(mfrow = c(2, 2))
plot(modelo_corregido)
Analizamos lo residuos \(r^{(w)}\):
rw = sqrt(diag(ws)) %*% (FEV - model.matrix(modelo_corregido) %*% coef(modelo_corregido))par(mfrow = c(1, 2))
plot(fitted.values(modelo_corregido), rw)
lines(lowess(rw ~ fitted.values(modelo_corregido)), col = "red")
qqnorm(rw)
qqline(rw)
Comprobemos normalidad de los residuos:
shapiro.test(rw)##
## Shapiro-Wilk normality test
##
## data: rw
## W = 0.99716, p-value = 0.3148Valor p: \(0.3148\). No podemos rechazar entonces \(H_0: \textrm{Los residuos} \ r^{(w)} \ \textrm{siguen una distribucion normal}\).
Realicemos un test estadístico para comprobar varianza constante.
library(lmtest)
bptest(modelo_corregido)##
## studentized Breusch-Pagan test
##
## data: modelo_corregido
## BP = 10944, df = 3, p-value < 2.2e-16Valor p: \(2.2 \cdot 10^{-16}\). Rechazamos \(H_0: \textrm{Los residuos} \ r^{(w)} \ \textrm{tienen la misma varianza}\).
Por lo tanto no se cumple el supuesto de homocedasticidad.
Pregunta 8. Utilice las funciones \(\texttt{plot3d()}\) y \(\texttt{persp3d()}\) de la librería \(\texttt{rgl}\) para graficar la superficie de respuesta del último modelo.
library(rgl)
regsurface <- function(x1, x2){
modelo_corregido$coef[1] + modelo_corregido$coef[2] * x1 +
modelo_corregido$coef[3] * x2 + modelo_corregido$coef[4] * x2^2
}
x1 <- seq(0, 20)
x2 <- seq(50, 80)
y <- matrix(0, length(x1), length(x2))
for (i in 1:length(x1)){
for (j in 1:length(x2)){
y[i, j] <- regsurface(x1[i], x2[j])
}
}
plot3d(age, height, FEV, type = "s", col = "red", size = 1)
persp3d(x1, x2, y, alpha = 0.3, add = TRUE)
lines3d(x1, 50, regsurface(x1, 50), add = TRUE, col = "red", lwd = 2)
lines3d(x1, 60, regsurface(x1, 60), add = TRUE, col = "red", lwd = 2)