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")

Vamos a volver a usar nuestros paquetes ggplot2, dplyr y sjlabelled.

library(ggplot2)
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

Gráficos de comparación cuanti y cuali

Los gráficos de comparaciones de cajas y distribuciones son muy útiles. Nos permiten ver cómo una variable cuantitativa se comporta entre la población dividida entre los grupos.

Boxplot

Empecemos con el “boxplot” univariado y con la función que viene en “graphics”, sin necesidad de instalar paquetería

boxplot(SDEMT219[SDEMT219$clase2==1,]$ing_x_hrs)

La función boxplot, también sirve para obtener el resumen de cinco números.

a<-boxplot(SDEMT219[SDEMT219$clase2==1,]$ing_x_hrs)

a$stats
##          [,1]
## [1,]  0.00000
## [2,]  0.00000
## [3,] 21.73913
## [4,] 37.50000
## [5,] 93.75000

Con formato dplyr

SDEMT219 %>% 
  filter(clase2==1) %>% 
    with(boxplot(ing_x_hrs))

Boxplot bivariado

Podemos comparar la distribución de una variable cuantitativa entre las categorías de una variable categórica

#otra forma
boxplot(ing_x_hrs ~ as_label(sex), data= SDEMT219[SDEMT219$clase2==1,])

Con dplyr:

SDEMT219 %>% 
  filter(clase2==1) %>% 
    with(boxplot(ing_x_hrs ~ as_label(sex)))

con ggplot

Para “producción”, quizás queremos un gráfico más acabado en ggplot. Como siempre, primero haremos nuestro lienzo, aprovecharemos para quitar el valor atípico para mejor visualización

lienzo <-SDEMT219 %>% 
          filter(clase2==1 & ing_x_hrs<4000) %>% 
          ggplot(aes(x=as_label(sex), y=ing_x_hrs))
lienzo

Hoy podemos agregar nuestra geometría

lienzo + geom_boxplot()

g<-lienzo +  geom_boxplot(fill='#A4A4A4', color="black")+
  theme_classic()
g

Para ponerlo de modeo

g  + coord_flip()

Una tercera variable

g + facet_wrap(~ as_label(t_loc))

g2<-lienzo +  geom_boxplot(aes(color=as_label(t_loc))) +
  theme_classic()
g2

Gráfico de densidad bivariado

Para poder editar la transparencia de nuestro relleno podemos usar la opción alpha y ponerle un valor menor a 1.

lienzo_bi <-SDEMT219 %>% 
          filter(clase2==1 & ing_x_hrs<4000) %>% 
          ggplot(aes(x=log(ing_x_hrs+1), fill=as_label(sex), 
           color=as_label(sex),
           alpha=I(0.5)))

lienzo_bi + geom_density()

Vamos -sin usar el ggTheme Assistant a ponerle otros elementos a nuestro gráfico.

lienzo_bi + geom_density() +
  ggtitle("Logaritmo de los ingresos por trabajo") +
    xlab("Log(Dólares mensuales)") + ylab("Densidad") +
    labs(fill = "sex", color="sex")

Análisis descriptivo

Un truco para hacer el análisis descriptivo es generar una sub-base con las variables que vamos analizar

mydata<- SDEMT219[SDEMT219$clase2==1, c("ing_x_hrs", "t_loc","eda", "sex", "anios_esc", "fac")]
tail(mydata)
## # A tibble: 6 x 6
##   ing_x_hrs                           t_loc   eda       sex anios_esc   fac
##       <dbl>                       <dbl+lbl> <dbl> <dbl+lbl>     <dbl> <dbl>
## 1      0    4 [Localidades menores de 2 50…    35 1 [Hombr…         9   239
## 2      0    4 [Localidades menores de 2 50…    14 1 [Hombr…         8   239
## 3      0    4 [Localidades menores de 2 50…    44 2 [Mujer]         6   186
## 4      8.31 4 [Localidades menores de 2 50…    69 2 [Mujer]         9   186
## 5     22.5  4 [Localidades menores de 2 50…    44 1 [Hombr…         9   259
## 6     19.0  4 [Localidades menores de 2 50…    35 2 [Mujer]         9   259

Si queremos obtener las medias de estas tres variables, usamos la función sapply. Esto nos permite repetir la función para todo el data frame Excluyendo missing values

sapply(as_label(mydata), mean, na.rm=TRUE) # checa los NA's
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## ing_x_hrs     t_loc       eda       sex anios_esc       fac 
##  28.68456        NA  39.46372        NA  10.60235 308.30198

Además de media (mean) tenemos las siguientes opciones: mean, sd, var, min, max, median, range, and quantile.

