rmarkdown_seattle.Rmd

---
title: "Airbnb"
author: "Sanjana Gupta"
date: "30/01/2020"
output:
  html_document:
    df_print: paged
  pdf_document: default
always_allow_html: yes
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

#### GitHub Repo - https://github.com/sanjanagupta16/Data-Science-Process

Importing Libraries and setting cwd
```{r libraries}
library(tidyverse)
library(plotly)
library(ggplot2)
library(naniar)
library(leaflet)
library(htmltools)
library(choroplethr)
library(dplyr)
library(choroplethrMaps)
library(ggplot2)
library(gridExtra)
library(ggmap)
library(htmlwidgets)
library(mapview)
library(dplyr)
library(plotly)
library(devtools)
library(kableExtra)
library(here)
library(janitor)
library(Hmisc)
library(corrplot)
library(jtools)
library(caret)
#setwd("~/Desktop/Winter 2020/Data 598C Data Science Process/Data Science Project")
#setwd(here("/Data Science Project"))
```

Importing data and creating a training dataset
```{r import-data}
set.seed(891)
data1 <- read.csv(file = "data/listings.csv", header = TRUE)
data2 <- read.csv(file = "data/calendar.csv")
data3 <- read.csv(file = "data/reviews.csv")

data_train <- data1
```

Viewing the data in different formats
```{r view-data}
#View(data_train)

#to view data kinda nicely - uses dplyr
#glimpse(data_train)
#str(data1)
#summary(data_train)
#head(data_train, 1000)
#tail(data_train, 50)
```

## Data Understanding & Data Cleaning

1. Data Cleaning - removing columns
```{r data-clean}
set.seed(891)
#cols_to_delete <- c("last_scraped", "thumbnail_url", "host_picture_url", "medium_url", "picture_url", "xl_picture_url", "host_thumbnail_url", "host_picture_url")
#select(data_train, -contains("url"))
data_train$last_scraped <- NULL
data_train$thumbnail_url <- NULL
data_train$host_picture_url <- NULL
data_train$medium_url <- NULL
data_train$picture_url <- NULL
data_train$xl_picture_url <- NULL
data_train$host_thumbnail_url <- NULL
data_train$host_picture_url <- NULL
data_train$host_url <- NULL
data_train$scrape_id <- NULL
data_train$experiences_offered <- NULL
data_train$neighborhood_overview <- NULL
data_train$host_about <- NULL
data_train$host_id <- NULL
data_train$host_verifications <- NULL
data_train$host_has_profile_pic <- NULL
data_train$host_identity_verified <- NULL
data_train$calendar_last_scraped <- NULL
data_train$license <- NULL
data_train$require_guest_profile_picture <- NULL
data_train$require_guest_phone_verification <- NULL
data_train$summary <- NULL
data_train$notes <- NULL
data_train$space <- NULL
data_train$description <- NULL
data_train$transit <- NULL
data_train$host_neighbourhood <- NULL
data_train$market <- NULL
data_train$country_code <- NULL
data_train$is_location_exact <- NULL

rownames(data_train) <- NULL

#Viewing the data
#glimpse(data_train)
#summary(data_train$is_location_exact)
```

2. Converting Data type for price, cleaning fee, extra_people, security deposit to numeric
```{r convert-num}
set.seed(891)
data_train$price <- sub("\\$","",data_train$price)
data_train$price <- sub(",","",data_train$price)
data_train$price <- as.integer(data_train$price)

data_train$cleaning_fee <- sub("\\$","",data_train$cleaning_fee)
data_train$cleaning_fee <- sub(",","",data_train$cleaning_fee)
data_train$cleaning_fee <- as.integer(data_train$cleaning_fee)

data_train$extra_people <- sub("\\$","",data_train$extra_people)
data_train$extra_people <- sub(",","",data_train$extra_people)
data_train$extra_people <- as.integer(data_train$extra_people)

data_train$security_deposit <- sub("\\$","",data_train$security_deposit)
data_train$security_deposit <- sub(",","",data_train$security_deposit)
data_train$security_deposit <- as.integer(data_train$security_deposit)
```

3. Changing data type of categorical variables to factor
```{r convert-factor}
set.seed(891)
data_train$host_response_time <- as.factor(data_train$host_response_time)
data_train$host_is_superhost <- as.factor(data_train$host_is_superhost)
data_train$neighbourhood_cleansed <- as.factor(data_train$neighbourhood_cleansed)
data_train$property_type <- as.factor(data_train$property_type)
data_train$room_type <- as.factor(data_train$room_type)
data_train$bed_type <- as.factor(data_train$bed_type)
data_train$calendar_updated <- as.factor(data_train$calendar_updated)
data_train$cancellation_policy <- as.factor(data_train$cancellation_policy)
data_train$instant_bookable <- as.factor(data_train$instant_bookable)
data_train$bed_type <- as.factor(data_train$bed_type)
```

4. Changing host response time, acceptance rate from factor to numeric
```{r convert-num2}
set.seed(891)
data_train$host_response_rate<- as.numeric(gsub("%", "", as.character(data_train$host_response_rate)))/100
data_train$host_acceptance_rate <- as.numeric(gsub("%", "", as.character(data_train$host_acceptance_rate)))/100
```

5. For Property Type: Replacing the categories with only a few listings to "Other" 
```{r}
set.seed(891)
#summary(data_train$property_type)
#data_train$property_type <- replace.value(data_train$property_type, c("Boat", "Bungalow", "Cabin", "Tent","Treehouse", "Yurt"),"Other")
data_train$property_type[data_train$property_type == "Boat"] <- "Other"
data_train$property_type[data_train$property_type == "Bungalow"] <- "Other"
data_train$property_type[data_train$property_type == "Dorm"] <- "Other"
data_train$property_type[data_train$property_type == "Chalet"] <- "Other"
data_train$property_type[data_train$property_type == "Treehouse"] <- "Other"
data_train$property_type[data_train$property_type == "Yurt"] <- "Other"
data_train$property_type[data_train$property_type == "Camper/RV"] <- "Other"

summary(data_train$property_type)
```


6. Checking if Price = 0 for any listings
```{r check-price}
set.seed(891)
0 %in% data_train$price
```
All Airbnb listings have a price associated with it

