TP 1: Análisis exploratorio e Introducción a Regresión lineal

Axel Mantalian (axel.manta@gmail.com)

1) Preparación de los datos (I)

Cargo las librerías necesarias

library(tidyverse)
library(corrr)
library(ggplot2)
library(ggthemes)
library(GGally)
library(data.table)
library(scales)

Leo el archivo csv

datos <- read_csv("./data/ar_properties.csv",progress = FALSE)

Estructura del dataset, conformado por 24 variables y 388891 observaciones de avisos inmobiliarios.

# Muestro la estructura 
head(datos)

glimpse(datos)

Observations: 388,891
Variables: 24
$ id              <chr> "S0we3z3V2JpHUJreqQ2t/w==", "kMxcmAS8NvrynGBVbMOEaQ==", "Ce3ojF+ZTOkB8d+LI9dpxg=="…
$ ad_type         <chr> "Propiedad", "Propiedad", "Propiedad", "Propiedad", "Propiedad", "Propiedad", "Pro…
$ start_date      <date> 2019-04-14, 2019-04-14, 2019-04-14, 2019-04-14, 2019-04-14, 2019-04-14, 2019-04-1…
$ end_date        <date> 2019-06-14, 2019-04-16, 9999-12-31, 9999-12-31, 2019-07-09, 2019-08-08, 2019-07-1…
$ created_on      <date> 2019-04-14, 2019-04-14, 2019-04-14, 2019-04-14, 2019-04-14, 2019-04-14, 2019-04-1…
$ lat             <dbl> -34.94331, -34.63181, NA, -34.65471, -34.65495, -32.93547, -34.65183, -34.91213, -…
$ lon             <dbl> -54.92966, -58.42060, NA, -58.79089, -58.78712, -60.68398, -58.65912, -54.84749, -…
$ l1              <chr> "Uruguay", "Argentina", "Argentina", "Argentina", "Argentina", "Argentina", "Argen…
$ l2              <chr> "Maldonado", "Capital Federal", "Bs.As. G.B.A. Zona Norte", "Bs.As. G.B.A. Zona Oe…
$ l3              <chr> "Punta del Este", "Boedo", NA, "Moreno", "Moreno", "Rosario", "Ituzaingó", "José I…
$ l4              <chr> NA, NA, NA, "Moreno", "Moreno", NA, "Ituzaingó", NA, NA, NA, NA, NA, NA, NA, NA, N…
$ l5              <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ l6              <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ rooms           <dbl> 2, NA, 2, 2, 2, 4, NA, 6, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1, 1, 1,…
$ bedrooms        <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ bathrooms       <dbl> 1, NA, 1, 2, 3, 1, 3, 3, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1, NA, 1, 1, 1, 1…
$ surface_total   <dbl> 45, NA, 200, 460, 660, NA, 70, NA, 1300, 405, 352, 373, 360, 1325, 250, 80142, 101…
$ surface_covered <dbl> 40, NA, NA, 100, 148, 89, 122, NA, NA, NA, NA, NA, NA, 2, NA, NA, NA, NA, 54, 180,…
$ price           <dbl> 13000, 0, NA, NA, NA, NA, NA, NA, 0, NA, 0, NA, NA, NA, NA, NA, 0, 0, 0, NA, 0, 0,…
$ currency        <chr> "UYU", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ price_period    <chr> "Mensual", "Mensual", NA, "Mensual", "Mensual", "Mensual", "Mensual", "Mensual", "…
$ title           <chr> "Departamento - Roosevelt", "PH - Boedo", "Ituzaingo  1100 - $ 1 - Casa Alquiler",…
$ property_type   <chr> "Departamento", "PH", "Casa", "Casa", "Casa", "Casa", "Casa", "Casa", "Lote", "Lot…
$ operation_type  <chr> "Alquiler", "Venta", "Alquiler", "Venta", "Venta", "Venta", "Venta", "Alquiler", "…

Filtrado y selección de variables

Aplicando diversos filtros, reduzco el dataset a un total de 61905 observaciones. También, de las 24 variables originales solamente conservo 9.

