Overview

Using Airbnb data, we created a model with the aim to analyse data about Airbnb listings and predict the total cost for two people staying four nights in Mexico city.

To create a successful model we had to:

• Apply data-cleaning techniques to extract data from Airbnb list and create a dataset.
• Use descriptive statistics to conduct data modification and relevant analysis. Apply multi-regression modeling to visualize and analyze the data that we have collected.

Our final model takes into consideration the:

• Neighbourhood
• Room type
• N. of bathrooms
• N. of bedrooms

The model predicts that an apartment with a private room in the neighbourhoood of Tlalpan, Mexico City, which has over 15 reviews and an average review score above 95/100 will cost, on average, 2,063 MXN$ with a 95% confidence interval ranging from 2,035 to 2,090.

listings <- vroom::vroom("http://data.insideairbnb.com/mexico/df/mexico-city/2020-06-20/data/listings.csv.gz")

Exploratory Data Analysis (EDA)

tidy_listings <- listings %>% 
dplyr::select(host_is_superhost, neighbourhood_cleansed, latitude, longitude, is_location_exact, property_type, room_type,
              bathrooms, bedrooms, beds, square_feet, price, cleaning_fee, extra_people, number_of_reviews, review_scores_rating,
              cancellation_policy, availability_365, amenities, host_identity_verified, minimum_nights, guests_included)

tidy_listings <- tidy_listings %>% 
dplyr::select(-square_feet)

### Changing characters variables to numerical ones
tidy_listings <- tidy_listings %>%
  mutate(price = parse_number(price), cleaning_fee = parse_number(cleaning_fee), 
       extra_people = parse_number(extra_people))

### Creating neighbourhood_cleased variable
tidy_listings$neighbourhood_cleansed <- as.factor(tidy_listings$neighbourhood_cleansed)

tidy_listings <- tidy_listings %>%
  mutate(neighbourhood_simplified = case_when(
    neighbourhood_cleansed %in% c("Cuauhtémoc", "Venustiano Carranza", "Gustavo A. Madero") ~ "Cuauhtémoc", 
    neighbourhood_cleansed %in% c("Miguel Hidalgo", "Azcapotzalaco") ~ "Miguel Hidalgo",
    neighbourhood_cleansed %in% c("Benito Juárez", "Iztacalco") ~ "Benito Juárez",
    neighbourhood_cleansed %in% c("Coyoacán", "Iztapalapa") ~ "Coyoacán",
    neighbourhood_cleansed %in% c("Álvaro Obregón", "Cuajimalpa de Morelos", "La Magdalena Contreras") ~ "Álvaro Obregón",
    TRUE ~ "Tlalpan"
  ))

### Creating amenities_clean variable 
tidy_listings <- tidy_listings %>% 
  mutate(amenities_clean = count.fields(textConnection(amenities), sep = ","))

### Sorting the data to show the four most common properties types
sort(table(tidy_listings$property_type),decreasing=TRUE)[1:4]

## 
##   Apartment       House Condominium        Loft 
##       13545        3289        1577        1120

### Sorting the data such that the top four property types have their own category and the remaining categories fall into 'other'.
tidy_listings <- tidy_listings %>%
  mutate(prop_type_simplified = case_when(
    property_type %in% c("Apartment","House", "Condominium","Loft") ~ property_type, 
    TRUE ~ "Other"
  ))

### Filtering minimum nights lower than or equal to 4 nights
tidy_listings <- filter(tidy_listings, minimum_nights <= 4)

### Change NA values to 0 in cleaning fees
tidy_listings <- tidy_listings %>%
  mutate(cleaning_fee = case_when(
    is.na(cleaning_fee) ~ 0, 
    TRUE ~ cleaning_fee
  ))

### Create a new variable price_4_nights
tidy_listings <- tidy_listings %>% 
  mutate(price_4_nights = case_when(
    guests_included >= 2 ~ 4*price + cleaning_fee,
    guests_included < 2 ~ 4*price + cleaning_fee + 4*extra_people))

ggplot(tidy_listings, aes(x = review_scores_rating,
                          y = price)) + 
  geom_point()

