Para ésta práctica, además de las funciones incluidas con R “base”, usaremos la librería tidyverse

library(tidyverse)

# Hacemos que R se absetnga de usar notación científica para representar números con muchos espacios decimales
options(scipen = 20)

Los datos

Trabajaremos con datos de airbnb. Se trata de un listado de alojamientos ofrecidos en la Ciudad de Buenos Aires, recopilado en enero 2018.

airbnBA <- read.csv("https://bitsandbricks.github.io/data/Airbnb_listings_Buenos_Aires_Jan_2018.csv")

Revisando el contenido:

summary(airbnBA)
##     room_id            host_id                    room_type   
##  Min.   :    7270   Min.   :     2616   Entire home/apt:8059  
##  1st Qu.: 5701972   1st Qu.:  6715024   Private room   :1826  
##  Median :14093556   Median : 27259276   Shared room    : 231  
##  Mean   :12857197   Mean   : 48017123                         
##  3rd Qu.:19876958   3rd Qu.: 78877917                         
##  Max.   :23219172   Max.   :172636284                         
##                                                               
##  country          city         neighborhood  
##  Mode:logical   Mode:logical   Mode:logical  
##  NA's:10116     NA's:10116     NA's:10116    
##                                              
##                                              
##                                              
##                                              
##                                              
##                                                      address    
##  Buenos Aires, Ciudad Autónoma de Buenos Aires, Argentina:4032  
##  Buenos Aires, Autonomous City of Buenos Aires, Argentina:1406  
##  Palermo, Buenos Aires, Argentina                        : 725  
##  Buenos Aires, Buenos Aires, Argentina                   : 685  
##  Buenos Aires, CABA, Argentina                           : 590  
##  Recoleta, Buenos Aires, Argentina                       : 413  
##  (Other)                                                 :2265  
##     reviews       overall_satisfaction  accommodates       bedrooms     
##  Min.   :  0.00   Min.   :2.500        Min.   : 1.000   Min.   : 0.000  
##  1st Qu.:  0.00   1st Qu.:4.500        1st Qu.: 2.000   1st Qu.: 1.000  
##  Median :  3.00   Median :5.000        Median : 2.000   Median : 1.000  
##  Mean   : 13.86   Mean   :4.745        Mean   : 2.758   Mean   : 1.144  
##  3rd Qu.: 15.00   3rd Qu.:5.000        3rd Qu.: 4.000   3rd Qu.: 1.000  
##  Max.   :360.00   Max.   :5.000        Max.   :16.000   Max.   :50.000  
##                   NA's   :4694                          NA's   :26      
##    bathrooms          price          deleted  minstay       
##  Min.   : 0.000   Min.   :    0   Min.   :0   Mode:logical  
##  1st Qu.: 1.000   1st Qu.:  525   1st Qu.:0   NA's:10116    
##  Median : 1.000   Median :  897   Median :0                 
##  Mean   : 1.276   Mean   : 1199   Mean   :0                 
##  3rd Qu.: 1.000   3rd Qu.: 1345   3rd Qu.:0                 
##  Max.   :45.000   Max.   :30570   Max.   :0                 
##  NA's   :64                                                 
##               last_modified      latitude        longitude     
##  2018-02-08T16:52:39Z:   18   Min.   :-34.70   Min.   :-58.53  
##  2018-02-10T00:44:52Z:   18   1st Qu.:-34.60   1st Qu.:-58.43  
##  2018-02-10T01:45:47Z:   18   Median :-34.59   Median :-58.42  
##  2018-02-09T07:56:56Z:   17   Mean   :-34.59   Mean   :-58.42  
##  2018-02-09T09:57:08Z:   17   3rd Qu.:-34.58   3rd Qu.:-58.39  
##  2018-02-09T18:58:45Z:   17   Max.   :-34.53   Max.   :-58.35  
##  (Other)             :10011                                    
##    survey_id                                               location    
##  Min.   :1   0101000020E61000006EA301BC05344DC0471FF301814A41C0:   11  
##  1st Qu.:1   0101000020E610000058552FBFD3324DC0897E6DFDF44B41C0:   10  
##  Median :1   0101000020E61000002BFC19DEAC334DC0F3C64961DE4B41C0:    8  
##  Mean   :1   0101000020E61000004F22C2BF08364DC098DEFE5C344C41C0:    8  
##  3rd Qu.:1   0101000020E610000006D671FC50354DC0DAACFA5C6D4F41C0:    7  
##  Max.   :1   0101000020E6100000DDEEE53E39364DC00C789961A34A41C0:    7  
##              (Other)                                           :10065  
##  coworker_hosted extra_host_languages
##  Mode:logical    Mode:logical        
##  NA's:10116      NA's:10116          
##                                      
##                                      
##                                      
##                                      
##                                      
##                                name           property_type  currency   
##  Departamento en Recoleta        :   15   Apartment  :8319   ARS:10116  
##                                  :   14   House      : 860              
##  Departamento en Palermo         :   10   Loft       : 218              
##  Luminoso departamento en Palermo:    6   Condominium: 205              
##  Apartment in Recoleta           :    5   Other      : 197              
##  Cozy Apt. Heart of Buenos Aires :    5   Dorm       :  78              
##  (Other)                         :10061   (Other)    : 239              
##    rate_type    
##  nightly:10116  
##                 
##                 
##                 
##                 
##                 
##

