Previo

Librerías

Hipótesis e intervalos de confianza

t-test

prop.test

Estimaciones bivariadas

Diferencias de medias por grupos
No tan bonito, las proporciones para dos poblaciones

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

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

Prueba chi-cuadrado chi-sq. Una aplicación más común

Análisis de varianza

Primero un gráfico

Comparación entre grupos

Supuestos de ANOVA
Kruskal-Wallis test

Ejercicio

Previo

Siempre: revisar el directorio

setwd("/Users/anaescoto/Dropbox/2019/Cursos ESA/UCA_R")
load("./datos/ehpm2017.RData")

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.0     ✔ purrr   0.3.2
## ✔ tibble  2.1.3     ✔ dplyr   0.8.1
## ✔ tidyr   0.8.3     ✔ 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

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(ehpm2017$money)

## 
##  One Sample t-test
## 
## data:  ehpm2017$money
## t = 116.8, df = 75132, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##   99.64088 103.04217
## sample estimates:
## mean of x 
##  101.3415

Univariado para hipótesis específica

t.test(ehpm2017$money, mu=200)

## 
##  One Sample t-test
## 
## data:  ehpm2017$money
## t = -113.7, df = 75132, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 200
## 95 percent confidence interval:
##   99.64088 103.04217
## sample estimates:
## mean of x 
##  101.3415

t.test(ehpm2017$money, mu=200, alternative = "two.sided") #default y de dos colas

## 
##  One Sample t-test
## 
## data:  ehpm2017$money
## t = -113.7, df = 75132, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 200
## 95 percent confidence interval:
##   99.64088 103.04217
## sample estimates:
## mean of x 
##  101.3415

t.test(ehpm2017$money, mu=200, alternative = "less") # cola izquierda

## 
##  One Sample t-test
## 
## data:  ehpm2017$money
## t = -113.7, df = 75132, p-value < 2.2e-16
## alternative hypothesis: true mean is less than 200
## 95 percent confidence interval:
##      -Inf 102.7687
## sample estimates:
## mean of x 
##  101.3415

t.test(ehpm2017$money, mu=200, alternative = "greater") #cola derecha

## 
##  One Sample t-test
## 
## data:  ehpm2017$money
## t = -113.7, df = 75132, p-value = 1
## alternative hypothesis: true mean is greater than 200
## 95 percent confidence interval:
##  99.9143     Inf
## sample estimates:
## mean of x 
##  101.3415

prop.test

Si recordamos un poco de la práctica 3, 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.asiste<-table(ehpm2017[ehpm2017$r106>4 & ehpm2017$r106<24,]$r203)
prop.table(freq.asiste)

## 
##         1         2 
## 0.6523062 0.3476938

freq.asiste

## 
##     1     2 
## 18003  9596

prop.test(freq.asiste)

## 
##  1-sample proportions test with continuity correction
## 
## data:  freq.asiste, null probability 0.5
## X-squared = 2560.3, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.6466486 0.6579213
## sample estimates:
##         p 
## 0.6523062

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(18003,(18003+9596))

## 
##  1-sample proportions test with continuity correction
## 
## data:  18003 out of (18003 + 9596), null probability 0.5
## X-squared = 2560.3, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.6466486 0.6579213
## sample estimates:
##         p 
## 0.6523062

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

prop.test(freq.asiste, p=0.75,alternative = "two.sided")

## 
##  1-sample proportions test with continuity correction
## 
## data:  freq.asiste, null probability 0.75
## X-squared = 1404.3, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.75
## 95 percent confidence interval:
##  0.6466486 0.6579213
## sample estimates:
##         p 
## 0.6523062

prop.test(freq.asiste, p=0.75, alternative = "greater")

## 
##  1-sample proportions test with continuity correction
## 
## data:  freq.asiste, null probability 0.75
## X-squared = 1404.3, df = 1, p-value = 1
## alternative hypothesis: true p is greater than 0.75
## 95 percent confidence interval:
##  0.6475581 1.0000000
## sample estimates:
##         p 
## 0.6523062

prop.test(freq.asiste, p=0.75, alternative = "less")

## 
##  1-sample proportions test with continuity correction
## 
## data:  freq.asiste, null probability 0.75
## X-squared = 1404.3, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is less than 0.75
## 95 percent confidence interval:
##  0.0000000 0.6570244
## sample estimates:
##         p 
## 0.6523062

La corrección ¿qué hace?

prop.test(18003,(18003+9596), alternative = "greater", p=0.5, correct = F)

## 
##  1-sample proportions test without continuity correction
## 
## data:  18003 out of (18003 + 9596), null probability 0.5
## X-squared = 2560.9, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is greater than 0.5
## 95 percent confidence interval:
##  0.6475763 1.0000000
## sample estimates:
##         p 
## 0.6523062

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?