cor(tidy_listings$review_scores_rating, tidy_listings$price,
    use = "complete.obs")

## [1] 0.00346

# tidy_listings %>%
  #dplyr::filter(!is.na(number_of_reviews),
         #!is.na(bathrooms),
         #!is.na(bedrooms),
         #!is.na(beds),
         #!is.na(review_scores_rating),
         #!is.na(cancellation_policy))

summary(tidy_listings$price)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0     408     747    1509    1358  349990

ggplot(tidy_listings, aes(x=price)) +
  geom_density()

#
tidy_listings <- tidy_listings %>% 
  filter(price <= 3000)

#
ggplot(tidy_listings, aes(x = price_4_nights)) +
  geom_density()

ggplot(tidy_listings, aes(x = price_4_nights)) +
  scale_x_log10() +
  geom_density()

#log_price_4_nights
tidy_listings <- tidy_listings %>% 
  mutate (log_price_4_nights = log(price_4_nights))

ggpairs1 <- tidy_listings %>%
 dplyr::select( price, cleaning_fee, extra_people, number_of_reviews,
                guests_included, room_type, host_is_superhost, prop_type_simplified) %>%
   rename("Price" = price,
          "Cleaning fee" = cleaning_fee,
          "Extra people" = extra_people,
          "N. Reviews" = number_of_reviews,
          "Guests included" = guests_included,
          "Room type" = room_type,
          "Superhost" = host_is_superhost) %>% 
  GGally::ggpairs() +
  theme_minimal() 

ggpairs1

ggplot(tidy_listings, aes(x = reorder(neighbourhood_simplified, price),
                          y = price, 
                          fill = neighbourhood_simplified)) +
  geom_boxplot(alpha = 0.6) +
  labs(title = "Central Mexico City is for the rich",
       subtitle = "Boxplot of prices across selected neighbourhoods",
       caption = "Source: insideairbnb.com",
       x= "Neighbourhood" ,
       y = "Price per night (MXN $)") +
  theme_bw() +
  scale_fill_brewer(palette = "Set3") +
  scale_fill_discrete(name = "Neighbourhood") +
  theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())

prop_plot <- tidy_listings  %>%
  filter(price <= 2000)

ggplot(prop_plot, aes(reorder(x = prop_type_simplified, price),
                          y = price,
                          fill = prop_type_simplified)) +
  geom_boxplot(alpha = 0.6) +
  labs(title = "Everyone dreams about a loft in Mexico City",
       subtitle = "Boxplot of prices across selected types of properties",
       caption = "Source: insideairbnb.com",
       x= "Property type" ,
       y = "Price per night (MXN $)") +
  theme_bw() +
  scale_fill_discrete(name = "Property type") +
  theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())

Bar_chart <- tidy_listings %>%
   mutate('room_type' = factor(`room_type`, levels = c("Hotel room",
                                                       "Shared room",
                                                       "Entire home/apt",
                                                       "Private room")))

ggplot(Bar_chart, aes(x= room_type,
                      fill = room_type)) +
  geom_bar(alpha = 0.6) +
  theme_bw() +
  labs(title = "Couples seek the right intimacy in Mexico City",
       subtitle = "Total number of rooms by most common types" ,
       caption = "Source: insideairbnb.com",
       x= "Room type",
       y = "Count" ) +
  scale_fill_discrete(name = "Room type") +
  theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())

ggplot(tidy_listings, aes(x = host_is_superhost,
                          y = price,
                          fill = host_is_superhost)) +
  geom_boxplot(alpha = 0.6) +
  labs(title = "Superhost can be super... expensive",
       subtitle = "Boxplot of prices for superhosts and normal hosts",
       caption = "Source: insideairbnb.com",
       x= "Superhost",
       y = "Price per night (MXN $)") +
   theme_bw() +
  theme(legend.position = "none") +
  scale_x_discrete(labels = c("No", "Yes"))

ggplot(tidy_listings, aes(x= number_of_reviews,
                          y= price)) +
  geom_point() + 
  labs(x = "Number of reviews (#)",
       y = "Price per night (MXN $)") +
  theme_bw()