Para obtener de una sola vez las funciones mean,median,25th and 75th quartiles,min,max, usamos la función “summary”

summary(as_label(mydata)) #total
##    ing_x_hrs                                             t_loc       
##  Min.   :   0.00   Localidades mayores de 100 000 habitantes:110024  
##  1st Qu.:   0.00   Localidades de 15 000 a 99 999 habitantes: 23661  
##  Median :  21.74   Localidades de 2 500 a 14 999 habitantes : 21726  
##  Mean   :  28.68   Localidades menores de 2 500 habitantes  : 24201  
##  3rd Qu.:  37.50                                                     
##  Max.   :6000.00                                                     
##       eda            sex           anios_esc         fac         
##  Min.   :12.00   Hombre:106938   Min.   : 0.0   Min.   :   12.0  
##  1st Qu.:28.00   Mujer : 72674   1st Qu.: 9.0   1st Qu.:  108.0  
##  Median :39.00                   Median :10.0   Median :  164.0  
##  Mean   :39.46                   Mean   :10.6   Mean   :  308.3  
##  3rd Qu.:50.00                   3rd Qu.:13.0   3rd Qu.:  396.0  
##  Max.   :98.00                   Max.   :99.0   Max.   :22449.0
summary(as_label(mydata[mydata$sex==1,])) #hombres
##    ing_x_hrs                                            t_loc      
##  Min.   :   0.0   Localidades mayores de 100 000 habitantes:63195  
##  1st Qu.:   0.0   Localidades de 15 000 a 99 999 habitantes:13961  
##  Median :  22.5   Localidades de 2 500 a 14 999 habitantes :13290  
##  Mean   :  29.1   Localidades menores de 2 500 habitantes  :16492  
##  3rd Qu.:  37.5                                                    
##  Max.   :3100.8                                                    
##       eda            sex           anios_esc          fac         
##  Min.   :12.00   Hombre:106938   Min.   : 0.00   Min.   :   12.0  
##  1st Qu.:27.00   Mujer :     0   1st Qu.: 8.00   1st Qu.:  111.0  
##  Median :38.00                   Median : 9.00   Median :  171.0  
##  Mean   :39.48                   Mean   :10.38   Mean   :  315.8  
##  3rd Qu.:50.00                   3rd Qu.:12.00   3rd Qu.:  406.0  
##  Max.   :98.00                   Max.   :99.00   Max.   :13564.0
summary(as_label(mydata[mydata$sex==2,])) #mujeres
##    ing_x_hrs                                             t_loc      
##  Min.   :   0.00   Localidades mayores de 100 000 habitantes:46829  
##  1st Qu.:   0.00   Localidades de 15 000 a 99 999 habitantes: 9700  
##  Median :  20.83   Localidades de 2 500 a 14 999 habitantes : 8436  
##  Mean   :  28.07   Localidades menores de 2 500 habitantes  : 7709  
##  3rd Qu.:  36.05                                                    
##  Max.   :6000.00                                                    
##       eda            sex          anios_esc          fac         
##  Min.   :12.00   Hombre:    0   Min.   : 0.00   Min.   :   12.0  
##  1st Qu.:28.00   Mujer :72674   1st Qu.: 9.00   1st Qu.:  105.0  
##  Median :39.00                  Median :11.00   Median :  155.0  
##  Mean   :39.44                  Mean   :10.93   Mean   :  297.3  
##  3rd Qu.:49.00                  3rd Qu.:15.00   3rd Qu.:  375.0  
##  Max.   :98.00                  Max.   :99.00   Max.   :22449.0

Tapply, se parece a sapply, nos permite hacer una función en un subconjunto que se deriva de un vector. Por ejemplo, a continuación podemos pedirle al programa que nos haga un summary del ingreso, tomando en cuenta un vector de factor

tapply(mydata$ing_x_hrs, as_label(mydata$sex), summary)
## $Hombre
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     0.0    22.5    29.1    37.5  3100.8 
## 
## $Mujer
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.00   20.83   28.07   36.05 6000.00
tapply(mydata$ing_x_hrs, as_label(mydata$sex), mean, na.rm=T)
##   Hombre    Mujer 
## 29.10184 28.07054

Con dplyr

Con dplyr a veces es más sencillo por la función “as group” y summarise (o summarize)

SDEMT219 %>% 
    filter(SDEMT219$clase2==1 )%>%
    summarise(avg_ing= mean(ing_x_hrs))
## # A tibble: 1 x 1
##   avg_ing
##     <dbl>
## 1    28.7

Con la función de grupos es más sencillo que tapply y sapply