ehpm2017 %>% 
    filter(ehpm2017$actpr2012==10) %>%
      group_by(as_label(r104)) %>%
      summarise(avg_ing = mean(money))

## # A tibble: 2 x 2
##   `as_label(r104)` avg_ing
##   <fct>              <dbl>
## 1 Hombre              243.
## 2 Mujer               234.

ehpm2017 %>% 
    filter(ehpm2017$actpr2012==10) %>%
      with(t.test(money~r104))

## 
##  Welch Two Sample t-test
## 
## data:  money by r104
## t = 2.7066, df = 31497, p-value = 0.006801
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   2.548946 15.933114
## sample estimates:
## mean in group 1 mean in group 2 
##        242.9812        233.7402

No tan bonito, las proporciones para dos poblaciones

Vamos a ver la proporción

#Prueba de proporciones
x<-ehpm2017 %>% 
    filter(r106>4 & r106<24) %>%
      group_by(as_label(r104), as_label(r203)) %>%
      tally()
x

## # A tibble: 4 x 3
## # Groups:   as_label(r104) [2]
##   `as_label(r104)` `as_label(r203)`     n
##   <fct>            <fct>            <int>
## 1 Hombre           Sí                9136
## 2 Hombre           No                4670
## 3 Mujer            Sí                8867
## 4 Mujer            No                4926

x$n

## [1] 9136 4670 8867 4926

x1<-9136
n1<-9136+4670
x2<-8867
n2<-8867+4926
x<-c(x1,x2)
n<-c(n1,n2)

prop.test(x,n)

## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  x out of n
## X-squared = 10.76, df = 1, p-value = 0.001037
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.007571487 0.030186413
## sample estimates:
##    prop 1    prop 2 
## 0.6617413 0.6428623

Y con ello también podemos ver el sentido de las pruebas de hipótesis

prop.test(x,n, alternative= "two.sided")

## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  x out of n
## X-squared = 10.76, df = 1, p-value = 0.001037
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.007571487 0.030186413
## sample estimates:
##    prop 1    prop 2 
## 0.6617413 0.6428623

prop.test(x,n, alternative= "greater")

## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  x out of n
## X-squared = 10.76, df = 1, p-value = 0.0005187
## alternative hypothesis: greater
## 95 percent confidence interval:
##  0.009377778 1.000000000
## sample estimates:
##    prop 1    prop 2 
## 0.6617413 0.6428623

prop.test(x,n, alternative= "less")

## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  x out of n
## X-squared = 10.76, df = 1, p-value = 0.9995
## alternative hypothesis: less
## 95 percent confidence interval:
##  -1.00000000  0.02838012
## sample estimates:
##    prop 1    prop 2 
## 0.6617413 0.6428623

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)

ehpm2017 %>% 
    filter(actpr2012==10) %>% 
      with(VarTest(money))

## 
##  One Sample Chi-Square test on variance
## 
## data:  money
## X-squared = 3199401940, df = 31817, p-value < 2.2e-16
## alternative hypothesis: true variance is not equal to 1
## 95 percent confidence interval:
##   99011.9 102137.4
## sample estimates:
## variance of x 
##      100556.4

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

ehpm2017 %>% 
    filter(actpr2012==10) %>% 
      with(VarTest(money, sigma.squared = 100000))

## 
##  One Sample Chi-Square test on variance
## 
## data:  money
## X-squared = 31994, df = 31817, p-value = 0.482
## alternative hypothesis: true variance is not equal to 1e+05
## 95 percent confidence interval:
##   99011.9 102137.4
## sample estimates:
## variance of x 
##      100556.4

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

test2<-ehpm2017 %>% 
    filter(actpr2012==10) %>% 
      with(VarTest(money))
test2$conf.int

## [1]  99011.9 102137.4
## attr(,"conf.level")
## [1] 0.95

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

## [1] 314.6616 319.5895
## 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.

ehpm2017 %>% 
    filter(actpr2012==10) %>% 
      with(var.test(x=money, y=r106, ratio=1))

## 
##  F test to compare two variances
## 
## data:  money and r106
## F = 423.11, num df = 31817, denom df = 31817, p-value < 2.2e-16
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  413.9178 432.5164
## sample estimates:
## ratio of variances 
##           423.1149

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

ehpm2017 %>% 
    filter(actpr2012==10) %>% 
      with(var.test(x=money, y=r106, ratio=1, conf.level = 0.98))

## 
##  F test to compare two variances
## 
## data:  money and r106
## F = 423.11, num df = 31817, denom df = 31817, p-value < 2.2e-16
## alternative hypothesis: true ratio of variances is not equal to 1
## 98 percent confidence interval:
##  412.2208 434.2969
## sample estimates:
## ratio of variances 
##           423.1149

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