7. Removing rows with NA
```{r remove-NA}
#is.null(data_train$zipcode)
data_train %>% filter_all(all_vars(!is.na(.)))
data_train %>% filter_all(all_vars(complete.cases(.)))  
```
8. Removing rows with null values
```{r remove-null}
data_train %>% filter_all(all_vars(!is.null(.)))
data_train %>% filter_all(all_vars(complete.cases(.)))  
```

9. Getting a summary of all fields
```{r summary}
summary_df <- data.frame(summary(data_train))
#write.csv(summary_df, "summary.csv")
```

10. Adding log price column and log accommodates column --- Transformations
```{r log-Price}
set.seed(891)
data_train$log_price <- log(data_train$price)
data_train$log_acc <- log(data_train$accommodates)

ggplot(data_train, aes(price)) + geom_bar() + theme_bw()
ggplot(data_train, aes(log_price)) + geom_bar() + theme_bw()

ggplot(data_train, aes(accommodates)) + geom_bar() + theme_bw()
ggplot(data_train, aes(log_acc)) + geom_bar() + theme_bw()

```


11. Checking for missing values and plotting the data
```{r missing-data}
set.seed(891)
check_missing<- function(x){
  if (is.character(x)) sum(x=="") else sum(is.na(x))
}
NMISS<-data.frame(nmiss=sapply(data_train,check_missing))
#write.csv(NMISS, "missing.csv")

#plotting missing data
missingdata <- data_train
missingdata[missingdata == ""] <- NA
missingdata <- missingdata %>% select(host_is_superhost, review_scores_rating,
                                      host_response_time, name, zipcode, latitude,longitude,
                                      host_location, host_response_rate,neighbourhood_cleansed,
                                      property_type, price,weekly_price, monthly_price, bedrooms, bathrooms)
gg_miss_var(missingdata) + labs(title = "Missing Values Plot",y = "Missing Data", x = "Attribute")  + theme_bw() + theme(text=element_text(color = "dark blue")) 
```

12. Replacing missing values
```{r replace-missing}
data_train$review_scores_rating <- ifelse(is.na(data_train$review_scores_rating)==T,97,data_train$review_scores_rating)
data_train$host_response_rate <- ifelse(is.na(data_train$host_response_rate)==T,97,data_train$host_response_rate)

#Using same plot as before to check the values
gg_miss_var(missingdata) + labs(title = "Missing Values Plot",y = "Missing Data", x = "Attribute")  + theme_bw() + theme(text=element_text(color = "dark blue")) 
```

## Data Understanding - plots

1. Growth of Airbnb - evolution of hosts over time
```{r}
set.seed(891)
#Removing NA 
new_hosts <- data_train %>% drop_na("host_since")
#str(new_hosts$host_since)

# Calculate the number of new hosts for each year 
new_hosts$host_since <- as.Date(new_hosts$host_since, '%Y-%m-%d')
new_hosts <- new_hosts[new_hosts$host_since < as.Date("2018-01-01"),]
new_hosts <- new_hosts[order(as.Date(new_hosts$host_since, format="%Y-%m-%d")),]
new_hosts$host_since <- format(as.Date(new_hosts$host_since, "%Y-%m-%d"), format="%Y-%m")
new_hosts_table <- table(new_hosts$host_since)

# Plot
plot(as.Date(paste(format(names(new_hosts_table), format="%Y-%m"),"-01", sep="")), as.vector(new_hosts_table), type = "l", xlab = "Time", ylab = "Number of new hosts", col = "Blue") 
axis(1, seq(0,2018, 100))
```
The number of new hosts have been steadily increasing since 2008, however, in the last couple of years there has been a decrease in this number.

2. Analyzing neighborhood grouping
```{r neighborhood}
set.seed(891)

#In neighborhood group cleansed, there are 794 rows with "Other Neighborhoods"
#Need to fix this - maybe look at neighorhood cleansed column instead?
#Looking at categorical data now
sorted <- sort(table(data_train$host_neighbourhood), decreasing = TRUE)
#head(subset(data_train, select = 'neighbourhood_group_cleansed'))
fact <- factor(data_train$neighbourhood_group_cleansed)
#count(data_train, 'neighbourhood_group_cleansed')
#count(data_train, 'neighbourhood_cleansed')

#Getting a count of listings per neighborhood
#write.csv(neigborhood_data, "neighborhoodcount.csv")
#count(data_train$neighbourhood_cleansed)
plot(data_train$neighbourhood_cleansed, xlab = 'Neighborhoods', ylab = 'Count')
#ggplot(data_train, aes(neighbourhood_cleansed)) + geom_bar() + theme_bw() 
neighborhood <- (data_train$neighbourhood_cleansed)
#x <- sort(count((neighborhood)))
#arrange(data_train, neighbourhood_cleansed)
neighborhood_data <- data_train[order(data_train[,21] ),]
#summary(neighborhood_data$neighbourhood_cleansed)

data_train %>% 
  tabyl(neighbourhood_cleansed) %>% 
  adorn_totals("row") %>% 
  adorn_pct_formatting() %>% 
  nrow()


data_train %>% 
  tabyl(neighbourhood_cleansed) %>% 
  adorn_pct_formatting() %>% 
  arrange(-n) %>% 
  filter(n > 10) %>% 
  adorn_totals("row") %>% 
  head()
```
There's a total of 88 neighborhoods.

```{r plot1}
set.seed(891)
#Trying to fix "Other Neighborhoods" value in neighborhood group cleaned column
#for(i in data_train)
 # if(i[neighbourhood_group_cleansed] == "Other neighborhoods")
  #  i[neighbourhood_group_cleansed] <-  i[neighborhood_cleansed]

#Plotting count of listings per neighborhood
neigborhood_data <- table(data_train$neighbourhood_cleansed)

barplot(sort(table(data_train$neighbourhood_group_cleansed), decreasing = TRUE), ylab = "Count",horiz = FALSE , main = "Distribution of Listings Per Neighborhood", beside = FALSE)
```

