Previo

¡Recuerda poner el directorio!

setwd("/Users/anaescoto/Dropbox/2020/R_invierno")

Vamos a importar la base completa del cuestionario sociodemográfico de la Encuesta Nacional de Ocupación y Empleo, trimestre III de 2019

library(haven)
SDEMT319 <- read_dta("./datos/SDEMT319.dta")

Librerías

#install.packages("DescTools", dependencies = T) #vamos a instalarlas
library(tidyverse) #porque usaremos dplyr y por cualquier cosa
## ── Attaching packages ──────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.1     ✔ purrr   0.3.3
## ✔ tibble  2.1.3     ✔ dplyr   0.8.3
## ✔ tidyr   1.0.0     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0
## ── Conflicts ─────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(sjlabelled)
## 
## Attaching package: 'sjlabelled'
## The following object is masked from 'package:forcats':
## 
##     as_factor
## The following object is masked from 'package:dplyr':
## 
##     as_label
## The following objects are masked from 'package:haven':
## 
##     as_factor, read_sas, read_spss, read_stata, write_sas,
##     zap_labels

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

t.test(SDEMT319$ing_x_hrs)
## 
##  One Sample t-test
## 
## data:  SDEMT319$ing_x_hrs
## t = 235.62, df = 405448, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  12.35883 12.56616
## sample estimates:
## mean of x 
##   12.4625

Univariado para hipótesis específica

t.test(SDEMT319$ing_x_hrs, mu=200)
## 
##  One Sample t-test
## 
## data:  SDEMT319$ing_x_hrs
## t = -3545.7, df = 405448, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 200
## 95 percent confidence interval:
##  12.35883 12.56616
## sample estimates:
## mean of x 
##   12.4625
t.test(SDEMT319$ing_x_hrs, mu=200, alternative = "two.sided") #default y de dos colas
## 
##  One Sample t-test
## 
## data:  SDEMT319$ing_x_hrs
## t = -3545.7, df = 405448, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 200
## 95 percent confidence interval:
##  12.35883 12.56616
## sample estimates:
## mean of x 
##   12.4625
t.test(SDEMT319$ing_x_hrs, mu=200, alternative = "less") # cola izquierda
## 
##  One Sample t-test
## 
## data:  SDEMT319$ing_x_hrs
## t = -3545.7, df = 405448, p-value < 2.2e-16
## alternative hypothesis: true mean is less than 200
## 95 percent confidence interval:
##     -Inf 12.5495
## sample estimates:
## mean of x 
##   12.4625
t.test(SDEMT319$ing_x_hrs, mu=200, alternative = "greater") #cola derecha 
## 
##  One Sample t-test
## 
## data:  SDEMT319$ing_x_hrs
## t = -3545.7, df = 405448, p-value = 1
## alternative hypothesis: true mean is greater than 200
## 95 percent confidence interval:
##  12.3755     Inf
## sample estimates:
## mean of x 
##   12.4625

