```{r, lib-read, results='hide', message=FALSE, include=FALSE, echo=FALSE}
setwd('.')
### Compuational Setup
# libraries
library(knitr)
library(kableExtra)
library(tidyverse)
library(olsrr)
library(gridExtra)
library(caret)
library(multcomp)
library(Hmisc)

# set a random seed for repodicibility
set.seed(123)

# helper code
source('./helper/visual.R')
source('./helper/data_munging.R')
source('./helper/performance.R')

load('./data/data.RData')

# filter to analysis data
train <- train %>% filter(GrLivArea < 4000)
```


# Introduction

What is the price of a home in Ames, Iowa? Our inaugural project for the Statistical Foundations for Data Science course in the Southern Methodist University Master of Science in Data Science (MSDS) program was to compete in an online Kaggle competition utilizing linear regression techniques we've learned in this course to date. Our team elected to use R as the preferred analysis platform under the consensus that it has more broad applicability for use in industry, including data gathering and wrangling, in addition to advanced visualization and analytic tools. The project objective was to apply various predictive models in order to assess the suitability of our parameter selections for determining the sales prices of homes in Ames. The measures of accuracy were applied in terms of the Root Mean Square Error, or RMSE, as well as other comparison models such as cross-validation, and the adjusted R-squared. Our approach outlined in this research document is limited in that we were not permitted to use more advanced algorithms we will be exposed to later in the MSDS program; rather, in conjunction with the aforementioned linear regression techniques, we were directed to apply the exploratory data analysis and data cleaning methods we have learned, which will surely be of use to us in our future personal and academic endeavors.

For more information see: https://github.com/sjmiller8182/RegressionHousingPrices

# Ames, Iowa Data Set