```{r plot2}
set.seed(891)
#barplot(table(data_train$neighbourhood_group_cleansed), density = 20, ylab = "Count",horiz = FALSE , main = "Distribution of Listings Per Neighborhood", beside = FALSE) + theme(axis.text.x = element_text(angle = 90, hjust = 1))
ggplot(data_train, aes(neighbourhood_group_cleansed)) + theme_bw()+ geom_bar(na.rm = TRUE, stat = "count") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) 
```

3. Looking at some more plots - Count of Host Listings
```{r count-listings}
set.seed(891)
#barplot(data_train$host_listings_count, main = "Count of Host Listings", xlab = "", ylab = "Count") 
ggplot(data_train, aes(host_listings_count)) + geom_bar(na.rm = TRUE) + xlim(0,20) + theme_bw()
```
Most hosts have a single Airbnb listing.

4. Looking at Types of Airbnb properties
```{r property_type}
set.seed(891)
#barchart(data_train$property_type, main = "Types of Airbnb Properties", xlab = "Count")
prop_type <- table(data_train$property_type)

ggplot(data_train, aes(property_type)) + geom_bar() + theme_bw()

ggplot(data_train, aes(property_type)) + geom_bar(na.rm = TRUE, stat = "count") + theme_bw() + theme(axis.text.x = element_text(angle = 90, hjust = 1))
#plot(data_train$host_is_superhost)
```

5. Prices per Neighborhood
```{r price-neighborhood}
set.seed(891)
z = table(data_train$room_type,data_train$neighbourhood_group)
kable(z)  %>%
  kable_styling(latex_options = "scale_down")
```


```{r boxplot1}
set.seed(891)
data_train %>% boxplot(price ~ neighbourhood_group_cleansed,data = ., main="Box Plot of price vs neighborhood", 
                     ylab="neighborhood", xlab="Price",horizontal=TRUE)
```

6. Number of listings per neighborhood
```{r listings}
set.seed(891)
ggplot(data_train, aes(x = fct_infreq(neighbourhood_cleansed), fill = room_type)) +
    geom_bar() +
    labs(title = "No. of listings by neighborhood",
         x = "Neighborhood", y = "No. of listings") +
    theme(legend.position = "bottom") + theme(axis.text.x = element_text(angle = 90, hjust = 1))
```
Looks messy, want to consider looking at only the top 10 neighborhoods

7. Top 10 neighborhoods - doesn't work - need to fix this
```{r top10}
set.seed(891)
#top_10 <- data_train %>% group_by(neighborhood_cleansed) %>% tally()
#top_10 <- data_train %>% top_n(n = 10, wt = neighbourhood_group_cleansed) %>% arrange(neighbourhood_group_cleansed) %>% #summarize(neighbourhood_group_cleansed = n())

top_10 <- data_train[data_train$neighbourhood_group_cleansed %in% names(sort(table(data_train$neighbourhood_group_cleansed), decreasing = TRUE)[0:1]),]
#top_10 <- data_train %>% top_n(10) %>% group_by(neighbourhood_group_cleansed) %>% tally() 
top_10.val <- table(top_10$neighbourhood_group_cleansed)
kable(top_10.val)
```

7. Looking at outliers for price
```{r outliers}
set.seed(891)
ggplot(data=data_train, aes(price)) + 
  geom_histogram(fill="red") + 
  labs(title="Histogram of Price") +
  labs(x="Price", y="Count")

# Percentile of price
quantile(data_train$price, c(.9, .95, .97, 0.975, 0.98, 0.99, 0.995, 0.999, 0.9999))
boxplot(data_train$price)
```
8. Room Type Analysis
# For a private room
```{r private-room}
set.seed(891)
Private_rooms <- data_train %>% filter(room_type == "Private room")
Private_rooms$price %>% summary()
```

# For an entire apartment

```{r entire-apt}
set.seed(891)
entire_house <- data_train %>% filter(room_type == "Entire home/apt")
entire_house$price %>% summary()
```

#Shared room
```{r shared-room}
set.seed(891)
shared_room <- data_train %>% filter(room_type == "Shared room")
shared_room$price %>% summary()
```

9. Average prices for each room type
```{r avg-price}
set.seed(891)
listingdf <- read.csv('data/listings.csv')
zipReviews <- listingdf %>% group_by(zipcode = zipcode) %>% summarise(avg_loc_review = mean(review_scores_location, na.rm = TRUE))

average_price <- data_train %>% group_by(room_type) %>% summarise(price = mean(price))
colnames(average_price)[1] <- "Room Type"
colnames(average_price)[2] <- "Average Price"
kable((average_price)) %>% kable_styling(latex_options = c("scale_down","condensed"))
```

10. Top 25 most expensive listings
```{r patial Data 1.3, echo=FALSE, message=FALSE, warning=FALSE}
set.seed(891)
top_df <- data_train %>% top_n(n = 25, wt = log_price)

# get background map
top_height <- max(top_df$latitude) - min(top_df$latitude)
top_width <- max(top_df$longitude) - min(top_df$longitude)
top_borders <- c(bottom  = min(top_df$latitude)  - 0.1 * top_height,
                 top     = max(top_df$latitude)  + 0.1 * top_height,
                 left    = min(top_df$longitude) - 0.1 * top_width,
                 right   = max(top_df$longitude) + 0.1 * top_width)

top_map <- get_stamenmap(top_borders, zoom = 12, maptype = "toner-lite")

# map of top 50 most expensive listings
ggmap(top_map) +
    geom_point(data = top_df, mapping = aes(x = longitude, y = latitude,
                                        col = price)) +
    scale_color_gradient(low = "blue", high = "red")
```

## Model Selection

#### Clustering
Used Clustering to compare which areas of Seattle have a higher concentration of Airbnb listings.

Plotting the listings on the map of Seattle
```{r clustering}
set.seed(891)

#Reading Listings Data
listingdf <- read.csv('data/listings.csv')
#Creating a link for each pop up

#Creating Listings across Seattle
map <- leaflet(listingdf) %>%
  addTiles() %>%
  addMarkers(~longitude, ~latitude,labelOptions = labelOptions(noHide = T),
             clusterOptions = markerClusterOptions(),
             popup = paste0("<b> Name: </b>", listingdf$name , "<br/><b> Host Name: </b>", 
                            listingdf$host_name, "<br> <b> Price/night: </b>", listingdf$price, "<br/><b> Room Type: </b>", 
                            listingdf$room_type, "<br/><b> Property Type: </b>", listingdf$property_type
             )) %>% 
  setView(-122.335167, 47.608013, zoom = 11) %>%
  addProviderTiles("CartoDB.Positron")
map
mapshot(map, url = paste0(getwd(),"/map.html"))
```
This shows listings in a clustered fashion. Clicking on a cluster zooms into the listings belonging to it.
Clicking on a pin gives you a pop up with some details about the property.
We can see that most of the listings are clustered around Downtown and Cap Hill.