Estimaciones 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 %>% 
    filter(SDEMT319$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             28.6
## 2 Mujer              27.6
SDEMT319 %>% 
    filter(SDEMT319$clase2==1) %>%
      with(t.test(ing_x_hrs~sex))
## 
##  Welch Two Sample t-test
## 
## data:  ing_x_hrs by sex
## t = 4.3816, df = 159468, p-value = 1.179e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.5329931 1.3957600
## sample estimates:
## mean in group 1 mean in group 2 
##        28.57272        27.60834

Estimación de varianzas y sus pruebas de hipótesis

Para poder hacer inferencia sobre la varianza utilizamos el comando varTest() del paquete “DescTools”

library(DescTools)

SDEMT319 %>% 
    filter(clase2==1) %>% 
      with(VarTest(ing_x_hrs))
## 
##  One Sample Chi-Square test on variance
## 
## data:  ing_x_hrs
## X-squared = 380450231, df = 179291, p-value < 2.2e-16
## alternative hypothesis: true variance is not equal to 1
## 95 percent confidence interval:
##  2108.148 2135.930
## sample estimates:
## variance of x 
##      2121.971

Podemos también decir algo sobre el valor objetivo de nuestra hipótesis

SDEMT319 %>% 
    filter(clase2==1) %>% 
      with(VarTest(ing_x_hrs, sigma.squared = 100000))
## 
##  One Sample Chi-Square test on variance
## 
## data:  ing_x_hrs
## X-squared = 3804.5, df = 179291, p-value < 2.2e-16
## alternative hypothesis: true variance is not equal to 1e+05
## 95 percent confidence interval:
##  2108.148 2135.930
## sample estimates:
## variance of x 
##      2121.971

Guardar como objeto nuestros resultados, siempres muy conveniente para pedir después o para realizar operaciones con ellos

test2<-SDEMT319 %>% 
    filter(clase2==1) %>% 
      with(VarTest(ing_x_hrs))
test2$conf.int
## [1] 2108.148 2135.930
## attr(,"conf.level")
## [1] 0.95
sqrt(test2$conf.int) ## sacamos la raíz cuadrada para tener las
## [1] 45.91457 46.21612
## attr(,"conf.level")
## [1] 0.95
#desviaciones estándar y sea más fácil de interpretar

Estimació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 %>% 
    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 = 10.12, num df = 179291, denom df = 179291, p-value < 2.2e-16
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  10.02671 10.21409
## sample estimates:
## ratio of variances 
##           10.11996

“x=” declara al numerador “y=” declara al denominador

SDEMT319 %>% 
    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 = 10.12, num df = 179291, denom df = 179291, p-value < 2.2e-16
## alternative hypothesis: true ratio of variances is not equal to 1
## 98 percent confidence interval:
##  10.00937 10.23178
## sample estimates:
## ratio of variances 
##           10.11996

Si lo que queremos es comparar la varianza entre dos grupos, usamos el signo ~

SDEMT319 %>% 
    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 = 1.0762, num df = 106754, denom df = 72536, p-value < 2.2e-16
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  1.061941 1.090649
## sample estimates:
## ratio of variances 
##           1.076208

Prueba 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$clase1, SDEMT319$sex)
##    
##          1      2
##   0  38488  36856
##   1 110989  75534
##   2  42647  93193
chisq.test(SDEMT319$clase1, SDEMT319$sex)
## 
##  Pearson's Chi-squared test
## 
## data:  SDEMT319$clase1 and SDEMT319$sex
## X-squared = 25156, df = 2, p-value < 2.2e-16

Aná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 %>% 
           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 %>% 
    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)      3   3017986 1005995   477.9 <2e-16 ***
## Residuals       179288 377432245    2105                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comparació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  -4.430829
## Localidades de 2 500 a 14 999 habitantes-Localidades mayores de 100 000 habitantes   -7.207890
## Localidades menores de 2 500 habitantes-Localidades mayores de 100 000 habitantes   -11.105366
## Localidades de 2 500 a 14 999 habitantes-Localidades de 15 000 a 99 999 habitantes   -2.777061
## Localidades menores de 2 500 habitantes-Localidades de 15 000 a 99 999 habitantes    -6.674537
## Localidades menores de 2 500 habitantes-Localidades de 2 500 a 14 999 habitantes     -3.897476
##                                                                                            lwr
## Localidades de 15 000 a 99 999 habitantes-Localidades mayores de 100 000 habitantes  -5.273939
## Localidades de 2 500 a 14 999 habitantes-Localidades mayores de 100 000 habitantes   -8.081144
## Localidades menores de 2 500 habitantes-Localidades mayores de 100 000 habitantes   -11.940902
## Localidades de 2 500 a 14 999 habitantes-Localidades de 15 000 a 99 999 habitantes   -3.881213
## Localidades menores de 2 500 habitantes-Localidades de 15 000 a 99 999 habitantes    -7.749107
## Localidades menores de 2 500 habitantes-Localidades de 2 500 a 14 999 habitantes     -4.995855
##                                                                                            upr
## Localidades de 15 000 a 99 999 habitantes-Localidades mayores de 100 000 habitantes  -3.587718
## Localidades de 2 500 a 14 999 habitantes-Localidades mayores de 100 000 habitantes   -6.334636
## Localidades menores de 2 500 habitantes-Localidades mayores de 100 000 habitantes   -10.269829
## Localidades de 2 500 a 14 999 habitantes-Localidades de 15 000 a 99 999 habitantes   -1.672910
## Localidades menores de 2 500 habitantes-Localidades de 15 000 a 99 999 habitantes    -5.599967
## Localidades menores de 2 500 habitantes-Localidades de 2 500 a 14 999 habitantes     -2.799096
##                                                                                     p adj
## Localidades de 15 000 a 99 999 habitantes-Localidades mayores de 100 000 habitantes     0
## Localidades de 2 500 a 14 999 habitantes-Localidades mayores de 100 000 habitantes      0
## Localidades menores de 2 500 habitantes-Localidades mayores de 100 000 habitantes       0
## Localidades de 2 500 a 14 999 habitantes-Localidades de 15 000 a 99 999 habitantes      0
## Localidades menores de 2 500 habitantes-Localidades de 15 000 a 99 999 habitantes       0
## Localidades menores de 2 500 habitantes-Localidades de 2 500 a 14 999 habitantes        0

