Práctica 3. Introducción a la inferencia
AE
26/06/2020
Previo
Paquetería
if (!require("pacman")) install.packages("pacman") # instala pacman si se requiere## Loading required package: pacman#cargan los paquetes necesarios para la práctica
pacman::p_load(tidyverse, haven, readxl, janitor, srvyr,
RColorBrewer, wesanderson, sjlabelled, esquisse, skimr,
GGally, DescTools)Directorio
setwd("/cloud/project/datos")Base de datos
Vamos a importar la base completa del cuestionario sociodemográfico de la Encuesta Nacional de Ocupación y Empleo, trimestre III de 2019
SDEMT319_10 <- read_dta("SDEMT319_10.dta", encoding="latin1")Hipótesis e intervalos de confianza
t-test
Este comando nos sirve para calcular diferentes tipos de test, que tienen como base la distribución t
Univariado para estimación
SDEMT319_10 %>%
filter(ing_x_hrs>0) %>% # filtro para quienes tiene ingreso
with(t.test(ing_x_hrs))##
## One Sample t-test
##
## data: ing_x_hrs
## t = 21.304, df = 924, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 39.59916 47.63523
## sample estimates:
## mean of x
## 43.61719Univariado para hipótesis específica
SDEMT319_10 %>%
filter(ing_x_hrs>0) %>% # filtro para quienes tiene ingreso
with(t.test(ing_x_hrs, mu=50))##
## One Sample t-test
##
## data: ing_x_hrs
## t = -3.1176, df = 924, p-value = 0.00188
## alternative hypothesis: true mean is not equal to 50
## 95 percent confidence interval:
## 39.59916 47.63523
## sample estimates:
## mean of x
## 43.61719SDEMT319_10 %>%
filter(ing_x_hrs>0) %>% # filtro para quienes tiene ingreso
with(t.test(ing_x_hrs, mu=50, alternative = "two.sided"))#default y de dos colas##
## One Sample t-test
##
## data: ing_x_hrs
## t = -3.1176, df = 924, p-value = 0.00188
## alternative hypothesis: true mean is not equal to 50
## 95 percent confidence interval:
## 39.59916 47.63523
## sample estimates:
## mean of x
## 43.61719SDEMT319_10 %>%
filter(ing_x_hrs>0) %>% # filtro para quienes tiene ingreso
with(t.test(ing_x_hrs, mu=50, alternative = "less")) # cola izquierda##
## One Sample t-test
##
## data: ing_x_hrs
## t = -3.1176, df = 924, p-value = 0.00094
## alternative hypothesis: true mean is less than 50
## 95 percent confidence interval:
## -Inf 46.9882
## sample estimates:
## mean of x
## 43.61719SDEMT319_10 %>%
filter(ing_x_hrs>0) %>% # filtro para quienes tiene ingreso
with(t.test(ing_x_hrs, mu=50, alternative = "greater")) #cola derecha ##
## One Sample t-test
##
## data: ing_x_hrs
## t = -3.1176, df = 924, p-value = 0.9991
## alternative hypothesis: true mean is greater than 50
## 95 percent confidence interval:
## 40.24619 Inf
## sample estimates:
## mean of x
## 43.61719Proporciones
Una opción es utilizar el comando prop.test, éste necesita un objeto de tipo “table”.
freq.sex<-table(SDEMT319_10$sex)
freq.sex##
## 1 2
## 1992 2067prop.table(freq.sex)##
## 1 2
## 0.4907613 0.5092387prop.test(freq.sex)##
## 1-sample proportions test with continuity correction
##
## data: freq.sex, null probability 0.5
## X-squared = 1.3491, df = 1, p-value = 0.2454
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.4752751 0.5062651
## sample estimates:
## p
## 0.4907613Por default, toma los valores de la primera categoría, también podemos hacer la estimación con los datos del total de “éxitos” y el total de “intentos”. Calculemos para las mujeres
prop.test(2067,(1992+2067))##
## 1-sample proportions test with continuity correction
##
## data: 2067 out of (1992 + 2067), null probability 0.5
## X-squared = 1.3491, df = 1, p-value = 0.2454
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.4937349 0.5247249
## sample estimates:
## p
## 0.5092387prop.test(freq.sex, p=0.5)##
## 1-sample proportions test with continuity correction
##
## data: freq.sex, null probability 0.5
## X-squared = 1.3491, df = 1, p-value = 0.2454
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.4752751 0.5062651
## sample estimates:
## p
## 0.4907613prop.test(freq.sex, p=0.5, alternative = "two.sided") #default y de dos colas##
## 1-sample proportions test with continuity correction
##
## data: freq.sex, null probability 0.5
## X-squared = 1.3491, df = 1, p-value = 0.2454
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.4752751 0.5062651
## sample estimates:
## p
## 0.4907613prop.test(freq.sex, p=0.5, alternative = "less") # cola izquierda##
## 1-sample proportions test with continuity correction
##
## data: freq.sex, null probability 0.5
## X-squared = 1.3491, df = 1, p-value = 0.1227
## alternative hypothesis: true p is less than 0.5
## 95 percent confidence interval:
## 0.0000000 0.5037929
## sample estimates:
## p
## 0.4907613prop.test(freq.sex, p=0.5, alternative = "greater") #cola derecha ##
## 1-sample proportions test with continuity correction
##
## data: freq.sex, null probability 0.5
## X-squared = 1.3491, df = 1, p-value = 0.8773
## alternative hypothesis: true p is greater than 0.5
## 95 percent confidence interval:
## 0.477742 1.000000
## sample estimates:
## p
## 0.4907613También no podemos cambiar \(\alpha\) o bien el nivel de confianza (IC=100%*(1-\(\alpha\)))
prop.test(freq.sex, p=0.5, alternative = "greater", conf.level = 0.98) #cola derecha ##
## 1-sample proportions test with continuity correction
##
## data: freq.sex, null probability 0.5
## X-squared = 1.3491, df = 1, p-value = 0.8773
## alternative hypothesis: true p is greater than 0.5
## 98 percent confidence interval:
## 0.4745411 1.0000000
## sample estimates:
## p
## 0.4907613Estimaciones bivariadas
Diferencias de medias por grupos
¿Podemos decir, con significancia estadística que los valores medios de una variable son diferentes entre los grupos?
SDEMT319_10 %>%
filter(SDEMT319_10$clase2==1) %>%
group_by(as_label(sex)) %>%
summarise(avg_ing = mean(ing_x_hrs))## # A tibble: 2 x 2
## `as_label(sex)` avg_ing
## <fct> <dbl>
## 1 Hombre 23.1
## 2 Mujer 18.7SDEMT319_10 %>%
filter(SDEMT319_10$clase2==1) %>%
with(t.test(ing_x_hrs~sex))##
## Welch Two Sample t-test
##
## data: ing_x_hrs by sex
## t = 2.1773, df = 1733.5, p-value = 0.02959
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.4347619 8.3314034
## sample estimates:
## mean in group 1 mean in group 2
## 23.06306 18.67997Estimación de varianzas y sus pruebas de hipótesis
Para poder hacer inferencia sobre la varianza utilizamos el comando varTest() del paquete “DescTools”
SDEMT319_10 %>%
filter(clase2==1) %>%
with(VarTest(ing_x_hrs))##
## One Sample Chi-Square test on variance
##
## data: ing_x_hrs
## X-squared = 4482994, df = 1893, p-value < 2.2e-16
## alternative hypothesis: true variance is not equal to 1
## 95 percent confidence interval:
## 2224.280 2526.616
## sample estimates:
## variance of x
## 2368.195Podemos también decir algo sobre el valor objetivo de nuestra hipótesis
SDEMT319_10 %>%
filter(clase2==1) %>%
with(VarTest(ing_x_hrs, sigma.squared = 100000))##
## One Sample Chi-Square test on variance
##
## data: ing_x_hrs
## X-squared = 44.83, df = 1893, p-value < 2.2e-16
## alternative hypothesis: true variance is not equal to 1e+05
## 95 percent confidence interval:
## 2224.280 2526.616
## sample estimates:
## variance of x
## 2368.195Guardar como objeto nuestros resultados, siempres muy conveniente para pedir después o para realizar operaciones con ellos
test2<-SDEMT319_10 %>%
filter(clase2==1) %>%
with(VarTest(ing_x_hrs))
test2$conf.int## [1] 2224.280 2526.616
## attr(,"conf.level")
## [1] 0.95sqrt(test2$conf.int) ## sacamos la raíz cuadrada para tener las## [1] 47.16228 50.26545
## attr(,"conf.level")
## [1] 0.95#desviaciones estándar y sea más fácil de interpretarEstimación de diferencias de varianzas y sus pruebas de hipótesis
Para comparar varianza, usamos su “ratio”, esto nos da un estadístico de prueba F, para comparar dos muestras de poblaciones normales.
SDEMT319_10 %>%
filter(clase2==1) %>%
with(var.test(x=ing_x_hrs, y=eda, ratio=1))##
## F test to compare two variances
##
## data: ing_x_hrs and eda
## F = 11.457, num df = 1893, denom df = 1893, p-value < 2.2e-16
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 10.46994 12.53787
## sample estimates:
## ratio of variances
## 11.45735“x=” declara al numerador “y=” declara al denominador
SDEMT319_10 %>%
filter(clase2==1) %>%
with(var.test(x=ing_x_hrs, y=eda, ratio=1, conf.level = 0.98))##
## F test to compare two variances
##
## data: ing_x_hrs and eda
## F = 11.457, num df = 1893, denom df = 1893, p-value < 2.2e-16
## alternative hypothesis: true ratio of variances is not equal to 1
## 98 percent confidence interval:
## 10.29496 12.75098
## sample estimates:
## ratio of variances
## 11.45735Si lo que queremos es comparar la varianza entre dos grupos, usamos el signo ~
SDEMT319_10 %>%
filter(clase2==1) %>%
with(var.test(ing_x_hrs ~ as_label(sex), ratio=1))##
## F test to compare two variances
##
## data: ing_x_hrs by as_label(sex)
## F = 4.3853, num df = 1132, denom df = 760, p-value < 2.2e-16
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 3.847204 4.990242
## sample estimates:
## ratio of variances
## 4.385301Prueba chi-cuadrado chi-sq. Una aplicación más común
Cuando tenemos dos variables cualitativas o nominales podemos hacer esta la prueba chi-cuadrado, o prueba de independencia. Esta tiene una lógica un poco diferente a las pruebas que hacemos, porque proviene de comparar la distribución de los datos dado que no hay independencia entre las variables y los datos que tenemos. La hipótesis nula postula una distribución de probabilidad totalmente especificada como el modelo matemático de la población que ha generado la muestra, por lo que si la rechazamos hemos encontrado evidencia estadística sobre la dependencia de las dos variables.
table(SDEMT319_10$clase1, SDEMT319_10$sex)##
## 1 2
## 0 325 313
## 1 1188 810
## 2 479 944chisq.test(SDEMT319_10$clase1, SDEMT319_10$sex)##
## Pearson's Chi-squared test
##
## data: SDEMT319_10$clase1 and SDEMT319_10$sex
## X-squared = 222.38, df = 2, p-value < 2.2e-16Análisis de varianza
Análisis de varianza. Haremos la versión más simple. Para ver el efecto de un factor sobre una variable cualitativa (oneway). Revisaremos si la región de residencia de los trabajadores tiene un efecto en la distribución de los ingresos por trabajo.
Primero un gráfico
la ANOVA se basa en que nuestra variable es normal. Quitaremos los outliers
lienzo_bi <-SDEMT319_10 %>%
filter(clase2==1 & !ing_x_hrs==0) %>%
ggplot(aes(x=log(ing_x_hrs), fill=as_label(t_loc),
color=as_label(t_loc),
alpha=I(0.5)))
lienzo_bi + geom_density()
La prueba ANOVA o análisis de varianza, nos dice cuánto de nuestra variable se ve explicado por un factor. En los modelos es mul útil guardar nuestros resultados como un objeto
anova<-SDEMT319_10 %>%
filter(clase2==1) %>%
with(aov(ing_x_hrs ~ as_label(t_loc)))
summary(anova)## Df Sum Sq Mean Sq F value Pr(>F)
## as_label(t_loc) 2 46501 23250 9.91 5.23e-05 ***
## Residuals 1891 4436493 2346
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Comparación entre grupos
¿si es significativo cuáles diferencias entre los grupos lo son?
TukeyHSD(anova)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = ing_x_hrs ~ as_label(t_loc))
##
## $`as_label(t_loc)`
## diff
## Localidades de 15 000 a 99 999 habitantes-Localidades mayores de 100 000 habitantes 21.291771
## Localidades de 2 500 a 14 999 habitantes-Localidades mayores de 100 000 habitantes 1.077048
## Localidades de 2 500 a 14 999 habitantes-Localidades de 15 000 a 99 999 habitantes -20.214724
## lwr
## Localidades de 15 000 a 99 999 habitantes-Localidades mayores de 100 000 habitantes 10.07399
## Localidades de 2 500 a 14 999 habitantes-Localidades mayores de 100 000 habitantes -15.89403
## Localidades de 2 500 a 14 999 habitantes-Localidades de 15 000 a 99 999 habitantes -40.19000
## upr
## Localidades de 15 000 a 99 999 habitantes-Localidades mayores de 100 000 habitantes 32.5095530
## Localidades de 2 500 a 14 999 habitantes-Localidades mayores de 100 000 habitantes 18.0481263
## Localidades de 2 500 a 14 999 habitantes-Localidades de 15 000 a 99 999 habitantes -0.2394504
## p adj
## Localidades de 15 000 a 99 999 habitantes-Localidades mayores de 100 000 habitantes 0.0000268
## Localidades de 2 500 a 14 999 habitantes-Localidades mayores de 100 000 habitantes 0.9878583
## Localidades de 2 500 a 14 999 habitantes-Localidades de 15 000 a 99 999 habitantes 0.0465547Supuestos de ANOVA
- Las observaciones se obtienen de forma independiente y aleatoria de la población definida por los niveles del factor
- Los datos de cada nivel de factor se distribuyen normalmente.
- Estas poblaciones normales tienen una varianza común.
#Prueba Bartlett para ver si las varianzas son iguales
SDEMT319_10 %>%
filter(clase2==1) %>%
with(bartlett.test(ing_x_hrs ~ as_label(t_loc)))##
## Bartlett test of homogeneity of variances
##
## data: ing_x_hrs by as_label(t_loc)
## Bartlett's K-squared = 912.01, df = 2, p-value < 2.2e-16La prueba tiene una Ho “Las varianzas son iguales”
#Test Normalidad # Shapiro-Wilk Normality Test
SDEMT319_10 %>%
filter(clase2==1) %>%
with(ks.test(ing_x_hrs, y='pnorm', alternative = "two.sided"))## Warning in ks.test(ing_x_hrs, y = "pnorm", alternative = "two.sided"): ties
## should not be present for the Kolmogorov-Smirnov test##
## One-sample Kolmogorov-Smirnov test
##
## data: ing_x_hrs
## D = 0.5, p-value < 2.2e-16
## alternative hypothesis: two-sidedLa prueba tiene una Ho “La variable es normal”
¿Qué hacer?
Kruskal-Wallis test
Hay una prueba muy parecida que se basa en el orden de las observaciones, y se lee muy parecida a la ANOVA
kruskal<-SDEMT319_10 %>%
filter(clase2==1) %>%
with(kruskal.test(ing_x_hrs ~ as_label(t_loc)))
kruskal##
## Kruskal-Wallis rank sum test
##
## data: ing_x_hrs by as_label(t_loc)
## Kruskal-Wallis chi-squared = 19.474, df = 2, p-value = 5.904e-05Para ver las comparaciones tenemos que usar el dunn.test(), del paquet DescTools
SDEMT319_10 %>%
filter(clase2==1) %>%
with(DunnTest(ing_x_hrs ~ as_label(t_loc)))##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff
## Localidades de 15 000 a 99 999 habitantes-Localidades mayores de 100 000 habitantes 216.78699
## Localidades de 2 500 a 14 999 habitantes-Localidades mayores de 100 000 habitantes 83.43737
## Localidades de 2 500 a 14 999 habitantes-Localidades de 15 000 a 99 999 habitantes -133.34962
## pval
## Localidades de 15 000 a 99 999 habitantes-Localidades mayores de 100 000 habitantes 4.8e-05
## Localidades de 2 500 a 14 999 habitantes-Localidades mayores de 100 000 habitantes 0.2724
## Localidades de 2 500 a 14 999 habitantes-Localidades de 15 000 a 99 999 habitantes 0.2724
##
## Localidades de 15 000 a 99 999 habitantes-Localidades mayores de 100 000 habitantes ***
## Localidades de 2 500 a 14 999 habitantes-Localidades mayores de 100 000 habitantes
## Localidades de 2 500 a 14 999 habitantes-Localidades de 15 000 a 99 999 habitantes
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1