Práctica 7. Inferencia

Previo

¡Recuerda, poner el directorio!

setwd("~/Dropbox/FCPyS-2020-i/EAIII/Prácticas_R")

A traer la base

library(haven)
SDEMT219 <- read_dta("SDEMT219.dta")

Librerías básicas

#install.packages("DescTools",  dependencies = T)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(sjlabelled)

## 
## Attaching package: 'sjlabelled'

## 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

Recuerda que para muchas de las variables de los ocupados tenemos que filtrar. Esto lo podemos hacer con dplyr y pipes y luego aplicar el t.test

SDEMT219 %>% 
  filter(clase2==1) %>% 
   with(t.test(ing_x_hrs))

## 
##  One Sample t-test
## 
## data:  ing_x_hrs
## t = 248.38, df = 179611, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  28.45821 28.91090
## sample estimates:
## mean of x 
##  28.68456

Para que todo sea más sencillo, vamos a hacer un subset de solo ocupados

ocupados<-SDEMT219 %>% 
  filter(clase2==1)

Univariado para hipótesis específica

t.test(ocupados$ing_x_hrs, mu=28)

## 
##  One Sample t-test
## 
## data:  ocupados$ing_x_hrs
## t = 5.9277, df = 179611, p-value = 3.078e-09
## alternative hypothesis: true mean is not equal to 28
## 95 percent confidence interval:
##  28.45821 28.91090
## sample estimates:
## mean of x 
##  28.68456

t.test(ocupados$ing_x_hrs, mu=28, alternative = "two.sided") #default y de dos colas

## 
##  One Sample t-test
## 
## data:  ocupados$ing_x_hrs
## t = 5.9277, df = 179611, p-value = 3.078e-09
## alternative hypothesis: true mean is not equal to 28
## 95 percent confidence interval:
##  28.45821 28.91090
## sample estimates:
## mean of x 
##  28.68456

t.test(ocupados$ing_x_hrs, mu=28, alternative = "less") # cola izquierda

## 
##  One Sample t-test
## 
## data:  ocupados$ing_x_hrs
## t = 5.9277, df = 179611, p-value = 1
## alternative hypothesis: true mean is less than 28
## 95 percent confidence interval:
##      -Inf 28.87451
## sample estimates:
## mean of x 
##  28.68456

t.test(ocupados$ing_x_hrs, mu=28, alternative = "greater") #cola derecha 

## 
##  One Sample t-test
## 
## data:  ocupados$ing_x_hrs
## t = 5.9277, df = 179611, p-value = 1.539e-09
## alternative hypothesis: true mean is greater than 28
## 95 percent confidence interval:
##  28.4946     Inf
## sample estimates:
## mean of x 
##  28.68456

Nivel de significancia

t.test(ocupados$ing_x_hrs, mu=28,
       alternative = "two.sided",
       conf.level=.99 ) #default y de dos colas + 99% de confianza

## 
##  One Sample t-test
## 
## data:  ocupados$ing_x_hrs
## t = 5.9277, df = 179611, p-value = 3.078e-09
## alternative hypothesis: true mean is not equal to 28
## 99 percent confidence interval:
##  28.38709 28.98203
## sample estimates:
## mean of x 
##  28.68456

prop.test

Si recordamos un poco prácticas pasadas, podemos guardar en un objeto la tabla de frecuencias. Podemos también obtener la estimación puntual de nuestra estimación y también podemos obtener de ahí nuestra prueba que incluye la inferencia.

freq.sex<-table(ocupados$sex)
prop.table(freq.sex)

## 
##         1         2 
## 0.5953834 0.4046166

freq.sex

## 
##      1      2 
## 106938  72674

prop.test(freq.sex)

## 
##  1-sample proportions test with continuity correction
## 
## data:  freq.sex, null probability 0.5
## X-squared = 6536.1, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.5931087 0.5976540
## sample estimates:
##         p 
## 0.5953834

Por 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(72674 , 179612)

