---
title: "House Price Prediction"
author: "Zhiyuan Xie"
output: 
  flexdashboard::flex_dashboard:
    theme: simplex
    orientation: columns
    social: 
    source_code: embed
---

```{r setup, include=FALSE}
# load necessary packages
library(ggplot2)
library(plotly)
library(plyr)
library(flexdashboard)  ## you need this package to create dashboard

# read the data set here, I use data: mtcars as an example
prodata <- read.csv("D:/ECE/MTH543/project/kc_house_data.csv")
```

Introduction
=======================================================================

Column {.tabset data-width=650}
-----------------------------------------------------------------------

### Get to know the data

```{r}
head(prodata)

prodata$condition=as.numeric(prodata$condition)

prodata$grade=as.numeric(prodata$grade)

prodata$yr_built=as.numeric(prodata$yr_built)

prodata$yr_renovated=as.numeric(prodata$yr_renovated)
  
#head(prodata)

index = sample(1:nrow(prodata),size=0.8*nrow(prodata))

train=prodata[index,]

test=prodata[-index,]

train$zipcat=rep(NA,nrow(train))

test$zipcat=rep(NA,nrow(test))

meanpr=mean(train$price)

zip=names(table(train$zipcode))

for(i in 1:length(table(train$zipcode)))
{
  
  index1=which(train$zipcode==zip[i])
  
  index2=which(test$zipcode==zip[i])
  
  if(mean(train$price[index1])>=meanpr)
  {
    train$zipcat[index1]=1
    test$zipcat[index2]=1
    
  }
  else
  {
    train$zipcat[index1]=0 
    test$zipcat[index2]=0
  }
  
}
```


### Zipcode region with regional median price 

![Caption for the picture.](D:/ECE/MTH543/project/median_price.jpg)


Column {data-width=350}
-----------------------------------------------------------------------

### Abstract

 

This project predicts the price of houses in King County, Washington by using the house selling data that is download from Kaggle between May 2014 and May 2015 of the county. We selected and processed the data carefully. After that, we built a multiple linear regression model. The model passed the diagnostic and hypothesis test, has a high R-square, and get a reasonable mean absolute error for predicting testing data.


 

### Introduction

 

In life, the house is an important part of living and growing, many people would face the price problem if they want to buy or sell a house. Therefore, the goal of this project is to predict the price of houses in King County, Washington. The dataset was provided by the King county, and download from Kaggle [1], it includes 21613 houses selling price between May 2014 and May 2015 of the county. According to the zip code region with regional median price, the different zip code region has large differences in price.

Available themes include:

- Id

- Date

- Price

- Bedrooms

- Sqft_living

- Sqft_lot

- Floors

- Waterfront

- View

- Condition

- Grade

- Sqft_basement

- Yr_built

- Yr_renovated

- Zipcode

- Lat

- Long

- Sqft_living15

- Sqft_lot15


 



Data Exploration & Method part 1
=======================================================================

Column {.tabset data-width=650}
-----------------------------------------------------------------------

### transform response

```{r}
par(mfrow=c(1,2))
hist(train$price,xlab="price", main="Histogram of price")
train$price <- log(train$price)
hist(train$price,xlab="log(price)", main="Histogram of log(price)")

```

### Price vs Floors

```{r}

boxplot(train$price~train$floors,xlab="floors",ylab="standardised price", col="#03f8fc")

```


### Price vs Grade

```{r}

boxplot(train$price~train$grade,col="#03f8fc",xlab="grade", ylab="standardised price")

```


### Price vs View

```{r}

boxplot(train$price~train$view,col="#03f8fc",xlab="view", ylab="standardised price")

```

### Price vs Condition

