Trabajar encuestas de hogares con R (usando dplyr)

Para ilustrar los siguientes ejemplos se utiliza las siguientes bases de miembros simuladas. Los códigos de simulación se anexan al final de la presente entrada. Son dos bases ejemplos: una con información sobre los miembros y otra de las viviendas. Como siempre, el primer paso es estudiar nuestra base de datos.

data_miembros
# A tibble: 14 x 6
   hogar_id  edad sexo  zona   empleo ing_lab_h
      <int> <int> <chr> <chr>   <int>     <dbl>
 1        1    33 M     Rural       1     286.
 2        1    34 F     Rural       1     199.
 3        1    34 M     Rural       0       0 
 4        1    17 M     Rural      NA       0 
 5        1    35 M     Rural       0       0 
 6        1    35 M     Rural       1      87.8
 7        2    13 F     Urbana     NA       0 
 8        2     7 F     Urbana     NA       0 
 9        2    33 F     Urbana      1      94.2
10        3    29 F     Urbana      1     257.
11        3    17 M     Urbana     NA       0 
12        4    18 M     Rural       1     477.
13        5    12 M     Urbana     NA       0 
14        5    22 F     Urbana      0       0 
 
data_vivienda
  hogar_id  pared estufa habitaciones
1        1  Block      0            3
2        2 Madera      0            1
3        3  Block      1            3
4        4  Block      1            2
5        5  Block      1            1

1.       Agregando información a nivel de miembros

Podemos contabilizar la cantidad de miembros por cada hogar combinando las funciones group_by+summarise. Summarise hace un collapse de la base de datos a nivel de los miembros identificados.

data_miembros |>
  group_by(hogar_id) |>
  summarise(n()) |>
  dplyr::select(`n()`)
# A tibble: 5 x 1
  `n()`
  <int>
1     6
2     3
3     2
4     1
5     2

Note que al final se usa la función n() para contabilizar la cantidad de miembros y solo aparece esta variable en la data final. En caso que quisiéramos poner esta variable delante de la base de miembros, usaríamos mutate.

data_miembros |>
  group_by(hogar_id) |>
  mutate(miembro = 1:n())
# A tibble: 14 x 7
# Groups:   hogar_id [5]
   hogar_id  edad sexo  zona   empleo ing_lab_h miembro
      <int> <int> <chr> <chr>   <int>     <dbl>   <int>
 1        1    33 M     Rural       1     286.        1
 2        1    34 F     Rural       1     199.        2
 3        1    34 M     Rural       0       0         3
 4        1    17 M     Rural      NA       0         4
 5        1    35 M     Rural       0       0         5
 6        1    35 M     Rural       1      87.8       6
 7        2    13 F     Urbana     NA       0         1
 8        2     7 F     Urbana     NA       0         2
 9        2    33 F     Urbana      1      94.2       3
10        3    29 F     Urbana      1     257.        1
11        3    17 M     Urbana     NA       0         2
12        4    18 M     Rural       1     477.        1
13        5    12 M     Urbana     NA       0         1
14        5    22 F     Urbana      0       0         2

Note que podemos contabilizar otras características de los miembros. Podemos contar el número de miembros y trabajadores de cada hogar, sumando las variables de empleo. Al ser binaria (0,1), la suma de esta indica la cantidad de observaciones que cumplen dichas características. Esto se cumple con cualquier condición lógica, siempre que sumemos estos vectores, obtendremos la cantidad de elementos que cumplen dichas condiciones. Finalmente, podemos verificar que podemos crear nuevas variables usando las variables previamente creadas. Puntualmente, se crea la variable percent_trab que indica que porcentaje de los miembros de cada hogar esta trabajando.

data_miembros |>
  group_by(hogar_id) |>
  summarise(miembro = max(1:n()),
            trabajo_hogar = sum(as.numeric(empleo), na.rm=TRUE),
            porcent_trab= trabajo_hogar/miembro)
# A tibble: 5 x 4
  hogar_id miembro trabajo_hogar porcent_trab
     <int>   <int>         <dbl>        <dbl>
1        1       6             3        0.5 
2        2       3             1        0.333
3        3       2             1        0.5 
4        4       1             1        1   
5        5       2             0        0 

Adicionalmente, podemos agregar la edad promedio y la edad mínima de cada hogar.

data_miembros |>
  group_by(hogar_id) |>
  mutate(edad_m = mean(edad),
         edad_min = min(edad))
