¡Recuerda, poner el directorio! Además vamos a cargar el ambiente que hemos estado trabajando desde la práctica 3
setwd("/Users/anaescoto/Dropbox/2019/Cursos ESA/UCA_R")
load("./datos/ehpm2017.RData")
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.
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
Empecemos con el “boxplot” univariado y con la función que viene en “graphics”, sin necesidad de instalar paquetería
boxplot(ehpm2017[ehpm2017$actpr2012==10,]$money)
La función boxplot, también sirve para obtener el resumen de cinco números.
a<-boxplot(ehpm2017[ehpm2017$actpr2012==10,]$money)
a$stats
## [,1]
## [1,] 0.0
## [2,] 80.0
## [3,] 195.0
## [4,] 325.0
## [5,] 692.4
Con formato dplyr
ehpm2017 %>%
filter(actpr2012==10) %>%
with(boxplot(money))
Podemos comparar la distribución de una variable cuantitativa entre las categorías de una variable categórica
#otra forma
boxplot(money ~ as_label(r104), data= ehpm2017[ehpm2017$actpr2012==10,])
Con dplyr:
ehpm2017 %>%
filter(actpr2012==10) %>%
with(boxplot(money ~ as_label(r104)))
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 <-ehpm2017 %>%
filter(actpr2012==10 & money<30000) %>%
ggplot(aes(x=as_label(r104), y=money))
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()
g + facet_wrap(~ as_label(region))
g2<-lienzo + geom_boxplot(aes(color=as_label(region))) +
theme_classic()
g2
Para poder editar la transparencia de nuestro relleno podemos usar la opción alpha y ponerle un valor menor a 1.
lienzo_bi <-ehpm2017 %>%
filter(actpr2012==10 & money<30000) %>%
ggplot(aes(x=log(money+1), fill=as_label(r104),
color=as_label(r104),
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 = "r104", color="r104")
Un truco para hacer el análisis descriptivo es generar una sub-base con las variables que vamos analizar
mydata<- ehpm2017[ehpm2017$actpr2012==10, c("money", "region","r106", "r104", "aproba1", "fac00")]
tail(mydata)
## # A tibble: 6 x 6
## money region r106 r104 aproba1 fac00
## <dbl> <dbl+lbl> <dbl> <dbl+lbl> <dbl> <dbl>
## 1 12 4 [Oriental] 84 2 [Mujer] 0 71
## 2 NA NA NA NA NA NA
## 3 0 4 [Oriental] 50 2 [Mujer] 3 71
## 4 0 4 [Oriental] 75 1 [Hombre] 0 71
## 5 52 4 [Oriental] 24 1 [Hombre] 7 71
## 6 3 4 [Oriental] 29 2 [Mujer] 2 71
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
## money region r106 r104 aproba1 fac00
## 239.301423 NA 38.039946 NA 7.546263 89.166950
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
## money region
## Min. : 0.0 Occidental :7475
## 1st Qu.: 80.0 Central I :7434
## Median : 195.0 Central II :5228
## Mean : 239.3 Oriental :6957
## 3rd Qu.: 325.0 Área metropolitana de San Salvador:4724
## Max. :32920.0 NA's :6286
## NA's :6286
## r106 r104 aproba1 fac00
## Min. : 6.00 Hombre:19148 Min. : 0.000 Min. : 13.00
## 1st Qu.:25.00 Mujer :12670 1st Qu.: 4.000 1st Qu.: 41.00
## Median :36.00 NA's : 6286 Median : 8.000 Median : 62.00
## Mean :38.04 Mean : 7.546 Mean : 89.17
## 3rd Qu.:48.00 3rd Qu.:12.000 3rd Qu.:104.00
## Max. :97.00 Max. :25.000 Max. :710.00
## NA's :6286 NA's :6286 NA's :6286
summary(as_label(mydata[mydata$r104==1,])) #hombres
## money region r106
## Min. : 0 Occidental :4673 Min. : 6.00
## 1st Qu.: 65 Central I :4427 1st Qu.:24.00
## Median : 208 Central II :3210 Median :35.00
## Mean : 243 Oriental :4297 Mean :37.65
## 3rd Qu.: 325 Área metropolitana de San Salvador:2541 3rd Qu.:49.00
## Max. :32920 NA's :6286 Max. :96.00
## NA's :6286 NA's :6286
## r104 aproba1 fac00
## Hombre:19148 Min. : 0.000 Min. : 13.00
## Mujer : 0 1st Qu.: 3.000 1st Qu.: 41.00
## NA's : 6286 Median : 7.000 Median : 61.00
## Mean : 7.326 Mean : 87.27
## 3rd Qu.:11.000 3rd Qu.:101.00
## Max. :25.000 Max. :710.00
## NA's :6286 NA's :6286
summary(as_label(mydata[mydata$r104==2,])) #mujeres
## money region
## Min. : 0.00 Occidental :2802
## 1st Qu.: 86.67 Central I :3007
## Median : 173.33 Central II :2018
## Mean : 233.74 Oriental :2660
## 3rd Qu.: 325.00 Área metropolitana de San Salvador:2183
## Max. :8012.00 NA's :6286
## NA's :6286
## r106 r104 aproba1 fac00
## Min. : 7.00 Hombre: 0 Min. : 0.000 Min. : 13.00
## 1st Qu.:27.00 Mujer :12670 1st Qu.: 4.000 1st Qu.: 42.00
## Median :37.00 NA's : 6286 Median : 8.000 Median : 63.00
## Mean :38.62 Mean : 7.879 Mean : 92.04
## 3rd Qu.:48.00 3rd Qu.:12.000 3rd Qu.:109.00
## Max. :97.00 Max. :25.000 Max. :710.00
## NA's :6286 NA's :6286 NA's :6286
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$money, as_label(mydata$r104), summary)
## $Hombre
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 65 208 243 325 32920
##
## $Mujer
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 86.67 173.33 233.74 325.00 8012.00
tapply(mydata$money, as_label(mydata$r104), mean, na.rm=T)
## Hombre Mujer
## 242.9812 233.7402
Con dplyr a veces es más sencillo por la función “as group” y summarise (o summarize)
ehpm2017 %>%
filter(ehpm2017$actpr2012==10 )%>%
summarise(avg_ing= mean(money))
## # A tibble: 1 x 1
## avg_ing
## <dbl>
## 1 239.
Con la función de grupos es más sencillo que tapply y sapply
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 )%>%
group_by(as_label(r104)) %>%
tally() # Nos da el conteo de casos
## # A tibble: 2 x 2
## `as_label(r104)` n
## <fct> <int>
## 1 Hombre 19148
## 2 Mujer 12670
Y podemos “summarizar” los indicadores que querramos. Recuerda que estos resultados se pueden guardar en un objeto.
resultados1<-ehpm2017 %>%
filter(ehpm2017$actpr2012==10 )%>%
group_by(as_label(r104)) %>%
summarise(avg_ing = mean(money),
avg_esc = mean(aproba1))
resultados1
## # A tibble: 2 x 3
## `as_label(r104)` avg_ing avg_esc
## <fct> <dbl> <dbl>
## 1 Hombre 243. 7.33
## 2 Mujer 234. 7.88
#Usando datos expandidos Un elemento importante, sobre todo cuando estamos usando una encuesta con representatividad nacional como la EHOM DE 2017, 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 “fac00”
Para tabulados podemos usar el comando tally
ehpm2017 %>%
filter(ehpm2017$actpr2012==10 )%>%
group_by(as_label(r104)) %>%
tally(fac00) #SUMA TODOS LOS RENGLONES DE LA VARIABLE FAC00
## # A tibble: 2 x 2
## `as_label(r104)` n
## <fct> <dbl>
## 1 Hombre 1670973
## 2 Mujer 1166141
Para sacar medias, podemos usar la función de media expandida
ehpm2017 %>%
filter(ehpm2017$actpr2012==10) %>%
group_by(as_label(r104)) %>%
summarise(avg_ing =weighted.mean(money, weights = fac00,na.rm = T))
## # A tibble: 2 x 2
## `as_label(r104)` avg_ing
## <fct> <dbl>
## 1 Hombre 243.
## 2 Mujer 234.
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$fac00)
## [1] "numeric"
Para variables cualitativas podemos hacer tablas con los factores de expansión
wtd.table(mydata$r104, weights = mydata$fac00)
## 1 2
## 1670973 1166141
prop.table(wtd.table(mydata$r104, weights = mydata$fac00))
## 1 2
## 0.5889693 0.4110307
addmargins(prop.table(wtd.table(mydata$r104, weights = mydata$fac00)))*100
## 1 2 Sum
## 58.89693 41.10307 100.00000
Para una tabla cruzada
wtd.table(mydata$r104,mydata$region, weights = mydata$fac00)
## 1 2 3 4 5
## 1 383636 371708 177382 324969 413278
## 2 249718 240018 116706 204680 355019
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 = ehpm2017$fac00)
}
Entonces puedo ya sólo poner la variable
expandir(ehpm2017$r106)
## [1] 31.3358
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
}
save.image("ambiente_p6.RData")
Envíe el código https://tinyurl.com/Ej6-ESA-R