## 
##  1-sample proportions test with continuity correction
## 
## data:  72674 out of 179612, null probability 0.5
## X-squared = 6536.1, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.4023460 0.4068913
## sample estimates:
##         p 
## 0.4046166

Para hacer la prueba con respecto a un nivel de proporción.

prop.test(freq.sex, p=0.5,alternative = "two.sided")

## 
##  1-sample proportions test with continuity correction
## 
## data:  freq.sex, null probability 0.5
## X-squared = 6536.1, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.5931087 0.5976540
## sample estimates:
##         p 
## 0.5953834

prop.test(freq.sex, p=0.5, alternative = "greater")

## 
##  1-sample proportions test with continuity correction
## 
## data:  freq.sex, null probability 0.5
## X-squared = 6536.1, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is greater than 0.5
## 95 percent confidence interval:
##  0.5934742 1.0000000
## sample estimates:
##         p 
## 0.5953834

prop.test(freq.sex, p=0.5, alternative = "less")

## 
##  1-sample proportions test with continuity correction
## 
## data:  freq.sex, null probability 0.5
## X-squared = 6536.1, df = 1, p-value = 1
## alternative hypothesis: true p is less than 0.5
## 95 percent confidence interval:
##  0.0000000 0.5972896
## sample estimates:
##         p 
## 0.5953834

La corrección ¿qué hace?

prop.test(72674,179612, alternative = "greater", p=0.5, correct = F)

## 
##  1-sample proportions test without continuity correction
## 
## data:  72674 out of 179612, null probability 0.5
## X-squared = 6536.4, df = 1, p-value = 1
## alternative hypothesis: true p is greater than 0.5
## 95 percent confidence interval:
##  0.4027131 1.0000000
## sample estimates:
##         p 
## 0.4046166

La función prop.test no realiza una prueba z, como casi todos los libros de estadística establece. ¡Hace una prueba de Chi cuadrado, basada en que hay una variable categórica con dos estados (éxito y fracaso)! Por ello vemos la línea que comienza con “X-cuadrado” La corrección de continuidad de Yates, que se ajusta a las diferencias que surgen al utilizar una aproximación normal a la distribución binomial, también se aplica automáticamente. Esto elimina 0.5 / n del límite inferior del intervalo de confianza y agrega 0.5 / n al límite superior. El intervalo de confianza dado por la prueba de propiedades no está la estimación de la muestra, p-hat. ¡Oh no! Pero, de nuevo, esto no es preocupante, ya que prop.test usa el intervalo de puntuación de Wilson para construir el intervalo de confianza. Esto da como resultado un intervalo de confianza asimétrico, pero presumiblemente más preciso (con respecto a la población real).

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?

ocupados %>% 
  group_by(sex) %>% 
  summarize(avg_ing=mean(ing_x_hrs))

## # A tibble: 2 x 2
##          sex avg_ing
##    <dbl+lbl>   <dbl>
## 1 1 [Hombre]    29.1
## 2 2 [Mujer]     28.1

t.test(ocupados$ing_x_hrs~ocupados$sex)

## 
##  Welch Two Sample t-test
## 
## data:  ocupados$ing_x_hrs by ocupados$sex
## t = 4.3423, df = 150862, p-value = 1.411e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.5658055 1.4967863
## sample estimates:
## mean in group 1 mean in group 2 
##        29.10184        28.07054

No tan bonito, las proporciones para dos poblaciones

Vamos a revisar si las mujeres tienen mayores participaciones en empleos informales que los hombres

table(as_label(ocupados$emp_ppal))

## 
##       No aplica Empleo informal   Empleo formal 
##               0           92701           86911

table(ocupados$emp_ppal)

## 
##     1     2 
## 92701 86911

Vamos a crear una variable que codifique para cuando se tienen un empleo informal(1) y formal(0)

ocupados<-ocupados %>% 
  mutate(dummy=ifelse(emp_ppal==2,T,0))

table(ocupados$dummy)

## 
##     0     1 
## 92701 86911