#### Linear Regression

1. Which variables influence price?
Predictor variables selected - bathrooms, price, property type, accommodates, bedrooms, room type, host is superhost, cancellation policy, availability_365, cleaning_fee, extra_people, security_deposit, instant_bookable, maximum_nights, minimum_nights, bed_type, latitude, longitude, number_of_reviews
```{r LR}
set.seed(891)
LM_model1 <- summary(lm(formula = (price) ~ room_type + accommodates + bedrooms + bathrooms + property_type + host_is_superhost + cancellation_policy + availability_365 +cleaning_fee + extra_people + security_deposit + instant_bookable + maximum_nights + minimum_nights + bed_type + latitude + longitude + number_of_reviews, data = data_train))

plot(price ~room_type + accommodates + bedrooms + bathrooms + property_type + factor(neighbourhood_cleansed) + host_is_superhost, data = data_train) + abline(lm(formula = price ~ room_type + accommodates + bedrooms + bathrooms, data = data_train))
LM_model1
qqnorm(LM_model1$residuals, main = "Normal qq plot of residuals")

LM1 <- lm(lm(formula = (price) ~ room_type + accommodates + bedrooms + bathrooms + property_type + host_is_superhost + cancellation_policy + availability_365 +cleaning_fee + extra_people + security_deposit + instant_bookable + maximum_nights + minimum_nights + bed_type + latitude + longitude + number_of_reviews, data = data_train))
plot(LM1)

summ(LM1)

#Calculate RMSE
RMSE <- sqrt(mean(residuals(LM1)^2))
cat("RMSE = ", RMSE)

```
Results:

The following features are significant to price:
- Private Room 
- Shared Room
- Accommodates
- Bedrooms
- Bathrooms
- Host is Superhost
- Cleaning Fee
- Security Deposit
- Instant bookable
- Location (Lat/Long)
- Number of reviews

R-squared is 61.99%, which is not bad.

RMSE = 63.25.

2. Taking log(price) with multiple linear regression:
```{r lR-with-logprice}
set.seed(891)
LM_model2 <- summary(lm(formula = log_price ~ room_type + accommodates + bedrooms + bathrooms + property_type + host_is_superhost +cancellation_policy + availability_365 + cleaning_fee + extra_people + security_deposit + instant_bookable + maximum_nights + minimum_nights + bed_type + latitude + longitude + number_of_reviews, data = data_train))

plot(price ~room_type + accommodates + bedrooms + bathrooms + factor(neighbourhood_cleansed) + host_is_superhost + cancellation_policy + availability_365 +  cleaning_fee + extra_people + security_deposit + instant_bookable + maximum_nights + minimum_nights + bed_type + latitude + longitude + number_of_reviews, data = data_train) + abline(lm(formula = log_price ~ room_type + accommodates + bedrooms + bathrooms, data = data_train))
LM_model2

LM2 <- lm(formula = log_price ~ room_type + accommodates + bedrooms + bathrooms + property_type + host_is_superhost +cancellation_policy + availability_365 + cleaning_fee + extra_people + security_deposit + instant_bookable + maximum_nights + minimum_nights + bed_type + latitude + longitude + number_of_reviews, data = data_train)
plot(LM2)

summ(LM2)

#Calculate RMSE
RMSE <- sqrt(mean(residuals(LM2)^2))
cat("RMSE = ", RMSE)
```
Results:
With transforming price on the log scale, the following features are found to be significant:
- Private Room Type
- Shared Room Type
- Accommodates
- Bedrooms
- Bathrooms
- Host is Superhost
- Availability 365
- Cleaning fee
- Security Deposit
- Instant Bookable
- Minimum Nights
- Latitude
- Longitude
- Number of Reviews

The R-squared metric is 68.68%.

RMSE = 0.322

3. Taking log for accommodates: log-accommodates
```{r log-acc}
set.seed(891)
LM_model3 <- summary(lm(formula = log_price ~ room_type + log_acc + bedrooms + bathrooms + property_type + host_is_superhost + cleaning_fee + cancellation_policy + extra_people + security_deposit + instant_bookable + maximum_nights + minimum_nights + bed_type + latitude + longitude + number_of_reviews + availability_365, data = data_train))

plot(price ~room_type + log_acc + bedrooms + bathrooms + factor(neighbourhood_cleansed) + host_is_superhost + cancellation_policy + availability_365 +  cleaning_fee + extra_people + security_deposit + instant_bookable + maximum_nights + minimum_nights + bed_type + latitude + longitude + number_of_reviews, data = data_train) + abline(lm(formula = log_price ~ room_type + accommodates + bedrooms + bathrooms, data = data_train))
LM_model3

LM3 <- lm(formula = log_price ~ room_type + log_acc + bedrooms + bathrooms + property_type + host_is_superhost + cleaning_fee + cancellation_policy + extra_people + security_deposit + instant_bookable + maximum_nights + minimum_nights + bed_type + latitude + longitude + number_of_reviews + availability_365, data = data_train)

plot(LM3)

summ(LM3)

#Calculate RMSE
RMSE <- sqrt(mean(residuals(LM3)^2))
cat("RMSE = ", RMSE)
```
Results:
With transforming price on the log scale, the following features are found to be significant (same as previous model):
- Private Room Type
- Shared Room Type
- Accommodates
- Bedrooms
- Bathrooms
- Host is Superhost
- Availability 365
- Cleaning fee
- Security Deposit
- Instant Bookable
- Minimum Nights
- Latitude
- Longitude
- Number of Reviews

The R-squared metric is 69.15%.

RMSE = 0.321


Since the R-squared value is slightly better with a log transformed accommodates feature, I will be using this going forward.