Los filtros aplicados son:

• Pertenecen a Argentina y Capital Federal
• Cuyo precio está en dólares (USD)
• El tipo de propiedad es: Departamento, PH o Casa
• El tipo de operacion es: Venta

Las variables seleccionadas son:

• id
• l3
• rooms
• bedrooms
• bathrooms
• surface_total
• surface_covered
• price
• property_type

datos <- datos %>% 
  filter(l1 == "Argentina" 
         & l2 == "Capital Federal"
         & currency == "USD"
         & property_type %in% c("Departamento","PH","Casa")
         & operation_type == "Venta") %>% 
  select(id,l3,rooms,bedrooms,bathrooms,surface_total,surface_covered,price,property_type)
datos

2) Análisis exploratorios (I)

Se comienza el análisis exploratorio verificando la cantidad de valores únicos y de valores faltantes (NAs) para cada una de las variables seleccionadas.

# Valores únicos
distintos <- map_dfc(datos,function(x) n_distinct(x))
# Valores faltantes
faltantes <- map_dfc(datos,function(x) sum(is.na(x)))
# Junto los datos y les doy formato de tabla
t_unicos_faltantes <- rbind(distintos,faltantes) %>% 
  transpose() 
  
rownames(t_unicos_faltantes) <- colnames(distintos)
colnames(t_unicos_faltantes) <- c('Valores únicos','Valores faltantes')
t_unicos_faltantes

Podemos observar de la tabla resultante, que los únicos campos que no presentan valores faltantes son price y property_type. Por otro lado, particularmente la variable bedrooms contiene una altísima proporción de faltantes (más del 40%). El resto de las variables presentan proporciones más aceptables de faltantes como para no considerarlas candidatas a ser eliminadas.

En cuanto a los valores únicos, se puede destacar que la variable categórica property_type posee 3 posibles valores (Casa, Departamento y PH) y variable l3 que corresponder al barrio posee 58 valores distintos. El resto de las variables, que son numéricas, serán analizadas a continuación mediante la matriz de correlación.

Matriz de correlación

# Formato tabla
datos %>% 
  select_if(.,is.numeric) %>% 
  correlate(use="complete.obs") %>% 
  shave() %>% 
  fashion()

# Formato gráfico
datos %>% 
  select_if(.,is.numeric) %>% 
  correlate(use="complete.obs") %>% 
  rplot() + labs(title= 'Gráfico de correlación')

De la matriz resultante se obtiene que existe una importante correlación entre las variables rooms y bedrooms. Esto resulta particularmente útil ya que, como se observó en el paso anterior, la variable bedrooms posee un alto porcentaje de valores faltantes. Así, se puede proceder a eliminar la variable bedrooms del dataset ya que en gran medida es explicada por la variable rooms. Esto es lo que se hará a continuación.

También cabe mencionar que se observa una correlación medianamente alta entre las variables relacionadas a la superficie de la propiedad (surface_total y surface_covered) lo cual tiene sentido dada la lógica y el conocimiento del dominio.

3) Preparación de los datos (II)

A continuación se elimina la variable bedrooms y adicionalmente todas las observaciones que posean datos faltantes.

El dataset resultante quedará con 51210 observaciones y 8 variables.

# Elimino todos los registros con datos faltantes
datos <- datos %>%
  select(.,-bedrooms) %>% 
  drop_na()
datos

4) Análisis exploratorios (II)

Se analiza la variable price utilizando estadísticas descriptivas (cuartiles, promedio, mínimo y máximo) junto con histogramas que muestren la distribución de la misma.

# Estadísticas descriptivas de la variable precio
datos %>% 
  summarise(min=min(price),
            max=max(price),
            q1 = quantile(price,0.25),
            q2 = quantile(price,0.5),
            q3 = quantile(price,0.75),
            promedio = mean(price)
  ) 

# Histograma con escala original
  ggplot(datos,aes(price))+
  geom_histogram(colour = "darkgreen", fill = "lightblue",bins=50) +
  labs(title = 'Distribución de la variable precio')

  
  # Histograma con escala logarítmica
  ggplot(datos,aes(price)) +
  geom_histogram(colour = "darkgreen", fill = "lightblue",bins=50) +
  labs(title = 'Distribución de la variable precio (escala logarítmica)')+
  scale_x_log10()