ggplot(tidy_listings, aes(x = number_of_reviews)) +
  geom_density() + 
  labs(x = "Number of reviews (#)",
       y = "Density") +
  theme_bw()

cor(tidy_listings$number_of_reviews, tidy_listings$price,
    use = "complete.obs")

## [1] 0.103

tidy_review <- tidy_listings %>% 
  filter(number_of_reviews <= 20)
ggplot(tidy_review, aes(x= number_of_reviews,
                        y= price)) +
  geom_point() + 
  labs(x = "Number of reviews (#)",
       y = "Price per night (MXN $)") +
  theme_bw()

cor(tidy_review$number_of_reviews, tidy_review$price,
    use = "complete.obs")

## [1] 0.0354

cleaning_hist <- tidy_listings %>%  filter(cleaning_fee <= 1000) 


ggplot(cleaning_hist, aes(x= cleaning_fee)) +
  geom_histogram(binwidth = 100) +
  theme_bw() +
  labs(title = "Most hosts don't charge extra for cleaning",
       subtitle = "Distribution of cleaning fee charge",
       caption = "Source: insideairbnb.com",
       x= "Cleaning fee (MXN $)",
       y = "Count") +
  scale_fill_discrete()

ggplot(tidy_listings, aes(x= guests_included,
                          y= price),
       fill = guests_included) +
  geom_bar(stat="identity") +
  theme_bw() +
  labs(title = "Sometimes not everyone is welcome",
       subtitle = "Number of guests included in the abitation",
       caption = "Source: insideairbnb.com",
       x= "Guests included (#)",
       y = "Price per night (MXN $)" )

ggplot(tidy_listings, aes(x = guests_included,
                          y = price)) +
  geom_point() +
  labs(x = "Guests included (#)",
       y = "Price per night (MXN $)") +
  theme_bw()

cor(tidy_listings$guests_included, tidy_listings$price,
    use = "complete.obs")

## [1] 0.267

ggplot(tidy_listings, aes(x = amenities_clean,
                            y = price)) +
  geom_point(aes(colour = host_is_superhost)) +
  theme_minimal() +
  scale_color_manual(values = c("darkorange", "cyan4"),
                     na.translate = FALSE) +
  labs(title = "Superhosts provide greater amenities",
       subtitle = "The effect of amenities on price across selected neighbourhoods",
       x = "Number of amenities (#)",
       y = "Price (MXN $)",
       caption = "Source: insideairbnb.com",
       colour = "Superhost") +
  theme(legend.position = "bottom",
        legend.background = element_rect(fill = "white", color = NA),
        plot.title.position = "plot",
        plot.caption = element_text(hjust = 0, face= "italic"),
        plot.caption.position = "plot") +
  scale_colour_discrete(labels = c("No", "Yes")) +
  facet_wrap(~neighbourhood_simplified)

Mapping

Another fun part: displaying the housing offer in Mexico City on a map!

leaflet(data = filter(listings, minimum_nights <= 4)) %>% 
  addProviderTiles("OpenStreetMap.Mapnik") %>% 
  addCircleMarkers(lng = ~longitude, 
                   lat = ~latitude, 
                   radius = 1, 
                   fillColor = "blue", 
                   fillOpacity = 0.4, 
                   popup = ~listing_url,
                   label = ~property_type)

Regression

Let us finally get now more analytical: we fit our first regression model using twelve different independent variables which resulted in a R 2 of 0.471 . However, this model had too many variables , some of which were not that significant, so we removed the one with the largest p value which is cancellation policy.

tidy_listings$cancellation_policy <- as.factor(tidy_listings$cancellation_policy)
 
final_1 <- lm(log_price_4_nights ~ guests_included + amenities_clean + host_is_superhost + neighbourhood_simplified + room_type + prop_type_simplified + bathrooms + bedrooms + beds + is_location_exact + cancellation_policy + host_identity_verified, data = tidy_listings)

huxreg(final_1)

car::vif(final_1)