#Prueba de proporciones

ocupados %>% 
  group_by(sex, dummy) %>%
  tally()

## # A tibble: 4 x 3
## # Groups:   sex [2]
##          sex dummy     n
##    <dbl+lbl> <dbl> <int>
## 1 1 [Hombre]     0 54148
## 2 1 [Hombre]     1 52790
## 3 2 [Mujer]      0 38553
## 4 2 [Mujer]      1 34121

ocupados %>% 
  group_by(sex) %>%
  tally()

## # A tibble: 2 x 2
##          sex      n
##    <dbl+lbl>  <int>
## 1 1 [Hombre] 106938
## 2 2 [Mujer]   72674

ocupados %>% 
  group_by(sex) %>%
  summarize(prop=mean(dummy))

## # A tibble: 2 x 2
##          sex  prop
##    <dbl+lbl> <dbl>
## 1 1 [Hombre] 0.494
## 2 2 [Mujer]  0.470

Un vector de éxitos y uno de intentos

prop.test(c(52790,34121 ),c(106938, 72674))

## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  c(52790, 34121) out of c(106938, 72674)
## X-squared = 100.89, df = 1, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.01942549 0.02886023
## sample estimates:
##    prop 1    prop 2 
## 0.4936505 0.4695077

prop.test(c(52790,34121 ),c(106938, 72674), alternative= "two.sided")

## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  c(52790, 34121) out of c(106938, 72674)
## X-squared = 100.89, df = 1, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.01942549 0.02886023
## sample estimates:
##    prop 1    prop 2 
## 0.4936505 0.4695077

prop.test(c(52790,34121 ),c(106938, 72674), alternative= "greater")

## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  c(52790, 34121) out of c(106938, 72674)
## X-squared = 100.89, df = 1, p-value < 2.2e-16
## alternative hypothesis: greater
## 95 percent confidence interval:
##  0.02018206 1.00000000
## sample estimates:
##    prop 1    prop 2 
## 0.4936505 0.4695077

prop.test(c(52790,34121 ),c(106938, 72674), alternative= "less")

## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  c(52790, 34121) out of c(106938, 72674)
## X-squared = 100.89, df = 1, p-value = 1
## alternative hypothesis: less
## 95 percent confidence interval:
##  -1.00000000  0.02810366
## sample estimates:
##    prop 1    prop 2 
## 0.4936505 0.4695077

Correlación

Veamos la relación que existe entre los ingresos por hora y los años de escolaridad.

Ojo que el 99 nos va a dar problemas

ocupados$anios_esc[ocupados$anios_esc==99]<-NA

Estimación puntual de correlación

cor(ocupados$ing_x_hrs, ocupados$anios_esc)

## [1] NA

cor(ocupados$ing_x_hrs, ocupados$anios_esc, use="pairwise")

## [1] 0.1191839

Otros tipos de correlación

cor(ocupados$ing_x_hrs,
    ocupados$anios_esc, 
    use="pairwise", method= "spearman")

## [1] 0.0573907

Prueba de hipotésis

cor.test(ocupados$ing_x_hrs, ocupados$anios_esc)

## 
##  Pearson's product-moment correlation
## 
## data:  ocupados$ing_x_hrs and ocupados$anios_esc
## t = 50.834, df = 179335, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1146189 0.1237438
## sample estimates:
##       cor 
## 0.1191839

cor.test(ocupados$ing_x_hrs, ocupados$anios_esc, conf.level = .99)

## 
##  Pearson's product-moment correlation
## 
## data:  ocupados$ing_x_hrs and ocupados$anios_esc
## t = 50.834, df = 179335, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 99 percent confidence interval:
##  0.1131835 0.1251756
## sample estimates:
##       cor 
## 0.1191839

Estimación de varianzas

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

Simple

library(DescTools)
VarTest(ocupados$ing_x_hrs)