NA

El histograma de la variable price, en su escala original (valores en USD), muestra claramente una asimetría a derecha. Esto, sumado a los resultados de las estadísticas descriptivas, nos indica que hay algunos valores extremos que podrían interferir en el análisis. Habrá que analizar con algún criterio si dichos valores extremos, sobre todo los mayores, pueden catalogarse como outliers y tratarlos en consecuencia. La media y la mediana no se encuentran muy cercanas, lo cual también es un indicador de posible presencia de outliers.

Adicionalmente al histograma en escala original se agregó uno escalado de manera logarítmica (log10) que permite observar de manera un poco más clara la distribución ya que no muestra valores absolutos sino relativos, en términos de proporciones. Este último histograma queda de una forma más normalizada.

A continuación calculo nuevamente las estadísticas descriptivas sobre la variable price pero agrupando por property_type. También se grafican boxplots con la misma configuración de variables y finalmente un correlograma general utilizando la librería GGAlly.

# Agrupo por property_type
datos %>% 
  group_by(property_type) %>% 
  summarise(min=min(price),
            max=max(price),
            q1 = quantile(price,0.25),
            q2 = quantile(price,0.5),
            q3 = quantile(price,0.75),
            promedio = mean(price)
  )

Boxplots

# Boxplot de precio, escala original
ggplot(datos,aes(x=property_type,y = price,fill=property_type)) +
  geom_boxplot()+
  labs(title='Boxplot de precio por tipo de propiedad')

# Boxplot de precio, escala logarítmica
ggplot(datos,aes(x=property_type,y = price,fill=property_type)) +
  geom_boxplot()+
  labs(title='Boxplot de precio por tipo de propiedad (escala log10)')+
  scale_y_log10()

Correlograma

# Correlograma
datos %>% 
  select(-id,-l3) %>% 
  ggpairs(., 
        title = "Matriz de correlaciones",
        mapping = aes(colour=property_type))

De estos últimos gráficos podemos observar que las casas son en promedio más caras que los otros tipos de propiedades, y que además presentan una mayor variabilidad en términos de precio. Adicionalmente podemos ver en uno de los scatter plots del correlograma, que hay una cierta linealidad entre las variables rooms y price lo cual es un dato útil para pensar en el posterior modelo lineal que se ajustará.

5) Outliers

Genero estadísticas agrupadas por tipo de propiedad, y agregando también una nueva variable calculada de precio por m2, la cual da una indicación un poco más precisa para detectar datos atípicos.

# Estadísticas descriptivas de la variable precio por m2 (pm2)
datos %>% 
  mutate(pm2 = price/surface_total) %>% 
  arrange(pm2) %>% 
  select(property_type,pm2) %>% 
  group_by(property_type) %>% 
  summarise(min=min(pm2),
            max=max(pm2),
            q1 = quantile(pm2,0.25),
            q2 = quantile(pm2,0.5),
            q3 = quantile(pm2,0.75),
            promedio = mean(pm2),
            iqr = IQR(pm2)
            )

NA

Como criterio para filtrar los outliers se tomarán en cuenta:

1. Los límites de la variable pm2, definidos por una magnitud mayor al criterio de Fisher. En este caso se consideraron 2.5 veces la distancia IQR, para tener un mayor margen y no eliminar tantos valores de precios altos que no necesariamente corresponden a outliers.
2. Como el criterio “a” no logra filtrar los valores atípicos muy bajos, se adicionó un filtro que descarta las propiedades cuyo valor pm2 fuera menor a 100 usd.
3. Por último también se consideró un filtro sobre la variable price que elimina 2 observaciones (una de Casa y otra de Departamento) cuyos valores superaban los 5M USD, monto que daba indicios junto a sus otras variables que se consideraban atípicos.

Con todos estos filtros aplicados, los valores resultantes fueron:

paste('dataset con outliers:',nrow(datos))