Comparando valores promedio entre dos categorías: el test de Student o t-test

El test de Student se utiliza para comparar el valor medio de los valores extraidos de dos muestras distintas. Por ejemplo, la concentración de alcohol en dos tandas distintas de producción de cerveza, las notas de examen entre los alumnos de dos docentes, o -en éste caso- el precio por noche de dos categorías distintas de alojamientos en Airbnb.

Se lo debemos a este señor, William Sealy Gosset…

… que lo desarrolló alrededor de 1907 mientras trabajaba para la cervercería Guinnes. Gosset usó sus conocimientos de estadística para desarrollar un método que permitiera identificar entre dos cepas de centeno aquella con mejor rinde para la cerveza. Publicó su método en un paper bajo el pseudónimo de “Student” para aplacar a Guinness, que temía que la competencia pudiera identificar uan técnica útil para sus propios procesos.

Nosotros usaremos el test de Student, también llamado t-test, para establecer si existe una diferencia estadísticamente significativa entre el valor por noche de habitaciones privadas versus compartidas.

¿Se ahorra alquilando una habitación compartida en lugar de una privada?

Primero, revisamos los valores medios para cada categoría:

airbnBA %>% 
    filter(room_type %in% c("Private room", "Shared room"), accommodates == 1) %>%  
    group_by(room_type) %>% 
    summarise(cantidad = n(),
              valor_promedio = mean(price))
## # A tibble: 2 x 3
##   room_type    cantidad valor_promedio
##   <fct>           <int>          <dbl>
## 1 Private room     1048           356.
## 2 Shared room       112           303.

Como era de esperarse, las habitaciones privadas tienen un valor más alto. Pero podemos asumir que en efecto hay una dinámica distinta de precios, o los datos no son concluyentes?

privadas_vs_compartidas <- airbnBA %>% 
    filter(room_type %in% c("Private room", "Shared room"), accommodates == 1)   


ggplot() +
    geom_density(data = privadas_vs_compartidas, 
                   aes(x = price)) +
    facet_wrap(~room_type, scales = "free_y")

TODO: Hablar de como verificar si la varianza es igual o no

T-TEST time!!

En líneas generales, un t test se realiza así: t.test(valor_de_cada_observacion ~ categoría_de_cada_observacion)

En código:

t.test(privadas_vs_compartidas$price ~ privadas_vs_compartidas$room_type)
## 
##  Welch Two Sample t-test
## 
## data:  privadas_vs_compartidas$price by privadas_vs_compartidas$room_type
## t = 5.3371, df = 134.06, p-value = 0.0000003919
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  33.04602 71.95861
## sample estimates:
## mean in group Private room  mean in group Shared room 
##                   355.9666                   303.4643

Podemos reportar los resultados así:

