script3.Rmd

---
output: html_document
---
---
title: "<center><h1>Detailed Data Analysis & Ensemble Modeling</h1></center>"
author: 
  "<center>Tanner Carbonati</center>"
date: "<center>April 5, 2017</center>"
output:
  pdf_document: default
  html_document:
    keep_md: yes
---

  There have been a lot of great scripts/notebooks for this dataset and I thought I would share my own. This was written for a course at UCSD so I did my best to go through intensive detail and include the theory used for cleaning the data and training models. The model used for training is an ensemble of XGBoost, Ridge, Lasso and Elastic-Net regularization. Many of ideas used in this notebook were influenced by a few other kernels including:
  
- https://www.kaggle.com/opanichev/house-prices-advanced-regression-techniques/ensemble-of-4-models-with-cv-lb-0-11489
- https://www.kaggle.com/humananalog/house-prices-advanced-regression-techniques/xgboost-lasso

  Our goal is to transform all of our data into numeric, while preserving as much information from the categoric variables as we're going to use a gradient boosting method (XGBoost). As XGBoost has trouble with extrapolation we will use ridge, lasso and elastic-net regularization, which will also account for a lot of the multicollinearity in the data. Once we have predictions from all 4 regressors we can average the results to get an accurate estimate.  We'll go more in depth on why we chose these methods and how we'll use them later on. 
  
Let the fun begin!

```{r warning=FALSE, message=FALSE}
require(ggplot2) # for data visualization
require(stringr) #extracting string patterns
require(Matrix) # matrix transformations
require(glmnet) # ridge, lasso & elastinet
require(xgboost) # gbm
require(randomForest)
require(Metrics) # rmse
require(dplyr) # load this in last so plyr doens't overlap it
require(caret) # one hot encoding
require(scales) # plotting $$
require(e1071) # skewness
require(corrplot) # correlation plot
```

The first thing we'll do is load in the training and testing data. The training data consists of 1460 rows and 81 columns while the testing has 1459 rows and 80 columns (excluding the `SalePrice` column), which in this dataset is the dependent variable we are trying to predict. For us to perform analysis at a more productive rate we can combine the 2 dataframes together and run analysis on both data sets at once, then split the data once we are ready to train a model. 

The `Id` feature is useless so we can toss it out of our dataset and we won't include `SalePrice` since it is our response variable. We also won't import string variables as factors since our ultimate goal is to tranform all our variables to numeric. 

```{r warning=FALSE, message=FALSE}
train <- read.csv('../input/train.csv', stringsAsFactors = FALSE)
test <- read.csv('../input/test.csv', stringsAsFactors = FALSE)


# combine the datasets
df.combined <- rbind(within(train, rm('Id','SalePrice')), within(test, rm('Id')))
dim(df.combined)
```

Our dataset is filled with many missing values, therefore, before we can build any predictive model we'll clean our data by filling in all NA's with appropriate values. For each column we'll try to replace NA's by using features that have a strong correlation, which will help us determine what values to fill in. 

```{r warning=FALSE, message=FALSE}
na.cols <- which(colSums(is.na(df.combined)) > 0)
sort(colSums(sapply(df.combined[na.cols], is.na)), decreasing = TRUE)
paste('There are', length(na.cols), 'columns with missing values')
```

```{r}
# helper function for plotting categoric data for easier data visualization
plot.categoric <- function(cols, df){
  for (col in cols) {
    order.cols <- names(sort(table(df.combined[,col]), decreasing = TRUE))
  
    num.plot <- qplot(df[,col]) +
      geom_bar(fill = 'cornflowerblue') +
      geom_text(aes(label = ..count..), stat='count', vjust=-0.5) +
      theme_minimal() +
      scale_y_continuous(limits = c(0,max(table(df[,col]))*1.1)) +
      scale_x_discrete(limits = order.cols) +
      xlab(col) +
      theme(axis.text.x = element_text(angle = 30, size=12))
  
    print(num.plot)
  }
}

```

**PoolQC: Pool quality**

There are 2909 missing values in PoolQC, which represents the quality of the houses pool. It is a reasonable assumption that houses with NA's for PoolQC most likely don't have a pool, however, for us to check this we can see if any of the houses that have a NA for PoolQC recorded a PoolArea > 0 (**PoolArea:** Pool area in square feet). it is safe to assume that houses with a PoolArea of 0 don't have a pool and for those greater than 0 do. 

```{r warning=FALSE, message=FALSE}
plot.categoric('PoolQC', df.combined)
df.combined[(df.combined$PoolArea > 0) & is.na(df.combined$PoolQC),c('PoolQC','PoolArea')]
```

Of the 2909 houses with NA's for PoolQC only 3 recorded a PoolArea greater than 0. To fill the missing values for these houses we can see what the mean pool area is for each quality of pool. To do this we'll group by each unique vlaue of `PoolQC` and compute the mean of PoolArea for that quality. Once we have the mean areas computed we can assign the quality that has the closest corresponding mean of PoolArea

```{r warning=FALSE, message=FALSE}
df.combined[,c('PoolQC','PoolArea')] %>%
  group_by(PoolQC) %>%
  summarise(mean = mean(PoolArea), counts = n()) 

```

Note that only 10 houses total have pools  The rest of the houses PoolQC can be filled with 'NoPool' since each has a corresponding PoolArea of 0. We'll fill in the 3 PoolQC's from above with the quality having the closest corresponding mean. 


```{r warning=FALSE, message=FALSE}
df.combined[2421,'PoolQC'] = 'Ex'
df.combined[2504,'PoolQC'] = 'Ex'
df.combined[2600,'PoolQC'] = 'Fa'
df.combined$PoolQC[is.na(df.combined$PoolQC)] = 'None'
```

**GarageType: Garage location**  
**GarageYrBlt: Year garage was built**  
**GarageFinish: Interior finish of the garage**  
**GarageCars: Size of garage in car capacity**  
**GarageArea: Size of garage in square feet**  
**GarageQual: Garage quality**  
**GarageCond: Garage condition**  

Lets see what the deal is with `GarageYrBlt`. It seems reasonable that most houses would build a garage when the house itself was built. We can check this by seeing how many houses were built the same year their garage was built. 


```{r}
length(which(df.combined$GarageYrBlt == df.combined$YearBuilt))
```

2216 of the 2919 houses have same year for for GarageYrBlt and YearBuilt. Lets replace any of the NA's for `GarageYrBlt` with the year from  `YearBuilt`. 

```{r}
idx <- which(is.na(df.combined$GarageYrBlt))
df.combined[idx, 'GarageYrBlt'] <- df.combined[idx, 'YearBuilt']
```

That leaves 6 garage features in our dataset and 4 of them have at least 157 missing values while `GarageArea` and `GarageCars` both only have 1, thus we can assume this particular house does not have a garage at all. For the rest of the houses we can check to see that if the NA's recorded also have 0 GarageArea and 0 GarageCars. If they do we can fill in their missing values with 'None' since having 0 area and 0 cars in their garage will imply that they do not have any at all. 

```{r}
garage.cols <- c('GarageArea', 'GarageCars', 'GarageQual', 'GarageFinish', 'GarageCond', 'GarageType')
df.combined[is.na(df.combined$GarageCond),garage.cols]
```

Only one house who had NA's in their garage columns had an area graeteer than 0. We can fill this house in manually and set the rest of the houses NA's to 0. 

For the house with *GarageArea = 360* and *GarageCars = 1*, but NA's in the other columns, we can use the most frequent values for each columns from houses with a similar area and car count. 

```{r warning=FALSE, message=FALSE}
idx <- which(((df.combined$GarageArea < 370) & (df.combined$GarageArea > 350)) & (df.combined$GarageCars == 1))

names(sapply(df.combined[idx, garage.cols], function(x) sort(table(x), decreasing=TRUE)[1]))
  
df.combined[2127,'GarageQual'] = 'TA'
df.combined[2127, 'GarageFinish'] = 'Unf'
df.combined[2127, 'GarageCond'] = 'TA'
```

Now we can fill in any missing numeric values with 0 and categoric with 'None' since these houses recorded having 0 area and 0 cars in their garage.