[1] "dataset con outliers: 51210"

paste('dataset sin outliers:',nrow(datos_clean))

[1] "dataset sin outliers: 50114"

paste('cantidad de outliers eliminados:',(nrow(datos) - nrow(datos_clean)))

[1] "cantidad de outliers eliminados: 1096"

paste('porcentaje de outliers:',percent((nrow(datos) - nrow(datos_clean)) / nrow(datos)))

[1] "porcentaje de outliers: 2.14%"

# Guardo un nuevo dataframe conservando el campo pm2 (precio por metro cuadrado)
datos_new <- datos %>% 
  mutate(pm2 = price/surface_total)
# Genero el dataframe sin outliers
  
datos_clean <- map_dfr(
  .x = unique(datos_new$property_type),
  .f = function(tprop){
    df <- datos_new %>%
      filter(property_type == tprop)
    q1  <- quantile(df$pm2,0.25)
    q3  <- quantile(df$pm2,0.75)
    iqr <- IQR(df$pm2)
    res <- df %>%
      filter( (pm2 >= q1 - 2.5 * iqr) & (pm2 <= q3 + 2.5 * iqr),
              pm2 > 100,
              price < 5000000
              )
    return(res)
  }
)
# Queda como resultado el nuevo dataframe sin las observaciones atípicas (50114 rows)
datos_clean

# 50114

6) Análisis exploratorios (III)

Repito nuevamente el análisis exploratorio luego de quitar los outliers del dataset, comenzando con estadisticas descriptivas para la variable precio, y continuando con el boxplot y correlograma.

# Estadísticas descriptivas de la variable precio
datos_clean %>% 
  summarise(min=min(price),
            max=max(price),
            q1 = quantile(price,0.25),
            q2 = quantile(price,0.5),
            q3 = quantile(price,0.75),
            promedio = mean(price)
  ) 

# Histograma con escala original
  ggplot(datos_clean,aes(price)) +
  geom_histogram(colour = "darkgreen", fill = "lightblue",bins=50) +
  labs(title = 'Distribución de la variable precio (sin outliers)')

# Agrupo por property_type
datos_clean %>% 
  group_by(property_type) %>% 
  summarise(min=min(price),
            max=max(price),
            q1 = quantile(price,0.25),
            q2 = quantile(price,0.5),
            q3 = quantile(price,0.75),
            promedio = mean(price)
  )

Se puede ver que las medias se acercaron a las medianas al quitar algunos valores muy extremos de precios en el paso anterior. La distribución reflejada en el histograma, si bien sigue presentando cola pesada a derecha, se ha normalizado un poco más.

Vuelvo a graficar el boxplot de la variable precio por tipo de propiedad y el correlograma completo.

Boxplot

# Boxplot de precio, escala original
ggplot(datos_clean,aes(x=property_type,y = price,fill=property_type)) +
  geom_boxplot()+
  labs(title='Boxplot de precio por tipo de propiedad (sin outliers)')

Correlograma

datos_clean %>% 
  select(-id,-l3) %>% 
  ggpairs(., 
        title = "Matriz de correlaciones (sin outliers)",
        mapping = aes(colour=property_type))

Se destaca del correlograma generado a partir del dataset sin outliers que, a diferencia del correlograma previo, ahora se visualiza en el scatter plot una relación lineal mucho más clara entre surface_total y price. Recordando que en el gráfico previo encontramos una relación similar entre rooms y price, tiene sentido que en el próximo paso evaluemos ambos modelos de regresión.

Para finalizar el análisis exploratorio agrego estadísticas descriptivas e histograma de la variable precio por metro cuadrado (pm2), en la cual se pueden ver distribuciones de mayor normalidad. También se ve una clara predominancia de la categoría Departamentos en cuanto a cantidad de observaciones, mayor mediana y media que las Casas y PHs. Esto es particularmente interesante de analizar ya que cuando se computaron las estadísticas descriptivas de la variable price se veía que las casas tienen en promedio un valor absoluto mayor que los Departamentos, pero no ocurre lo mismo en el precio por metro cuadrado.