# A tibble: 14 x 8
# Groups:   hogar_id [5]
   hogar_id  edad sexo  zona   empleo ing_lab_h edad_m edad_min
      <int> <int> <chr> <chr>   <int>     <dbl>  <dbl>    <int>
 1        1    33 M     Rural       1     286.    31.3       17
 2        1    34 F     Rural       1     199.    31.3       17
 3        1    34 M     Rural       0       0     31.3       17
 4        1    17 M     Rural      NA       0     31.3       17
 5        1    35 M     Rural       0       0     31.3       17
 6        1    35 M     Rural       1      87.8   31.3       17
 7        2    13 F     Urbana     NA       0     17.7        7
 8        2     7 F     Urbana     NA       0     17.7        7
 9        2    33 F     Urbana      1      94.2   17.7        7
10        3    29 F     Urbana      1     257.    23         17
11        3    17 M     Urbana     NA       0     23         17
12        4    18 M     Rural       1     477.    18         18
13        5    12 M     Urbana     NA       0     17         12
14        5    22 F     Urbana      0       0     17         12

2.       Estimar la incidencia de la pobreza monetaria

Ahora se muestra un ejemplo de como estimar la pobreza monetaria. Asumiendo una línea de pobreza de 5,000RD$, como tenemos el ingreso por hora, necesitamos estimar los salarios el ingreso mensual ing_t_h y luego dividiendo ese ingreso entre el total de miembros, obtenemos el ingreso per cápita del hogar para finalmente obtener una variable binaria que verifica si este ingreso per cápita está por encima o por debajo de la línea de pobreza.

data_pobre <- data_miembros |>
  group_by(hogar_id) |>
  mutate(miembro = max(1:n()),
   ing_t_h = sum(ing_lab_h)*8*23.8,
  ing_pc_h = ing_t_h/max(1:n()),
  pobres = (ing_pc_h < 5000)*1)
 
data_pobre
# A tibble: 12 x 9
# Groups:   hogar_id [4]
   hogar_id  edad sexo  empleo ing_lab_h miembro ing_t_h ing_pc_h pobres
      <int> <int> <chr>  <int>     <dbl>   <int>   <dbl>    <dbl>  <dbl>
 1        1    30 M          1      428.       2  81554.   40777.      0
 2        1    18 M          0        0        2  81554.   40777.      0
 3        2    16 F         NA        0        6 206198.   34366.      0
 4        2    27 F          1      409.       6 206198.   34366.      0
 5        2    26 M          1      452.       6 206198.   34366.      0
 6        2    10 F         NA        0        6 206198.   34366.      0
 7        2     8 M         NA        0        6 206198.   34366.      0
 8        2    27 M          1      222.       6 206198.   34366.      0
 9        3    10 F         NA        0        3      0        0       1
10        3     6 F         NA        0        3      0        0       1
11        3    11 M         NA        0        3      0        0       1
12        4    11 F         NA        0        1      0        0       1

Ahora, como tenemos la variable pobre, podemos obtener la incidencia de la pobreza y la cantidad de personas pobres.

data_pobre |>
  ungroup() |>
  summarise(num_pob = sum(pobres),
            por_pobr = mean(pobres))
# A tibble: 1 x 2
  num_pob por_pobr
    <dbl>    <dbl>
1       4    0.333

Ahora, podemos segmentar la según una categoría determinada. Esto lo hacemos en los casos donde necesitamos obtener incidencia de la pobreza por grupos específicos de mi población. Según los datos el 50% de las mujeres está en condición de pobreza, frente a un 17% en caso de los hombres.

data_pobre |>
  group_by(sexo) |>
  summarise(num_pob = sum(pobres),
            por_pobr = mean(pobres))
# A tibble: 2 x 3
  sexo  num_pob por_pobr
  <chr>   <dbl>    <dbl>
1 F           3    0.5 
2 M           1    0.167

3.       Combinando data

Finalmente, dado que tememos una basa de vivienda donde obtenemos información sobre las características de esta, podemos estar interesados en estimar que porcentaje de mi población vive en casa con paredes de zinc. En tal sentido necesitamos llevar esta información a la base de miembros, o por el contrario llevar el numero de miembros de cada vivienda a la base de viviendas. Ambas opciones pueden hacerse mediante un join usando el id del hogar como llave. Esto es especialmente útil cuando necesitamos combinar información contenida en distintas bases de datos.

Alt. 1. Llevo la información de la data de vivienda a la data de miembros y luego contabilizo mis datos.

left_join(data_1, data_2, by="id")

Noten como opera la función, vamos a pesar la información referida en la data_2, sobre la data_1 (que esta a la izquierda), usando el id como identificador. Por ejemplo, para el id 1 le vamos a colocar el material de la vivienda.