```{r warning=FALSE, message=FALSE}
for (col in garage.cols){
  if (sapply(df.combined[col], is.numeric) == TRUE){
    df.combined[sapply(df.combined[col], is.na), col] = 0
  }
  else{
    df.combined[sapply(df.combined[col], is.na), col] = 'None'
  }
}
```

**KitchenQual: Kitchen quality**  
**Electrical: Electrical system**  

With only 1 missing value for `KitchenQual` and `Electrical` each we can fill in the missing value with the most frequent value from each column. 

```{r warning=FALSE, message=FALSE}
plot.categoric('KitchenQual', df.combined)
df.combined$KitchenQual[is.na(df.combined$KitchenQual)] = 'TA'

plot.categoric('Electrical', df.combined)
df.combined$Electrical[is.na(df.combined$Electrical)] = 'SBrkr'
```

- **BsmtQual: Height of the basement**  
- **BsmtCond: General condition of the basement**  
- **BsmtExposure: Walkout or garden level basement walls**  
- **BsmtFinType1: Quality of basement finished area**  
- **BsmtFinSF1: Type 1 finished square feet**  
- **BsmtFinType2: Quality of second finished area (if present)**  
- **BsmtFinSF2: Type 2 finished square feet**  
- **BsmtUnfSF: Unfinished square feet of basement area**  
- **TotalBsmtSF: Total square feet of basement area**  
- **BsmtFullBath: Basement full bathrooms**  
- **BsmtHalfBath: Basement half bathrooms**  

There are 11 basement features each with at least 1 missing value. We can take a look at the subset of just these columns from our data.

```{r warning=FALSE, message=FALSE}

bsmt.cols <- names(df.combined)[sapply(names(df.combined), function(x) str_detect(x, 'Bsmt'))]

df.combined[is.na(df.combined$BsmtExposure),bsmt.cols]
plot.categoric('BsmtExposure', df.combined)
```

Almost all of the missing values for  each categoric basement feature comes from houses with 0 on each features corresponding to area. We can fill in these values with 'None' since these houses certainly don't have basements. Rows 949, 1488 and 2349 are the only missing values from BsmtExposure, we can fill this with *No* as that is the most frequent value and these houses most likely don't have any exposure for their basements. The rest of the basement columns corresponding to area will be filled with 0 since they likely don't have a basement and the categoric missing values will be filled with *NoBsmt*. 

```{r warning=FALSE, message=FALSE}
df.combined[c(949, 1488, 2349), 'BsmtExposure'] = 'No'

for (col in bsmt.cols){
  if (sapply(df.combined[col], is.numeric) == TRUE){
    df.combined[sapply(df.combined[col], is.na),col] = 0
  }
  else{
    df.combined[sapply(df.combined[col],is.na),col] = 'None'
  }
}
```

- **Exterior1st: Exterior covering on house**  
- **Exterior2nd: Exterior covering on house (if more than one material)** 

```{r warning = FALSE, message=FALSE}
#plot.categoric(c('Exterior1st', 'Exterior2nd'), df.combined)
idx <- which(is.na(df.combined$Exterior1st) | is.na(df.combined$Exterior2nd))
df.combined[idx,c('Exterior1st', 'Exterior2nd')]
```

There is only 1 missing value for `Exterior1st` and `Exterior2nd` coming from the same hosue and there aren't any other features that can help us predict what variable should be filled so we can fill this with *'Other'* since the NA is likely due to having an exterior cover that is not listed. 

```{r warning=FALSE, message=FALSE}
df.combined$Exterior1st[is.na(df.combined$Exterior1st)] = 'Other'
df.combined$Exterior2nd[is.na(df.combined$Exterior2nd)] = 'Other'

```

- **SaleType: Type of sale**  
- **Functional: Home functionality rating**  
- **Utilities: Type of utilities available**  

`SaleType`, `Functional` and `Utilities` have less than 3 missing values. For SaleType we can see what the SaleCondition of the house was and use a contingency table to see which SaleType and SaleCondition overlap together the most. 

```{r warning=FALSE, message=FALSE}
plot.categoric('SaleType', df.combined)

df.combined[is.na(df.combined$SaleType),c('SaleCondition')]

table(df.combined$SaleCondition, df.combined$SaleType)
```

Most houses with a SaleCondition of 'Normal' almost all have a SaleType of 'WD'. We'll replace the missing value accordingly. 

```{r warning=FALSE, message=FALSE}
df.combined$SaleType[is.na(df.combined$SaleType)] = 'WD'

plot.categoric('Functional', df.combined)
df.combined$Functional[is.na(df.combined$Functional)] = 'Typ'
```

Utilities only has 1 value for NoSeWa and the rest AllPub. We can drop this feature from our dataset as the house with 'NoSeWa' is from our training set and will have won't help with any predictive modelling

```{r warning=FALSE, message=FALSE}
plot.categoric('Utilities', df.combined)
which(df.combined$Utilities == 'NoSeWa') # in the training data set


col.drops <- c('Utilities')
df.combined <- df.combined[,!names(df.combined) %in% c('Utilities')]
```

- **MSZoning: The general zoning classification**  
- **MSSubClass: The building class**  


There are only 4 missing values for MSZoning. We can see what the subclass is for the houses with missing values for Zoning.

```{r warning=FALSE, message=FALSE}
df.combined[is.na(df.combined$MSZoning),c('MSZoning','MSSubClass')]
plot.categoric('MSZoning', df.combined)

table(df.combined$MSZoning, df.combined$MSSubClass)
```

For Subclasses with 20 'RL' has the largest frequency, however, for Subclasses with 30 and 70 'RM' has the most frequency. We will fill the missing values accordingly. 

```{r warning=FALSE, message=FALSE}
df.combined$MSZoning[c(2217, 2905)] = 'RL'
df.combined$MSZoning[c(1916, 2251)] = 'RM'
```

- **MasVnrType: Masonry veneer type**  
- **MasVnrArea: Masonry veneer area in square feet**

There are 23 missing values for `MasVnrArea` and 24 for `MasVnrType`. We can see if both missing values come from the same houses

```{r warning=FALSE, message=FALSE}
df.combined[(is.na(df.combined$MasVnrType)) | (is.na(df.combined$MasVnrArea)), c('MasVnrType', 'MasVnrArea')]
```

All but one house has missing values for both columns. For houses with NA's on both columns we can fill 0 for the area and None for the type since they likely do not have a masonry veneer. For the house with a MasVnrArea of 198 but NA for MasVnrType we can record the median areas for each type and see which type is closest to 198

```{r warning=FALSE, message=FALSE}
na.omit(df.combined[,c('MasVnrType','MasVnrArea')]) %>%
  group_by(na.omit(MasVnrType)) %>%
  summarise(MedianArea = median(MasVnrArea,na.rm = TRUE), counts = n()) %>%
  arrange(MedianArea)

#plot.categoric('MasVnrType', df.combined)

df.combined[2611, 'MasVnrType'] = 'BrkFace'
```

The areas we can replace with 0 and types can be replaced with 'None'

```{r warning=FALSE, message=FALSE}
df.combined$MasVnrType[is.na(df.combined$MasVnrType)] = 'None'
df.combined$MasVnrArea[is.na(df.combined$MasVnrArea)] = 0
```

- **LotFrontage: Linear feet of street connected to property**  

There are 486 missing values for LotFrontage, which is quite a lot of values to fill and we can't just replace these with 0. We're given that "LotFrontage: Linear feet of street connected to property." The area of each street connected to the house property is most likely going to have a similar area to other houses in its neighborhood. We can group by each neighborhood and take the median of each LotFrontage and fill the missing values of each LotFrontage based on what neighborhood the house comes from. 

```{r warning=FALSE, message=FALSE}

df.combined['Nbrh.factor'] <- factor(df.combined$Neighborhood, levels = unique(df.combined$Neighborhood))

lot.by.nbrh <- df.combined[,c('Neighborhood','LotFrontage')] %>%
  group_by(Neighborhood) %>%
  summarise(median = median(LotFrontage, na.rm = TRUE))
lot.by.nbrh

idx = which(is.na(df.combined$LotFrontage))

for (i in idx){
  lot.median <- lot.by.nbrh[lot.by.nbrh == df.combined$Neighborhood[i],'median']
  df.combined[i,'LotFrontage'] <- lot.median[[1]]
}
```