```{r}
train$floors<- ifelse(train$floors>3,3,train$floors)
train$floors<-as.factor(train$floors)
train$zipcat<-as.factor(train$zipcat)
train$sqft_basement<-as.factor(ifelse(train$sqft_basement>0,0,1))
train$waterfront<-as.factor(train$waterfront)
train$yr_renovated<-ifelse((train$yr_renovated-train$yr_built)<0,0,train$yr_renovated-train$yr_built)

#train$renovated<-as.factor(ifelse(train$yr_renovated>0,0,1))

varind=c(3,4,5,6,7,11,12,15,16)

trainmean=apply(train[varind],2,mean)

trainvar=apply(train[varind],2,var)

test$floors<- ifelse(test$floors>3,3,test$floors)
test$price <- log(test$price)
#test$sqft_living <- log(test$sqft_living)
test$floors<-as.factor(test$floors)
test$zipcat<-as.factor(test$zipcat)
test$sqft_basement<-as.factor(ifelse(test$sqft_basement>0,0,1))
test$waterfront<-as.factor(test$waterfront)
test$yr_renovated<-ifelse((test$yr_renovated-test$yr_built)<0,0,test$yr_renovated-test$yr_built)

par(mfrow=c(1,2))

boxplot(exp(train$price)~train$condition,col="#03f8fc",xlab="condition", ylab="price",ylim=c(0,2e+6))

boxplot(prodata$price~prodata$condition,col="#03f8fc",xlab="condition", ylab="price",ylim=c(0,2e+6))

```


Column {data-width=350}
-----------------------------------------------------------------------
 

The dataset is divided into training and testing, 80% for training (17290 observations), 20% for testing (4323 observations). We standardize the training and testing data, the testing data is standardized by using the mean and standard deviation of training data.

After observing the price distribution from the histogram, we decide to transform the price by using log() function. As we can see, the histogram of price close to the normal distribution after transformation.

Selected regressors for multiple regression include (blue fonts represent indicator variables):

  - Bedrooms

  - Sqft_living

  - Sqft_lot

  - $\color{blue}{\text{Floors}}$

  - $\color{blue}{\text{Waterfront}}$

  - View

  - Condition

  - Grade

  - $\color{blue}{\text{Sqft_basement}}$

  - Yr_built

  - Yr_renovated

  - $\color{blue}{\text{Zipcat}}$

In the above variables, the sqrt_living is the summation of sqrt_above and sqrt_basement, therefore, we decide to skip sqrt_above, and transform the sqrt_basement to the sign to represent that the house has a basement or not. The waterfront is a sign to represent the house near the water or not, so we also use it as an indicator variable. 

To further study the data, we observe the relationships between price and some variables. According to these box plots. For the condition, we found the condition 1 has fewer samples, to avoid the effect of extremums, so we plot both training data and all data. We found that there is a significant positive correlation between price and these variables except floors. Therefore, the floors also becomes an indicator variable.

Then, the Yr_renovated is a variable representing the year of renovation for the house, the Yr_renovated is equal to zero if the house was not renovated. We change it to the variable represent the house that was renovated how many years after the house was built. We also found that there have no enough samples for 3.5 floors, so we transform all samples that have 3.5 floors to 3 floors.

Furthermore, we add a new variable named zipcat, we will talk about it in next part of Data Exploration & Method section





Data Exploration & Method part 2
=======================================================================

```{r}

#test$renovated<-as.factor(ifelse(test$yr_renovated>0,0,1))

for(i in c(1:length(varind)))
{
  train[varind[i]]=(train[varind[i]]-trainmean[i])/sqrt(trainvar[i])
  test[varind[i]]=(test[varind[i]]-trainmean[i])/sqrt(trainvar[i])
  
}
```

Column {.tabset data-width=650}
-----------------------------------------------------------------------

### model summary for using Zipcode as an indicator variable

```{r}

train$zipcode<-as.factor(train$zipcode)
fit_zipcode=lm(formula = price~bedrooms+bathrooms+sqft_living+floors+condition+grade
+yr_built+yr_renovated+sqft_basement+sqft_lot+waterfront+zipcode+0,data=train)

summary(fit_zipcode)

```

### Linearity