Para el alojamiento de un huésped, el test de medias halló un precio promedio significativamente* menor en las habitaciones compartidas respecto a las privadas (356 vs 303 pesos por noche), con un p-valor menor a 0,001.

*: aquí hablamos de significancia estadística, que es independiente de que la diferencia hallada sea “significativa” en el sentido de “contundente”

Comparando valores promedio entre múltiples categorías: el test ANOVA

Un test ANOVA (“Analysis of Variance”) sirve cuando queremos comparar atribtuos de múltiples grupos par determinar si existe alguna diferencia entre ellos. Por ejemplo,

Buscando un ejemplo relacionado con nuestros datos, podríamos preguntarnos…

Cómo impacta el tipo de alojamiento en el precio por noche?

Evaluemos ahora el precio por alojamientos en propiedades de distinta clase. Podemos echar un vistazo a las categorías encontradas en la base, la frecuencia con la que aparecen, y su valor promedio:

airbnBA %>% 
    group_by(property_type) %>% 
    summarise(cantidad = n(),
              valor_promedio = mean(price)) %>% 
    arrange(desc(valor_promedio))
## # A tibble: 23 x 3
##    property_type      cantidad valor_promedio
##    <fct>                 <int>          <dbl>
##  1 Boat                      3          6771.
##  2 Cabin                     1          5156 
##  3 Serviced apartment       18          1820.
##  4 Loft                    218          1594.
##  5 Villa                     1          1325 
##  6 House                   860          1294.
##  7 Boutique hotel           11          1238.
##  8 Timeshare                 1          1223 
##  9 Condominium             205          1205.
## 10 Apartment              8319          1203.
## # … with 13 more rows

Podemos aislar sólo los alojamientos con una cama, para no comparar capacidades distintos entre si (es decir, no mezclar en el promedio alojamientos para una pareja con otros donde pueden dormir muchas mas personas):

airbnBA %>% 
    filter(bedrooms == 1) %>% 
    group_by(property_type) %>% 
    summarise(cantidad = n(),
              valor_promedio = mean(price)) %>% 
    arrange(desc(valor_promedio))
## # A tibble: 20 x 3
##    property_type       cantidad valor_promedio
##    <fct>                  <int>          <dbl>
##  1 Boat                       3          6771.
##  2 Boutique hotel             4          2498 
##  3 Serviced apartment         6          1874 
##  4 Loft                     135          1361.
##  5 Condominium              135           992.
##  6 Apartment               5207           985.
##  7 Chalet                     1           836 
##  8 Other                    141           719.
##  9 Townhouse                  9           706.
## 10 Casa particular            5           534.
## 11 Vacation home              1           525 
## 12 Camper/RV                  2           505 
## 13 In-law                    12           494.
## 14 Pension (Korea)            2           453 
## 15 House                    599           423.
## 16 Guest suite               44           383.
## 17 Dorm                      72           376.
## 18 Bed &amp; Breakfast       52           374.
## 19 Guesthouse                22           359.
## 20 Hostel                    26           307.

Excepto para los botes, (hay 3 botes!) el promedio desciende al concentrarnos en alojamientos de una sola cama.

También deberíamos retirar las categorías que aparecen pocas veces, dado que con una muestra tan chica podrían no ser representativas. Así podemos elegir sólo aquellas categorías que tienen más de 100 representantes.

Identificamos los tipos frecuentes:

tipos_frecuentes <- airbnBA %>% 
    filter(bedrooms == 1) %>% 
    group_by(property_type) %>% 
    summarise(cantidad = n()) %>% 
    filter(cantidad > 100) %>% 
    pull(property_type)
tipos_frecuentes
## [1] Apartment   Condominium House       Loft        Other      
## 23 Levels: Apartment Bed &amp; Breakfast Boat Boutique hotel ... Villa

Y los usamos para filtrar el dataset, extrayendo sólo los alojamientos que tienen un tipo dentro de la lista identificada:

propiedades_1dorm <- airbnBA %>% 
    filter(property_type %in% tipos_frecuentes, bedrooms == 1)

Cómo están distribuidos?