- **Fence: Fence quality**

We can replace any missing vlues for Fence and MiscFeature with 'None' as they probably don't have this feature with their property. 

```{r warning=FALSE, message=FALSE}
plot.categoric('Fence', df.combined)

df.combined$Fence[is.na(df.combined$Fence)] = 'None'


table(df.combined$MiscFeature)
df.combined$MiscFeature[is.na(df.combined$MiscFeature)] = 'None'
```

- **Fireplaces: Number of fireplaces**  
- **FireplaceQu: Fireplace quality **

FireplaceQu denotes the fireplace quality. We can check to see if any of the missing values for FireplaceQu come from houses that recorded having at least 1 fireplace.

```{r warning=FALSE, message=FALSE}
plot.categoric('FireplaceQu', df.combined)
which((df.combined$Fireplaces > 0) & (is.na(df.combined$FireplaceQu)))
```

All the houses that have missing values did not record having any fireplaces. We can replace the NA's with 'None' since these houses don't have any fireplaces at all. 

```{r warning=FALSE, message=FALSE}
df.combined$FireplaceQu[is.na(df.combined$FireplaceQu)] = 'None'
```


- **Alley: Type of alley access**

There are 2721 missing values for Alley and only 2 potential options - Grvl and Pave. We can fill 'None' for any of the houses with NA's as these houses must not have any type of alley access. 

```{r warning=FALSE, message=FALSE}
plot.categoric('Alley', df.combined)
df.combined$Alley[is.na(df.combined$Alley)] = 'None'
```

```{r warning=FALSE, message=FALSE}
paste('There are', sum(sapply(df.combined, is.na)), 'missing values left')
```

##Adding custom numeric features

Now that we've filled in all missing values for both the training and testing sets we can split our data into a numeric set and a categoric set. Since we want our dataset to be strictly numeric values while retaining as much information as possible, we will have to transform any categoric into a binary feature using one-hot encoding. First off, we can start by adding all our numeric features to a new dataframe along with adding any extra custom made features that will help with predictive modeling.

We can start by creating a new data frame of just the numeric features

```{r warning=FALSE, message=FALSE}
num_features <- names(which(sapply(df.combined, is.numeric)))
cat_features <- names(which(sapply(df.combined, is.character)))

df.numeric <- df.combined[num_features]
```

Next we can transform any of the ordinal variables (variables that can be scaled) into numeric values. We can do this by determining which order the categories follow and assigning the values an order from 1,2,..,n. We'll group each feature by its possible values and return the median SalePrice and mean OverallQual for each unique value, then map the higher priced/better quality homes larger numeric values and the lower priced/lower quality homes smaller numeric values. 

```{r warning=FALSE, message=FALSE}
group.df <- df.combined[1:1460,]
group.df$SalePrice <- train$SalePrice
      
# function that groups a column by its features and returns the mdedian saleprice for each unique feature. 
group.prices <- function(col) {
  group.table <- group.df[,c(col, 'SalePrice', 'OverallQual')] %>%
    group_by_(col) %>%
    summarise(mean.Quality = round(mean(OverallQual),2),
      mean.Price = mean(SalePrice), n = n()) %>%
    arrange(mean.Quality)
    
  print(qplot(x=reorder(group.table[[col]], -group.table[['mean.Price']]), y=group.table[['mean.Price']]) +
    geom_bar(stat='identity', fill='cornflowerblue') +
    theme_minimal() +
    scale_y_continuous(labels = dollar) +
    labs(x=col, y='Mean SalePrice') +
    theme(axis.text.x = element_text(angle = 45)))
  
  return(data.frame(group.table))
}

## functional to compute the mean overall quality for each quality
quality.mean <- function(col) {
  group.table <- df.combined[,c(col, 'OverallQual')] %>%
    group_by_(col) %>%
    summarise(mean.qual = mean(OverallQual)) %>%
    arrange(mean.qual)
  
  return(data.frame(group.table))
}


# function that maps a categoric value to its corresponding numeric value and returns that column to the data frame
map.fcn <- function(cols, map.list, df){
  for (col in cols){
    df[col] <- as.numeric(map.list[df.combined[,col]])
  }
  return(df)
}
```

Any of the columns with the suffix 'Qual' or 'Cond' denote the quality or condition of that specific feature. Each of these columns have the potential values: `r unique(df.combined$GarageQual)`. We'll compute the mean house prices for these unique values to get a better sense of what their abbreviations mean. 

```{r warning=FALSE, message=FALSE}
qual.cols <- c('ExterQual', 'ExterCond', 'GarageQual', 'GarageCond', 'FireplaceQu', 'KitchenQual', 'HeatingQC', 'BsmtQual')

group.prices('FireplaceQu')
group.prices('BsmtQual')
group.prices('KitchenQual')

```

From seeing the mean saleprices from a few of the quality and condition features we can infer that the abbreviations mean **poor**, **fair**, **typical/average**, **good** and **excelent**. We'll map numeric values from 0-5 to their corresponding categoric values (including 0 for **None**) and combine that to our dataframe.

**Note:** we will set 'None' = 0 for all categories as *None* signifies that the house does not have that particular quality/condition to rank and regardless of the houses overall quality or sale price we will keep 'None' = 0 for consistency. 

```{r warning=FALSE, mpessage=FALSE}
qual.list <- c('None' = 0, 'Po' = 1, 'Fa' = 2, 'TA' = 3, 'Gd' = 4, 'Ex' = 5)

df.numeric <- map.fcn(qual.cols, qual.list, df.numeric)
```


```{r warning=FALSE, message=FALSE}
group.prices('BsmtExposure')
bsmt.list <- c('None' = 0, 'No' = 1, 'Mn' = 2, 'Av' = 3, 'Gd' = 4)

df.numeric = map.fcn(c('BsmtExposure'), bsmt.list, df.numeric)

```

`BsmtFinType1` and `BsmtFineType2` represent the quality of the 1st and 2nd finished areas. There should be some order between the different types of qualities, we'll see if we can find a similarity between the 2 features.

```{r warning=FALSE, message=FALSE}
group.prices('BsmtFinType1')
```

- **BsmtFinType1: Quality of basement finished area**  
- **BsmtFinSF1: Type 1 finished square feet**  

Here returning the mean sale prices might not be as helpful as computing the median basement areas for both columns to determine which quality is better than the other. 

```{r warning=FALSE, message=FALSE}
# visualization for BsmtFinTyp2 instead of another table
df.combined[,c('BsmtFinType1', 'BsmtFinSF1')] %>%
  group_by(BsmtFinType1) %>%
  summarise(medianArea = median(BsmtFinSF1), counts = n()) %>%
  arrange(medianArea) %>%
  ggplot(aes(x=reorder(BsmtFinType1,-medianArea), y=medianArea)) +
  geom_bar(stat = 'identity', fill='cornflowerblue') +
  labs(x='BsmtFinType2', y='Median of BsmtFinSF2') +
  geom_text(aes(label = sort(medianArea)), vjust = -0.5) +
  scale_y_continuous(limits = c(0,850)) +
  theme_minimal()
```

Through investigating the relationships between the basement quality and areas we an see the true order of qualities of each basement to be *'None' < 'Unf' < 'LwQ' < 'BLQ' < 'Rec' < 'ALQ' < 'GLQ'*. 

```{r warning=FALSE, message=FALSE}
bsmt.fin.list <- c('None' = 0, 'Unf' = 1, 'LwQ' = 2,'Rec'= 3, 'BLQ' = 4, 'ALQ' = 5, 'GLQ' = 6)
df.numeric <- map.fcn(c('BsmtFinType1','BsmtFinType2'), bsmt.fin.list, df.numeric)
```

- **Functional: Home functionality rating**

This feature doesn't really tell us much and functionality is very vague to tie which other features have a correlation with it. We can compute the mean sale prices for each functional category and assign numeric values accordingly.