left_join(data_miembros, data_vivienda, by="hogar_id")
# A tibble: 14 x 7
   hogar_id  edad sexo  zona   empleo ing_lab_h pared
      <int> <int> <chr> <chr>   <int>     <dbl> <chr>
 1        1    33 M     Rural       1     286.  Block
 2        1    34 F     Rural       1     199.  Block
 3        1    34 M     Rural       0       0   Block
 4        1    17 M     Rural      NA       0   Block
 5        1    35 M     Rural       0       0   Block
 6        1    35 M     Rural       1      87.8 Block
 7        2    13 F     Urbana     NA       0   Madera
 8        2     7 F     Urbana     NA       0   Madera
 9        2    33 F     Urbana      1      94.2 Madera
10        3    29 F     Urbana      1     257.  Block
11        3    17 M     Urbana     NA       0   Block
12        4    18 M     Rural       1     477.  Block
13        5    12 M     Urbana     NA       0   Block
14        5    22 F     Urbana      0       0   Block

Ahora, como tememos la información de la pared en la base de miembros, podemos obtener el porcentaje de individuos que viven en casas con un determinado tipo de material (si quisiéramos saber el porcentaje de vivienda, repetimos este proceso en la base de viviendas). [Nota: tenga pendiente que ahora con esta base podemos obtener cualquier información adicional sobre el perfil de la población según el material de la vivienda].

left_join(data_miembros, data_vivienda, by="hogar_id") |>
  group_by(pared) |>
  summarise(num_v = n()) |>
  mutate(percent = num_v/sum(num_v))
 
# A tibble: 2 x 3
  pared  num_v percent
  <chr>  <int>   <dbl>
1 Block     11   0.786
2 Madera     3   0.214

Alt. 2. Llevo información de la data de miembros a la data de vivienda. Por tanto, debo antes de, agregar la data de miembros para poder pasar la data. Primero hacemos esta agregación y obtenemos el numero de miembros por hogar.

data_miembros |>
    group_by(hogar_id) |>
    summarise(n())
# A tibble: 5 x 2
  hogar_id `n()`
     <int> <int>
1        1     6
2        2     3
3        3     2
4        4     1
5        5     2

Luego anexamos estas cantidades en al base de viviendas usando el mismo left join.

left_join(data_vivienda,
          data_miembros |> group_by(hogar_id) |> summarise(n()),
          by="hogar_id")
 
  hogar_id  pared n()
1        1  Block   6
2        2 Madera   3
3        3  Block   2
4        4  Block   1
5        5  Block   2

Ahora con esta base podemos calcular los porcentajes anteriores, solo debemos sumar las frecuencias de los miembros del hogar, contenida en la data anterior.

left_join(data_vivienda,
          data_miembros |> group_by(hogar_id) |> summarise(n()),
          by="hogar_id") |>
  group_by(pared) |>
  summarise(num_v = sum(`n()`))
# A tibble: 2 x 2
  pared  num_v
  <chr>  <int>
1 Block     11
2 Madera     3

4.       Porcentaje de personas por sexo según rango de edad 

De las herramientas vistas hasta aquí, primero se recodifica la variable edad por rangos usando la función cut. Posteriormente es cuestión de contabilizar categorías y obtener las frecuencias relativas.

 left_join(data_miembros, data_vivienda, by="hogar_id") |>
    mutate(r_edad = cut(edad, seq(0,100,20))) |>
    group_by(r_edad, sexo) |>
    summarise(num_v = n()) |>
    mutate(percent = num_v/sum(num_v))

# A tibble: 4 x 4
# Groups:   r_edad [2]
  r_edad  sexo  num_v percent
  <fct>   <chr> <int>   <dbl>
1 (0,20]  F         2   0.333
2 (0,20]  M         4   0.667
3 (20,40] F         4   0.5 
4 (20,40] M         4   0.5

Ahora se agrega la opción pivot_inder que se utiliza para convertir una data de formato long a wide. Ahora se lleva cada sexo a una columna para obtener mejor representación de la información.

left_join(data_miembros, data_vivienda, by="hogar_id") |>
    mutate(r_edad = cut(edad, seq(0,100,10))) |>
    group_by(r_edad, sexo) |>
    summarise(num_v = n()) |>
    mutate(percent = num_v/sum(num_v)) |>
    dplyr::select(-num_v) |>
    pivot_wider(names_from = sexo, values_from = percent)
 
# A tibble: 4 x 3
# Groups:   r_edad [4]
  r_edad      F      M
  <fct>   <dbl>  <dbl>