SDEMT219 %>% 
    filter(SDEMT219$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             29.1
## 2 Mujer              28.1
SDEMT219 %>%
    filter(SDEMT219$clase2==1 )%>%
      group_by(as_label(sex)) %>%
      tally() # Nos da el conteo de casos
## # A tibble: 2 x 2
##   `as_label(sex)`      n
##   <fct>            <int>
## 1 Hombre          106938
## 2 Mujer            72674

Y podemos “summarizar” los indicadores que querramos. Recuerda que estos resultados se pueden guardar en un objeto.

resultados1<-SDEMT219 %>% 
    filter(SDEMT219$clase2==1 )%>%
      group_by(as_label(sex)) %>%
      summarise(avg_ing = mean(ing_x_hrs),
                avg_esc = mean(anios_esc)) 
resultados1
## # A tibble: 2 x 3
##   `as_label(sex)` avg_ing avg_esc
##   <fct>             <dbl>   <dbl>
## 1 Hombre             29.1    10.4
## 2 Mujer              28.1    10.9

#Usando datos expandidos Un elemento importante, sobre todo cuando estamos usando una encuesta con representatividad nacional es que podemos establecer los valores poblacionales. Para ello, los descriptivos se presentan con los datos de la población y no de la muestra. Esto se hace expandiendo los datos muestrales según el peso o factor de espansión.

En el caso de la enigh este es factor cambia de acuerdo de los módulos “fac”

Para tabulados podemos usar el comando tally

SDEMT219 %>%
    filter(SDEMT219$clase2==1 )%>%
      group_by(as_label(sex)) %>%
      tally(fac) #SUMA TODOS LOS RENGLONES DE LA VARIABLE fac
## # A tibble: 2 x 2
##   `as_label(sex)`        n
##   <fct>              <dbl>
## 1 Hombre          33771831
## 2 Mujer           21602905

Para sacar medias, podemos usar la función de media expandida

SDEMT219 %>% 
    filter(SDEMT219$clase2==1) %>%
      group_by(as_label(sex)) %>%
      summarise(avg_ing =weighted.mean(ing_x_hrs, weights = fac,na.rm = T))
## # A tibble: 2 x 2
##   `as_label(sex)` avg_ing
##   <fct>             <dbl>
## 1 Hombre             29.1
## 2 Mujer              28.1

También existen librerías para hacer tabulados de manera más tradicional. Es más fácil cuando queremos calcular porcentajes y proporciones. Usaremos la librería “questionr”, instalada en prácticas anteriores.

library(questionr)

Para usar nuestro factor de expansión primero tenemos que verificar qué tipo de valor es en nuestro objeto.

class(mydata$fac)
## [1] "numeric"

Para variables cualitativas podemos hacer tablas con los factores de expansión

wtd.table(mydata$sex, weights = mydata$fac)
##        1        2 
## 33771831 21602905
prop.table(wtd.table(mydata$sex, weights = mydata$fac))
##        1        2 
## 0.609878 0.390122
addmargins(prop.table(wtd.table(mydata$sex, weights = mydata$fac)))*100
##        1        2      Sum 
##  60.9878  39.0122 100.0000

Para una tabla cruzada

wtd.table(mydata$sex,mydata$t_loc, weights = mydata$fac)
##          1        2        3        4
## 1 16064549  4921153  4853674  7932455
## 2 11697396  3382344  2929406  3593759

Mi primera función

Unos de los elementos más poderosos de R es hacer nuestra propias funciones. Y cuando trabajamos bases de datos a veces hacemos un procedimiento que lo queremos repetir. Tal vez para una base de datos, tal vez para una variable, tal vez para grupos de casos.

Para ello haremos una función sencilla. Para sumarle un valor un 1

mi_funcion<-function(x) {
    x<-x+1
    x
}

mi_funcion(5)
## [1] 6

Qué tal hoy una una función donde podamos establecer qué valor sumamos

mi_funcion<-function(x,a) {
    x<-x+a
    x
}

mi_funcion(5,2)
## [1] 7

Mi función para expandir medias

expandir<- function(x){
  weighted.mean(x, w = SDEMT219$fac)
}

Entonces puedo ya sólo poner la variable

expandir(SDEMT219$eda)
## [1] NA

Otras funciones

percentage1<- function(x) {
  addmargins(prop.table(table(x)))*100
}
percentage2<- function(x, y) {
  addmargins(prop.table(table(x,y)))*100
}
percent_col<- function(x, y) {
  addmargins(prop.table(table(x,y),2))*100}
percent_row<- function(x, y) {
  addmargins(prop.table(table(x,y),1))*100
}

Guardando todo nuestro Environment

save.image("ambiente_p6.RData")