```{r warning=FALSE, message=FALSE}
group.prices('Functional')
functional.list <- c('None' = 0, 'Sal' = 1, 'Sev' = 2, 'Maj2' = 3, 'Maj1' = 4, 'Mod' = 5, 'Min2' = 6, 'Min1' = 7, 'Typ'= 8)

df.numeric['Functional'] <- as.numeric(functional.list[df.combined$Functional])
```

We'll continue to do the same for GarageFinish and Fence:

```{r warning=FALSE, message=FALSE}
group.prices('GarageFinish')
garage.fin.list <- c('None' = 0,'Unf' = 1, 'RFn' = 1, 'Fin' = 2)

df.numeric['GarageFinish'] <- as.numeric(garage.fin.list[df.combined$GarageFinish])


group.prices('Fence')
fence.list <- c('None' = 0, 'MnWw' = 1, 'GdWo' = 1, 'MnPrv' = 2, 'GdPrv' = 4)

df.numeric['Fence'] <- as.numeric(fence.list[df.combined$Fence])

MSdwelling.list <- c('20' = 1, '30'= 0, '40' = 0, '45' = 0,'50' = 0, '60' = 1, '70' = 0, '75' = 0, '80' = 0, '85' = 0, '90' = 0, '120' = 1, '150' = 0, '160' = 0, '180' = 0, '190' = 0)

df.numeric['NewerDwelling'] <- as.numeric(MSdwelling.list[as.character(df.combined$MSSubClass)])

      
```

We've transformed all the categoric features with an ordianl scale into a numeric columns. Let's see which variables have the strongest effect on a houses sale price. 