1 (0,10]  1     NA   
2 (10,20] 0.2    0.8 
3 (20,30] 1     NA   
4 (30,40] 0.333  0.667

5.       Data simulada

# simulando datos de miembros
 
set.seed(20)
library(tidyverse)
#simula data 
 
# Define number of households and members per household
num_hogares <- 5
max_miembros <- 6
 
num_miembro_hogar <- sample(1:max_miembros, num_hogares, replace = TRUE)
hogar_id <- rep(1:num_hogares, num_miembro_hogar)
 
data_miembros <- data.frame(hogar_id=hogar_id)
 
n<-nrow(data_miembros)
data_miembros$edad <- sample(5:35, n, replace = TRUE)
data_miembros$sexo <- sample(c("M", "F"), n, replace = TRUE)
data_miembros$zona <- ifelse(hogar_id %in% sample(1:num_hogares, round(num_hogares*0.3)), "Rural", "Urbana")
 
# Simulate employment status for each adult
 
data_miembros <- data_miembros |>
  mutate(
    ramdon_var = runif(n, min = 0, max = 1),
    empleo = case_when(
      edad >= 18 & ramdon_var>0.4~ "1",
      edad >= 18 & ramdon_var<0.4~ "0",
      TRUE ~ NA_character_
    ),
    empleo = as.integer(empleo),
    ing_lab_h = case_when(
      empleo == 1 ~ runif(n, min = 85, max = 500),
      empleo == 0 ~  0,
      TRUE ~ 0
    )
  ) |>
  as.tibble() |>
  dplyr::select(-ramdon_var)
 
# simulando datos de viviendas
 
data_vivienda <- data.frame(hogar_id=1:num_hogares)
data_vivienda$pared <- sample(c("Madera", "Concreto", "Block"), num_hogares, replace = TRUE)
data_vivienda$estufa <-  sample(c(0, 1), num_hogares, replace = TRUE)
data_vivienda$habitaciones <-  sample(c(1:3), num_hogares, replace = TRUE)

Test de independencia estadística en R

 1. Análisis preliminar

En la siguiente entrada se muestran algunos ejemplos de cómo realizar test de independencia estadística en R, para responder a uno de los ejercicios de la asignatura de Seminario de Investigación (clase 6 de la parte de R). Los ejercicios se realizan con la base utown del libro Principles of Econometrics 4e, by Carter Hill, William Griffiths, and Guay Lim. Los mismos se encuentran disponibles en el paquete PoEdata.

library(PoEdata)

library(tidyverse)
             
data(utown)        
head(utown,5)        

##     price  sqft age utown pool fplace
## 1 205.452 23.46   6     0    0      1
## 2 185.328 20.03   5     0    0      1
## 3 248.422 27.77   6     0    0      0
## 4 154.690 20.17   1     0    0      0
## 5 221.801 26.45   0     0    0      1

1. Análisis de independencia

1.1. Usando ggplo2, presente un histograma de precios y muestre un summary (min, max,median, q1, q5) condicionado a si la casa tiene o no piscina (pool).

 theme_set(theme_minimal())

#histograma
utown%>%
ggplot(aes(price)) +
geom_histogram(col = "black",fill = "darkblue", alpha = 0.6)

# Summary de precios por grupo
tapply(utown$price, utown$pool, summary)

## $`0`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##   134.3   214.8   244.3   246.5   276.8   345.2
##
## $`1`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##   160.3   219.5   249.1   252.0   285.4   340.7

1.2. Obtenga un boxplot donde relacione el tener piscina o no, con el precio de la casa, y establezca la independencia entre las variables.

utown%>%
  ggplot(aes(y=price, x=factor(pool), fill=factor(pool))) +
  geom_boxplot(col = "black", alpha = 0.6)

1.3. A partir de la comparación de promedios, que conclusión podemos obtener y que limitaciones puede enfrentar esta conclusión sobre la posible obtención de un efecto causal.

Al considerar los intervalos del gráfico anterior se puede observar que los promedios son relativamente parecidos, sin embargo, esta comparación no se puede interpretar como un efecto causal, dado que esta simple comparación omite otros factores relevantes que inciden sobre el precio de las casas.

1.4. Presente una tabla de frecuencia, donde los porcentajes se presenten en relación a las tablas, interprete sus resultados y posteriormente realice el test chi2 para determinar la posible independencia entre las variables.