```{r}

Fitted.Values <- fit_zipcode$fitted.values

# Extract residuals
Residuals <- fit_zipcode$residuals

# Calculate standardized residuals 
Standardized.Residuals <- scale(fit_zipcode$residuals)

# Extract fitted values for lm() object
Theoretical.Quantiles <- qqnorm(Residuals, plot.it = F)$x

# find Square root of abs(residuals)
Root.Residuals <- sqrt(abs(Standardized.Residuals))

# Calculate Leverage
Leverage <- lm.influence(fit_zipcode)$hat

# Create data frame 
# Will be used as input to plot_ly

diagnostics <- data.frame(Fitted.Values, 
                     Residuals, 
                     Standardized.Residuals, 
                     Theoretical.Quantiles,
                     Root.Residuals,
                     Leverage)

m <- list(
  l = 100,
  r = 100,
  b = 100,
  t = 100,
  pad = 4
)

# Fitted vs Residuals
p1 <- plot_ly(diagnostics, x = Fitted.Values, y = Residuals, 
          type = "scatter", mode = "markers", hoverinfo = "x+y", name = "Data", 
          marker = list(size = 10, opacity = 0.5))%>%
  layout(title = "Residuals vs Fitted Values", 
       xaxis = list(title="Fitted Values", font=list(size=14)), 
       yaxis = list(title="Residuals", font=list(size=14)), 
       plot_bgcolor = "#e6e6e6", 
       font=list(size=14), margin=m)

ggplotly(p1)

```

### Normality
```{r}
# QQ Pot
p2 <- plot_ly(diagnostics, x = Theoretical.Quantiles, y = Standardized.Residuals, type = "scatter", mode = "markers", hoverinfo = "x+y", name = "Data", marker = list(size = 10, opacity = 0.5), showlegend = F)%>%
    add_trace(x = Theoretical.Quantiles, y = Theoretical.Quantiles, type = "scatter", mode = "line", name = "", line = list(width = 2))%>%
  layout(title = "Q-Q Plot", plot_bgcolor = "#e6e6e6",
              xaxis = list(title="Theoretical Quantiles", font=list(size=14)), 
       yaxis = list(title="Standardized Residuals", font=list(size=14)), font=list(size=14), margin=m)

ggplotly(p2)
```

### Equality Variance
```{r}
# Scale Location
p3 <- plot_ly(diagnostics, x = Fitted.Values, y = Root.Residuals, 
          type = "scatter", mode = "markers", hoverinfo = "x+y", name = "Data",
          marker = list(size = 10, opacity = 0.5), showlegend = F)%>%
  layout(title = "Scale-Location", plot_bgcolor = "#e6e6e6", xaxis = list(title="Fitted Values", font=list(size=14)), 
       yaxis = list(title=expression(sqrt("|Standardized Residuals|")), font=list(size=14)), font=list(size=14), margin=m)

ggplotly(p3)
```

### Residuals vs Leverage
```{r}
s <- loess.smooth(Leverage, Residuals)
p4 <- plot_ly(diagnostics, x = Leverage, y = Residuals, 
            type = "scatter", mode = "markers", hoverinfo = "x+y", name = "Data", marker = list(size = 10, opacity = 0.5), showlegend = F) %>% 
    add_trace(x = s$x, y = s$y, type = "scatter", mode = "line", name = "Smooth", line = list(width = 2)) %>% 
    layout(title = "Leverage vs Residuals", plot_bgcolor = "#e6e6e6", xaxis = list(title="Leverage", font=list(size=14)), 
       yaxis = list(title="Residuals", font=list(size=14)), font=list(size=14), margin=m)
ggplotly(p4)  
```



Column {.tabset data-width=350}
-----------------------------------------------------------------------


We use the Zipcode as the indicator variable for the linear regreesion. Although the R-square of the model is not low, according to the t-test, some zip codes are insignificant at 0.05 level, and the Q-Q plot shows that the points do not follow the 45-degree line.