To determine which features have the strongest relationship with **SalePrice** we can compute the sample correlation coefficient between the 2 variables, $$r_{xy} = \frac{s_{xy}}{s_xs_y}$$ where $s_{xy} = \mathbf{E}[(X-\mathbf{E}[X])(Y-\mathbf{E}[Y])$ is the sample covariance and $s_x, s_y$ are the sample standard deviations. The correlation coefficient measures the how linearly related the 2 variables are. A coefficient of 0 means that the 2 variables show no linear relationship, a coefficient between (0,1] shows that they have positive relationship and a coefficient between [-1,0) means they have a negative relationship. We're particularly interested in variables that show strong relationship with **SalePrice** so we will focus primarily on features that have a coefficient > .5 or < -.5. 

```{r}
# need the SalePrice column
corr.df <- cbind(df.numeric[1:1460,], train['SalePrice'])

# only using the first 1460 rows - training data
correlations <- cor(corr.df)
# only want the columns that show strong correlations with SalePrice
corr.SalePrice <- as.matrix(sort(correlations[,'SalePrice'], decreasing = TRUE))

corr.idx <- names(which(apply(corr.SalePrice, 1, function(x) (x > 0.5 | x < -0.5))))

corrplot(as.matrix(correlations[corr.idx,corr.idx]), type = 'upper', method='color', addCoef.col = 'black', tl.cex = .7,cl.cex = .7, number.cex=.7)
```

From the correlation plot we can see the 10 features with the strongest effect on SalePrice.   

 - **OverallQual: Overall material and finish quality**  
 - **GrLivArea: Above grade (ground) living area square feet**  
These 2 features have the highest correlation and very reasonably make sense.

- **GarageCars: Size of garage in car capacity**  
- **GarageArea: Size of garage in square feet**  
These 2 features from our plot are the next best thing, however, this might not have been something we expected. We can also see that the correlation between the 2 is extremely high, 0.88, which makes sense as the area of the garage is a constraint on how many cars can fit. 

- **TotalBsmtSF: Total square feet of basement area**  
- **1stFlrSF: First Floor square feet**  
These features also make intuitive sense as anything corresponding to area/size of the house will have an effect on the price. We can also see that these 2 features have a strong linear relationship with one another, which makes sense as the size of the basement can certainly depend on the size of the houses 1st floor. 

**FullBath: Full bathrooms above grade**  
**TotRmsAbvGrd: Total rooms above grade (does not include bathrooms)**  
Something we might have expected as well is that the number of rooms in the house having a relationship with the houses price. We can further see that the strong correlation between **GrLivArea** and **TotRmsAbvGrd** since the size of the livign area will most likely be a contstraint on the number of rooms above ground. An interesting relationship we might want to look into are the relationship between houses with a small living area but many rooms, which will result in smaller rooms and vise versa. We may want to check to see if this is an indicator on higher/lower house prices. 

**YearBuilt: Original construction date**  
**YearRemodAdd: Remodel date**  
A positive correlation for both of these is not something we should be surprised by. The newer homes will likely give higher listings compared to older models. 

We can print a matrix of scatter plots to see what these relationships look like under the hood to get a better sense of whats going on. 

```{r warning=FALSE, message=FALSE}
require(GGally)
lm.plt <- function(data, mapping, ...){
   plt <- ggplot(data = data, mapping = mapping) + 
    geom_point(shape = 20, alpha = 0.7, color = 'darkseagreen') +
    geom_smooth(method=loess, fill="red", color="red") +
    geom_smooth(method=lm, fill="blue", color="blue") +
    theme_minimal()
  return(plt)
}

ggpairs(corr.df, corr.idx[1:6], lower = list(continuous = lm.plt))
```

The blue lines in the scatter plots represent a simple linear regression fit while the red lines represent a local polynomial fit. We can see both OverallQual and GrLivArea and TotalBsmtSF follow a linear model, but have some outliers we may want to look into. For instance, there are multiple houses with an overall quality of 10, but have suspisciously low prices. We can see similar behavior in GrLivArea and TotalBsmtSF. GarageCars and GarageArea both follow more of a quadratic fit. It seems that having a 4 car garage does not result in a higher house price and same with an extremely large area. 

```{r warning=FALSE, message=FALSE}
ggpairs(corr.df, corr.idx[c(1,7:11)], lower = list(continuous = lm.plt))
```

More of the same with the remaining features with 1stFlrSF, FullBath, TotRmsAbvGrd following linear model while YearBuilt and YearRemodAdd both having non-linear/quadratic models. Taking a closer look to YearBuilt and YearRemodAdd we can see that the most expensive homes are the most recently built and remodelled. 


Now for some of the nominal variables we can take one of the categories that is distict from the others and create a binary feature that returns 1 if the house has that specific value and 0 if it does not.

```{r warning=FALSE, message=FALSE}
plot.categoric('LotShape', df.combined)
```

LotShape has 3 values for having an irregular shape and only 1 for regular. We can create a binary column that returns 1 for houses with a regular lot shape and 0 for houses with any of the 3 irregular lot shapes. Using this method of turning a categoric feature into a binary column will ultimately  help our data train better through boosted models without using numeric placeholders on nominal data. 

```{r warning=FALSE, message=FALSE}
df.numeric['RegularLotShape'] <- (df.combined$LotShape == 'Reg') * 1
```

We'll use this exact same method for:

**LandContour: Flatness of the property**  
**LandSlope: Slope of property**  
**Electrical: Electrical system**  
**GarageType: Garage location**  
**PavedDrive: Paved driveway**  
**MiscFeature: Miscellaneous feature not covered in other categories** 

```{r warning=FALSE, message=FALSE}
plot.categoric('LandContour', df.combined)

df.numeric['LandLeveled'] <- (df.combined$LandContour == 'Lvl') * 1

plot.categoric('LandSlope', df.combined)

df.numeric['LandSlopeGentle'] <- (df.combined$LandSlope == 'Gtl') * 1

plot.categoric('Electrical', df.combined)

df.numeric['ElectricalSB'] <- (df.combined$Electrical == 'SBrkr') * 1

plot.categoric('GarageType', df.combined)

df.numeric['GarageDetchd'] <- (df.combined$GarageType == 'Detchd') * 1

plot.categoric('PavedDrive', df.combined)

df.numeric['HasPavedDrive'] <- (df.combined$PavedDrive == 'Y') * 1

df.numeric['HasWoodDeck'] <- (df.combined$WoodDeckSF > 0) * 1

df.numeric['Has2ndFlr'] <- (df.combined$X2ndFlrSF > 0) * 1

df.numeric['HasMasVnr'] <- (df.combined$MasVnrArea > 0) * 1
```

```{r warning=FALSE, message=FALSE}
plot.categoric('MiscFeature', df.combined)
```

For `MiscFeature` the only feature with a significant amount of houses having it is *Shed*. We can one-hot encode houses that have Sheds vs those who do not. 

```{r warning=FALSE, message=FALSE}
df.numeric['HasShed'] <- (df.combined$MiscFeature == 'Shed') * 1
```

**YearBuilt: Original construction date**  
**YearRemodAdd: Remodel date**  

Many of the houses recorded the same year for `YearBuilt` and `YearRemodAdd`. We can create a new column that records that a house was remodelled if the year it was built is different than the remodel year. This 

```{r warning=FALSE, message=FALSE}
df.numeric['Remodeled'] <- (df.combined$YearBuilt != df.combined$YearRemodAdd) * 1
```

We can also create a column that seperates which houses have been recently remodelled vs those who are not. Houses that have been remodelled after the year they were sold will fall into this category. 

```{r warning=FALSE, message=FALSE}
df.numeric['RecentRemodel'] <- (df.combined$YearRemodAdd >= df.combined$YrSold) * 1
```

There can be potential value to homes who were sold the same year they were built as this could be an indicator that these houses were hot in the market. 

```{r warning=FALSE, message=FALSE}
df.numeric['NewHouse'] <- (df.combined$YearBuilt == df.combined$YrSold) * 1
```

What about the houses with area based features equal to 0? Houses with 0 square footage for a columnshows that the house does not have that feature at all. We add a one-hot encoded column for returning 1 for any house with an area greater than 0 since this means that the house *does* have this feature and 0 for those who do not! 

```{r warning=FALSE, message=FALSE}
#cols.binary <- c('X2ndFlrSF', 'MasVnrArea', 'WoodDeckSF')
cols.binary <- c('X2ndFlrSF', 'MasVnrArea', 'WoodDeckSF', 'OpenPorchSF', 'EnclosedPorch', 'X3SsnPorch', 'ScreenPorch')

for (col in cols.binary){
  df.numeric[str_c('Has',col)] <- (df.combined[,col] != 0) * 1
}
```

We know how important the year a house was built and sold but what about what the specific month it was sold? How do houses sold during summer compare to the other seasons? 

```{r warning=FALSE, message=FALSE}
ggplot(df.combined, aes(x=MoSold)) +
  geom_bar(fill = 'cornflowerblue') +
  geom_text(aes(label=..count..), stat='count', vjust = -.5) +
  theme_minimal() +
  scale_x_continuous(breaks = 1:12)
```

The largest proportion of houses sold is during the summer months: May, June, July. Let's add a column that seperates the the summer houses from the rest. 

```{r warning=FALSE, message=FALSE}
df.numeric['HighSeason'] <- (df.combined$MoSold %in% c(5,6,7)) * 1
```

What about which Neighborhoods are more expensive than others? 

```{r warning=FALSE, message=FALSE, results=FALSE}
train[,c('Neighborhood','SalePrice')] %>%
  group_by(Neighborhood) %>%
  summarise(median.price = median(SalePrice, na.rm = TRUE)) %>%
  arrange(median.price) %>%
  mutate(nhbr.sorted = factor(Neighborhood, levels=Neighborhood)) %>%
  ggplot(aes(x=nhbr.sorted, y=median.price)) +
  geom_point() +
  geom_text(aes(label = median.price, angle = 45), vjust = 2) +
  theme_minimal() +
  labs(x='Neighborhood', y='Median price') +
  theme(text = element_text(size=12),
        axis.text.x = element_text(angle=45))
```


StoneBr, NoRidge, NridgHt have a large gap between them versus the rest of the median prices from any of the other neighborhods. It would be wise of us to check if this is from outliers or if these houses are much pricier as a whole. 

```{r warning=FALSE, message=FALSE}
other.nbrh <- unique(df.combined$Neighborhood)[!unique(df.combined$Neighborhood) %in% c('StoneBr', 'NoRidge','NridgHt')]

ggplot(train, aes(x=SalePrice, y=GrLivArea, colour=Neighborhood)) +
  geom_point(shape=16, alpha=.8, size=4) +
  scale_color_manual(limits = c(other.nbrh, 'StoneBr', 'NoRidge', 'NridgHt'), values = c(rep('black', length(other.nbrh)), 'indianred',
                                    'cornflowerblue', 'darkseagreen')) +
  theme_minimal() +
  scale_x_continuous(label=dollar)
```

3 houses from StoneBr, NoRidge and NridgHt with house prices over $500,000 which no other Neighborhood is in the same range. What pops out even more in this plot is the 2 houses with an enormous GrLivArea of over 4500 square ft yet very low sale prices. We certianly have an opportunity to explore these outliers in depth, we'll take a look at this later. 

In the mean time lets one-hot encode the more expensive neighborhoods and add that to our dataframe

```{r warnings=FALSE, message=FALSE}
nbrh.rich <- c('Crawfor', 'Somerst, Timber', 'StoneBr', 'NoRidge', 'NridgeHt')
df.numeric['NbrhRich'] <- (df.combined$Neighborhood %in% nbrh.rich) *1
```

How about a numeric mapping to the neighborhoods who have higher quality homes and run for larger sale prices?

```{r}
group.prices('Neighborhood')

nbrh.map <- c('MeadowV' = 0, 'IDOTRR' = 1, 'Sawyer' = 1, 'BrDale' = 1, 'OldTown' = 1, 'Edwards' = 1, 
             'BrkSide' = 1, 'Blueste' = 1, 'SWISU' = 2, 'NAmes' = 2, 'NPkVill' = 2, 'Mitchel' = 2,
             'SawyerW' = 2, 'Gilbert' = 2, 'NWAmes' = 2, 'Blmngtn' = 2, 'CollgCr' = 2, 'ClearCr' = 3, 
             'Crawfor' = 3, 'Veenker' = 3, 'Somerst' = 3, 'Timber' = 3, 'StoneBr' = 4, 'NoRidge' = 4, 
             'NridgHt' = 4)

df.numeric['NeighborhoodBin'] <- as.numeric(nbrh.map[df.combined$Neighborhood])
```

**SaleCondition:** Condition of sale

```{r warning=FALSE, message=FALSE}
group.prices('SaleCondition')
```

Tried to do some research to get a sense of what the terms mean for SaleCondition but couldn't find much. The houses with 'Partial' for SaleCondition have both the highest quality and highest sale price well above the other qualities. Instead of mapping each quality to its own numeric value we can one hot encode for houses with *Partial*

```{r warning=FALSE, message=FALSE}
df.numeric['PartialPlan'] <- (df.combined$SaleCondition == 'Partial') * 1

group.prices('HeatingQC')

heating.list <- c('Po' = 0, 'Fa' = 1, 'TA' = 2, 'Gd' = 3, 'Ex' = 4)

df.numeric['HeatingScale'] <- as.numeric(heating.list[df.combined$HeatingQC])
```


**1stFlrSF:** First Floor square feet  
**2ndFlrSF:** Second floor square feet  
**LowQualFinSF:** Low quality finished square feet (all floors)  
1stFlrSF + 2nFlrSF + LowQualFinSF = GrLivArea. 

```{r warning=FALSE, message=FALSE}
area.cols <- c('LotFrontage', 'LotArea', 'MasVnrArea', 'BsmtFinSF1', 'BsmtFinSF2', 'BsmtUnfSF',
             'TotalBsmtSF', 'X1stFlrSF', 'X2ndFlrSF', 'GrLivArea', 'GarageArea', 'WoodDeckSF', 
             'OpenPorchSF', 'EnclosedPorch', 'X3SsnPorch', 'ScreenPorch', 'LowQualFinSF', 'PoolArea')

df.numeric['TotalArea'] <- as.numeric(rowSums(df.combined[,area.cols]))

#paste('There are', sum((df.combined$X1stFlrSF + df.combined$X2ndFlrSF) != df.combined$GrLivArea), 'houses with LowQualFinSF > 0')
```


```{r warning=FALSE, message=FALSE}
df.numeric['AreaInside'] <- as.numeric(df.combined$X1stFlrSF + df.combined$X2ndFlrSF)
```

We've seen how strong of an effect the year of a house built has on the house price, therefore, as this dataset collects houses up until 2010 we can determine how old a house is and how long ago the house was sold:

```{r warning=FALSE, message=FALSE}
df.numeric['Age'] <- as.numeric(2010 - df.combined$YearBuilt)

df.numeric['TimeSinceSold'] <- as.numeric(2010 - df.combined$YrSold)

# how many years since the house was remodelled and sold 
df.numeric['YearSinceRemodel'] <- as.numeric(df.combined$YrSold - df.combined$YearRemodAdd)
```

Correlation plot with OverallQual
```{r warnings=FALSE, message=FALSE}
corr.OverallQual <- as.matrix(sort(correlations[,'OverallQual'], decreasing = TRUE))

corr.idx <- names(which(apply(corr.OverallQual, 1, function(x) (x > 0.5 | x < -0.5))))

corrplot(as.matrix(correlations[corr.idx, corr.idx]), type = 'upper',
         method = 'color', addCoef.col = 'black', tl.cex =.7, cl.cex = .7,
         number.cex = .7)
```

We may want to use PCA and go more in depth, if we have time we'll come back to this

(we didn't have time ): )

##Outliers

Earlier we saw some suspicious houses with abnormally large GrLivArea's, 2 of which had very low SalePrices. These outliers may disrupt our ability to accurately predict. Lets take a closer look at these values and see if we may want to remove them

```{r warning=FALSE, message=FALSE}
train.test.df <- rbind(dplyr::select(train,-SalePrice), test)
train.test.df$type <- c(rep('train',1460),rep('test',1459))

ggplot(train, aes(x=GrLivArea)) +
  geom_histogram(fill='lightblue',color='white') +
  theme_minimal()

outlier_values <- boxplot.stats(train$GrLivArea)$out  # outlier values.
boxplot(train$GrLivArea, main="GrLivArea", boxwex=0.1)
mtext(paste("Outliers: ", paste(outlier_values[outlier_values>4000], collapse=", ")), cex=0.6)

ggplot(train.test.df, aes(x=type, y=GrLivArea, fill=type)) +
  geom_boxplot() +
  theme_minimal() +
  scale_fill_manual(breaks = c("test", "train"), values = c("indianred", "lightblue"))
```

For the training data we can see 4 houses whose GrLivArea is greater than 4000 yet there is one in the testing set. While it is always a great option to remove outliers from our data as they are usually telling us something more than meets the eye about what is going on with the data, these houses in the trainng set are obnoxiously large and ultimately do not add much value and are causing heavy skewness in both the SalePrice and GrLivArea and particularly the 2 values that have above a 4000 GrLivArea but low SalePrice are putting a constraint on the correlation between the 2 variables. Lets throw these houses out of training set begin preprocessing our data now that there are no heavy outliers. 

```{r warning=FALSE, message=FALSE}
idx.outliers <- which(train$GrLivArea > 4000)
df.numeric <- df.numeric[!1:nrow(df.numeric) %in% idx.outliers,]
df.combined <- df.combined[!1:nrow(df.combined) %in% idx.outliers,]
dim(df.numeric)
```


##PCA

I unfortunately did not make use of this, if you have a strong familiarity this could certainly improve your score. 

```{r warning=FALSE, message=FALSE}
require(factoextra)
pmatrix <- prcomp(df.numeric, center = TRUE, scale. = TRUE)

pcaVar <- as.data.frame(c(get_pca_var(pmatrix)))

# lets
pcaVarNew <- pcaVar[, 1:10]
```

##Preprocessing

If we want to use any type of linear regression model an important assumption we need to check is for normality in any of the dependant variables. We can use a Kolmogorov-Smirnof test or compute the skewness/kurtosis in each column to verify normality. A Kolmogorov-Smirnof test compares the sample distribution to a normal and returns a p-value determining if the the 2 distributions are similar. Skewness is a measure of symmetry where distributions with 0 skew follow a normal shape. Kurtosis measures the taildness of the distribution. For skewnesses outside the range of -0.8 to 0.8 and kurtosises outside the of -3.0 to 3.0 do not satisfy the assumption of normality  For any features that are not normally distributed we can make a non-linear transformation like a log-transformation such that $f(x)=log(x+1)$ when a column has 0 and $f(x)=log(x)$ when there are no 0's in the column. We do this because $log(0)=-\infty$. We will also scale all of the numeric data by standardizing the data (opposed to normalizing) as we know our data has potential outliers we don't want to bound each column. To standardize the observations at indivisual column we compute $$z_{ij} = \frac{x_{ij}-\bar{x_j}}{s_j}$$ where $\bar{x_j}$ is the sample mean at column *j* and $s_j$ is the sample deviation at column *j*. 

```{r warning=FALSE, message=FALSE}
require(psych)
# linear models assume normality from dependant variables 
# transform any skewed data into normal
skewed <- apply(df.numeric, 2, skewness)
skewed <- skewed[(skewed > 0.8) | (skewed < -0.8)]

kurtosis <- apply(df.numeric, 2, kurtosi)
kurtosis <- kurtosis[(kurtosis > 3.0) | (kurtosis < -3.0)]

# not very useful in our case
ks.p.val <- NULL
for (i in 1:length(df.numeric)) {
  test.stat <- ks.test(df.numeric[i], rnorm(1000))
  ks.p.val[i] <- test.stat$p.value
}

for(col in names(skewed)){
  if(0 %in% df.numeric[, col]) {
    df.numeric[,col] <- log(1+df.numeric[,col])
  }
  else {
    df.numeric[,col] <- log(df.numeric[,col])
  }
}

# normalize the data
scaler <- preProcess(df.numeric)
df.numeric <- predict(scaler, df.numeric)

```

For the rest of the categoric features we can one-hot encode each value to get as many splits in the data as possible  

```{r warning=FALSE, message=FALSE}
# one hot encoding for categorical data
# sparse data performs better for trees/xgboost
dummy <- dummyVars(" ~ .",data=df.combined[,cat_features])
df.categoric <- data.frame(predict(dummy,newdata=df.combined[,cat_features]))
```

What about YearBuilt, GarageYrBlt and YearRemodAdd?

YearBuilt had a correlation coefficient of 0.57 with SalePrice, GarageYrBlt had a coefficiient of 0.56 and YearRemodAdd had scored a 0.55. Knowing how important these features are we want to make as much use with them as possible. The houses in our data start from 1871 and were built up until 2010. We can bin houses into sequences of 20, which will give us 7 different bins for a feature built in a year to fall into. This will allow us to differentiate the different times/era houses, garages and remodeling took place. 

```{r warning=FALSE, message=FALSE}
# every 20 years create a new bin
# 7 total bins
# min year is 1871, max year is 2010!
year.map = function(col.combined, col.name) {
  for (i in 1:7) {
    year.seq = seq(1871+(i-1)*20, 1871+i*20-1)
    idx = which(df.combined[,col.combined] %in% year.seq)
    df.categoric[idx,col.name] = i
  }
  return(df.categoric)
}

```

```{r warning=FALSE, message=FALSE}
# we'll c
df.categoric['GarageYrBltBin'] = 0
df.categoric <- year.map('GarageYrBlt', 'GarageYrBltBin')
df.categoric['YearBuiltBin'] = 0
df.categoric <- year.map('YearBuilt','YearBuiltBin')
df.categoric['YearRemodAddBin'] = 0
df.categoric <- year.map('YearRemodAdd', 'YearRemodAddBin')
```

Now that we 3 new columns that generalize what year a house, garage and remodeling took place we'll need to one-hot encode these columns so that each bin is turned into a binary column. The new columns are ranged from 1-7, but we don't know what specific order they follow so it will benefit us more to give each bin its own column. 

```{r warning=FALSE, message=FALSE}
bin.cols <- c('GarageYrBltBin', 'YearBuiltBin', 'YearRemodAddBin')

for (col in bin.cols) {
  df.categoric <- cbind(df.categoric, model.matrix(~.-1, df.categoric[col]))
}

# lets drop the orginal 'GarageYrBltBin', 'YearBuiltBin', 'YearRemodAddBin' from our dataframe
df.categoric <- df.categoric[,!names(df.categoric) %in% bin.cols]
```

We're finally done clearning, manipulating and adding features to our data. We'll combine the numeric and categoric dataframes into one, which we will use to build or models on. 

```{r warning=FALSE, message=FALSE}
df <- cbind(df.numeric, df.categoric)
```

What does our distribution of housing prices look like?

```{r warning=FALSE, message=FALSE}
require(WVPlots)
y.true <- train$SalePrice[which(!1:1460 %in% idx.outliers)]

qplot(y.true, geom='density') +# +(train, aes(x=SalePrice)) +
  geom_histogram(aes(y=..density..), color='white', 
                 fill='lightblue', alpha=.5, bins = 60) +
  geom_line(aes(y=..density..), color='cornflowerblue', lwd = 1, stat = 'density') + 
  stat_function(fun = dnorm, colour = 'indianred', lwd = 1, args = 
                  list(mean(train$SalePrice), sd(train$SalePrice))) +
  scale_x_continuous(breaks = seq(0,800000,100000), labels = dollar) +
  scale_y_continuous(labels = comma) +
  theme_minimal() +
  annotate('text', label = paste('skewness =', signif(skewness(train$SalePrice),4)),
           x=500000,y=7.5e-06)

qqnorm(train$SalePrice)
qqline(train$SalePrice)
```

We can see from the histogram and the quantile-quantile plot that the distribution of sale prices is right-skewed and does not follow a normal distribution. Lets make a log-transformation and see how our data looks

```{r warning=FALSE, message=FALSE}
y_train <- log(y.true+1)

qplot(y_train, geom = 'density') +
  geom_histogram(aes(y=..density..), color = 'white', fill = 'lightblue', alpha = .5, bins = 60) +
  scale_x_continuous(breaks = seq(0,800000,100000), labels = comma) +
  geom_line(aes(y=..density..), color='dodgerblue4', lwd = 1, stat = 'density') + 
  stat_function(fun = dnorm, colour = 'indianred', lwd = 1, args = 
                  list(mean(y_train), sd(y_train))) +
  #scale_x_continuous(breaks = seq(0,800000,100000), labels = dollar) +
  scale_y_continuous(labels = comma) +
  theme_minimal() +
  annotate('text', label = paste('skewness =', signif(skewness(y_train),4)),
           x=13,y=1) +
  labs(x = 'log(SalePrice + 1)') 

qqnorm(y_train)
qqline(y_train)
```

```{r}
paste('The dataframe has', dim(df)[1], 'rows and', dim(df)[2], 'columns')
```

We've manipulated a lot of variables and added many features, 345 to be exact, in our dataset and their could potentially be some variables that won't give our data any value when we're modeling. Some of these features may have become zero-variance predictors, such that a few samples may have an insignificant influence on the model. These *near-zero-variance* may cause overfitting or will prevent our model from generalizing over the data at a more sufficient rate. The package `caret` offers a function `nearZeroVar`, which checks the frequency of the most common value over the second most frequent value, which would be closer to 1 for well-behaved predictors and very large for highly-unbalanced features. It also checks the number of unique values divided by the *n* number of samples which will approach zero as the level of detail in the feature increases. We can remove all of the near-zero-variance variables from out dataframe.

```{r warning=FALSE, message=FALSE}
nzv.data <- nearZeroVar(df, saveMetrics = TRUE)

# take any of the near-zero-variance perdictors
drop.cols <- rownames(nzv.data)[nzv.data$nzv == TRUE]

df <- df[,!names(df) %in% drop.cols]

paste('The dataframe now has', dim(df)[1], 'rows and', dim(df)[2], 'columns')
```


Our new distribution with the transformed data follows much closer to a normal distribution and we can verify this from the quantile-quantile plot. Now that our independent and dependent variables satisfy the assumption of normality we can begin building our model.

**XGBoost** fits shallow regression trees to our data and then additional trees to the residuals, we'll repeat this process for 30000 rounds so that our model has learned from the data as much as possible without overfitting. XGBoost is a **Gradient Boosted Method (GBM)**, which is an ensemble learning method that uses a very large number of decision trees, which are typically weak learners and combines them into one final prediction. For gradient boosted trees, the new trees added to the model are the weak learnings where $$F_0=0, \quad F_t(x)=F_{t-1}(x)+h(x)$$ suchat that $F(x)$ is the entire model after time $t-1$ and $h(x)$ is the new weak learner (new tree) that will be added to our model. In particular for XGBoost we will minimize the following objective function at time *t* such that $$Obj(F_t)=L(F_{t-1}+F_t) + \Omega(F_t)$$ such that $L(F_t)$ is our loss function and $\Omega(F_t)$ is the regularization function. In the XGBoost package the regularization function is computed such that $$\Omega(F_t) = \gamma T+\frac{1}{2}\sum_{j=1}^Tw_j^2$$ where $w_j$ is the score of the $j^{th}$ leaf, $\lambda$ is the leaf weight penalty, $\gamma$ is the tree size penalty and *T* represents the number of leaves in the tree. When we call `xgb.train` from the `XGBoost` package we will choose which parameters minimize the loss in our model using the `xgb.csv` by searching through a grid of parameters. The parameters we'll be interested in *tuning* is:

$\eta$ - Shrinkage term  
$\gamma$ - Tree size penalty  
$\lambda$ - $L^2$ leaf node weight penalty  
*Max Depth* -  The maximum depth of each tree  
*Subsample* -  proportion of data for ”bagging” i.e. randomly sampling with replacement  
*Min Child Weight* -  The minimum weight that a node can have


One limitation to using GBM's and XGBoost in particular is its inability to extrapolate and because of this we can use linear model to better predict any sale prices outside the range of prices given in our training set.

We know there exists multicollinearity in our dependent variables as area features are a determinant for a house having a certain amount of rooms, we know that the garage variables have heavy dependency and the list goes on. Due to this a simple linear regression model will not be of much help to predict accurate sale prices, which is why we can make use of both **ridge** and **lasso** regression. The Ridge penalty is known to shrink the coefficients of correlated predictors towards each other with the use of the $\ell^2$-norm, while the lasso tends to pick one of them and discard the others by construct of the $\ell^1$-norm. Adding the $\ell^1$ and $\ell^2$ penalties give a nice contstraint on the coefficients of our model and solve our problem of having collinear variables. Another useful property of these regression models is their ability to *extrapolate*, i.e. they can predict house prices that our outside of the range of prices we were given in our training set. Along with Ridge and Lasso we will use  elastic-net regularization which can be derived such that, $$min_{\beta_0,\beta}\frac{1}{N}\sum_{n=1}^Nw_il(y_i,\beta_0+\beta^Tx_i) + \lambda[(1-\alpha)||\beta||^2_2/2+\alpha||\beta||_1]$$
where $l(y,n)$ is the negative log-likelihood contribution at observation *i* and $\lambda$ is a shrinkage parameter, when $\lambda=0$, no shrinkage will be performed and and as $\lambda$ increases the coefficients are shrunken at a larger rate. Elastic-net has a great property such that when we call `glmnet`, the equation is controlled by $\alpha$, which leverages the gap between ridge and lasso. When $\alpha = 0$ we have $$min_{\beta_0,\beta}\frac{1}{N}\sum_{n=1}^Nw_il(y_i,\beta_0+\beta^Tx_i) + \lambda[(1-\alpha)||\beta||^2_2/2]$$ which is our equation for ridge and when $\alpha=1$ we have, $$min_{\beta_0,\beta}\frac{1}{N}\sum_{n=1}^Nw_il(y_i,\beta_0+\beta^Tx_i) + \lambda[\alpha||\beta||_1]$$ which is our requation for lasso. 

Once we have made predictions with all of our models we can take the mean of the 4 predictions and let that be our final prediction. 

To better estimate the test error of a predictive model we'll use **K-fold cross-validation**. For K-fold cross-validaiton the observations in our data are split into K partitions, then our model is trained on the K – 1 partitions, where the test error is predicted on the left out partition k. We repeat this process for k = 1,2…K times and we will average the results. To score how well our model predicts output we will take the **root mean squared error (RMSE)**. The RMSE can be computed such that $$RMSE = \sqrt{\frac{1}{N}\sum_{i=1}^N(y_i-\hat{y}_i)^2}$$ The RMSE is one of the most used error metrics used in practice and gives us a good score to determine how well our model trains on new data, specifically through the use of cross-validation. 



```{r warning=FALSE, message=FALSE}
x_train <- df[1:1456,]

x_test <- df[1457:nrow(df),]

dtrain <- xgb.DMatrix(as.matrix(x_train), label = y_train)
dtest <- xgb.DMatrix(as.matrix(x_test))
```


`expand.grid` is your best friend when it comes to parameter tuning. The code below is a baseline approach, however, you should have a set of values you would like to try and run yourself. I won't run this in the kernel as this would take too long.
`xgb.cv` is also a very useful function that will show you how well your model generalizes to testing data. To avoid overfit `early_stopping_rounds` is a great way to find at which point your model begins to overfit, which will change based on the paramters used.
```{r warning=FALSE, message=FALSE}

cv.ctrl <- trainControl(method = "repeatedcv", repeats = 1,number = 4, 
                        allowParallel=T)

    xgb.grid <- expand.grid(nrounds = 750,
        eta = c(0.01,0.005,0.001),
        max_depth = c(4,6,8),
        colsample_bytree=c(0,1,10),
        min_child_weight = 2,
        subsample=c(0,0.2,0.4,0.6),
        gamma=0.01)
set.seed(45)
    #xgb_tune <- train(as.matrix(x_train),
    #        y_train,
    #        method="xgbTree",
    #        trControl=cv.ctrl,
    #        tuneGrid=xgb.grid,
    #        verbose=T,
    #        metric="RMSE",
    #        nthread =3)

xgb_params <- list(
  booster = 'gbtree',
  objective = 'reg:linear',
  colsample_bytree=1,
  eta=0.005,
  max_depth=4,
  min_child_weight=3,
  alpha=0.3,
  lambda=0.4,
  gamma=0.01, # less overfit
  subsample=0.6,
  seed=5,
  silent=TRUE)

#xgb.cv(xgb_params, dtrain, nrounds = 5000, nfold = 4, early_stopping_rounds = 500)
```

```{r warning=FALSE, message=FALSE}
bst <- xgb.train(xgb_params,dtrain, nrounds = 1000)#, early_stopping_rounds = 300, watchlist = list(train=dtrain))
```

**Running XGBoost in a Kaggle Kernel is computationally heavy with large datasets so here I am just using `nrounds = 1000` as an example, but determining the right amount of `nrounds` for you model relies heavily on $\eta$. I would recammend `nrounds = 10000` and `eta = 0.01` as a starter.
A very effective way of choosing optimal parameters is with the use of `expand.grid`, however, this can take quite a while to let many models train over your data, I would strongly recammend reading up on the XGBoost and GBM's in general.**

This is a very good reference: https://www.slideshare.net/ShangxuanZhang/kaggle-winning-solution-xgboost-algorithm-let-us-learn-from-its-author

XGBoost has an extremely useful function called `xgb.plot.importance`, which plots the **feature importance** of the model.  

Feature importance is computed by averaging the gain of each feature for all split and all trees in the model. We can take a look at the 10 most important features used in our model. 

```{r warning=FALSE, message=FALSE}
#model.names <- dimnames(dtrain)[[2]]

#importance_matrix <- xgb.importance(model.names, model = bst)

#xgb.plot.importance(importance_matrix[1:10])
```

Having the opportunity to see which features are most important allows us to give an emphasis on making use of specific features to improve our model even more. We can also look deeper in the model to see which features are adding little to no value for predictions. Certainly somethig for us to consider. 

## How well did our model do?

```{r warning=FALSE, message=FALSE}
rmse_eval <- function(y.true, y.pred) {
  mse_eval <- sum((y.true - exp(y.pred)-1)^2) / length(y.true)
  return(sqrt(mse_eval))
}

#y_pred.xgb <- predict(bst, dtrain)
#rmse_eval(y.true, y_pred.xgb)
```

How accurate is the RMSE that our XGBoost produced? Lets see how well the linear models train on the data for comparison. 


Lucky for us the `glmnet` package offers a `cv.glmnet` function which will run cross validation for us and pick the $\lambda$ the produces the smallest error rate. We can call on the same glmnet function and use a ridge regression when $\alpha=0$, lasso when $\alpha=1$ and elastic-net when $\alpha \in (0,1)$. 

```{r warning=FALSE, message=FALSE}
# ridge, lasso, elasticnet

# these won't run on kaggle uncomment these in your own R-env
#glm.cv.ridge <- cv.glmnet(as.matrix(x_train), y_train, alpha = 0)
#glm.cv.lasso <- cv.glmnet(as.matrix(x_train), y_train, alpha = 1)
#glm.cv.net <- cv.glmnet(data.matrix(x_train), y_train, alpha = 0.001)

# use the lamdba that minimizes the error
#penalty.ridge <- glm.cv.ridge$lambda.min
#penalty.lasso <- glm.cv.lasso$lambda.min
#penalty.net <- glm.cv.net$lambda.min
```

Lets train our model using our chosen $\lambda$'s

```{r warning=FALSE, message=FALSE}
#glm.ridge <- glmnet(x = as.matrix(x_train), y = y_train, alpha = 0, lambda = penalty.ridge )
#glm.lasso <- glmnet(x = as.matrix(x_train), y = y_train, alpha = 1, lambda = penalty.lasso)
#glm.net <- glmnet(x = as.matrix(x_train), y = y_train, alpha = 0.001, lambda = penalty.net)

#y_pred.ridge <- as.numeric(predict(glm.ridge, as.matrix(x_train)))
#y_pred.lasso <- as.numeric(predict(glm.lasso, as.matrix(x_train)))
#y_pred.net <- as.numeric(predict(glm.net, as.matrix(x_train)))
```

```{r warning=FALSE, message=FALSE}
#rmse_eval(y.true, y_pred.ridge)
#rmse_eval(y.true, y_pred.lasso)
#rmse_eval(y.true, y_pred.net)
```

Here we have the RMSE for elastic-net, lasso and ridge are all very similar yet not as low as the RMSE we got from XGBoost, *however*, ridge, lasso, and elastic-nets ability to extrapolate saleprices beyond the range we were given in the training set will be crucial for predicting more accurately towards unseen data (future houses). If we combine all 4 of these predictions on our test data by averaging all of our predictions we can get what we hope to be an accurate submission. **Note**: our predictions are in logarithmic form since we made a transformation due to the skewness in the Saleprices so we will need to tranform our predictions back into their original form by taking $f(x) = e^x-1$. 


```{r warning=FALSE, message=FALSE}
#y_pred.ridge <- as.double(predict(glm.ridge, as.matrix(x_test)))
#y_pred.lasso <- as.double(predict(glm.lasso, as.matrix(x_test)))
#y_pred.net <- as.double(predict(glm.net, as.matrix(x_test)))

#y_pred.ridge <- as.double(exp(y_pred.ridge) - 1)
#y_pred.lasso <- as.double(exp(y_pred.lasso) - 1)
#y_pred.net <- as.double(exp(y_pred.net) - 1)

#y_pred.xgb <- as.double(predict(bst, dtest))
#y_pred.xgb <- as.double(exp(y_pred.xgb) - 1)

# take the average of our predictions for our ensemble
#y_pred <- (y_pred.xgb + y_pred.ridge + y_pred.lasso + y_pred.net)/4.0
#head(y_pred)
```

Uncomment everything above for modeling, this ensemble placed top 20%, but can certainly use improvements. Please give a thumbs up if this was helpful in any way it will encourage me to make more kernels!