Necesitamos una tabla de frecuencia de precio con pool, pero price es una variable continua, por lo que, debemos discretizarla, antes de realizar la tabla de frecuencia. Una alternativa es crear rangos de precios (cut(price, breaks = 5)), utilizar los quintiles (ejemplo debajo), o separarla en grupos alrededor de alguna medida como la media (price>mean(price)).

# discretizar la variable precio
utown <- utown %>%
  mutate(quintiprecio = ntile(price, 5))

utown %>%
select(quintiprecio, pool) %>%
table() %>%
prop.table(.,2) 

##             pool
## quintiprecio         0         1
##            1 0.2085427 0.1666667
##            2 0.2010050 0.1960784
##            3 0.1959799 0.2156863
##            4 0.2022613 0.1911765
##            5 0.1922111 0.2303922

Los valores son bastante diferentes, eso es evidencia de independencia, pero se debe realizar el test chi2 de independencia entre variables discretas (p.valor=0.5508, no se rechaza ho de independencia):

# discretizar la variable precio
utown %>%
select(quintiprecio, pool) %>%
table() %>%
chisq.test()

##
##  Pearson's Chi-squared test
##
## data:  .
## X-squared = 3.0422, df = 4, p-value = 0.5508

2. Análisis de independencia: Correlación

2.1. Sabiendo que tenemos una medida del size del terreno de la casa (sqft), crear los quintiles de esta variable y muestre una tabla resumen donde indique el precio max, min, mediano, q1 y q3 del precio de las casas para cada quintil.

utown %>%
  mutate(quintiprecio = ntile(price, 5)) %>%
  select(price,quintiprecio) %>%
  group_by(quintiprecio) %>%
  summarise(min = min(price),
    q1 = quantile(price,0.2),
    mean = mean(price),
    q3 = quantile(price,0.6),
    max = max(price))

## # A tibble: 5 x 6
##   quintiprecio   min    q1  mean    q3   max
##          <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1            1  134.  178.  190.  197.  208.
## 2            2  209.  214.  222.  225.  235.
## 3            3  235.  239.  246.  249.  258.
## 4            4  258.  263.  272.  274.  287.
## 5            5  287.  293.  309.  311.  345.

2.2. Obtenga un scatter donde relacione el size de la casa con el precio, interprete los resultados, coloreando los puntos según la casa tenga o no piscina (pool).

utown %>%
ggplot(aes(sqft, price, color =factor(pool)))+
geom_point()+
geom_smooth()

2.3. Realice el test de correlación de pearson y determine si rechaza o no ho.

Se rechaza ho.

cor.test(utown$price,utown$sqft, method = "pearson")

##
##  Pearson's product-moment correlation
##
## data:  utown$price and utown$sqft
## t = 23.367, df = 998, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5530750 0.6333235
## sample estimates:
##       cor
## 0.5946785

2.4. Utilice el qqplot para realizar un test de normalidad sobre el price y sqft.

car::qqPlot(utown$price)

2.5. En base al test de normalidad, que podemos concluir sobre las medidas de resumen y el análisis de correlación realizado.

Como para las colas se rechaza que la variable se comporte como una normal, se entiende que es preferible usar alguna estimación no parámetrica para obtener el coeficiente de correlación, tal como el método de Spearman.

3. Test de media

3.1. Realice un test de media para testear si las casas con chimeneas son estadísticamente más grandes que aquellas casas que no tienen.

t.test(sqft~fplace,
       data=utown,
       var.equal=F,
       alternative = "greater",
       conf.level=0.95)

##
##  Welch Two Sample t-test
##
## data:  sqft by fplace
## t = -3.1236, df = 989.19, p-value = 0.9991
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  -0.8778771        Inf
## sample estimates:
## mean in group 0 mean in group 1
##        24.91187        25.48674

3.2. Realice boxplot donde compare el precio de la casa según el size esté por encima o por debajo de la mediana.

utown%>%
  mutate(sqftmedian = sqft>median(sqft)) %>%
  ggplot(aes(y=price, x=factor(sqftmedian), fill=factor(sqftmedian))) +
  geom_boxplot(col = "black", alpha = 0.6)

3.3. Realice un test de media para testear si las casas con edad por encima de la mediana, tienen a tener un precio menor al resto de casa.

utown <- utown %>%
  mutate(medianEdad = as.numeric(age>median(age)))

t.test(price~medianEdad,
       data=utown,
       var.equal=F,
       alternative = "less",
       conf.level=0.95)

##
##  Welch Two Sample t-test
##
## data:  price by medianEdad
## t = 1.2334, df = 994.12, p-value = 0.8911
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##      -Inf 7.685081
## sample estimates:
## mean in group 0 mean in group 1
##        249.2916        246.0001