##                          GVIF Df GVIF^(1/(2*Df))
## guests_included          1.24  1            1.11
## amenities_clean          1.38  1            1.17
## host_is_superhost        1.18  1            1.09
## neighbourhood_simplified 1.19  5            1.02
## room_type                1.81  3            1.10
## prop_type_simplified     1.54  4            1.06
## bathrooms                2.45  1            1.57
## bedrooms                 3.12  1            1.77
## beds                     2.12  1            1.46
## is_location_exact        1.02  1            1.01
## cancellation_policy      1.17  4            1.02
## host_identity_verified   1.06  1            1.03

Seen from the statistics, R2 of the model is great with a value of 0.423 but this model has too many variables and some of which are not that significant so we will delete the one with the largest p value which is cancellation policy.

final_2 <- lm(log_price_4_nights ~ guests_included + amenities_clean + host_is_superhost + neighbourhood_simplified + room_type + prop_type_simplified + bathrooms + bedrooms + beds + is_location_exact + host_identity_verified, data = tidy_listings)

huxreg(final_2)

car::vif(final_2)

##                          GVIF Df GVIF^(1/(2*Df))
## guests_included          1.23  1            1.11
## amenities_clean          1.35  1            1.16
## host_is_superhost        1.17  1            1.08
## neighbourhood_simplified 1.18  5            1.02
## room_type                1.77  3            1.10
## prop_type_simplified     1.54  4            1.06
## bathrooms                2.45  1            1.56
## bedrooms                 3.12  1            1.77
## beds                     2.12  1            1.46
## is_location_exact        1.02  1            1.01
## host_identity_verified   1.05  1            1.02

We then fit the model with other 11 variables and this model gives us a R2 of 0.422 which is pretty close to before. Thus, we are pretty confident that cancellation policy does not have a significant effect on the price_4_night. Therefore, we continue to delete variables guests_included, amenities_clean, and beds which has largest p value.

We ran the model again with 8 variables but realised we had to further eliminate some independent variables to bring the AIC down.

final_3 <- lm(log_price_4_nights ~ host_is_superhost + neighbourhood_simplified + room_type + prop_type_simplified + bathrooms + bedrooms + is_location_exact + host_identity_verified, data = tidy_listings)

huxreg(final_3)

car::vif(final_3)

##                          GVIF Df GVIF^(1/(2*Df))
## host_is_superhost        1.05  1            1.02
## neighbourhood_simplified 1.17  5            1.02
## room_type                1.53  3            1.07
## prop_type_simplified     1.51  4            1.05
## bathrooms                2.28  1            1.51
## bedrooms                 2.42  1            1.56
## is_location_exact        1.02  1            1.01
## host_identity_verified   1.03  1            1.01

We once again fit the model with other 8 variables and this model gives us a R2 of 0.407 and as we can tell, there are still too many variables in the model.

final_4 <- lm(log_price_4_nights ~ neighbourhood_simplified + room_type + bathrooms + bedrooms, data = tidy_listings)

huxreg(final_4)

car::vif(final_4)

##                          GVIF Df GVIF^(1/(2*Df))
## neighbourhood_simplified 1.04  5            1.00
## room_type                1.13  3            1.02
## bathrooms                2.25  1            1.50
## bedrooms                 2.37  1            1.54

autoplot(final_4)

summary(final_4)

## 
## Call:
## lm(formula = log_price_4_nights ~ neighbourhood_simplified + 
##     room_type + bathrooms + bedrooms, data = tidy_listings)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.690 -0.334 -0.014  0.323  2.539 
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                             8.30001    0.01591  521.58  < 2e-16 ***
## neighbourhood_simplifiedBenito Juárez  -0.16428    0.01655   -9.93  < 2e-16 ***
## neighbourhood_simplifiedCoyoacán       -0.18016    0.01832   -9.83  < 2e-16 ***
## neighbourhood_simplifiedCuauhtémoc      0.03013    0.01507    2.00    0.046 *  
## neighbourhood_simplifiedMiguel Hidalgo  0.16214    0.01703    9.52  < 2e-16 ***
## neighbourhood_simplifiedTlalpan        -0.33312    0.02068  -16.11  < 2e-16 ***
## room_typeHotel room                    -0.13207    0.03391   -3.89  9.9e-05 ***
## room_typePrivate room                  -0.68779    0.00820  -83.88  < 2e-16 ***
## room_typeShared room                   -0.90030    0.02868  -31.39  < 2e-16 ***
## bathrooms                               0.04560    0.00605    7.53  5.2e-14 ***
## bedrooms                                0.05559    0.00549   10.12  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.514 on 18058 degrees of freedom
##   (44 observations deleted due to missingness)
## Multiple R-squared:  0.402,  Adjusted R-squared:  0.401 
## F-statistic: 1.21e+03 on 10 and 18058 DF,  p-value: <2e-16