ggplot() +
    geom_histogram(data = propiedades_1dorm, 
                   aes(x = price, fill = property_type)) +
    facet_wrap(~property_type, scales = "free_y") +
    guides(fill = FALSE)

Aparecen claros valores extremos, “outliers” que empujan la distribución hacia la derecha (hacia valores michísimo más altos que el promedio).

En este caso particular, una forma efectiva de retirar del análisis aquellos alojamientos que figuran con precios poco creíbles es retener solamente aquellos que tengan reviews; es decir, que sepamos que hayan sido alquilados por alguien alguna vez, dado que han dejado la opinión post estadía.

propiedades_1dorm <- airbnBA %>% 
    filter(property_type %in% tipos_frecuentes, bedrooms == 1, reviews > 0)

ggplot() +
    geom_histogram(data = propiedades_1dorm, 
                   aes(x = price, fill = property_type)) +
    facet_wrap(~property_type, scales = "free_y") +
    guides(fill = FALSE)

ANOVA time!

aov(data = misdatos, variable_dependiente ~ variable_independiente)

analisis_tipo_propiedad <- aov(data = propiedades_1dorm, price ~ property_type)

analisis_tipo_propiedad
## Call:
##    aov(formula = price ~ property_type, data = propiedades_1dorm)
## 
## Terms:
##                 property_type  Residuals
## Sum of Squares      109381676 1707874538
## Deg. of Freedom             4       4375
## 
## Residual standard error: 624.797
## Estimated effects may be unbalanced
summary(analisis_tipo_propiedad)
##                 Df     Sum Sq  Mean Sq F value              Pr(>F)    
## property_type    4  109381676 27345419   70.05 <0.0000000000000002 ***
## Residuals     4375 1707874538   390371                                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

De acuedo a nuestro análisis, no caben dudas de que el tipo de alojamiento afecta su precio? Pero… ¿cuánto? ¿Qué tipos de alojamiento exhiben diferencias mayores?

El test ANOVA no nos permite extraer conclusiones espécificas; sólo indica que existe una diferencia significativa entre los grupos.

Para poder extraer más detalles, podemos recurrr a…

“El Test Honesto de Diferencias Significativas” de Tukey, o HSD

Desarrollado por el estadístico John Tukey, que tenía todo un talento para nombrar cosas -se lo acredita también como inventor de la palabra “bit” para expresar la mínima unidad de datos en un sistema de información.

Al hacer el test HSD, asumimos tres cosas:

  • Las observaciones son independientes a nivel intra y extra grupos
  • El atributo medido tiene una distribución normal a nivel grupo
  • La varianza de cada grupo es similar entre todos.

Para realizarlo en R podemos usar la función TukeyHSD(). Como primer argumento usamos un objeto “aov”, el resultado de un análisis ANOVA como el que hicimos en el paso anterior. El segundo argumento para la función es la variable independiente que define los grupos, en nuestro caso el tipo de propiedad:

analisis_tukey_tipo_propiedad <- TukeyHSD(analisis_tipo_propiedad, "property_type")

Veamos los resultados

analisis_tukey_tipo_propiedad
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = price ~ property_type, data = propiedades_1dorm)
## 
## $property_type
##                             diff        lwr        upr     p adj
## Condominium-Apartment  -85.94412 -264.93784   93.04959 0.6848290
## House-Apartment       -529.64574 -622.35245 -436.93902 0.0000000
## Loft-Apartment         279.20131  107.26296  451.13965 0.0000938
## Other-Apartment       -252.00926 -437.97454  -66.04399 0.0020504
## House-Condominium     -443.70161 -641.37228 -246.03095 0.0000001
## Loft-Condominium       365.14543  120.11084  610.18002 0.0004661
## Other-Condominium     -166.06514 -421.13815   89.00787 0.3873750
## Loft-House             808.84704  617.54170 1000.15239 0.0000000
## Other-House            277.63647   73.63151  481.64144 0.0019299
## Other-Loft            -531.21057 -781.38308 -281.03806 0.0000001