To solve the problem, we decided not to use the zip code as an indicator variable directly. The zip code includes 70 different Zip codes for houses, each Zip code is a category. To reduce the number of categories, we transform the Zip codes to two categories by comparing the mean price of each category to the overall mean price, if the mean price of the category is smaller than the overall mean price, it belongs to zip category zero, otherwise, it belongs to zip category one. In all data, 9933 samples belong to zip category one, and 11680 samples belong to zip category zero.

Furthermore, Both latitude, longitude and zip code are location information for the house. In this project, we decided to use zip code but not latitude and longitude. Because the zip code is already a well-grouped category, but the latitude and longitude need some group strategy before use.






Diagnostics
=======================================================================
Column {.tabset data-width=650}
-----------------------------------------------------------------------

```{r}


#train[varind]<-as.data.frame(apply(train[varind],2,scale))

fit=lm(formula = price~bedrooms+bathrooms+sqft_living+floors+condition+grade
+yr_built+yr_renovated+sqft_basement+sqft_lot+waterfront+zipcat+0,data=train)

#fit=lm(formula = price~bedrooms+bathrooms+sqft_living+floors+condition+grade
#      +yr_built+renovated+yr_renovated*renovated+sqft_basement+sqft_lot+waterfront+zipcat+0,data=train)

testprice <- exp(test$price*sqrt(trainvar[1])+trainmean[1])
predictprice <- exp(predict(fit,newdata=test)*sqrt(trainvar[1])+trainmean[1])

MAE=mean(abs(testprice - predictprice))



Fitted.Values <- fit$fitted.values

# Extract residuals
Residuals <- fit$residuals

# Calculate standardized residuals 
Standardized.Residuals <- scale(fit$residuals)

# Extract fitted values for lm() object
Theoretical.Quantiles <- qqnorm(Residuals, plot.it = F)$x

# find Square root of abs(residuals)
Root.Residuals <- sqrt(abs(Standardized.Residuals))

# Calculate Leverage
Leverage <- lm.influence(fit)$hat

# Create data frame 
# Will be used as input to plot_ly

diagnostics <- data.frame(Fitted.Values, 
                     Residuals, 
                     Standardized.Residuals, 
                     Theoretical.Quantiles,
                     Root.Residuals,
                     Leverage)
```

### Linearity

```{r}
m <- list(
  l = 100,
  r = 100,
  b = 100,
  t = 100,
  pad = 4
)

# Fitted vs Residuals
p1 <- plot_ly(diagnostics, x = Fitted.Values, y = Residuals, 
          type = "scatter", mode = "markers", hoverinfo = "x+y", name = "Data", 
          marker = list(size = 10, opacity = 0.5))%>%
  layout(title = "Residuals vs Fitted Values", 
       xaxis = list(title="Fitted Values", font=list(size=14)), 
       yaxis = list(title="Residuals", font=list(size=14)), 
       plot_bgcolor = "#e6e6e6", 
       font=list(size=14), margin=m)

ggplotly(p1)
```

### Normality
```{r}
# QQ Pot
p2 <- plot_ly(diagnostics, x = Theoretical.Quantiles, y = Standardized.Residuals, type = "scatter", mode = "markers", hoverinfo = "x+y", name = "Data", marker = list(size = 10, opacity = 0.5), showlegend = F)%>%
    add_trace(x = Theoretical.Quantiles, y = Theoretical.Quantiles, type = "scatter", mode = "line", name = "", line = list(width = 2))%>%
  layout(title = "Q-Q Plot", plot_bgcolor = "#e6e6e6",
              xaxis = list(title="Theoretical Quantiles", font=list(size=14)), 
       yaxis = list(title="Standardized Residuals", font=list(size=14)), font=list(size=14), margin=m)

ggplotly(p2)
```

### Equality Variance
```{r}
# Scale Location
p3 <- plot_ly(diagnostics, x = Fitted.Values, y = Root.Residuals, 
          type = "scatter", mode = "markers", hoverinfo = "x+y", name = "Data",
          marker = list(size = 10, opacity = 0.5), showlegend = F)%>%
  layout(title = "Scale-Location", plot_bgcolor = "#e6e6e6", xaxis = list(title="Fitted Values", font=list(size=14)), 
       yaxis = list(title=expression(sqrt("|Standardized Residuals|")), font=list(size=14)), font=list(size=14), margin=m)

ggplotly(p3)
```