4. Pairwise comparison
```{r pairwise}
set.seed(891)
# price vs bathrooms
p1 <- ggplot(data_train, aes(x=bathrooms, y=log_price)) +
  geom_point(colour = "orange", size = 1.5) +
  geom_smooth(method='lm', color='red') +
  ylab("price")

# price vs accommodates
p2 <- ggplot(data_train, aes(x=log_acc, y=log_price)) +
  geom_point(colour = "orange", size = 1.5) +
  geom_smooth(method='lm', color='red') +
  ylab("price")

# price vs bedrooms
p3 <- ggplot(data_train, aes(x=bedrooms, y=log_price)) +
  geom_point(colour = "orange", size = 1.5) +
  geom_smooth(method='lm', color='red') +
  ylab("price")

# price vs room_type
p4 <- ggplot(data_train, aes(x=room_type, y=log_price)) +
  geom_point(colour = "orange", size = 1.5) +
  geom_smooth(method='lm', color='red') +
  ylab("price")

#price vs host_is_superhost
p5 <- ggplot(data_train, aes(x=host_is_superhost, y=log_price)) +
  geom_point(colour = "orange", size = 1.5) +
  geom_smooth(method='lm', color='red') +
  ylab("price")

# price vs cancellation policy
p6 <- ggplot(data_train, aes(x=cancellation_policy, y=log_price)) +
  geom_point(colour = "orange", size = 1.5) +
  geom_smooth(method='lm', color='red') +
  ylab("price")

# price vs cleaning fee
p7 <- ggplot(data_train, aes(x=cleaning_fee, y=log_price)) +
  geom_point(colour = "orange", size = 1.5) +
  geom_smooth(method='lm', color='red') +
  ylab("price")

#price vs extra people
p8 <- ggplot(data_train, aes(x=extra_people, y=log_price)) +
  geom_point(colour = "orange", size = 1.5) +
  geom_smooth(method='lm', color='red') +
  ylab("price")

#price vs security deposit
p9 <- ggplot(data_train, aes(x=security_deposit, y=log_price)) +
  geom_point(colour = "orange", size = 1.5) +
  geom_smooth(method='lm', color='red') +
  ylab("price")

#price vs availability 365
p10 <- ggplot(data_train, aes(x=availability_365, y=log_price)) +
  geom_point(colour = "orange", size = 1.5) +
  geom_smooth(method='lm', color='red') +
  ylab("price")

#price vs instant_bookable
p11 <- ggplot(data_train, aes(x=instant_bookable, y=log_price)) +
  geom_point(colour = "orange", size = 1.5) +
  geom_smooth(method='lm', color='red') +
  ylab("price")

# combining all plots
grid.arrange(p1, p2, p3, p4,p5,p6, p7,p8, p9, p10, p11, ncol = 4)
```

5. Preparing test and train data 
```{r}
set.seed(891)
library(caTools)
set.seed(1)
data_train <- data_train %>% drop_na(log_price) %>% drop_na(log_acc) %>% drop_na(bedrooms) %>% drop_na(bathrooms) %>% drop_na(room_type) %>% drop_na(host_is_superhost) %>% drop_na(cancellation_policy) %>% drop_na(availability_365) %>% drop_na(cleaning_fee) %>% drop_na(extra_people) %>% drop_na(security_deposit) %>% drop_na(instant_bookable) %>% drop_na(maximum_nights) %>% drop_na(minimum_nights) %>% drop_na(property_type) %>% drop_na(bed_type) %>% drop_na(latitude) %>% drop_na(longitude) %>% drop_na(number_of_reviews) %>% drop_na(availability_365)

sample_size <- floor(0.75 * nrow(data_train))

set.seed(1)
train_ind <- sample(seq_len(nrow(data_train)), size = sample_size)
train <- data_train[train_ind,]
test <- data_train[-train_ind,]
```


6. Predicting Price using Linear Regression Model
Using log transformed Price and Accommodates - LM_model3
```{r predict-price}
set.seed(891)
#trainingRowIndex <- sample(1:nrow(data_train), 0.8 *nrow(data_train))
#trainingData <- data_train[trainingRowIndex,]
#testData <- data_train[-trainingRowIndex,]

lm_model <- lm(formula = log_price ~ room_type + log_acc + bedrooms + bathrooms + property_type + host_is_superhost + cleaning_fee + cancellation_policy + extra_people + security_deposit + instant_bookable + maximum_nights + minimum_nights + bed_type + latitude + longitude + number_of_reviews + availability_365, data = train)

price_pred <- predict(lm_model, test)

summary(lm_model)

LM_pred <- data.frame(cbind(Actual = test$log_price, Predicted= price_pred))
correlation_accuracy <- cor(LM_pred)

head(LM_pred)
correlation_accuracy
plot(lm_model)
summ(lm_model)

#Calculate prediction error - RMSE
RMSE(price_pred, test$log_price)

```
R-squared value is 69.9%.

RMSE = 0.335

7. Correlation matrix
```{r corr}
set.seed(891)


data.corr <- data_train %>% filter(!is.na(log_acc)) %>% filter(!is.na(bathrooms)) %>% filter(!is.na(bedrooms)) %>% filter(!is.na(availability_365)) %>% filter(!is.na(cleaning_fee)) %>% filter(!is.na(extra_people)) %>% filter(!is.na(security_deposit)) %>% filter(!is.na(maximum_nights)) %>% filter(!is.na(minimum_nights)) %>% filter(!is.na(latitude)) %>% filter(!is.na(longitude)) %>% filter(!is.na(number_of_reviews)) %>%
  select(log_price, price, bedrooms, bathrooms, log_acc, availability_365, cleaning_fee, security_deposit, maximum_nights, minimum_nights, extra_people, latitude, longitude, number_of_reviews)
kable(cor(data.corr)) %>% kable_styling()

#rcorr(data.corr, type = c("pearson", "spearman"))
#corrplot(data.corr)
```
There is some correlation between the variables, but not strong enough to worry about multicollinearity.

#### Decision Trees
1. Create a CART model - with log price and bedroom
```{r Cart1}
set.seed(891)
# Implement CART model
library(rpart)
library(rpart.plot)
CARTmodel1 = rpart(log_price ~ bedrooms, data = train, cp =0.001)

prp(CARTmodel1)
printcp(CARTmodel1)
plotcp(CARTmodel1)
summary(CARTmodel1)

plot(CARTmodel1, uniform=TRUE,
   main="Regression Tree for Price ")
text(CARTmodel1, use.n=TRUE, all=TRUE, cex=.8)
```