For our final model, we take out is_location_exact, host_is_superhost,and host_identity-verified and the R2 is 0.402.

target_listings <- tidy_listings %>%
dplyr::filter(number_of_reviews > 15,
         review_scores_rating > 95,
         room_type == "Private room",
         prop_type_simplified == "Apartment",
         neighbourhood_simplified == "Tlalpan")

We are essentially facing a clear trade off whilst building our model and had to weigh the benefits of additional variables (higher R 2) against the amount of error in our model ( higher AIC ).

ci <- predict(final_4, newdata = target_listings, interval = "confidence")

ci = exp(ci)

summary(ci)

##       fit            lwr            upr      
##  Min.   :1517   Min.   :1469   Min.   :1567  
##  1st Qu.:1604   1st Qu.:1556   1st Qu.:1654  
##  Median :1604   Median :1556   Median :1654  
##  Mean   :1606   Mean   :1557   Mean   :1656  
##  3rd Qu.:1604   3rd Qu.:1556   3rd Qu.:1654  
##  Max.   :1679   Max.   :1626   Max.   :1734

For prediction and confidence interval, we get a median of 2063 as our lower bound, 2035 as our prediction, and 2090 as our higher bound.

model1 <- lm(price_4_nights ~ prop_type_simplified + number_of_reviews + review_scores_rating, data = tidy_listings)
model1

## 
## Call:
## lm(formula = price_4_nights ~ prop_type_simplified + number_of_reviews + 
##     review_scores_rating, data = tidy_listings)
## 
## Coefficients:
##                     (Intercept)  prop_type_simplifiedCondominium  
##                         3003.87                           240.01  
##       prop_type_simplifiedHouse         prop_type_simplifiedLoft  
##                        -1318.25                           178.82  
##       prop_type_simplifiedOther                number_of_reviews  
##                         -422.57                             4.53  
##            review_scores_rating  
##                           10.55

model2 <- lm(price_4_nights ~ prop_type_simplified + number_of_reviews + review_scores_rating + room_type, data = tidy_listings)
model2

## 
## Call:
## lm(formula = price_4_nights ~ prop_type_simplified + number_of_reviews + 
##     review_scores_rating + room_type, data = tidy_listings)
## 
## Coefficients:
##                     (Intercept)  prop_type_simplifiedCondominium  
##                        4367.016                          222.816  
##       prop_type_simplifiedHouse         prop_type_simplifiedLoft  
##                          21.968                         -522.028  
##       prop_type_simplifiedOther                number_of_reviews  
##                          55.893                           -0.375  
##            review_scores_rating              room_typeHotel room  
##                           9.637                         -565.933  
##           room_typePrivate room             room_typeShared room  
##                       -2855.795                        -3492.220

For each increase in review_scores_rating, the price of two people for 4 nights will increase by 10.55.

Furthermore, we notice the following change in prices according to the property type:

• if the property type is condominium, the price of two people for 4 nights will increase by 240.0
• if the property type is house, the price of two people for 4 nights will decrease by 1318.3
• if the property type is loft, the price of two people for 4 nights will increase by 178.8
• if the property type is other, the price of two people for 4 nights will decrease by 422.6

Outcomes and concerns

• To check that our results were valid we searched Airbnb to find out the rough cost of two people staying four nights in Tlalpan (Mexico City) the results showed an average price of 2,100 MXN$ which is comparable our findings .

• Our model is rather limited in its use as it focuses on Mexico City alone to make this a general model we would need Airbnb data from around the world.