### Residuals vs Leverage
```{r}
s <- loess.smooth(Leverage, Residuals)
p4 <- plot_ly(diagnostics, x = Leverage, y = Residuals, 
            type = "scatter", mode = "markers", hoverinfo = "x+y", name = "Data", marker = list(size = 10, opacity = 0.5), showlegend = F) %>% 
    add_trace(x = s$x, y = s$y, type = "scatter", mode = "line", name = "Smooth", line = list(width = 2)) %>% 
    layout(title = "Leverage vs Residuals", plot_bgcolor = "#e6e6e6", xaxis = list(title="Leverage", font=list(size=14)), 
       yaxis = list(title="Residuals", font=list(size=14)), font=list(size=14), margin=m)
ggplotly(p4)  
```


Column {data-width=350}
-----------------------------------------------------------------------
 

The final model is built after data selection and transformation. So We need to check the model by diagnostics process.

**1. Linearity Assumption**

We check the linearity assumption: according to the residual vs fitted value plot, there is a linear relationship between the price and all used variable.

**2. Normality Assumption**

We check the normality assumption: the residuals almost follow the normal distribution. We exam this assumption by checking the Q-Q plot.

**3. Equal Variance Assumption**

We check the equal variance assumption: the variability of the errors should be about the same for all values of each predictor. We exam this assumption by checking the Scale-Location plot.

**4. leverage or influencial points**
Since the Q-Q plot looks good, there is no need to check leverage or influencial points. But we can still check it by exam the Residuals vs Leverage plot. There are no prominence samples in the plot.



Results & discussion
=======================================================================

Column {data-width=600}
-----------------------------------------------------------------------

### model summary

```{r}
summary(fit)
```

Column {data-width=400}
-----------------------------------------------------------------------


### The Analysis of Variance Table
```{r}
options(knitr.kable.NA = '')
variance_table <- anova(fit)
V_table <- data.frame(matrix(NA, nrow=3, ncol=5))
rownames(V_table) <- c("Regression", "Residual Error", "Total")
colnames(V_table) <- c("DF", "Sum Sq", "Mean Sq", "F Value", "Pr(>F)")
V_table$DF <- c(15, 17274, 17289)
V_table$`Sum Sq` <- c(sum(variance_table$`Sum Sq`[1:3]), variance_table$`Sum Sq`[4], sum(variance_table$`Sum Sq`[1:4]))
V_table$`Mean Sq`[1:2] <- c(mean(variance_table$`Sum Sq`[1:3]), variance_table$`Mean Sq`[4])
V_table$`F Value`[1] <- V_table$`Mean Sq`[1]/V_table$`Mean Sq`[2]
V_table$`Pr(>F)`[1] <- pf(V_table$`F Value`[1], V_table$DF[1], V_table$DF[2], lower.tail = F)
library(knitr)
kable(V_table)
```


### Mean absolute error & standard deviation of testing price

```{r}
cat("Mean absolute error:", MAE)
cat("Standard deviation of testing price:", sd(testprice))

```

### discussion

 

According to the model summary, all variables are significant at 0.01 level of t-test, and more than 77% of the total variability in the price can be explained by the used variables.

Then, we can reject the $H_0$ at a significance level of 0.01 by examining the p-value (9.429656e-11<0.01) from the Analysis of Variance Table.

Furthermore, the predicted mean absolute error for the testing price seems not low, but it is reasonable because of the high standard deviation of testing price.



### conclusion

 

In this project, we predict the price of houses in King County, Washington. The model passed the diagnostic and hypothesis test, has a high R-square, and get a reasonable mean absolute error for predicting testing data. Overall, the project is succeeded.



### reference

[1] “House Sales in King County, USA.” [Online]. Available: https://kaggle.com/harlfoxem/housesalesprediction. [Accessed: 04-Dec-2019].