COMPONENTE APLICADA 4
BF 2021-I - UNAL
Procedimiento 1
Usando la tabla natalidad.csv realice el diagrama de dispersión y realice alguna prueba de independencia para las variables “talla o estatura del recién nacido” (talla_nac) y “peso del recién nacido” (peso_nac). Explique descriptivamente e inferencialmente si se debería rechazar que las variables son independientes o no. Si alguna de las variables depende de la otra ¿el peso depende de la talla o viceversa? Justifique.
1.1 Preparación de los datos
Cargamos los datos desde el archivo natalidad.csv con la función read.csv
natalidad <- read.csv("C:/Users/johan/Documents/Johann/UN/2021-I/BIOESTADISTICA/Componente_aplicada4/natalidad.csv",
header = TRUE, sep = ";")Creamos un data.frame con las variables de interés. En este documento se tomará (en caso de existir dependencia) “Peso nacimiento” en función de “Talla nacimiento”.
datos_nat <- data.frame(talla_nac, peso_nac)## [1] "table table-condensed table-bordered"1.2 Diagrama de dispersión
1.2.1 Uso de SQL para filtrar datos
Se considera atípico que un recién nacido pese más de 6000 g y talle más de 75 cm, es por ello que se eliminan dichos datos mediante SQL en R.
nat <- sqldf("SELECT * FROM datos_nat WHERE talla_nac < 75 AND peso_nac < 6000")
1.3 Prueba de independencia
1.3.1 Determinar normalidad en las variables
Para determinar por cuál método deberíamos hacer la prueba para hallar el coeficiente de correlación, primero hacemos uso de la función ad.test() y verificamos si los datos siguen una distribución Normal. Recordemos que detrás de este test de Normalidad hay una prueba de hipótesis de la siguente forma:
\[ \begin{equation*} \begin{cases} H_0: X \sim N(\mu, \sigma^{2}) \\ H_1: X \nsim N(\mu, \sigma^{2}) \end{cases} \end{equation*} \]
ad.test(nat$talla_nac)##
## Anderson-Darling normality test
##
## data: nat$talla_nac
## A = 4548.8, p-value < 2.2e-16ad.test(nat$peso_nac)##
## Anderson-Darling normality test
##
## data: nat$peso_nac
## A = 1295.1, p-value < 2.2e-16De esta manera, ninguna de las variables de interés sigue una distribución normal, pues en ambos casos se tiene un p-value menor al nivel de significancia por defecto (0.05)
1.3.2 Correlación
Mediante la función cor.test() hallamos el coeficiente de correlación por el método de Spearman, ya que los datos no siguen una distribución normal.
cor.test(nat$talla_nac, nat$peso_nac, method = "spearman")##
## Spearman's rank correlation rho
##
## data: nat$talla_nac and nat$peso_nac
## S = 3.6781e+14, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.7099799Notemos que se tiene un rho = 0.7099799, es decir el peso de un recién nacido depende sustancialmente de su talla (estatura). Sin embargo se tiene un p-valor inferior a 2.2 x 10^-16, lo cual sugiere que el coeficiente de correlación es no significativo.
También podríamos estar interesados en usar el coeficiente tau de Kendall. Para ello tenemos el siguiente código
# cor.test(nat$talla_nac, nat$peso_nac, method = "kendall")
# Kendall's rank correlation tau
# data: nat$talla_nac and nat$peso_nac
# z = 348.87, p-value < 2.2e-16
# alternative hypothesis: true tau is not equal to 0
# sample estimates:
# tau
# 0.5576303 Nótese que en este caso se tiene una dependencia moderada. De todos modos, el coeficiente tau de Kendall es también no significativo.
Procedimiento 2
Verifique o pruebe si se satisfacen o no cada una de las hipótesis del modelo de regresión lineal simple: (1) Independencia de los residuales, (2) Homocedasticidad y (3) Normalidad de los residuales ¿es adecuado el modelo de regresión lineal simple para los datos?
2.1 Regresión lineal
regresion1 <- lm(nat$peso_nac ~ nat$talla_nac, data = nat)
summary(regresion1)##
## Call:
## lm(formula = nat$peso_nac ~ nat$talla_nac, data = nat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2627.1 -200.8 -9.3 191.2 3854.4
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3891.0462 12.1556 -320.1 <2e-16 ***
## nat$talla_nac 141.4657 0.2471 572.4 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 304.9 on 196688 degrees of freedom
## Multiple R-squared: 0.6249, Adjusted R-squared: 0.6249
## F-statistic: 3.277e+05 on 1 and 196688 DF, p-value: < 2.2e-16Es de notar el resultado de Coefficients, dado que este determina la recta de regresión. Para nuestro caso la recta está dada por:
\[ \begin{align*} y&= 141.4657x + (-3891.0462) + \epsilon \\ &= 141.4657x -3891.0462 + \epsilon \end{align*} \]
Por otro lado también hay que tener en cuenta el resultado de Multiple R-squared, este es igual a 0.6249 y nos resultará útil más adelante.
En adición el gráfico de regresión es:
2.2 Independencia residuales
2.2.1 Descriptivamente
Es de resaltar, que los datos parecen no seguir un patrón y descriptivamente, esto implica que se no rechaze la independencia de los residuales. Sin embargo debemos verificar lo que ocurre si observamos la independencia de los residuales de una forma más minuciosa a través de una prueba de hipótesis (inferencialmente).
2.2.2 Inferencialmente
Aplicamos el test de Durbin-Watson. En adición recordemos que lo que hay detrás del siguiente test, es una prueba de hipótesis de la siguiente forma:
\[ \begin{equation*} \begin{cases} H_0: \epsilon_{1}, \epsilon_{2}, \cdots, \epsilon_{n} \text{ son independientes} \\ H_1: \epsilon_{1}, \epsilon_{2}, \cdots, \epsilon_{n} \text{ NO son independientes} \end{cases} \end{equation*} \]
donde
\[ \epsilon_{i} \]
corresponde al iésimo residual.
dwtest(regresion1)##
## Durbin-Watson test
##
## data: regresion1
## DW = 1.9737, p-value = 2.739e-09
## alternative hypothesis: true autocorrelation is greater than 0Como el p-valor es menor a 0.05 rechazamos la hipótesis nula y con ello que haya independencia en los residuales.
2.3 Homocedasticidad
Esto implica que
\[ \epsilon_{1}, \epsilon_{2}, \cdots, \epsilon_{n} \text{ tienen una misma varianza } \sigma^{2} \]
2.3.1 Descriptivamente
Al observar el anterior gráfico, notamos que los datos no se podrían encerrar en una especie de tubo horizontal de manera óptima incluyendo la gran mayoría posible de datos. Por lo cual descriptivamente rechazamos que se cumpla la homocedasticidad de los residuales.
2.3.2 Inferencialmente
Usamos el test de Breusch-Pagan para determinar si se puede asumir homocedasticidad. Se sabe que lo que hay detrás de dicho test, es una prueba de hipótesis de la siguiente forma:
\[ \begin{equation*} \begin{cases} H_0: \epsilon_{1}, \epsilon_{2}, \cdots, \epsilon_{n} \text{ tienen una misma varianza } \sigma^{2} \\ H_1: \epsilon_{1}, \epsilon_{2}, \cdots, \epsilon_{n}\text{ NO tienen una misma varianza } \sigma^{2} \end{cases} \end{equation*} \]
bptest(regresion1)##
## studentized Breusch-Pagan test
##
## data: regresion1
## BP = 1146.7, df = 1, p-value < 2.2e-16Sin embargo, como p-valor es menor a 0.05 NO se puede asumir homocedasticidad en los residuales.
2.4 Normalidad de los residuales
2.4.1 Descriptivamente
Notemos que, descriptivamente a partir del QQ-Plot graficado, no hay un buen ajuste para los valores iniciales y tampoco para los finales, luego se decide que NO hay normalidad en los residuales.
2.4.2 Inferencialmente
Aplicamos el test de Anderson-Darling sobre los residuales de la regresión. Esto a su vez implica que se considere una prueba hipótesis de la siguiente forma:
\[ \begin{equation*} \begin{cases} H_0: \epsilon_{1}, \epsilon_{2}, \cdots, \epsilon_{n}, \epsilon \sim N(0, \sigma^{2}) \\ H_1: \epsilon_{1}, \epsilon_{2}, \cdots, \epsilon_{n}, \epsilon \nsim N(0, \sigma^{2}) \end{cases} \end{equation*} \]
ad.test(regresion1$residuals)##
## Anderson-Darling normality test
##
## data: regresion1$residuals
## A = 130.83, p-value < 2.2e-16Al tener un p-valor menor a 0.05 se decide rechazar la hipótesis nula, es decir se decide que NO hay Normalidad en los residuales.
2.5 ¿El modelo se ajusta a los datos?
2.5.1 Descriptivamente
Para verificarlo, usamos el coeficiente de determinación (el cual es el coef. de correlación usado, elevado al cuadrado). Este mide la bondad de ajuste del MRLS, es decir de la recta de regresión.
resumen_reg1 <- summary(regresion1)
resumen_reg1$r.squared## [1] 0.6249046Sabiendo que un coef. de determinación indica un ajuste adecuado cuando este está entre 0.8 y 0.95, decidimos que nuestro MRLS NO es adecuado para los datos.
2.5.2 Inferencialmente
Lo comprobaremos mediante Bootstrap en la siguiente sección.
Procedimiento 3
Determine usando el coeficiente de determinación si el modelo de regresión lineal simple es adecuado o no para los datos. Haga la siguiente prueba de hipótesis usando bootstrap
\[ \begin{equation*} H_0: R^{2} \leq 0.8 \\ H_1: R^{2} > 0.8 \end{equation*} \]
donde \[ R^{2} \] es el coeficiente de determinación ¿se rechaza o no que el modelo de regresión lineal simple es adecuado para los datos?
Bootstrap
Se usa la técnica no paramétrica de re-muestreo conocida como Bootstrap dado que no contamos con un estadístico de prueba asociado al coeficiente de determinación.
Considerando la prueba de hipótesis dada:
\[ \begin{equation*} \begin{cases} H_0: R^{2} \leq 0.8 \\ H_1: R^{2} > 0.8 \end{cases} \end{equation*} \]
Definimos la siguiente función:
R2 <- function(data, indices){
regresion <- lm(data[,1][indices] ~ data[,2][indices])
resumen <- summary(regresion)
return(resumen$r.squared)
}Aplicamos la función anteriormente definida:
R2(nat, c(1,2,3))## [1] 0.9156823Haciendo uso de la librería boot realizamos las n simulaciones para el Bootstrap
library(boot)
muestras <- boot(data = nat, statistic = R2, R= 1000)
muestras##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = nat, statistic = R2, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.6249046 -1.562752e-05 0.00213862A continuación se muestran los resultados de las simulaciones, es decir el coef. de determinación hechas pruebas.
## [1] "table table-condensed table-bordered"Además se presenta el histograma donde se puede observar (descriptivamente) la distribución que sigue el coeficiente de determinación.
Se determina a continuación dónde cae nuestro estadístico de prueba t, y se establece cuál es la región crítica para nuestra prueba de hipótesis
DECISIÓN: Como nuestro estadístico de prueba t NO cae fuera de la región de confianza, entonces no se rechaza la hipótesis nula, lo cual sugiere que no aceptemos la hipótesis alternativa (que era la que nos interesaba aceptar). Por lo cual el coeficiente de determinación NO es mayor que 0.8 y descartamos que este sea sea adecuado para nuestros datos.
Procedimiento 4
Haga un diagrama de dispersión donde se muestren las bandas de pronóstico (bandas de confianza para datos observados) y las bandas de predicción (bandas de confianza para datos NO observados). ¿Hay muchos datos por fuera de las bandas de confianza de predicción?
Recordemos que los parámetros de la recta de regresión de nuestro modelo lineal pueden ser hallados mediante:
regresion1$coefficients## (Intercept) nat$talla_nac
## -3891.0462 141.4657De modo que esta es:
\[ \begin{align*} y= 141.4657x -3891.0462 + \epsilon \end{align*} \]
Determinamos un intervalo de confianza para los parámetros de la recta:
confint(regresion1)## 2.5 % 97.5 %
## (Intercept) -3914.8709 -3867.221
## nat$talla_nac 140.9813 141.9504.1 Bandas de confianza para datos observados
A través de la función predict definimos la banda de pronóstico (banda de datos observados) y especificamos el argumento interval = "confidence"
bandas_pron <- predict(regresion1, data.frame(nat$peso_nac),
interval = "confidence")
head(bandas_pron)## fit lwr upr
## 1 1484.650 1479.103 1490.196
## 2 2757.841 2756.150 2759.531
## 3 2899.306 2897.856 2900.757
## 4 3182.238 3180.823 3183.653
## 5 3606.635 3604.318 3608.952
## 6 2899.306 2897.856 2900.757
4.2 Bandas de confianza para datos NO observados
Por medio de la función predict definimos la banda de predicción (banda de datos NO observados) y especificamos el argumento interval = "prediction"
bandas_pred <- predict(regresion1, data.frame(nat$peso_nac),
interval = "prediction")## fit lwr upr
## 1 1484.650 887.0468 2082.252
## 2 2757.841 2160.2612 3355.420
## 3 2899.306 2301.7275 3496.885
## 4 3182.238 2584.6590 3779.817
## 5 3606.635 3009.0532 4204.216
## 6 2899.306 2301.7275 3496.885
Basados en este último gráfico afirmamos que sí hay una gran cantidad de datos atípicos fuera de las bandas de pronóstico y predicción.
Procedimiento 5
Haga un diagrama de dispersión donde se muestren las bandas de pronóstico y las bandas de predicción pero esta vez utilizando el modelo de regresión lineal ponderada, utilizando como pesos estimados
\[ w_{i}= \dfrac{1}{\epsilon^{2}_{i}} \]
con
\[ |\epsilon| = aY+b+\epsilon \text{ tal que } \epsilon \sim N(0, \delta^{2}) \]
es decir, usando los pesos dados que como se toman en el siguiente código de ejemplo
# pesos <- 1/fitted(lm(abs(residuals(regresion_lineal)) ~ fitted(regresion_lineal)))^25.1 Regresión lineal ponderada
Definimos los pesos o ponderaciones mediante:
w_i <- 1/(fitted(lm(abs(residuals(regresion1)) ~ fitted(regresion1)))^{2})Luego, hacemos la regresión lineal ponderada especificando el argumento weights :
regresion_ponderada <- lm(nat$peso_nac ~ nat$talla_nac, weights = w_i)De esta manera, el gráfico de regresión lineal basada en pesos, es:
En adición mediante summary(regresion_ponderada) observamos el resumen de la regresión realizada
summary(regresion_ponderada)##
## Call:
## lm(formula = nat$peso_nac ~ nat$talla_nac, weights = w_i)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -11.4653 -0.8442 -0.0378 0.7987 13.8427
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3902.473 12.569 -310.5 <2e-16 ***
## nat$talla_nac 141.698 0.255 555.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.281 on 196688 degrees of freedom
## Multiple R-squared: 0.6109, Adjusted R-squared: 0.6109
## F-statistic: 3.088e+05 on 1 and 196688 DF, p-value: < 2.2e-165.2 Bandas de confianza para datos observados
Especificando el argumento interval="confidence" hallamos las bandas de pronóstico bajo regresión ponderada haciendo uso de la función predict()
bandas_pron_pond <- predict(regresion_ponderada, data.frame(nat$peso_nac),
interval = "confidence")
head(bandas_pron_pond)## fit lwr upr
## 1 1482.061 1476.297 1487.825
## 2 2757.345 2755.604 2759.087
## 3 2899.044 2897.569 2900.518
## 4 3182.440 3181.040 3183.841
## 5 3607.535 3605.214 3609.856
## 6 2899.044 2897.569 2900.518De modo que, el diagrama de dispersión solamente con las bandas de pronóstico es:
5.3 Bandas de confianza para datos NO observados
Especificando el argumento interval="prediction" hallamos las bandas de predicción bajo regresión ponderada haciendo uso de la función predict()
bandas_pred_pond <- predict(regresion_ponderada, data.frame(nat$peso_nac),
interval = "prediction")
head(bandas_pred_pond)## fit lwr upr
## 1 1482.061 1475.774 1488.348
## 2 2757.345 2754.290 2760.401
## 3 2899.044 2896.132 2901.955
## 4 3182.440 3179.565 3185.315
## 5 3607.535 3604.116 3610.954
## 6 2899.044 2896.132 2901.955De manera que, el diagrama de dispersión es:
Nótese que las bandas, tanto la de prónostico como la de predicción se superponen, esto dado que los extremos de ambas bandas del iésimo dato tienen valores muy cercanos. Si pudieramos hacer zoom, veríamos de mejor manera la diferencia; es relevante mencionar que se podría lograr con la librería plotly sin embargo son muchos datos y ralentiza el computador desde donde se está originando este documento, por lo cual resulta imposible ver con más detalle la diferencia en las bandas.
Como notamos, son demasiados datos los que quedan por fuera de la región que encierran las bandas, es por ello que el modelo de regresión lineal ponderado, para este caso también resulta impreciso y no se ajusta a los datos (también es importante mencionar que hay un coeficiente de determinación igual a 0.6109).