## 
##  One Sample Chi-Square test on variance
## 
## data:  ocupados$ing_x_hrs
## X-squared = 430245888, df = 179611, p-value < 2.2e-16
## alternative hypothesis: true variance is not equal to 1
## 95 percent confidence interval:
##  2379.842 2411.176
## sample estimates:
## variance of x 
##      2395.432

Con valor objetivo

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

VarTest(ocupados$ing_x_hrs, sigma.squared = 28*28)

## 
##  One Sample Chi-Square test on variance
## 
## data:  ocupados$ing_x_hrs
## X-squared = 548783, df = 179611, p-value < 2.2e-16
## alternative hypothesis: true variance is not equal to 784
## 95 percent confidence interval:
##  2379.842 2411.176
## sample estimates:
## variance of x 
##      2395.432

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

test2<-VarTest(ocupados$ing_x_hrs)
test2$conf.int

## [1] 2379.842 2411.176
## attr(,"conf.level")
## [1] 0.95

sqrt(test2$conf.int) ## sacamos la raíz cuadrada para tener las

## [1] 48.78362 49.10373
## attr(,"conf.level")
## [1] 0.95

#desviaciones estándar y sea más fácil de interpretar

Estimación de diferencias de varianzas

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.

var.test(ocupados$ing_x_hrs, ocupados$eda, ratio=1)

## 
##  F test to compare two variances
## 
## data:  ocupados$ing_x_hrs and ocupados$eda
## F = 11.569, num df = 179611, denom df = 179611, p-value < 2.2e-16
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  11.46251 11.67653
## sample estimates:
## ratio of variances 
##           11.56903

var.test(ocupados$ing_x_hrs, ocupados$eda, ratio=0.5, conf.level = 0.99)

## 
##  F test to compare two variances
## 
## data:  ocupados$ing_x_hrs and ocupados$eda
## F = 23.138, num df = 179611, denom df = 179611, p-value < 2.2e-16
## alternative hypothesis: true ratio of variances is not equal to 0.5
## 99 percent confidence interval:
##  11.42925 11.71052
## sample estimates:
## ratio of variances 
##           11.56903

¿Tendrá sentido comparar la varianza del ingreso y de la edad? Piénsalo. ## Tiene más sentido comparar grupos Si lo que queremos es comparar la varianza entre dos grupos, usamos el signo ~

var.test(ocupados$ing_x_hrs ~ocupados$sex, ratio=1, conf.level = 0.99)

## 
##  F test to compare two variances
## 
## data:  ocupados$ing_x_hrs by ocupados$sex
## F = 0.9062, num df = 106937, denom df = 72673, p-value < 2.2e-16
## alternative hypothesis: true ratio of variances is not equal to 1
## 99 percent confidence interval:
##  0.8904522 0.9221931
## sample estimates:
## ratio of variances 
##          0.9061952

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.

addmargins(
  table(as_label(ocupados$c_ocu11c), 
      as_label(ocupados$sex)))

##                                                  
##                                                   Hombre  Mujer    Sum
##   No aplica                                            0      0      0
##   Profesionales, técnicos y trabajadores del arte  11291   8279  19570
##   Trabajadores de la educación                      2732   4267   6999
##   Funcionarios y directivos                         2074   1176   3250
##   Oficinistas                                       7119   9656  16775
##   Trabajadores industriales artesanos y ayudantes  35463  13459  48922
##   Comerciantes                                     14417  17361  31778
##   Operadores de transporte                          8709    141   8850
##   Trabajadores en servicios personales             10788  16421  27209
##   Trabajadores en protección y vigilancia           1284    253   1537
##   Trabajadores agropecuarios                       12994   1630  14624
##   No especificado                                     67     31     98
##   Sum                                             106938  72674 179612

chisq.test(as_label(ocupados$c_ocu11c), 
      as_label(ocupados$sex))

## 
##  Pearson's Chi-squared test
## 
## data:  as_label(ocupados$c_ocu11c) and as_label(ocupados$sex)
## X-squared = 24971, df = 10, p-value < 2.2e-16