Supuestos 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 %>% 
    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 = 12718, df = 3, p-value < 2.2e-16

La prueba tiene una Ho “Las varianzas son iguales”

#Test Normalidad # Shapiro-Wilk Normality Test
SDEMT319 %>% 
    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.6745, p-value < 2.2e-16
## alternative hypothesis: two-sided

La 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 %>% 
    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 = 2003.5, df = 3, p-value < 2.2e-16

Para ver las comparaciones tenemos que usar el dunn.test(), del paquet DescTools

SDEMT319 %>% 
    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      -4891.943
## Localidades de 2 500 a 14 999 habitantes-Localidades mayores de 100 000 habitantes       -8266.128
## Localidades menores de 2 500 habitantes-Localidades mayores de 100 000 habitantes       -15175.431
## Localidades de 2 500 a 14 999 habitantes-Localidades de 15 000 a 99 999 habitantes       -3374.184
## Localidades menores de 2 500 habitantes-Localidades de 15 000 a 99 999 habitantes       -10283.487
## Localidades menores de 2 500 habitantes-Localidades de 2 500 a 14 999 habitantes         -6909.303
##                                                                                        pval
## Localidades de 15 000 a 99 999 habitantes-Localidades mayores de 100 000 habitantes < 2e-16
## Localidades de 2 500 a 14 999 habitantes-Localidades mayores de 100 000 habitantes  < 2e-16
## Localidades menores de 2 500 habitantes-Localidades mayores de 100 000 habitantes   < 2e-16
## Localidades de 2 500 a 14 999 habitantes-Localidades de 15 000 a 99 999 habitantes  1.5e-12
## Localidades menores de 2 500 habitantes-Localidades de 15 000 a 99 999 habitantes   < 2e-16
## Localidades menores de 2 500 habitantes-Localidades de 2 500 a 14 999 habitantes    < 2e-16
##                                                                                        
## 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 menores de 2 500 habitantes-Localidades mayores de 100 000 habitantes   ***
## Localidades de 2 500 a 14 999 habitantes-Localidades de 15 000 a 99 999 habitantes  ***
## Localidades menores de 2 500 habitantes-Localidades de 15 000 a 99 999 habitantes   ***
## Localidades menores de 2 500 habitantes-Localidades de 2 500 a 14 999 habitantes    ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ejercicio

Este es un ejercicio más libre y para que revisen más la base de datos. Presente dos pruebas de hipótesis revisadas en esta práctica con otras variables de su elección. Revise la ayuda de las pruebas para mejor interpretación de sus resultados.

Envíelo a la siguiente liga. https://www.dropbox.com/request/uParoat6jNX25S5UyA2c