TP3_Estadistica
Daniela Rattia
28/9/2021
Instalamos / invocamos las librerías a utilizar
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.0.3 v dplyr 1.0.2
## v tidyr 1.1.2 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0## Warning: package 'ggplot2' was built under R version 4.0.5## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()# install.packages("ggpubr")
# install.packages("mosaic")
# install.packages("bootstrap")
library(dplyr)
library(ggpubr)## Warning: package 'ggpubr' was built under R version 4.0.5library(mosaic)## Warning: package 'mosaic' was built under R version 4.0.5## Registered S3 method overwritten by 'mosaic':
## method from
## fortify.SpatialPolygonsDataFrame ggplot2##
## The 'mosaic' package masks several functions from core packages in order to add
## additional features. The original behavior of these functions should not be affected by this.##
## Attaching package: 'mosaic'## The following object is masked from 'package:Matrix':
##
## mean## The following objects are masked from 'package:dplyr':
##
## count, do, tally## The following object is masked from 'package:purrr':
##
## cross## The following object is masked from 'package:ggplot2':
##
## stat## The following objects are masked from 'package:stats':
##
## binom.test, cor, cor.test, cov, fivenum, IQR, median, prop.test,
## quantile, sd, t.test, var## The following objects are masked from 'package:base':
##
## max, mean, min, prod, range, sample, sumlibrary(bootstrap)## Warning: package 'bootstrap' was built under R version 4.0.3Importar datasets
A los efectos de este trabajo práctico, utilizaremos datos de las localidades de Venado Tuerto y Avellaneda de la Provincia de Santa Fe. Los datos han sido obtenidos de la página poblaciones.org
#setwd("C:/Users/u621723/Desktop/TP_PEC")
venado_tuerto <- read.csv("VenadoTuerto.csv", encoding = "UTF-8")
reconquista_avellaneda <- read.csv("ReconquistaAvellaneda.csv", encoding = "UTF-8")CONSIGNA 1
Seleccion de atributos
Por jefe de hogar, seleccionamos la variable: con primaria completa o secundario incompleto. El interés es explorar si el hecho de que el o la jefa de hogar sin secundario completo (nivel educativo) es útil para estratificar muestras de radios censales y hogares
#Venado Tuerto
jefes_hogar_vt <- venado_tuerto %>%
select(1,2,34,prim_secincompleto="Hogares.con.jefes.con.primaria.completa.o.secundaria.incompleto")
#Reconquista Avellaneda
jefes_hogar_ra <- reconquista_avellaneda %>%
select(1,2,34,prim_secincompleto="Hogares.con.jefes.con.primaria.completa.o.secundaria.incompleto")Cálculo de medidas
para conocer las distribuciones, calculo primero la media y la desviación de ambos dataset
# media o promedio
media_vt <- mean(jefes_hogar_vt$prim_secincompleto)
media_ra <- mean(jefes_hogar_ra$prim_secincompleto)
# desviación estandar
desv_vt <- sd(jefes_hogar_vt$prim_secincompleto)
desv_ra <- sd(jefes_hogar_ra$prim_secincompleto) Ya de entrada observamos que es mayor el primedio de jefes con este atributo en Reconquista - Avellaneda, pero esto aún no nos indica mucho más. Observemos un gráfico de distribución para ver mejor el comportamiento.
Grafico de distribución
Hacemos un gráfico de distribución normal de la población de cada localidad, para comprender cómo se distribuye
# Venado Tuerto
gvt1 <- ggplot(jefes_hogar_vt,
aes(prim_secincompleto))+
labs(title = "Distribuciones (normal y empírica) - Venado Tuerto",
subtitle = "Jefes y jefas de hogar con Secundario Incompleto")+
stat_ecdf(geom = "step", pad = F)+
stat_function(fun = pnorm, color="purple", args = list(mean = media_vt,
sd= desv_vt))
# Reconquista Avellaneda
gra1 <- ggplot(jefes_hogar_ra,
aes(prim_secincompleto))+
labs(title = "Distribuciones (normal y empírica) - Reconquista_Avenllaneda",
subtitle = "Jefes y jefas de hogar con Secundario Incompleto")+
stat_ecdf(geom = "step", pad = F)+
stat_function(fun = pnorm, color="purple", args = list(mean = media_ra,
sd= desv_ra))
gvt1
gra1
En negro se observa, para ambos casos, la distribución empírica, es decir, que proviene de los datos, mientras que en violeta (construido con la función pnorm) la curva de la distribución normal. Hasta el momento, se observa que no son símiles en alguna medida.
Veamos los histogramas
para comprender la distribución normal
# Venado Tuerto
gvt2 <- ggplot(jefes_hogar_vt, aes(x= prim_secincompleto))+
geom_histogram(aes(y = ..density..),
fill = "purple",
alpha = 0.5,
col = "black", bins= 10) +
labs(title = "Distribución Jefes/as Hogar por Sec. Incompleto, Venado Tuerto",
x = "x", y = "")+
stat_function(fun= dnorm,
color="blue",
size= 1,
args=list(mean=media_vt,sd=desv_vt))
# Reconquista Avellaneda
gra2 <- ggplot(jefes_hogar_ra, aes(x= prim_secincompleto))+
geom_histogram(aes(y = ..density..),
fill = "purple",
alpha= 0.5,
col = "black", bins= 10) +
labs(title = "Distribución Jefes/as Hogar por Sec. Incompleto, Reconquista Avellaneda",
x = "x", y = "")+
stat_function(fun= dnorm,
color="blue",
size= 1,
args=list(mean=media_ra,sd=desv_ra))
gvt2
gra2
Se observa desde ahora, que en el caso de Venado Tuerto, la distribución es asimétrica a la derecha ligeramente, mientras que la distribución en Reconquista-Avellenada tiene mayor variabilidad y que destaca en meyor medida la asimetría.
Gráficos de cuantil-cuantil
veamos gráficos qq para entender en qué medida está cerca está la distribución de un conjunto de datos a alguna distribución ideal
# venado Tuerto
gvt3 <- ggplot(jefes_hogar_vt,
aes(sample = prim_secincompleto)) +
stat_qq(size = 0.8, color="purple") +
stat_qq_line(size=0.5,
alpha= 0.5) +
labs(title = "qq - Venado Tuerto: cuánto se acerca a la distribución ideal",
x = "q-norm")
# Reconquista - Avellaneda
gra3 <- ggplot(jefes_hogar_ra,
aes(sample = prim_secincompleto)) +
stat_qq(size = 0.8, color="purple") +
stat_qq_line(size=0.5,
alpha= 0.5) +
labs(title = "qq - Reconquista Avellaneda: cuánto se acerca a la distribución ideal",
x = "q-norm")
gvt3
gra3
Observamos ahora con mayor claridad que la distribución (a nivel comparativo) no alcanza similitud, y se ve una disperción aún mayor en el caso de Reconquista Avellaneda (los puntos en púrpura, como distribución propia de los datos con respecto a la distribución ideal teórica, la línea negra)
Estimación suavizada
Vemos ahora los gráficos de los dos casos utilizando la técnica de suavizado, para poder aprenciar (a simple vista) cómo se comporta la normal en ambas series de forma alisada y minimizando el efecto de los valores atípicos.
# Venado Tuerto
gvt4 <- ggplot(jefes_hogar_vt,
aes(x = prim_secincompleto)) +
stat_density(fill = "purple",
alpha= 0.5)+
labs(title="Venado Tuerto - Est. suavizada",
x = "Q", y = " ")+
stat_function(fun = dnorm, #declaramos la distribución normal contra la cual se hará el suavizado
color= "black",
args=list(mean = media_vt,
sd=desv_vt))
gra4 <- ggplot(jefes_hogar_ra,
aes(x = prim_secincompleto)) +
stat_density(fill = "purple",
alpha= 0.5)+
labs(title="Reconquista_Avellaneda - Est. suavizada",
x = "Q", y = " ")+
stat_function(fun = dnorm, #declaramos la distribución normal contra la cual se hará el suavizado
color= "black",
args=list(mean = media_ra,
sd=desv_ra))
gvt4
gra4
Notemos que, como ya se había apreciado en el gráfico anterior, las distribuciones (empíricas y teóricas) son disímiles, pudiendose evidenciar mejor en éste último gráfico.
Esto nos indica, que probablemente la variable de Jefes y Jefas de hogar con secundario incompleto no sea el mejor estimador, o bien, utilizando estas técnicas de exploración se evidencia que los datos (observaciones empíricas) no se ajustan a los modelos normales.
En este sentido, es una posibilidad que, para obtener un buen criterio de estratificación de muestras, se consideren otras variables socioconómicas, con ciertas condiciones materiales de la viviendas o, más adelante, entender si se puede hacer una definición territorial, es decir, que la distribución espacial de los hogares impacte de alguna forma en los resultados de los modelos o la generación de un índice de estratificación espacial. (quizá para la otra materia)
CONSIGNA 2
En este punto, veremos si la misma variable de nivel educativo de los y las jefas de hogar es un buen predictor para saber si un hogar cuenta con cocina por garrafa o leña.
selección de las variables
#Venado Tuerto
garrafa_vt <- venado_tuerto %>%
select(1,2,34,prim_secincompleto="Hogares.con.jefes.con.primaria.completa.o.secundaria.incompleto", garrafa= "Hogares.con.garrafa.o.leña.como.combustible.usado.principalmente.para.cocinar")
# Reconquista - Avellaneda
garrafa_ra <- reconquista_avellaneda %>%
select(1,2,34,prim_secincompleto="Hogares.con.jefes.con.primaria.completa.o.secundaria.incompleto", garrafa= "Hogares.con.garrafa.o.leña.como.combustible.usado.principalmente.para.cocinar")Corelacionan ambas variables?
?cor## starting httpd help server ... donecor(garrafa_vt$prim_secincompleto, garrafa_vt$garrafa)## [1] 0.8346137# 0.8346137
cor(garrafa_ra$prim_secincompleto, garrafa_ra$garrafa)## [1] 0.9646976# 0.9646976Observamos en primer lugar que las variables correlacionan con alto grado de dependencia: +0.83 para Venado Tuerto y +0.96 para Reconquista Avellaneda. Esto quiere decir que ambas variables están fuertemente correlacionadas, sin embargo, no implica causalidad. Por lo cual, deberemos observar otras medidas para asumir a la variable X como buen predictor.
Graficamente:
Elaboro un gráfico tipo scatterplot para observar el resultado anterior
#Venado Tuerto
gvt5 <- ggplot(garrafa_vt,
aes(x = prim_secincompleto, y = garrafa)) +
geom_point(colour="purple")+
labs(title = "Venado Tuerto - corr: sec.incompleto ~ cocina a garrafa o leña")
# Reconquista Avellaneda
gra5 <- ggplot(garrafa_ra,
aes(x = prim_secincompleto, y = garrafa)) +
geom_point(colour="purple")+
labs(title = "Reconquista Avellaneda - corr: sec.incompleto ~ cocina a garrafa o leña")
gvt5
gra5
Se observa claramente en el gráfico que, para el caso de Reconquista Avellaneda, la relación es positiva y muy fuerte, como se habia detallado en el apartado anterior. Para Venado Tuerto, es tambiepn positiva, sin embargo el comportamiento es más difuso con respecto a la otra localidad.
Generamos un modelo de regresión lineal
Dado que el coeficiente R es cercano a 1 y nos indica que las dos variables están correlacionadas, conviene que generemos un modelo de regresión lineal simple para ver en qué medida lo están.
#Venado Tuerto
lm_vt <- lm(data = garrafa_vt, garrafa~prim_secincompleto)
# Reconquista
lm_ra <- lm(data = garrafa_ra, garrafa~prim_secincompleto)
lm_vt##
## Call:
## lm(formula = garrafa ~ prim_secincompleto, data = garrafa_vt)
##
## Coefficients:
## (Intercept) prim_secincompleto
## -92.461 0.963lm_ra##
## Call:
## lm(formula = garrafa ~ prim_secincompleto, data = garrafa_ra)
##
## Coefficients:
## (Intercept) prim_secincompleto
## 88.526 1.083En las medidas detalladas, Intercepto nos muestra dónde corta la recta en el modelo, pero vemos con más interés lo que nos dice el coficiente de variación.
En el modelo, se intenta estimar el efecto de la variable “garrafa” cuando se incrementa la cantidad de “jefes de hogares con sec. incompleto”. Dicho de manera abstracta: cada vez que aumenta una unidad de “jefe de hogar etc”, la variable “garrafa” aumenta 0.96 y 1.083 veces para Venado Tuerto y Reconquista-Avellaneda respectivamente.
Ambos valores rondan la unidad.
R cuadrado o Coeficiente de determinacion
Nos permitirá saber el nivel de ajuste (a la variable que se desea explicar) del modelo elaborado en el apartado anterior.
summary(lm_vt)$r.squared## [1] 0.69658summary(lm_ra)$r.squared## [1] 0.9306414PAra el caso de Venado tuerto, R2 es de 0.69 mientras que para Reconquista es de 0.93. En el primer caso, R2 es menor que R, y se considera bajo. Por otra parte, en Reconquiste, el valor de R2 es mayor, y esto nos indica, grosso modo, que es un modelo cuyas estimaciones se ajustan bastante bien a la variable real, o mejor, podría decir que el modelo explica en un 93% a la variable real.
Graficar la recta del modelo de regresion
#Venado Tuerto
gvt6 <- ggplot(garrafa_vt,
aes(x = prim_secincompleto, y = garrafa)) +
geom_point( color="purple") +
geom_smooth(method=lm , color="black") +
labs(title = "LM de Venado Tuerto - secundario incompleto vs hogar con cocina a garrafa o leña")
# Reconquista - Avellaneda
gra6 <- ggplot(garrafa_ra,
aes(x = prim_secincompleto, y = garrafa)) +
geom_point( color="purple") +
geom_smooth(method=lm , color="black") +
labs(title = "LM de Reconquista Avellaneda - secundario incompleto vs hogar con cocina a garrafa o leña")
gvt6## `geom_smooth()` using formula 'y ~ x'
gra6## `geom_smooth()` using formula 'y ~ x'
Lo que obsevamos es que el modelo se ajusta mejor (sobre la recta y los intervalos de confianza) en los valores más bajos (de hogares) para el caso de Venado Tuerto, mientras que para Reconquista, los valores se ajustan en gran medida a lo largo de toda la recta. Esto nos indica que, el hecho de que el nivel educativo de jefes y jefas en un hogar, prediga el tipo de combustion en la cocina, es más óptimo (o ajusta mejor) en el caso de Reconquista que de Venado Tuerto.
Revisamos el resultado del modelo lineal
# Venado Tuerto
summary(lm_vt)##
## Call:
## lm(formula = garrafa ~ prim_secincompleto, data = garrafa_vt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -155.68 -29.46 -10.02 31.64 169.31
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -92.46136 15.95406 -5.795 1.46e-07 ***
## prim_secincompleto 0.96304 0.07291 13.209 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 63.08 on 76 degrees of freedom
## Multiple R-squared: 0.6966, Adjusted R-squared: 0.6926
## F-statistic: 174.5 on 1 and 76 DF, p-value: < 2.2e-16# Reconquista
summary(lm_ra)##
## Call:
## lm(formula = garrafa ~ prim_secincompleto, data = garrafa_ra)
##
## Residuals:
## Min 1Q Median 3Q Max
## -88.038 -20.369 -1.815 23.227 132.020
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 88.52620 9.39591 9.422 5.12e-14 ***
## prim_secincompleto 1.08264 0.03558 30.427 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 38.89 on 69 degrees of freedom
## Multiple R-squared: 0.9306, Adjusted R-squared: 0.9296
## F-statistic: 925.8 on 1 and 69 DF, p-value: < 2.2e-16en ambos casos existe un nivel de significancia muy alto es decir, con menos de 0.01% de error. Sin embargo, para el caso de Reconquista, el modelo se ajusta en mayor medida, como fue mencionado anteriormente.
Haciendo un doble clic, lo que podemos rechazar es la Ho (Hipótesis Nula), y considerar, en efecto, y que el fenomeno corrleativo entre el nivel educativo el tipo de combustión para cocinar en el hogar no son producto de la aletoreidad, al menos para este modelo. Lo que acabamos de analizar, no explica las razones de causalidad entre una variable y la ocurrencia de la otra, pero si nos indica con certeza la relación entre ambas.
Podemos concluir que el modelo lineal se ajusta y es óptimo para el análisis.
Residuos
residuos_vt <- residuals(lm_vt)
residuos_ra <- residuals(lm_ra)
residuos_vt## 1 2 3 4 5 6
## 85.8634628 89.5656163 94.4177945 147.7719281 81.4525470 21.2197994
## 7 8 9 10 11 12
## -90.1216115 80.6241064 147.4696757 39.1392796 37.4762846 -40.7779976
## 13 14 15 16 17 18
## 21.8326212 -8.9649778 62.2959133 -140.4172550 -37.0410917 -117.3780966
## 19 20 21 22 23 24
## -31.7757946 -2.9953243 1.7437846 19.1155420 -84.9280223 -28.1889134
## 25 26 27 28 29 30
## -76.5932204 -128.5672798 -14.4387897 -70.2236659 -13.6627253 -12.1039877
## 31 32 33 34 35 36
## -41.0322798 -11.0670322 -13.7344332 -155.6759431 34.7525965 -18.1039877
## 37 38 39 40 41 42
## -51.2192600 -16.7498540 3.3786361 -24.0300768 -11.9561659 39.5308638
## 43 44 45 46 47 48
## -11.3996312 33.8634628 -12.8822550 6.0090816 -8.9561659 -19.8452996
## 49 50 51 52 53 54
## 32.1960618 -9.5104976 -12.2931709 -29.8844580 19.4547499 -19.1801015
## 55 56 57 58 59 60
## -14.8083441 -47.8844580 16.6025717 -54.1889134 -1.1497550 -1.2909679
## 61 62 63 64 65 66
## 0.6003688 -61.0714382 169.3086360 -28.0736412 4.4828935 -67.5171065
## 67 68 69 70 71 72
## -25.9605718 78.3913588 18.0416311 -10.5193095 29.9892549 -5.3692847
## 73 74 75 76 77 78
## -33.5628738 -66.1150025 135.8392301 63.9522994 72.7415816 94.0873985residuos_ra## 1 2 3 4 5 6
## -9.7407638 22.4748959 -8.5336093 6.3667375 -10.7241915 38.8294497
## 7 8 9 10 11 12
## 6.4001012 60.2675224 31.2592362 17.8460221 -4.3275497 -54.1864657
## 13 14 15 16 17 18
## -9.1371866 36.8625944 -24.4016874 -1.8151205 -0.3109774 -17.3441221
## 19 20 21 22 23 24
## -8.9143356 79.5485955 132.0195939 43.6973088 -0.1622641 -4.7904809
## 25 26 27 28 29 30
## -18.3194825 5.8623755 -17.5342663 -24.4104115 -47.4762630 -28.2285535
## 31 32 33 34 35 36
## 9.8211635 -4.8482652 71.4247408 13.0608058 44.0690920 -88.0381904
## 37 38 39 40 41 42
## -51.4021254 -73.3521893 38.6063798 36.8874530 -20.2697655 -55.7732516
## 43 44 45 46 47 48
## 6.2838757 33.1683072 -42.1369676 -63.4847682 -28.2863378 34.2511690
## 49 50 51 52 53 54
## -20.4673199 28.4413132 -44.9636148 -10.2531931 -14.1705503 -6.8234066
## 55 56 57 58 59 60
## 2.9700958 -59.1042608 2.7228243 65.9700958 1.3089533 9.6973088
## 61 62 63 64 65 66
## -2.6993329 34.4576665 15.3338118 -5.3853340 57.0778162 -1.3358359
## 67 68 69 70 71
## 23.9786010 -23.7900429 1.4083874 -40.6825416 -64.7492690Incorporo los resultados a los df
Usando mutate, incorporo el resutlado de los residuos para cada observación (cada unidad de observación o RC) para tener en el mismo dataset los datos reales y el cálculo de errores a partir del modelo
#Venado tuerto
garrafa_vt <- garrafa_vt %>%
mutate(res=residuos_vt)
# Reconquista Avellaneda
garrafa_ra <- garrafa_ra %>%
mutate(res=residuos_ra)Graficar los residuos y ver cómo se comportan: ## Homocedasticidad de los residuos
#Venado Tuerto
gvt7 <- ggplot(garrafa_vt,
aes(x= prim_secincompleto, y=res)) +
geom_point() +
geom_hline(yintercept = 0, col = "purple") +
labs(title="Venado Tuerto - Residuos del modelo lineal")
#Reconquista Avellaneda
gra7 <- ggplot(garrafa_ra,
aes(x= prim_secincompleto, y=res)) +
geom_point() +
geom_hline(yintercept = 0, col = "purple") +
labs(title="Reconquista Avellaneda - Residuos del modelo lineal")
gvt7
gra7
Observamos errores más dispersos en el caso de Reconquista, mientras que en Venado Tuerto vemos un comportamiento a nivel visual que invita a explorar con más profundidad.
Distribucion normal de los errores
gvt8 <- ggplot(garrafa_vt,
aes(sample=res)) +
stat_qq(size = 0.5) +
stat_qq_line(color="purple") +
labs(title = "Venado Tuerto - Residuos por qq normal")
gra8 <- ggplot(garrafa_ra,
aes(sample=res)) +
stat_qq(size = 0.5) +
stat_qq_line(color="purple") +
labs(title = "Reconquista Avellaneda - Residuos por qq normal")
gvt8
gra8
Si bien se observa un comportamiento similar en la escala de este gráfico, lo que podemos observar, usando quintil-quintil de dist. normal, es que los errores no se distribuyen de forma normal.
Valores atípicos
Finalmente exploramos los valores atípicos, considerando +100 t -100 para venado tuerto
filtro_outlier_vt <- garrafa_vt %>%
filter(res >=100 | res <=-100)
filtro_outlier_ra <- garrafa_ra %>%
filter(res>=100 | res <=-60)Graficamos
gvt9 <- ggplot(filtro_outlier_vt)+
geom_point(aes(x=garrafa, y=res), color="purple")+
geom_hline(yintercept = 0, col = "black") +
labs(title = "Errores Atípicos en Venado Tuerto")
gra9 <- ggplot(filtro_outlier_ra)+
geom_point(aes(x=garrafa, y=res), color="purple")+
geom_hline(yintercept = 0, col = "black") +
labs(title = "Errores Atípicos en Reconquista Avellaneda")
gvt9
gra9
Vistos los gráficos anteriores, se puede concluir que, si bien las variables están relacionadas fuertemente, y el modelo ajusta de forma óptima, los errores no presentan una distribución normal, presentan valores atípicos y además tiene un comportamiento que se distancia de la línea teorética y no son independientes.
De esta forma, podemos saber que los errores distorsionan el modelo en alguna medida, y que la presencia de outliers (propios de los datos o por error real) trae problemas a futuro si se desea usar estas variables como predictoras.
Comos e mencionó anteriormente, convendría ver si existe algun tipo de dependencia espacial que pueda explicar estos errores (como la presencia de ruralidad o no, nivel socioecnómico complejo o en la definición de otras características del hogar)
Veremos en la materia de Geoestadística si esto nos puede resolver algo más.
ADRM