2. Making predictions
```{r pred1}
set.seed(891)
# Make predictions
predDT = predict(CARTmodel1, newdata = test)

#SSE
SSE = sum((predDT - test$log_price)^2)
cat("SSE = ",SSE)

# RMSE
#Calculate RMSE
RMSE(predDT, test$log_price)

# Baseline
baseline = mean(train$log_price)
cat("\nbaseline = ",baseline)

# SSE of baseline model on testing set
SSEb = sum((baseline - test$log_price)^2)
cat("\nSSEb = ", SSEb)

# R^2
Rsquared = 1 - SSE/SSEb
cat("\nR squared = ",Rsquared)

#Actual vs Predicted
DT_pred <- data.frame(cbind(Actual = test$log_price, Predicted = predDT))
head(DT_pred)

```
R-squared value of 44% with 1 predictor - bedrooms.

RMSE = 0.426


3. Creating CART model to predict log(price) using variables: room_type, log_acc, bedrooms, bathrooms, property_type, host_is_superhost, cleaning_fee, cancellation_policy, extra_people, security_deposit, instant_bookable, maximum_nights,  minimum_nights, bed_type, latitude, longitude, number_of_reviews

```{r}
set.seed(891)
CARTmodel2 = rpart(log_price ~ room_type + log_acc + bedrooms + bathrooms + property_type + host_is_superhost + cleaning_fee + cancellation_policy + extra_people + security_deposit + instant_bookable + maximum_nights + minimum_nights + bed_type + latitude + longitude + number_of_reviews, data = train, cp =0.001)
prp(CARTmodel2)

rpart.plot(CARTmodel2)
printcp(CARTmodel2)
plotcp(CARTmodel2)
#summary(CARTmodel2)

plot(CARTmodel2, uniform=TRUE,
   main="Regression Tree for Price ")
text(CARTmodel2, use.n=TRUE, all=TRUE, cex=.8)
```

4. Make predictions:
```{r}
set.seed(891)
# Make predictions
predDT2 = predict(CARTmodel2, newdata = test)

#head(predTest)

#SSE
SSE = sum((predDT2 - test$log_price)^2)
cat("SSE = ",SSE)

#Calculate RMSE
RMSE(price_pred, test$log_price)

# Baseline
baseline = mean(train$log_price)
cat("\nbaseline = ",baseline)

# SSE of baseline model on testing set
SSEb = sum((baseline - test$log_price)^2)
cat("\nSSEb = ", SSEb)

# R^2
Rsquared = 1 - SSE/SSEb
cat("\nR squared = ",Rsquared)

#Actual vs Predicted
DT_pred2 <- data.frame(cbind(Actual = test$log_price, Predicted = predDT2))
head(DT_pred2)
plot(DT_pred2)

```
Using Decision Trees, R-squared value is 61.69%.

RMSE = 0.335

Predictor variables found to be significant are:
- bathrooms 
- bedrooms
- cancellation_policy
- cleaning_fee
- extra_people 
- latitude
- log_accommodates
- longitude,
- maximum_nights 
- minimum_nights 
- number_of_reviews 
- property_type
- room_type  
- security_deposit.

#### Random Forest

1. Building a Random Forest

mtry = # of variables: ranges from 1 to p (# of predictor variables)
mtry = p -> bagging
mtry = 1 -> split variable is completely random
sampsize = the number of samples to train on. The default value is 63.25% of the training set since this is the expected value of unique observations in the bootstrap sample. Lower sample sizes can reduce the training time but may introduce more bias than necessary. Increasing the sample size can increase performance but at the risk of overfitting because it introduces more variance. Typically, when tuning this parameter we stay near the 60-80% range.

```{r RTF}
set.seed(891)
library(randomForest)
library(party)
library(inTrees)
library(randomForestExplainer)
library(sandwich)
set.seed(1)

RFmodel1 <- randomForest(formula = log_price ~room_type + log_acc + bedrooms + bathrooms + property_type + host_is_superhost + cleaning_fee + cancellation_policy + extra_people + security_deposit + instant_bookable + maximum_nights + minimum_nights + bed_type + latitude + longitude + number_of_reviews + availability_365, data=train, nodesize = 100, ntree = 1000, mtry = 8, importance = TRUE, na.action = na.omit)

RFmodel1
#getTree(RFmodel1, 1)
plot(RFmodel1)

round(importance(RFmodel1),2)
```

2. Using Random Forest to make predictions
```{r predict2}
# Make predictions
predTest <- predict(RFmodel1, newdata = test, type = "response")

RT_pred <- data.frame(cbind(Actual = test$log_price, Predicted= predTest))

RT_pred
```

3. Estimating Performance for Random Forest 
````{r}
set.seed(891)
# SSE
SSE = sum((predTest - test$log_price)^2)
cat("SSE = ",SSE)

#Calculate RMSE
RMSE(predTest, test$log_price)

# Baseline
baseline = mean(test$log_price)
cat("\nBaseline = ",baseline)

# SSE of baseline model on testing set
SSEb = sum((baseline - test$log_price)^2)
cat("\nSSEb = ",SSEb)

# R^2
Rsquared = 1 - SSE/SSEb
cat("\nRsquared = ",Rsquared)
```
Random Forest model gives an R-squared value of 66.55%.

RMSE = 0.329

4. Parameter Tuning - Grid Search
https://uc-r.github.io/random_forests
```{r}
set.seed(891)
library(mlbench)
library(caret)
library(e1071)

set.seed(1)

trControl <- trainControl(method = "cv", number = 10, search = "grid", savePredictions = TRUE)
rf_default <- train(log_price ~room_type + log_acc + bedrooms + bathrooms + property_type + host_is_superhost + cleaning_fee + cancellation_policy + extra_people + security_deposit + instant_bookable + maximum_nights + minimum_nights + bed_type + latitude + longitude + number_of_reviews, data = train, method = "rf", metric = "Rsquared", trControl = trControl)
print(rf_default)