El test compara los grupos de a pares, y estima la diferencia entre sus medias, y el grado de significancia estádistica. Por ejemplo, aquí ha establecido que la diferencia entre condominios y apartamentos es de unos 86 dólares a favor de la segunda categoría (la media de los condos es 86 dólares menor) pero con un p-valor muy por encima de 0.05 debemos aceptar la hipótesis nula: no hay una diferencia significativa. Por otro lado, diferencias entre pares como House y Apartment, o entre Loft y Condominium si muestran significancia.

Prediciendo un resultado binario (“Sí” o “No”): La regresión logística

En base a las preferencias y características demográficas de una persona, ¿podemos predecir si votará a un partido político u otro?

Para el ejercicio usaremos datos recopilados por el proyecto ANES, American National Election Studies. ANES realiza en forma periódica extensas encuestas individuales en EEUU, para proveer información que permita a Cientistas Sociales teorizar y explicar resultados electorales. Se pueden encontrar en https://electionstudies.org/

Comencemos por cargar un dataset con atributos a nivel persona, que incluye el partido por el cual votaron en las elecciones presidenciales de 2008 en EEUU.

votantes <- read.csv("https://bitsandbricks.github.io/data/anes_2008tr.csv")
names(votantes)
## [1] "edad"                   "raza_blanco_a"         
## [3] "genero_varon"           "educacion"             
## [5] "ingresos"               "ideologia_conservadora"
## [7] "ha_votado_republicano"  "ha_votado"
summary(votantes)
##       edad       raza_blanco_a     genero_varon      educacion    
##  Min.   : 0.00   Min.   :0.0000   Min.   :0.0000   Min.   :1.000  
##  1st Qu.:32.00   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:3.000  
##  Median :46.00   Median :1.0000   Median :0.0000   Median :5.000  
##  Mean   :46.45   Mean   :0.5043   Mean   :0.4302   Mean   :4.127  
##  3rd Qu.:59.00   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:5.000  
##  Max.   :93.00   Max.   :1.0000   Max.   :1.0000   Max.   :7.000  
##                                                                   
##     ingresos     ideologia_conservadora ha_votado_republicano
##  Min.   :1.000   Min.   :1.000          Min.   :0.000        
##  1st Qu.:2.000   1st Qu.:4.000          1st Qu.:0.000        
##  Median :3.000   Median :4.000          Median :0.000        
##  Mean   :2.705   Mean   :4.098          Mean   :0.334        
##  3rd Qu.:3.000   3rd Qu.:5.000          3rd Qu.:1.000        
##  Max.   :5.000   Max.   :7.000          Max.   :1.000        
##                                         NA's   :783          
##    ha_votado     
##  Min.   :0.0000  
##  1st Qu.:0.0000  
##  Median :1.0000  
##  Mean   :0.6904  
##  3rd Qu.:1.0000  
##  Max.   :1.0000  
##

Exploremos algunas de las variables

ggplot() +
    geom_histogram(data = votantes, aes(x = edad))

ggplot() +
    geom_histogram(data = votantes, aes(x = edad)) +
    facet_wrap(~ha_votado_republicano)

ggplot() +
    geom_bar(data = votantes, aes(x = ideologia_conservadora))

ggplot() +
    geom_bar(data = votantes, aes(x = ideologia_conservadora)) +
    facet_wrap(~ha_votado_republicano)

ggplot() +
    geom_bar(data = votantes, aes(x = as.factor(educacion))) +
    facet_wrap(~raza_blanco_a)

¡Listos para hacer una regresión logística! En R se resuelve de manera similar a la regresión lineal, con una fórmula que expresa la variable a predecir versus las predictoras. La función a la que recurriremos es glm(), que permite utilizar varias técnicas de regresión; con family="binomial" le aclaramos que queremos una regresión logística.

Por ejemplo, para establecer el grado en que la edad, la raza, el nivel de ingresos y la inclinación ideológica predicen el voto republicano, usamos:

vota_republicano <- glm(data = votantes, 
                        ha_votado_republicano ~ edad + raza_blanco_a + ingresos + ideologia_conservadora,
                        family="binomial")
