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

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
#library(DescTools) ## la instalamos en la guía pasada

ANOVA

Análisis de varianza. Haramos la versión más simple. Para ver el efecto de un factor sobre una variable cualitativa. Recordemos que los ingresos por hora son diferentes entre la residencia según tamaño de localidad.

Lo primero que haremos es checar la distribución de los ingresos para observar si se parece a una normal.

[Y recuerda que no tienes que instalar el paquete si ya lo instalaste]

#install.packages("ggplot2",  dependencies = TRUE)
library(ggplot2)

Hoy graficamos como lo hiciemos en prácticas pasadas.

lienzo_bi <-SDEMT219 %>% 
          filter(clase2==1) %>% 
          ggplot(aes(x=ing_x_hrs, fill=as_label(t_loc), 
           color=as_label(t_loc),
           alpha=I(0.5)))

lienzo_bi + geom_density()

¿Esto es normal? Si nos quedamos con los que tienen ingresos mayores que cero y transformamos en logaritmo…

lienzo_bi2 <-SDEMT219 %>% 
          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_bi2 + geom_density()

Para hacer más fácil nuestro análisis, vamos a quedarnos con un subset…

dt_anova<-SDEMT219 %>% 
          filter(clase2==1 & ing_x_hrs>0) %>% 
          mutate(log_ing=log(ing_x_hrs))

La prueba ANOVA o análisis de varianza, nos dice cuánto de nuestra variable se ve explicado por un factor

anova <- aov(log_ing ~ as.factor(t_loc), data=dt_anova)
summary(anova)
##                      Df Sum Sq Mean Sq F value Pr(>F)    
## as.factor(t_loc)      3   3612  1204.1    2433 <2e-16 ***
## Residuals        122781  60776     0.5                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Extra

¿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 = log_ing ~ as.factor(t_loc), data = dt_anova)
## 
## $`as.factor(t_loc)`
##           diff        lwr        upr p adj
## 2-1 -0.1721923 -0.1878343 -0.1565504     0
## 3-1 -0.2961286 -0.3121458 -0.2801115     0
## 4-1 -0.4673519 -0.4828074 -0.4518963     0
## 3-2 -0.1239363 -0.1442683 -0.1036043     0
## 4-2 -0.2951595 -0.3150521 -0.2752670     0
## 4-3 -0.1712232 -0.1914121 -0.1510343     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
bartlett.test(log_ing ~ as.factor(t_loc), data=dt_anova)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  log_ing by as.factor(t_loc)
## Bartlett's K-squared = 963.88, df = 3, p-value < 2.2e-16

La prueba tiene una Ho “Las varianzas son iguales”

#Test Normalidad # Shapiro-Wilk Normality Test
ks.test(dt_anova$log_ing, y='pnorm')
## Warning in ks.test(dt_anova$log_ing, y = "pnorm"): ties should not be
## present for the Kolmogorov-Smirnov test
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  dt_anova$log_ing
## D = 0.95347, p-value < 2.2e-16
## alternative hypothesis: two-sided

La prueba tiene una Ho “La variable es normal”

¿Qué hacer?

Hay una prueba muy parecida que se basa en el orden de las observaciones, y se lee muy parecida a la ANOVA

kruskal<- kruskal.test(log_ing ~ as_label(t_loc), data=dt_anova)
kruskal
## 
##  Kruskal-Wallis rank sum test
## 
## data:  log_ing by as_label(t_loc)
## Kruskal-Wallis chi-squared = 6573.3, df = 3, p-value < 2.2e-16

Para ver las comparaciones tenemos que usar el dunn.test()

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

DunnTest(log_ing ~ as_label(t_loc), data=dt_anova,
         method="bonferroni")
## 
##  Dunn's test of multiple comparisons using rank sums : bonferroni  
## 
##                                                                                     mean.rank.diff
## Localidades de 15 000 a 99 999 habitantes-Localidades mayores de 100 000 habitantes      -9411.742
## Localidades de 2 500 a 14 999 habitantes-Localidades mayores de 100 000 habitantes      -15018.278
## Localidades menores de 2 500 habitantes-Localidades mayores de 100 000 habitantes       -21719.728
## Localidades de 2 500 a 14 999 habitantes-Localidades de 15 000 a 99 999 habitantes       -5606.536
## Localidades menores de 2 500 habitantes-Localidades de 15 000 a 99 999 habitantes       -12307.986
## Localidades menores de 2 500 habitantes-Localidades de 2 500 a 14 999 habitantes         -6701.450
##                                                                                       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  <2e-16
## 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

Regresión lineal

Regresión lineal simple

La regresión lineal nos ayuda a describir esta relación a través de una línea recta. Trabajaremos con la clásica “ingresos y escolaridad”

Quitamos los 99 de los anios_esc y sacamos una correlación preliminar.

dt_reg<-dt_anova %>% 
        filter(anios_esc<99) 

cor(dt_reg$log_ing, dt_reg$anios_esc, use = "pairwise")
## [1] 0.3712612

El modelo se ajusta así:

modelo <-lm(log_ing ~anios_esc, data=dt_reg, na.action=na.exclude)

summary(modelo) # show results
## 
## Call:
## lm(formula = log_ing ~ anios_esc, data = dt_reg, na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1866 -0.3736 -0.0085  0.3709  5.8956 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.8038968  0.0050062   560.1   <2e-16 ***
## anios_esc   0.0641848  0.0004584   140.0   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6725 on 122632 degrees of freedom
## Multiple R-squared:  0.1378, Adjusted R-squared:  0.1378 
## F-statistic: 1.961e+04 on 1 and 122632 DF,  p-value: < 2.2e-16
plot(modelo)

##Diagnósticos

Instalar paquete “car”

#install.packages("car", dependencies = TRUE)
library(carData)
library(car)
## 
## Attaching package: 'car'
## The following object is masked from 'package:DescTools':
## 
##     Recode
## The following object is masked from 'package:dplyr':
## 
##     recode
  1. Outliers
# Assessing Outliers
outlierTest(modelo) # Bonferonni p-value for most extreme obs
##         rstudent unadjusted p-value Bonferroni p
## 122082  8.769555         1.8161e-18   2.2271e-13
## 98353  -7.714226         1.2263e-14   1.5039e-09
## 101155 -7.301888         2.8547e-13   3.5008e-08
## 101466 -7.218372         5.2917e-13   6.4894e-08
## 51772   7.016179         2.2919e-12   2.8107e-07
## 68210   6.832785         8.3668e-12   1.0261e-06
## 56893   6.434264         1.2453e-10   1.5271e-05
## 92568  -6.383935         1.7321e-10   2.1241e-05
## 101671 -6.383935         1.7321e-10   2.1241e-05
## 101674 -6.326837         2.5110e-10   3.0793e-05
qqPlot(modelo, main="QQ Plot") #qq plot for studentized resid  (también sirve para normalidad)

## [1]  98353 122082
  1. Homocedasticidad
# non-constant error variance test
ncvTest(modelo)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 166.2335, Df = 1, p = < 2.22e-16
# plot studentized residuals vs. fitted values 
spreadLevelPlot(modelo)

## 
## Suggested power transformation:  0.6716695