```

5. Building a new Random Forest with parameter - mtry = 17
```{r RTF2}
set.seed(891)
library(randomForest)
library(party)
library(inTrees)
library(randomForestExplainer)
set.seed(1)

RFmodel2 <- randomForest(formula = log_price ~room_type + log_acc + bedrooms + bathrooms + property_type + host_is_superhost + cleaning_fee + cancellation_policy + extra_people + security_deposit + instant_bookable + maximum_nights + minimum_nights + bed_type + latitude + longitude + number_of_reviews +availability_365, data=train, nodesize = 100, ntree = 1000, mtry = 17, importance = TRUE, na.action = na.omit)

RFmodel2
#getTree(RFmodel1, 1)
plot(RFmodel2)

round(importance(RFmodel2),2)
```

6. Interpreting New Model
```{r RTF-results3}
set.seed(891)
explain_forest(RFmodel2, interactions = TRUE, data = train)

#Variable Importance Measure
importance_frame <- measure_importance(RFmodel2)
save(importance_frame, file = "importance_frame.rda")
load("importance_frame.rda")
importance_frame

#Multi Way Importance Plot
plot_multi_way_importance(RFmodel2, size_measure = "no_of_nodes") # gives the same result as below but takes longer
plot_multi_way_importance(importance_frame, size_measure = "no_of_nodes")

#Comparing Measures Using ggpairs
plot_importance_ggpairs(RFmodel2) 
plot_importance_ggpairs(importance_frame)


```
Based on the Random Forest Model, the top 10 predictor variables were found to be:
- bedrooms
- room type
- bathrooms
- cleaning fee
- longitude
- latitude
- log accommodates
- security deposit
- minimum nights
- availability 365


7. Using new Random Forest to make predictions
```{r predict}
set.seed(891)
# Make predictions
predTest2 <- predict(RFmodel2, newdata = test, type = "response")

RT_pred2 <- data.frame(cbind(Actual = test$log_price, Predicted= predTest2))

RT_pred2
```

8. Parameters for Random Forest 
````{r}
set.seed(891)
# SSE
new_SSE = sum((predTest2 - test$log_price)^2)
cat("SSE = ",new_SSE)

#Calculate RMSE
RMSE(predTest2, test$log_price)

# Baseline
new_baseline = mean(test$log_price)
cat("\nBaseline = ",new_baseline)

# SSE of baseline model on testing set
new_SSEb = sum((new_baseline - test$log_price)^2)
cat("\nSSEb = ",new_SSEb)

# R^2
new_Rsquared = 1 - (new_SSE/new_SSEb)
cat("\nRsquared = ",new_Rsquared)
```
Random Forest model gives an R-squared value of 66.94%.

RMSE = 0.327

9. Measuring Predictive Accuracy - Not complete
```{r}
# create training and validation data 
library(rsample)
set.seed(123)
valid_split <- initial_split(train, .8)

# training data
ames_train_v2 <- analysis(valid_split)

# validation data
ames_valid <- assessment(valid_split)
x_test <- ames_valid[setdiff(names(ames_valid), c("log_acc","bedrooms", "bathrooms","property_type","host_is_superhost", "cleaning_fee", "cancellation_policy", "extra_people", "security_deposit","instant_bookable","maximum_nights","minimum_nights","bed_type","latitude","longitude","number_of_reviews"))]
y_test <- ames_valid$log_price

#rf_oob_comp <- randomForest(formula = log_price ~room_type + log_acc + bedrooms + bathrooms + property_type + host_is_superhost + cleaning_fee + cancellation_policy + extra_people + security_deposit + instant_bookable + maximum_nights + minimum_nights + bed_type + latitude + longitude + number_of_reviews,data    = ames_train_v2,xtest   = x_test,ytest   = y_test)

# extract OOB & validation errors
oob <- sqrt(rf_oob_comp$mse)
validation <- sqrt(rf_oob_comp$test$mse)

# compare error rates
tibble::tibble(
  `Out of Bag Error` = oob,
  `Test error` = validation,
  ntrees = 1:rf_oob_comp$ntree
) %>%
  gather(Metric, RMSE, -ntrees) %>%
  ggplot(aes(ntrees, RMSE, color = Metric)) +
  geom_line() +
  scale_y_continuous(labels = scales::dollar) +
  xlab("Number of trees")
```


### ANOVA - to compare categorical data

```{r anova}
AOV_res <- summary(aov(log_price ~ factor(property_type) + factor(neighbourhood_cleansed), data = data_train))
AOV_res
```
Results:
As we can see, property type and neighborhood do have a significant impact on price since the p-value is much smaller than the significance level.

### GLM

1. Fitting model

```{r glm}
set.seed(891)
library(jtools)
library(kableExtra)

glm_1 <- glm(log_price ~ room_type + log_acc + bedrooms + bathrooms + property_type + host_is_superhost + cleaning_fee + cancellation_policy + extra_people + security_deposit + instant_bookable + maximum_nights + minimum_nights + bed_type + latitude + longitude + number_of_reviews + availability_365, data = data_train, family = gaussian(link  = "log"))
summary(glm_1)

## review the model summary
#summary(glm_1)
#hist(glm_1$fitted.values)
summ(glm_1,robust = 'HC1', scale = TRUE, confint = TRUE)

plot(glm_1)

## check out fitted values against the originals
plot(data_train$log_price, glm_1$fitted.values)
abline(a=0, b=1)

#plot(data_train$log_price, glm_1$fitted.values, xlim = c(0, 37000), ylim = c(0, 37000))
#abline(a=0, b=1)


## MSE and RMSE
MSE <-  mean((glm_1$fitted.values - data_train$log_price)^2)
RMSE <- sqrt(mean((glm_1$fitted.values - data_train$log_price)^2))

cat("MSE = ", MSE)
cat("\n RMSE = ", RMSE)
```
Results:
With transforming price on the log scale, the following features are found to be significant (same as previous model):
- Private Room Type
- Shared Room Type
- Accommodates
- Bedrooms
- Bathrooms
- Property Type
- Host is Superhost
- Availability 365
- Cleaning fee
- Security Deposit
- Instant Bookable
- Minimum Nights
- Latitude
- Longitude
- Number of Reviews

The R-squared metric is 67%.