The Ames, Iowa Data Set describes the sale of individual residential properities from 2006-2010 in Ames, Iowa \cite{Cock}. The data was retreved from the dataset hosting site Kaggle, where it is listed under a machine learning competition named \href{https://www.kaggle.com/c/house-prices-advanced-regression-techniques/overview}{\textit{House Prices: Advanced Regression Techniques}} \cite{Kaggle2016}. The data is comprised of 37 numeric features, 43 non-numeric features and an observation index split between a training set and a testing set, which contain 1460 and 1459 observations, respectively. The response variable (`SalePrice`) is only provided for the training set. The output of a model on the test set can be submitted to the Kaggle competition for scoring the performance of the model in terms of RMSE. The first analysis models property sale prices (`SalePrice`) as the response of living room area (`GrLivArea`) of the property and neighborhood (`Neighborhood`) where it is located. In the second analysis, variable selection techniques are used to determine which explanatory varaibles are associated with `SalePrice` to find a predictive model.

# Analysis Question I

## Question of Interest

Century 21 has commissioned an analysis of this data to determine how the sale price of property is related to living room area of the property in the Edwards, Northwest Ames, and Brookside neighborhoods of Ames, IA.

## Modeling

Linear regression will be used to model sale price as a response of the living room area. From the initial exploratory data analysis, it was determined that sale prices should be log-transformed to meet the model assumptions for linearity (see section \ref{appendix:linearity}), thus improving our models fit and reducing standard error. Additionally, two observations were removed as they appeared to be from a different population than the other observations in the dataset (see section \ref{appendix:infleu-points}); therefore, analysis only considers properties with living rooms less than 3500 sq. ft. in area.

We will use extra sums of square (ESS) tests to determine if a neighborhood should be added in the model.
We start with the logarithm of sale price as the response of living room area and build up a model with neighborhood.
The equations (1-3) below show the models considered. 
The Edwards neightborhood is used for reference.

**Base Model**

\begin{equation}
\mu \lbrace log(SalePrice) \rbrace = \hat{\beta_0} + \hat{\beta_1} (LivingRoomArea) (\#eq:Base)
\end{equation}

**Additive Model**

\begin{equation}
\mu \lbrace log(SalePrice) \rbrace = \hat{\beta_0} + \hat{\beta_1} (LivingRoomArea) +  \hat{\beta_2} (Brookside) + \hat{\beta_3} (NorthwestAmes) (\#eq:Additive)
\end{equation}

**Interaction Model**

\begin{align}
\mu \lbrace log(SalePrice) \rbrace = \hat{\beta_0} + \hat{\beta_1} (LivingRoomArea) +  \hat{\beta_2} (Brookside) + \hat{\beta_3} (NorthwestAmes) + \nonumber\\
\hat{\beta_3} (Brookside)(LivingRoomArea) + \hat{\beta_4} (NorthwestAmes)(LivingRoomArea) (\#eq:Interaction)
\end{align}

\newpage

### ESS Tests Between Models

```{r, ESS-setup, results='hide', echo=FALSE}
# filter data for analysis 1
train.mod <- train %>% 
  filter(Neighborhood %in% c("Edwards", "BrkSide", "NAmes"))
train.mod$Neighborhood <- factor(train.mod$Neighborhood, levels = c("Edwards", "BrkSide", "NAmes"))


# create dummy variables with Neighborhood == 'Edwards' as reference
train.mod <- get.dummies(train.mod, "Neighborhood", reference = 'Edwards')

# interaction model formula
model.interaction.formula = log(SalePrice) ~ (GrLivArea) + 
     Neighborhood_BrkSide + 
     Neighborhood_NAmes +
     (GrLivArea) * Neighborhood_BrkSide + 
     (GrLivArea) * Neighborhood_NAmes
# additive model formula
model.add.formula = log(SalePrice) ~ (GrLivArea) + 
     Neighborhood_BrkSide + 
     Neighborhood_NAmes 
# base model formula
model.base.formula = log(SalePrice) ~ (GrLivArea)
```

**Comparison: Base Model to Additive Model**

The following ESS test provides convincing evidence that the addition of the additive neighborhood terms
is an improvement over the base model (p-value < 0.0001).

```{r, ESS1, echo=FALSE}
# fit models
model <- lm(formula = model.add.formula, data = train.mod)
model.reduced <- lm(formula = model.base.formula, data = train.mod)
# ESS test on models
anova(model.reduced, model)
```

**Comparison: Additive Model to Interaction Model**

The following ESS test provides convincing evidence that the addition of the interaction neighborhood terms
is an improvement over the additive model (p-value < 0.0001).

```{r, ESS2, echo=FALSE}
# fit models
model <- lm(formula = model.interaction.formula, data = train.mod)
model.reduced <- lm(formula = model.add.formula, data = train.mod)
# ESS test on models
anova(model.reduced, model)
```

Based on the two ESS tests, the interaction model appears to be significant;
thus, we will continue with the interaction model.

\newpage

## Model Assumptions Assessment

The following assessments for model assumptions are made based on Figure \@ref(fig:diag-plots) and Figure \@ref(fig:scatter-plots):

* The residuals of the model appear to be approximately normally distrubited based on the QQ plot of the residuals and histogram of the residuals, suggesting the assumption of normality is met.
* No patterns are evident in the scatter plots of residuals and studentized residuals vs predicted value, suggesting the assumption of constant variance is met.
* While some observations appear to be influential and have high leverage, removing these observations does not have a significant impact on the result of the model fit.
* Based on the scatter plot of the log transform of `SalePrice` vs `GrLivArea`, it appears that a linear model is reasonable (see section \ref{appendix:linearity}).

The sampling procedure is not known. We will assume the independence assumption is met.

```{r, diag-plots, echo=FALSE, fig.width=5, fig.height=5, fig.align='center', fig.cap='Diagnostic Plots', out.width = '45%', fig.pos="htbp", fig.show = 'hold'}
# option 'htbp' is used to lock the position of the image
# option 'hold' is used to arrange images side-by-side

# create plots of residuals
basic.fit.plots(train.mod, model)
# create leverage / outlier plot
ols_plot_resid_lev(model)
```


## Comparing Competing Models

The three models were trained and validated on the training dataset using 10-fold cross validation. The table below summerizes the performance of the models with RMSE, adjusted $R^2$, and PRESS. These results show that the interaction model produces the best performance, which is consistent with the result of the ESS test.

```{r, cross-validation, results='hide', message=FALSE, include=FALSE, echo=FALSE}
## cross validate the full model

# Set up k-fold cross-validation
train.control <- trainControl(method = "cv", number = 10)
# Train the model
model.cv <- train(model.interaction.formula, 
                    data = train.mod,
                    method = 'lm',
                    trControl = train.control)
# print model summary
model.cv

# get the CV results
res <- model.cv$results

# get cross-validated PRESS statistic
PCV <- PRESS.cv(model.cv)



## cross validate the reduced model

# Set up k-fold cross-validation
train.control <- trainControl(method = "cv", number = 10)
# Train the model
model.reduced.cv <- train(model.add.formula, 
                    data = train.mod,
                    method = 'lm',
                    trControl = train.control)
# print model summary
model.reduced.cv

# get the CV results
res.red <- model.reduced.cv$results

# get cross-validated PRESS statistic
PCV.red <- PRESS.cv(model.reduced.cv)

## cross validate the base model

# Set up k-fold cross-validation
train.control <- trainControl(method = "cv", number = 10)
# Train the model
model.base.cv <- train(model.base.formula, 
                    data = train.mod,
                    method = 'lm',
                    trControl = train.control)
# print model summary
model.base.cv

# get the CV results
base.red <- model.base.cv$results

# get cross-validated PRESS statistic
PCV.base <- PRESS.cv(model.base.cv)
```

```{r, echo=FALSE}
# print accuracy metrics to md table
kable(data.frame('Model' = c('Interaction Model', 'Additive Model', 'Base Model'), 
                 'RMSE'=c(res$RMSE, res.red$RMSE, base.red$RMSE),
                 'CV Press'=c(PCV, PCV.red, PCV.base),
                 'Adjused R Squared'=c(res$Rsquared, res.red$Rsquared, base.red$Rsquared)),
      "latex", booktabs = T)  %>%
  kable_styling(position = "center")
```

\newpage

## Parameters

The following table summerizes the parameter estimates for the interaction model.

```{r, echo=FALSE}
# extract the model estimates from the model summary
sm <- summary(model)
sm.coe <- sm$coefficients
# get the CIs for the coefficients
model.conf <- confint(model)

# print model estimates to md / latex table
kable(data.frame('Parameter' = c('Intercept', 'GrLivArea', 
                                 'Neighborhood_BrkSide', 'Neighborhood_NAmes', 
                                 'GrLivArea:Neighborhood_BrkSide', 'GrLivArea:Neighborhood_NAmes '), 
                 'Estimate'=c(sm.coe[[1]],sm.coe[[2]],sm.coe[[3]],sm.coe[[4]],sm.coe[[5]],sm.coe[[6]]),
                 'CI Lower' = c(model.conf[[1]],model.conf[[2]],model.conf[[3]],
                                model.conf[[4]],model.conf[[5]],model.conf[[6]]),
                 'CI Upper' = c(model.conf[[1,2]],model.conf[[1,2]],model.conf[[3,2]],
                                model.conf[[4,2]],model.conf[[5,2]],model.conf[[6,2]])),
      "latex", booktabs = T)  %>%
  kable_styling(position = "center")
```

Where `Intercept` is $\beta_0$, `GrLivArea` is $\beta_1$, `Neighborhood_BrkSide` is $\beta_2$, `Neighborhood_NAmes` is $\beta_3$, `GrLivArea:Neighborhood_BrkSide` is $\beta_4$, and `GrLivArea:Neighborhood_NAmes` is $\beta_5$

## Model Interpretation

We estimate that for increase in 100 sq. ft., there is associated multiplicative increase in median price of

* 1.055 for the Edwards neighborhood with a 95% confidence interval of [1.044 , 1.066]
* 1.033 for the Northwest Ames neighorhood with a 95% confidence interval of [1.026 , 1.040]
* 1.077 for the Brookside neighorhood with a 95% confidence interval of [1.063 , 1.090]

Since the sampling procedure is not known and this is an observational study, the results only apply to this data.

## Conclusion

In response to the analysis commissioned by Century 21, the log transform of property sale price was modeled as a linear response to the property living room area for residential properties in Ames, IA. It was determined that it was necessary to include interaction terms to allow for the influence of neighborhood on sale price. Based on the model, there is strong evidence of an associated multiplicative increase in median sale price for an increase in living room area (p-vlue < 0.0001, overall F-test).

\newpage

# Analysis Question II

## Question of Interest

Century 21 has commissioned a second analysis using the same dataset for the creation of a very predictive model of `SalePrice`. The analysis will be expanded to include as many of the 80 total features as required to determine the sale price of residential properties across all neighborhoods of Ames, Iowa, beyond only the three - Edwards, Northwest Ames, and Brookside - previously commissioned for analysis.

## Modeling

Through analyzing our variable selection and cross-validation processes - along with our nascant domain knowledge of residential real estate - we ultimately arrived at a multiple linear regression model featuring 11 linear predictor variables and two interaction terms. Specifically, our variable selection process included direct analysis of a correlation plot and a correlation matrix as well as performing forward selection, backward elimination, and stepwise regression. 

Regarding missing data, we imputed NA values for 19 variables using a combination of the data dictionary provided by Century 21 as well as our domain knowledge. After building models with and without transformations applied to variables, we noted no significaznt difference in variable selection from our selection process so elected to use non-transformed predictor variables. We did, however, use the log-transformed `SalePrice` applied in the first analysis.

**Forward Selection**

Forward selection is a variable selection methodology that begins with a constant mean and adds explanatory variables one-by-one until no further additonal predictor variables significantly improve the model's fit. This employess the "F-to-enter" method from the extra-sum-of-squares F-statistic. This was the first method we employed. For this process, we provided the test a starting model with no predictor variables and a model from which terms can be selected, which included all predictor variables available. The process worked forward with selecting one parameter. The suggested model shown in section \ref{appendix:forSelection}.


```{r, echo=F, results='hide'}
# Initial Model - Forward Selection
model2.forward.Start <- lm(log(SalePrice)~1,data = train)

# All Variables Model - Forward Selection
model2.Allvar <- lm(log(SalePrice) ~ OverallQual + GrLivArea + Neighborhood + 
                      BsmtFinSF1 + MSSubClass + OverallCond + YearBuilt + GarageCars + 
                      TotalBsmtSF + SaleCondition + LotArea + MSZoning + Functional + 
                      CentralAir + KitchenQual + Condition1 + FireplaceQu + BsmtExposure + 
                      BsmtFullBath + ScreenPorch + Exterior1st + YearRemodAdd + 
                      GarageQual + WoodDeckSF + OpenPorchSF + Street + LotConfig + 
                      LotFrontage + Foundation + Heating + KitchenAbvGr + EnclosedPorch + 
                      HalfBath + FullBath + MasVnrType + BsmtFinSF2 + HeatingQC + 
                      GarageArea + SaleType + ExterCond + PoolArea + BsmtFinType1 + 
                      GarageYrBlt + Electrical + `3SsnPorch` + LowQualFinSF, data = train
                    )

#### Forward Selection, first pass
model2.Forward <- stepAIC(model2.forward.Start, direction = "forward", trace = F, scope = formula(model2.Allvar))

summary(model2.Forward)
model2.Forward$anova

# Suggested output from Forward Selection:
final.Forward.Model <- lm(log(SalePrice) ~ OverallQual + GrLivArea + Neighborhood + 
                            BsmtFinSF1 + MSSubClass + OverallCond + YearBuilt + GarageCars + 
                            TotalBsmtSF + SaleCondition + LotArea + MSZoning + Functional + 
                            CentralAir + KitchenQual + Condition1 + FireplaceQu + BsmtExposure + 
                            BsmtFullBath + ScreenPorch + Exterior1st + YearRemodAdd + 
                            GarageQual + WoodDeckSF + OpenPorchSF + Street + LotConfig + 
                            LotFrontage + Foundation + Heating + KitchenAbvGr + EnclosedPorch + 
                            HalfBath + FullBath + MasVnrType + BsmtFinSF2 + HeatingQC + 
                            GarageArea + SaleType + ExterCond + PoolArea + BsmtFinType1 + 
                            GarageYrBlt + Electrical + `3SsnPorch` + LowQualFinSF, data = train
                          )
```


**Backward Elimination**

Backward elimination is a variable selection methodology that begins with all possible predictor variables and works backward, eliminating variables using all possible combinations until only the best for the fit are provided. This employess the "F-to-remove" method from the extra-sum-of-squares F-statistic. For this process, we provided the test a model with all available predictor variables from which insignificant variables were eliminated. The suggested model shown in section \ref{appendix:backSelection}.

```{r, echo=F, results='hide'}
# Backward Elimination
model2.Backward <- stepAIC(model2.Allvar, direction = "backward", trace = F, scope = formula(model2.forward.Start))

summary(model2.Backward)
model2.Backward$anova

# Suggested output from Backward Elimination:
final.Backward.Model <- lm(log(SalePrice) ~ MSSubClass + MSZoning + LotFrontage +
                              LotArea + Street + LotConfig + LandSlope + Neighborhood +
                              Condition1 + OverallQual + OverallCond + YearBuilt + YearRemodAdd +
                              RoofMatl + Exterior1st + MasVnrType + ExterCond + Foundation +
                              BsmtExposure + BsmtFinType1 + BsmtFinSF1 + BsmtFinSF2 + BsmtUnfSF +
                              Heating + HeatingQC + CentralAir + Electrical + `1stFlrSF` +
                              `2ndFlrSF` + LowQualFinSF + BsmtFullBath + FullBath + HalfBath +
                              KitchenAbvGr + KitchenQual + Functional + FireplaceQu + GarageYrBlt +
                              GarageCars + GarageArea + GarageQual + WoodDeckSF + OpenPorchSF +
                              EnclosedPorch + ScreenPorch + PoolArea + SaleType + SaleCondition, data = train
                           )
```

**Stepwise Regression**

Stepwise regression is a variable selection methodology that performs one step of forward selection for each step of backward elimination. The steps are repeated, concurrently, until no further predictor variables can be added or removed. This is the third model approach we used. The suggested model shown in section \ref{appendix:stepWSelection}.

```{r, echo=F, results='hide'}
# Stepwise Selection
model2.Stepwise <- stepAIC(model2.Allvar, direction = "both", trace = F)

summary(model2.Stepwise)
model2.Stepwise$anova

# Suggested output from Stepwise Regression:
final.Stepwise.Model <- lm(log(SalePrice) ~ OverallQual + GrLivArea + Neighborhood + 
                              BsmtFinSF1 + MSSubClass + OverallCond + YearBuilt + GarageCars + 
                              TotalBsmtSF + SaleCondition + LotArea + MSZoning + Functional + 
                              CentralAir + KitchenQual + Condition1 + FireplaceQu + BsmtExposure + 
                              BsmtFullBath + ScreenPorch + Exterior1st + YearRemodAdd + 
                              GarageQual + WoodDeckSF + OpenPorchSF + Street + LotConfig + 
                              LotFrontage + Foundation + Heating + KitchenAbvGr + EnclosedPorch + 
                              HalfBath + FullBath + MasVnrType + BsmtFinSF2 + HeatingQC + 
                              GarageArea + SaleType + ExterCond + PoolArea + BsmtFinType1 + 
                              GarageYrBlt + Electrical + `3SsnPorch` + LowQualFinSF, data = train
                           )

```

\newpage

**Custom Variable Selection** 

To develop the custom model, we employed a combination of a correlation matrix for quantitative data, analysis of the summarization of the suggested model from stepwise selection, and through direct analysis of the pairs plots. As previously mentioned, our final model included 11 linear terms and two interaction terms. We removed all variables suggested to be removed by the stepwise regression and backward elimination tests, then reprocessed the updated models until forward selection, backward elimination, and stepwise regression were in agreeance with respect to the linear terms. Once this trial-and-error process was completed, we added interaction terms based on domain knowledge and re-applied the forward selection, backward elimination, and stepwise regression methods until only significant terms - both linear and interactive - remained. We then used graphical analysis to visually confirm interaction between the interactive terms remaining. The custom model shown in section \ref{appendix:customSelection}.

```{r, echo=F, results='hide'}

fitFull.all4.Final <- lm(log(SalePrice) ~ MSSubClass + MSZoning + LotFrontage + 
                            LotArea + Street + Neighborhood + 
                            Condition1 + OverallQual + OverallCond + YearBuilt + YearRemodAdd + 
                            ExterCond + Foundation + 
                            BsmtExposure + BsmtFinType1 + BsmtFinSF1 + BsmtFinSF2 + BsmtUnfSF + 
                            HeatingQC + CentralAir + Electrical + `1stFlrSF` + 
                            `2ndFlrSF` + LowQualFinSF + BsmtFullBath + FullBath + HalfBath + 
                            KitchenAbvGr + KitchenQual + Functional + FireplaceQu + GarageYrBlt + 
                            GarageCars + GarageArea + GarageQual + WoodDeckSF + OpenPorchSF + 
                            EnclosedPorch + ScreenPorch + SaleType + SaleCondition, data = train)

```

## Model Assumption Assessment

The assumption assessment plots were similar for all four models. The assumption assessment plots and discussion for the custom model are provided here with Figure \@ref(fig:custom-assumptions). The assumption assessment plots for the other three models are provided for reference in section \ref{appendix:assessmentPlots}.

Based on the diagnostic plots below, the custom model appears to reasonably meet the assumptions of linear regression. The standardized residuals do not appear to exhibit a discernible pattern, indicating constant variance along the regression, or homoscedasticity. While there are some outliers, this does not appear to be an egregious violation. Based on the QQ plot, there is a small level of deviation on the ends of the distribution of the errors, but for the most part, the errors adhere to normality. The sample size should be sufficient to protect against this non-normality. Based on the standardized residuals vs. leverage plot, only a few values have high leverage and are outlying. However, these violations do not appear to be egregious.


```{r, custom-assumptions, echo=F, warning=F, fig.width=5, fig.height=5, fig.align='center', fig.cap='Custom Assumption Assessment Plots', out.width = '45%', fig.pos="htbp", fig.show = 'hold'}

# create plots of residuals
basic.fit.plots(train, fitFull.all4.Final)
# create leverage / outlier plot
ols_plot_resid_lev(fitFull.all4.Final)
```

## Comparing Competing Models

While the models from forward, backward, and stepwise selection produce higher adjusted $R^2$ values on the training data, the yield much higher errors when applied to the Kaggle test set. These selection methods appear to be overfitting to the training data, thus fail to generalize to the Kaggle test set. Undisputedly, the custom model outperformed the model built strictly on the output of the forward selection, backward elimination, and stepwise regression variable selection procedures when applied to a new dataset.

```{r, echo=F, results='hide', warning=F}
#The final adjusted $R^2$ of our custom model is 0.9303:
R2.custom <- summary(fitFull.all4.Final)$adj.r.squared ##Custom Model
```

```{r, echo=F, results='hide', warning=F}
#The final adjusted $R^2$ of our Forward Selection model is 0.9327:
R2.forward <- summary(model2.Forward)$adj.r.squared
```


```{r, echo=F, results='hide', warning=F}
#The final adjusted $R^2$ of our Backward Elimination model is 0.9329:
R2.backward <- summary(model2.Backward)$adj.r.squared
```


```{r, echo=F, results='hide', warning=F}
#The final adjusted $R^2$ of our Stepwise Regression model is 0.9327:
R2.step <- summary(model2.Stepwise)$adj.r.squared
```

```{r, echo=F, results='hide', warning=F}
#Our cross-validated PRESS statistic output for our custom model is below.

# Set up repeated k-fold cross-validation
train.control2 <- trainControl(method = "cv", number = 10)
# Train the model
model.cv2 <- train(log(SalePrice) ~ MSSubClass + MSZoning + LotFrontage + 
                      LotArea + Street + Neighborhood + 
                      Condition1 + OverallQual + OverallCond + YearBuilt + YearRemodAdd + 
                      ExterCond + Foundation + 
                      BsmtExposure + BsmtFinType1 + BsmtFinSF1 + BsmtFinSF2 + BsmtUnfSF + 
                      HeatingQC + CentralAir + Electrical + `1stFlrSF` + 
                      `2ndFlrSF` + LowQualFinSF + BsmtFullBath + FullBath + HalfBath + 
                      KitchenAbvGr + KitchenQual + Functional + FireplaceQu + GarageYrBlt + 
                      GarageCars + GarageArea + GarageQual + WoodDeckSF + OpenPorchSF + 
                      EnclosedPorch + ScreenPorch + SaleType + SaleCondition,
                   data = train,
                   method = 'lm',
                   trControl = train.control2)

# print model summary
model.cv2

# get the CV results
res <- model.cv2$results

# get cross-validated PRESS statistic
PCV.custom <- PRESS.cv(model.cv2)
```

```{r, echo=F, results='hide', warning=F}
#Our cross-validated PRESS statistic output for our Forward Selection model is below.

# Set up repeated k-fold cross-validation
train.control2 <- trainControl(method = "cv", number = 10)
# Train the model
model.cv2 <- train(log(SalePrice) ~ OverallQual + GrLivArea + Neighborhood + 
                      BsmtFinSF1 + MSSubClass + OverallCond + YearBuilt + GarageCars + 
                      TotalBsmtSF + SaleCondition + LotArea + MSZoning + Functional + 
                      CentralAir + KitchenQual + Condition1 + FireplaceQu + BsmtExposure + 
                      BsmtFullBath + ScreenPorch + Exterior1st + YearRemodAdd + 
                      GarageQual + WoodDeckSF + OpenPorchSF + Street + LotConfig + 
                      LotFrontage + Foundation + Heating + KitchenAbvGr + EnclosedPorch + 
                      HalfBath + FullBath + MasVnrType + BsmtFinSF2 + HeatingQC + 
                      GarageArea + SaleType + ExterCond + PoolArea + BsmtFinType1 + 
                      GarageYrBlt + Electrical + `3SsnPorch` + LowQualFinSF,
                   data = train,
                   method = 'lm',
                   trControl = train.control2)

# print model summary
model.cv2

# get the CV results
res <- model.cv2$results

# get cross-validated PRESS statistic
PCV.forward <- PRESS.cv(model.cv2)
```

```{r, echo=F, results='hide', warning=F}
#Our cross-validated PRESS statistic output for our Backward Elimination model is below.

# Set up repeated k-fold cross-validation
train.control2 <- trainControl(method = "cv", number = 10)
# Train the model
model.cv2 <- train(log(SalePrice) ~ MSSubClass + MSZoning + LotFrontage +
                      LotArea + Street + LotConfig + LandSlope + Neighborhood +
                      Condition1 + OverallQual + OverallCond + YearBuilt + YearRemodAdd +
                      RoofMatl + Exterior1st + MasVnrType + ExterCond + Foundation +
                      BsmtExposure + BsmtFinType1 + BsmtFinSF1 + BsmtFinSF2 + BsmtUnfSF +
                      Heating + HeatingQC + CentralAir + Electrical + `1stFlrSF` +
                      `2ndFlrSF` + LowQualFinSF + BsmtFullBath + FullBath + HalfBath +
                      KitchenAbvGr + KitchenQual + Functional + FireplaceQu + GarageYrBlt +
                      GarageCars + GarageArea + GarageQual + WoodDeckSF + OpenPorchSF +
                      EnclosedPorch + ScreenPorch + PoolArea + SaleType + SaleCondition,
                   data = train,
                   method = 'lm',
                   trControl = train.control2)

# print model summary
model.cv2

# get the CV results
res <- model.cv2$results

# get cross-validated PRESS statistic
PCV.backward <- PRESS.cv(model.cv2)
```

```{r, echo=F, results='hide', warning=F}

#Our cross-validated PRESS statistic output for our Stepwise Regression model is below.
# Set up repeated k-fold cross-validation
train.control2 <- trainControl(method = "cv", number = 10)
# Train the model
model.cv2 <- train(log(SalePrice) ~ OverallQual + GrLivArea + Neighborhood + 
                      BsmtFinSF1 + MSSubClass + OverallCond + YearBuilt + GarageCars + 
                      TotalBsmtSF + SaleCondition + LotArea + MSZoning + Functional + 
                      CentralAir + KitchenQual + Condition1 + FireplaceQu + BsmtExposure + 
                      BsmtFullBath + ScreenPorch + Exterior1st + YearRemodAdd + 
                      GarageQual + WoodDeckSF + OpenPorchSF + Street + LotConfig + 
                      LotFrontage + Foundation + Heating + KitchenAbvGr + EnclosedPorch + 
                      HalfBath + FullBath + MasVnrType + BsmtFinSF2 + HeatingQC + 
                      GarageArea + SaleType + ExterCond + PoolArea + BsmtFinType1 + 
                      GarageYrBlt + Electrical + `3SsnPorch` + LowQualFinSF,
                   data = train,
                   method = 'lm',
                   trControl = train.control2)

# print model summary
model.cv2

# get the CV results
res <- model.cv2$results

# get cross-validated PRESS statistic
PCV.stepwise <- PRESS.cv(model.cv2)
```

```{r,echo=F, results='hide', warning=F}
#Our custom model achieved a Kaggle score of 0.13290. This was our best score.

## To test in Kaggle, submit the produced "submit" file
test$predicted.log.price <- predict.lm(fitFull.all4.Final, test)
test$predicted.log.price[is.na(test$predicted.log.price)] <- mean(test$predicted.log.price, na.rm = T)

submit <- test %>% mutate(SalePrice = exp(predicted.log.price)) %>% subset(select=c(Id, SalePrice))

write.csv(submit, file = "./kaggle_submission.csv", row.names = F)
```

```{r,echo=F, results='hide', warning=F}
#Our Forward Selection model achieved a Kaggle score of 0.13476

## To test in Kaggle, submit the produced "submit" file
test$predicted.log.price <- predict.lm(final.Forward.Model, test)
test$predicted.log.price[is.na(test$predicted.log.price)] <- mean(test$predicted.log.price, na.rm = T)

submit <- test %>% mutate(SalePrice = exp(predicted.log.price)) %>% subset(select=c(Id, SalePrice))

# write.csv(submit, file = "./kaggle_submission.csv", row.names = F)
```

```{r,echo=F, results='hide', warning=F}
#Our Backward Selection model achieved a Kaggle score of 0.13475

## To test in Kaggle, submit the produced "submit" file
test$predicted.log.price <- predict.lm(final.Backward.Model, test)
test$predicted.log.price[is.na(test$predicted.log.price)] <- mean(test$predicted.log.price, na.rm = T)

submit <- test %>% mutate(SalePrice = exp(predicted.log.price)) %>% subset(select=c(Id, SalePrice))

# write.csv(submit, file = "./kaggle_submission.csv", row.names = F)
```

```{r,echo=F, results='hide', warning=F}
#Our Stepwise Regression model also achieved a Kaggle score of 0.13476

## To test in Kaggle, submit the produced "submit" file
test$predicted.log.price <- predict.lm(final.Stepwise.Model, test)
test$predicted.log.price[is.na(test$predicted.log.price)] <- mean(test$predicted.log.price, na.rm = T)

submit <- test %>% mutate(SalePrice = exp(predicted.log.price)) %>% subset(select=c(Id, SalePrice))

# write.csv(submit, file = "./kaggle_submission.csv", row.names = F)
```


```{r, echo=FALSE}
# print accuracy metrics to md table

kable(data.frame('Model' = c('Custom', 'Forward Selection', 'Backward Selection', 'Stepwise Regression'), 
                 'Kaggle Score'=c(0.13290,0.13476,0.13475,0.13476),
                 'CV Press'=c(PCV.custom, PCV.forward, PCV.backward, PCV.stepwise),
                 'Adjused R Squared'=c(R2.custom, R2.forward, R2.backward, R2.step)),
      "latex", booktabs = T)  %>%
  kable_styling(position = "center")
```

## Conclusion

In an effort to produce a highly accurate and repeatable predictive model using linear regression, all explanatory variables were considered with three types of variable selection techniques: forward selection, backward elimination, and stepwise regression. Additionally, a custom model was initially produced by eliminating variables suggested by the automatic selection processes, visually exploring the data with pairwise scatter plots, and adding interaction terms based on graphical analysis and domain knowledge. Automatic selection was reapplied to suggest terms from the initial custom model, which was then again adjusted for final inspection by the automatic techniques. The final models suggested strictly by the automatic techniques produced high $R^2$ values, but performed poorly on the Kaggle test set. This suggests the automatic techniques alone were overfitting to the training data. The final custom model, however, produced a high $R^2$ value and performed remarkably well on the Kaggle test set (see section \ref{appendix:kaggle}). This suggests that the custom model is not overfitting to the training data and generalizes well to an unseen dataset. Ultimately, we determined the best approach is a combination of automatic selection, visual and analytic inspection, and the application of domain knowledge.

\newpage

# Appendix

## Checking for Linearity in `SalePrice` vs `GrLivArea`

\label{appendix:linearity}

The scatter plot in Figure \@ref(fig:scatter-plot) shows relationship of `SalePrice` vs `GrLivArea` for all three neighborhoods of interest to Century 21. Based on this plot, it does not appear that this relationship meets the assumptions of linear regression, specifically the constant varaince assumption. The response will be transformed to attempt to handle the changing variance.

```{r, scatter-plot, echo=FALSE, fig.width=5, fig.height=5, fig.align='center', fig.cap='Scatter Plot of Sale Price vs Living Room Area', out.width = '45%', fig.pos="htbp"}

# scatter plot of observations from all three neighborhoods
train %>% 
  ggplot(aes(x = (GrLivArea), y = (SalePrice))) +
  geom_point(alpha = 0.3) +
  labs(title = 'Sale Price vs Living Room Area', 
       y = 'Sale Price', x = 'Living Room Area (sq. ft.)')
```

The images below show the scatter plots of log sale price vs living room area (Figure \@ref(fig:scatter-plots)). In the image on the right, the scatter plot is shown for each neighborhood. In the image on the left the observations for all three neighborhoods are included. In all cases, a linear model appears to be reasonable to model this data.

```{r, scatter-plots, echo=FALSE, fig.width=5, fig.height=5, fig.align='center', fig.cap='Scatter Plots of Log of Sale Price vs Living Room Area', out.width = '45%', fig.pos="htbp", fig.show='hold'}

# plots of log of sale price ~ living room area

# create scatter plot for northwest ames
regplot.names <- train %>%
  filter(Neighborhood == 'NAmes') %>%
  ggplot(aes(x = (GrLivArea), y = log(SalePrice))) +
  geom_point(alpha = 0.3) +
  ylim(10, 13) +
  xlim(0, 3500) +
  labs(subtitle = 'Northwest Ames', 
       y = 'Log of Sale Price', x = 'Living Room Area (sq. ft.)')

# create scatter plot for edwards
regplot.ed <- train %>%
  filter(Neighborhood == 'Edwards') %>%
  ggplot(aes(x = (GrLivArea), y = log(SalePrice))) +
  geom_point(alpha = 0.3) +
  ylim(10, 13) +
  xlim(0, 3500) +
  labs(subtitle = 'Edwards', 
       y = 'Log of Sale Price', x = 'Living Room Area (sq. ft.)')

# create regression plot for brookside
regplot.brk <- train %>%
  filter(Neighborhood == 'BrkSide') %>%
  ggplot(aes(x = (GrLivArea), y = log(SalePrice))) +
  geom_point(alpha = 0.3) +
  ylim(10, 13) +
  xlim(0, 3500) +
  labs(subtitle = 'Brook Side', 
       y = 'Log of Sale Price', x = 'Living Room Area (sq. ft.)')

# add the scatter plots for the neighborhood into a single plot
grid.arrange(regplot.names,regplot.ed,regplot.brk, nrow = 2,
             top = 'Regression Plots for Neighborhoods')

# scatter plot of observations from all three neighborhoods
train %>% 
  ggplot(aes(x = (GrLivArea), y = log(SalePrice))) +
  geom_point(alpha = 0.3) +
  labs(title = 'Log of Sale Price vs Living Room Area', 
       y = 'Log of Sale Price', x = 'Living Room Area (sq. ft.)')
```


## Analysis of Influential points

\label{appendix:infleu-points}

The two outlying observations with living room areas greater than 4000 sq. ft. appear to be from a different distribution than the main dataset. Since these are partial sales, it is possible that the sale prices do not reflect market value. For this reason, we will limit the analysis to properities with less than 3500 sq. ft. (Figure \@ref(fig:infleu-points))

```{r, infleu-points, echo=FALSE, fig.width=5, fig.height=5, fig.align='center', fig.cap='Influential Points', out.width = '50%', fig.pos="htbp"}

# scatter plot of observations from all three neighborhoods with labeling by `SaleCondition`
train.mod %>% 
  filter(Neighborhood %in% c("Edwards", "BrkSide", "NAmes")) %>%
  ggplot(aes(x = (GrLivArea), y = log(SalePrice))) +
  geom_point(alpha = 0.3) +
  labs(title = 'Log of Sale Price vs Living Room Area', 
       y = 'Log of Sale Price', x = 'Living Room Area (sq. ft.)') +
  geom_text(aes(label = ifelse((log(GrLivArea) > 7.75 & log(SalePrice) > 11) |
                                 (log(SalePrice) > 12.45),
                               SaleCondition, '')), hjust=0, vjust=0)
```


\newpage

## Models Suggested by Automated Selection

### Forward Selection

\label{appendix:forSelection}

The model suggested by forward selection.

```r
log(SalePrice) ~ 
              OverallQual + GrLivArea + Neighborhood + BsmtFinSF1 +
              OverallCond + YearBuilt + TotalBsmtSF + GarageCars + MSZoning +
              SaleCondition + BldgType + Functional + LotArea + KitchenQual +
              BsmtExposure + CentralAir + Condition1 + ScreenPorch + BsmtFullBath +
              Heating + Fireplaces + YearRemodAdd + Exterior1st + GarageQual +
              WoodDeckSF + SaleType + OpenPorchSF + HeatingQC + LotConfig +
              EnclosedPorch + ExterCond + PoolQC + Foundation + LandSlope +
              RoofMatl + GarageArea + MasVnrType + HalfBath + PoolArea +
              `3SsnPorch` + Street + KitchenAbvGr + GarageCond + FullBath +
              BsmtQual + BsmtFinSF2
```

### Backward Selection

\label{appendix:backSelection}

The model suggested by backward selection.

```r
log(SalePrice) ~ 
           MSZoning + LotArea + Street + LotConfig + LandSlope +
           Neighborhood + Condition1 + Condition2 + BldgType + OverallQual +
           OverallCond + YearBuilt + YearRemodAdd + RoofStyle + RoofMatl +
           Exterior1st + MasVnrType + ExterCond + Foundation + BsmtQual +
           BsmtCond + BsmtExposure + BsmtFinSF1 + BsmtFinSF2 + BsmtUnfSF +
           Heating + HeatingQC + CentralAir + `1stFlrSF` + `2ndFlrSF` +
           LowQualFinSF + BsmtFullBath + FullBath + HalfBath + KitchenAbvGr +
           KitchenQual + Functional + Fireplaces + GarageCars + GarageArea +
           GarageQual + GarageCond + WoodDeckSF + OpenPorchSF + EnclosedPorch +
           `3SsnPorch` + ScreenPorch + PoolArea + PoolQC + SaleType +
           SaleCondition
```

\newpage

### Stepwise Selection

\label{appendix:stepWSelection}

The model suggested by stepwise selection.

```r
log(SalePrice) ~ 
           MSZoning + LotArea + Street + LotConfig + LandSlope +
           Neighborhood + Condition1 + Condition2 + BldgType + OverallQual +
           OverallCond + YearBuilt + YearRemodAdd + RoofStyle + RoofMatl +
           Exterior1st + MasVnrType + ExterCond + Foundation + BsmtQual +
           BsmtCond + BsmtExposure + BsmtFinSF1 + BsmtFinSF2 + BsmtUnfSF +
           Heating + HeatingQC + CentralAir + `1stFlrSF` + `2ndFlrSF` +
           LowQualFinSF + BsmtFullBath + FullBath + HalfBath + KitchenAbvGr +
           KitchenQual + Functional + Fireplaces + GarageCars + GarageArea +
           GarageQual + GarageCond + WoodDeckSF + OpenPorchSF + EnclosedPorch +
           `3SsnPorch` + ScreenPorch + PoolArea + PoolQC + SaleType +
           SaleCondition
```

### Custom Model

\label{appendix:customSelection}

The model constructed by hand.

```r
log(SalePrice) ~ 
            BsmtUnfSF + CentralAir + HalfBath + KitchenQual + Neighborhood +
            OverallCond + OverallQual + RoofMatl + `1stFlrSF` + `2ndFlrSF` + 
            YearBuilt + MSZoning:Neighborhood + YearBuilt:Neighborhood
```

## Kaggle Score

\label{appendix:kaggle}

The following image shows the result on Kaggle for the custom model.

```{r, kaggle, echo=F, out.width='90%', fig.align='center'}
knitr::include_graphics('./Evidence for Model Score on Kaggle.PNG')
```

\newpage

## Assumption Assessment Plots for Automatic Selection Models

\label{appendix:assessmentPlots}

The following discussion applies to the assumption assessment for the three models produced by automatic selection.

Generally, based on the diagonstic plots, these models appear to reasonably meet the assumptions for linear regression. The standardized residuals do not appear to exhibit a discernible pattern, indicating constant variance along the regression, or homoscedasticity. However, there are a small number of observations - relative to the overall sample size - with unusually high residuals. Nonetheless, this is not enough to add detrimental impact to the model. Based on the QQ plot, there is a small level of deviation on the ends of the distribution of the errors, but for the most part, the errors adhere to normality. The sample size should be sufficient to protect against this non-normality. Based on the standardized residuals vs. leverage plot, only a few values have high leverage and are outlying. Compared to the custom model (Figure \@ref(fig:custom-assumptions)), these diagnostic plots for these models show a few more influential observations with high leverage. However, these observations cannot be excluded from the model.

```{r, forward-assumptions, echo=F, warning=F, fig.width=5, fig.height=5, fig.align='center', fig.cap='Forward Selection Assumption Assessment Plots', out.width = '45%', fig.pos="htbp", fig.show='hold'}
# create plots of residuals
basic.fit.plots(train, model2.Forward)
# create leverage / outlier plot
ols_plot_resid_lev(model2.Forward)
```


```{r, backward-selection, echo=F, warning=F, fig.width=5, fig.height=5, fig.align='center', fig.cap='Backward Selection Assumption Assessment Plots', out.width = '45%', fig.pos="htbp", fig.show='hold'}
# create plots of residuals
basic.fit.plots(train, model2.Backward)
# create leverage / outlier plot
ols_plot_resid_lev(model2.Backward)
```


```{r, stepwise-assimptions, echo=F, warning=F, fig.width=5, fig.height=5, fig.align='center', fig.cap='Stepwise Selection Assumption Assessment Plots', out.width = '45%', fig.pos="htbp", fig.show='hold'}
# create plots of residuals
basic.fit.plots(train, model2.Stepwise)
# create leverage / outlier plot
ols_plot_resid_lev(model2.Stepwise)
```

\newpage

# References