# Agrego además estadísticas de la variable pm2 (precio por metro cuadrado) y muestro su distribución mediante un histograma agrupado por tipo de propiedad.
datos_clean %>% 
  select(property_type,pm2) %>% 
  group_by(property_type) %>% 
  summarise(min=min(pm2),
            max=max(pm2),
            q1 = quantile(pm2,0.25),
            q2 = quantile(pm2,0.5),
            q3 = quantile(pm2,0.75),
            promedio = mean(pm2),
            iqr = IQR(pm2)
            )

  ggplot(datos_clean,aes(x=pm2,color=property_type)) + 
  geom_histogram()+
  labs(title = 'Distribución de la variable precio/m2')+
  facet_grid(. ~property_type)

NA

7) Modelo lineal

El primer modelo lineal que se muestra a continuación intenta explicar el precio en función de las habitaciones (rooms).

# Modelo lineal 1:
  # target: price
  # predictor: rooms
modelo_rooms <- lm(price ~ rooms, data = datos_clean)
summary(modelo_rooms)

Call:
lm(formula = price ~ rooms, data = datos_clean)

Residuals:
     Min       1Q   Median       3Q      Max 
-2198233   -75799   -18967    38033  3224481 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -15977.4     1639.7  -9.744   <2e-16 ***
rooms        87944.1      527.2 166.818   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 163600 on 50112 degrees of freedom
Multiple R-squared:  0.357, Adjusted R-squared:  0.357 
F-statistic: 2.783e+04 on 1 and 50112 DF,  p-value: < 2.2e-16

Los resultados obtenidos mediante el primer modelo arrojan coeficientes estimados de B0 = -15977.4 y B1 = 87944.1, siendo el primero de ellos el valor correspondiente a la ordenada al origen (en este caso particular sería una propiedad con 0 habitaciones, lo cual no tendría sentido); y el segundo de ellos indica que por cada habitación extra, el valor estimado del inmueble aumenta en 87944.1 USD.

ggplot(datos_clean,aes(x=rooms,y=price))+
  geom_point(aes(color=property_type)) +
  geom_smooth(method = "lm") +
  labs(x="Cantidad de habitaciones",y="Precio (usd)",title="Modelo lineal basado en cantidad de habitaciones")

El segundo modelo lineal intenta explicar el precio (rooms) en función de la superficie total (surface_total) (surface_total)

# Modelo lineal 2:
  # target: price
  # predictor: surface_total
modelo_surface_total <- lm(price ~ surface_total, data = datos_clean)
summary(modelo_surface_total)

Call:
lm(formula = price ~ surface_total, data = datos_clean)

Residuals:
     Min       1Q   Median       3Q      Max 
-3048394   -44973   -18310    23255  2077479 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   61141.078    938.550   65.14   <2e-16 ***
surface_total  1868.627      7.819  238.97   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 139500 on 50112 degrees of freedom
Multiple R-squared:  0.5326,    Adjusted R-squared:  0.5326 
F-statistic: 5.711e+04 on 1 and 50112 DF,  p-value: < 2.2e-16

Los resultados obtenidos mediante el segundo modelo arrojan coeficientes estimados de B0 = 61141.078 y B1 = 1868.627, siendo el primero de ellos el valor correspondiente a la ordenada al origen (en este caso particular sería una propiedad con 0 metros cuadrados de superficie, lo cual no tendría sentido); y el segundo de ellos indica que por cada metro cuadrado extra de superficie, el valor estimado del inmueble aumenta en 1868.627 USD.

ggplot(datos_clean,aes(x=surface_total,y=price))+
  geom_point(aes(color=property_type)) +
  geom_smooth(method = "lm") +
  labs(x="Superficie total (m2)",y="Precio (usd)",title="Modelo lineal basado en superficie total")

Habiendo analizado los resultados de ambos modelos y comparándolos mediante la métrica resultante R2, se concluye que el segundo modelo es más preciso a la hora de estimar el precio del inmueble, alcanzando un valor de variabilidad explicada del 53.26% frente al 35.7% del primer modelo.