RMSE = 0.323

2. Price Prediction with GLM
```{r glm_predict}
set.seed(891)

glm_model <- glm(log_price ~ room_type + log_acc + bedrooms + bathrooms + property_type + host_is_superhost + cleaning_fee + cancellation_policy + extra_people + security_deposit + instant_bookable + maximum_nights + minimum_nights + bed_type + latitude + longitude + number_of_reviews + availability_365, family = gaussian(link  = "log"), data = train)

price_glm <- predict(glm_model, test)

summary(glm_model)
summ(glm_model, robust = 'HC1')

GLM_pred <- data.frame(cbind(Actual = test$log_price, Predicted= exp(price_glm)))
correlation_accuracy <- cor(GLM_pred)

head(GLM_pred)
correlation_accuracy
#plot(glm_model)

#effect_plot(glm_model, pred = bedrooms, interval = TRUE, plot.points = TRUE)
plot_summs(glm_model)
plot_summs(glm_model, scale = TRUE, plot.distributions = TRUE, inner_ci_level = .9)

#Calculate RMSE
RMSE(price_glm, test$log_price)
```
Using the McFadden Pseudo R-squared parameter, we get an R-squared value of 68%.

RMSE = 3.264

## Model Evaluation

1. Actual vs Predicted Plot

```{r actual-vs-predicted}
set.seed(891)

ggplot(LM_pred, aes(x = Predicted, y = Actual)) + geom_point() + geom_smooth(method = "lm") + ggtitle("Linear Regression Prediction")

ggplot(RT_pred, aes(x = Predicted, y = Actual)) + geom_point() + ggtitle("Random Forest Prediction") + geom_smooth()

ggplot(GLM_pred, aes(x = Predicted, y = Actual)) + geom_point() + geom_smooth(method = "glm") + ggtitle("Generalized Linear Model Prediction")

ggplot(DT_pred2, aes(x = Predicted, y = Actual)) + geom_point() + ggtitle("Decision Trees Prediction")  + geom_smooth()

#hist(LM_pred$Predicted, main = "Linear Regression")
#hist(RT_pred$Predicted, main = "Random Forest")
#hist(GLM_pred$Predicted, main = "Generalized Linear Regression")
#hist(DT_pred2$Predicted, main = "Decision Trees")
```
Random Forest?

2. $R^2$ Value

Model is evaluated based on R-squared value i.e. % of variability explained by the model.
```{r plot-R2}
error_val <- data.frame("Rsquared" = c(61.7, 67, 68.0, 70), "Model" = c("Decision Trees", "Random Forest", "Generalized Linear Regression", "Linear Regression"))
#sort(error_val$Rsquared)
ggplot(error_val, aes(Model,Rsquared)) + geom_bar(stat = "identity", position = "dodge") + theme_bw() + xlab("Model") + ylab("R-squared value") + geom_text(aes(label = Rsquared, vjust = -0.25, hjust = 0)) 
```

Predicting price on:
- Linear Regression has an R-squared value of 69.82%.
- Decision Trees have an R-squared value of  61.69%
- Random Forest has an R-square value of 66.94%.
- GLM has an R-squared value of 68%.

Based on this, Linear Regression wins.
But, $R^2$ value is not sufficient to assess a model's accuracy.

3. RMSE - Prediction Error
```{r plot-RMSE}

error_val2 <- data.frame("RMSE" = c(0.335, 0.355, 0.327, 3.264), "Model" = c("Linear Regression", "Decision Trees", "Random Forest", "GLM"))
#sort(error_val$Rsquared)
ggplot(error_val2, aes(Model,RMSE)) + geom_bar(stat = "identity", position = "dodge") + theme_bw() + xlab("Model") + ylab("RMSE") + geom_text(aes(label = RMSE, vjust = -0.25, hjust = 0))
```

- Linear Regression has an RMSE value of 0.335
- Decision Trees have an RMSE value of 0.335
- Random Forest has an RMSE value of 0.327
- GLM has an RMSE value of 3.264

Based on this, Random Forest model has the least prediction error, and hence is the best model.


## Utility Function

Utility Function created  calculates the profit predicted using each model. This utility function sums up the difference between Actual and Predicted values.
Since price was log transformed, both Actual and Predicted values are first exponentiated to find the sum of all differences.

1. For Linear Regression
```{r util-LR}
set.seed(891)
library(formattable)
LM_pred$diff <- exp(LM_pred$Actual) - exp(LM_pred$Predicted)
diff1 <- sum(LM_pred$diff)
head(LM_pred)
cat("Sum of diff = ", diff1)
```

A profit of $1978 is predicted with Linear Regression.

2. Decision Trees
```{r util-DT}
set.seed(891)
DT_pred2$diff <- exp(DT_pred2$Actual) - exp(DT_pred2$Predicted)
diff2 <- sum(DT_pred2$diff)
head(DT_pred2)
cat("\n Diff = ", (diff2))
```

A profit of $1249 is predicted with Decision Trees.

3. Random Forest
```{r util-RTF}
set.seed(891)
RT_pred2$diff <- exp(RT_pred2$Actual) - exp(RT_pred2$Predicted)
diff3 <- sum(RT_pred2$diff)

head(RT_pred2)
cat("\n Diff = ",(diff3))
```

A profit of $2815 predicted with a random forest.

4. GLM
```{r util-GLM}
set.seed(891)
GLM_pred$diff <- exp(GLM_pred$Actual) - exp(GLM_pred$Predicted)
diff4 <- sum(GLM_pred$diff)
head(GLM_pred)
cat("Sum of diff = ", (diff4))
```

A profit of $1886 is predicted with GLM model.

```{r plot-util}
library(formattable)
utility <- data.frame("Profit" = c(currency(ceiling(diff1)), currency(ceiling(diff2)), currency(ceiling(diff3)), currency(ceiling(diff4))), "Model" = c("Linear Regression", "Decision Trees", "Random Forest", "GLM"))

ggplot(utility, aes(Model,Profit)) + geom_bar(stat = "identity", position = "dodge") + theme_bw() + xlab("Model") + ylab("Profit") + geom_text(aes(label = Profit, vjust = -0.25, hjust = 0))
```

Based on this, Random Forest model gives the highest profit.