ehpm2017 %>% 
    filter(actpr2012==10) %>% 
      with(var.test(money ~ as_label(r104), ratio=1))

## 
##  F test to compare two variances
## 
## data:  money by as_label(r104)
## F = 1.8662, num df = 19147, denom df = 12669, p-value < 2.2e-16
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  1.807835 1.926338
## sample estimates:
## ratio of variances 
##           1.866244

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(ehpm2017$r203, ehpm2017$r104)

##    
##         1     2
##   1  9679  9439
##   2 23316 27745

chisq.test(ehpm2017$r203, ehpm2017$r104)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  ehpm2017$r203 and ehpm2017$r104
## X-squared = 137.43, df = 1, 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 <-ehpm2017 %>% 
           filter(actpr2012==10 & money<30000 & !money==0) %>% 
           ggplot(aes(x=log(money), fill=as_label(region), 
           color=as_label(region),
           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<-ehpm2017 %>% 
    filter(actpr2012==10) %>% 
      with(aov(money ~ as_label(region)))

summary(anova)

##                     Df    Sum Sq  Mean Sq F value Pr(>F)    
## as_label(region)     4 7.402e+07 18505167   188.4 <2e-16 ***
## Residuals        31813 3.125e+09    98242                   
## ---
## 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 = money ~ as_label(region))
## 
## $`as_label(region)`
##                                                     diff        lwr
## Central I-Occidental                           19.342574   5.338158
## Central II-Occidental                          -2.649243 -18.064026
## Oriental-Occidental                            -4.759991 -19.003083
## Área metropolitana de San Salvador-Occidental 136.792959 120.901667
## Central II-Central I                          -21.991817 -37.424084
## Oriental-Central I                            -24.102565 -38.364579
## Área metropolitana de San Salvador-Central I  117.450385 101.542132
## Oriental-Central II                            -2.110748 -17.759929
## Área metropolitana de San Salvador-Central II 139.442202 122.279325
## Área metropolitana de San Salvador-Oriental   141.552950 125.434187
##                                                      upr     p adj
## Central I-Occidental                           33.346991 0.0015483
## Central II-Occidental                          12.765539 0.9901151
## Oriental-Occidental                             9.483102 0.8925685
## Área metropolitana de San Salvador-Occidental 152.684251 0.0000000
## Central II-Central I                           -6.559551 0.0009619
## Oriental-Central I                             -9.840552 0.0000396
## Área metropolitana de San Salvador-Central I  133.358638 0.0000000
## Oriental-Central II                            13.538433 0.9961015
## Área metropolitana de San Salvador-Central II 156.605079 0.0000000
## Área metropolitana de San Salvador-Oriental   157.671713 0.0000000

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

ehpm2017 %>% 
    filter(actpr2012==10) %>% 
      with(bartlett.test(money ~ as_label(region)))

## 
##  Bartlett test of homogeneity of variances
## 
## data:  money by as_label(region)
## Bartlett's K-squared = 4944.5, df = 4, p-value < 2.2e-16

La prueba tiene una Ho “Las varianzas son iguales”

#Test Normalidad # Shapiro-Wilk Normality Test
ehpm2017 %>% 
    filter(actpr2012==10) %>% 
      with(ks.test(money, y='pnorm', alternative = "two.sided"))

## Warning in ks.test(money, y = "pnorm", alternative = "two.sided"): ties
## should not be present for the Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  money
## D = 0.82098, 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<-ehpm2017 %>% 
    filter(actpr2012==10) %>% 
      with(kruskal.test(money ~ as_label(region)))

kruskal

## 
##  Kruskal-Wallis rank sum test
## 
## data:  money by as_label(region)
## Kruskal-Wallis chi-squared = 1437.6, df = 4, p-value < 2.2e-16

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

ehpm2017 %>% 
    filter(actpr2012==10) %>% 
      with(DunnTest(money ~ as_label(region)))

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                                               mean.rank.diff    pval    
## Central I-Occidental                                858.1691 3.2e-08 ***
## Central II-Occidental                              -517.2855  0.0035 ** 
## Oriental-Occidental                                -953.8508 1.6e-09 ***
## Área metropolitana de San Salvador-Occidental      5064.4950 < 2e-16 ***
## Central II-Central I                              -1375.4547 4.3e-16 ***
## Oriental-Central I                                -1812.0199 < 2e-16 ***
## Área metropolitana de San Salvador-Central I       4206.3258 < 2e-16 ***
## Oriental-Central II                                -436.5652  0.0092 ** 
## Área metropolitana de San Salvador-Central II      5581.7805 < 2e-16 ***
## Área metropolitana de San Salvador-Oriental        6018.3457 < 2e-16 ***
## ---
## 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://tinyurl.com/Ej7-ESA-R

Práctica 7. Introducción a la inferencia

AE

09/07/2019