summary(vota_republicano)
## 
## Call:
## glm(formula = ha_votado_republicano ~ edad + raza_blanco_a + 
##     ingresos + ideologia_conservadora, family = "binomial", data = votantes)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3038  -0.4875  -0.2900   0.4763   3.3359  
## 
## Coefficients:
##                         Estimate Std. Error z value             Pr(>|z|)
## (Intercept)            -8.076556   0.462932 -17.447 < 0.0000000000000002
## edad                    0.004738   0.004190   1.131                0.258
## raza_blanco_a           2.436605   0.167899  14.512 < 0.0000000000000002
## ingresos                0.404587   0.074955   5.398         0.0000000675
## ideologia_conservadora  1.041977   0.066277  15.721 < 0.0000000000000002
##                           
## (Intercept)            ***
## edad                      
## raza_blanco_a          ***
## ingresos               ***
## ideologia_conservadora ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1960.6  on 1538  degrees of freedom
## Residual deviance: 1170.3  on 1534  degrees of freedom
##   (783 observations deleted due to missingness)
## AIC: 1180.3
## 
## Number of Fisher Scoring iterations: 5

Que tan bien predice nuestra regresión logística?

La regresión logística calcula un porcentaje de probabilidad (del 0 al 100%) de que cada votante haya votado republicano de acuerdo a los atributos utilizados como predictores. Una vez realizado un modelo, como el que hemos guardado en la variable vota_republicano, se puede usar la funcion predict() para obtener los valores que el modelo “predice” para cualquier conjunto de datos, por ejemplo con el resultado de futuras encuestas.

Si usamos predict() con los mismos datos con los que hicimos el modelo, lo que obtenemos es la probabilidad de votar republicano asignada por el modelo a cada uno de los votantes en el dataset.

Verifiquemos el porcentaje de votantes que es clasificado en forma correcta utilizando

Veamos esas probabilidadaes:

predicciones <- predict(vota_republicano, votantes) 

head(predicciones)
##          1          2          3          4          5          6 
## -2.1244712  2.1004745  2.4150400  0.3832037  2.5713932  2.1620682

Aquí hay un problema. En lugar de valores entre 0 y 1, expresando del 0 al 10o% de pobabilidad, hemos obtenido números mayores a 1, y menores a 0… negativos! ¿Qué ocurre?

Las probabilidades están expresadas en “logits”. Son el logaritmo de las chances u odds (por ejemplo “una en mil”) de que cada votante sea republicano.

La buena noticia es que es fácil pasar de logits a chances: usamos los logits como exponente de e, y obtenemos las chances. En R es tan simple como exp(predicciones)

Y luego es facil pasar de chances a probabilidad porcentual; la cuenta es:

\[ \Leftrightarrow \dfrac{chances}{1 + chances} = P(Y = 1)\]

Donde P(Y = 1) es, de acuerdo a las chances estilo 1 en 10000, la probabilidad de que la variable Y sea de la categoría que estamos analizando… en éste caso, que la persona haya votado por el partido republicano.

Calculemos entonces las probabilidades de nuestra predicción, y agreguémoslas a los datos:

votantes <- votantes %>% 
  mutate(prediccion.logit = predict(vota_republicano, votantes),
         prediccion.odds = exp(prediccion.logit),
         prediccion.probabilidad = prediccion.odds / (1 + prediccion.odds))

Y dado que sabemos si el votante optó o no por el partido Republicano, hagamos una tabla de cruce y veamos que tan buenas resultan las predicciones.

Si la probilidad es mayor a 50% (mayor a 0.5) consideramos que se predice un votante republicano. de lo contrario, que no vota por ese partido:

votantes <- votantes %>% 
  mutate(predice_republicano = ifelse(prediccion.probabilidad > 0.5, 1, 0))

y ahora, la tabla cruzada:

xtabs(data = votantes, ~ predice_republicano + ha_votado_republicano)
##                    ha_votado_republicano
## predice_republicano   0   1
##                   0 921 160
##                   1 104 354

Total de casos: 921 + 104 + 160 + 354 = 1539

Aciertos: 921 + 354 = 1275

Tasa de acierto: 1275 / 1539 = 0.82846

Casi 83% de acierto… resulta mejor que predecir tirando una moneda al aire.