# Assignment 4 
# Chapter 6: Multiple Linear Regression 


> (c) 2019-2020 Galit Shmueli, Peter C. Bruce, Peter Gedeck 
>
> _Data Mining for Business Analytics: Concepts, Techniques, and Applications in Python_ (First Edition) 
> Galit Shmueli, Peter C. Bruce, Peter Gedeck, and Nitin R. Patel. 2019.
>
> Date: 2020-03-08
>
> Python Version: 3.8.2
> Jupyter Notebook Version: 5.6.1
>
> Packages:
>   - dmba: 0.0.12
>   - matplotlib: 3.2.0
>   - numpy: 1.18.1
>   - pandas: 1.0.1
>   - seaborn: 0.10.0
>   - scikit-learn: 0.22.2
>   - statsmodels: 0.11.1
>
> The assistance from Mr. Kuber Deokar and Ms. Anuja Kulkarni in preparing these solutions is gratefully acknowledged.


# import required functionality for this chapter
from pathlib import Path

import pandas as pd
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
import statsmodels.formula.api as sm

import matplotlib.pylab as plt
import seaborn as sns

from dmba import regressionSummary, exhaustive_search
from dmba import backward_elimination, forward_selection, stepwise_selection
from dmba import adjusted_r2_score, AIC_score #, BIC_score

%matplotlib inline

!pip install dmba

# Problem 6.1 Predicting Boston Housing Prices 

The file _BostonHousing.csv_ contains information collected by the US Bureau of the Census concerning housing in the area of
Boston, Massachusetts. The dataset includes information on 506 census housing tracts in the Boston area. The goal is to predict the median house price in new tracts based on information such as crime rate, pollution, and number of rooms. The dataset contains 13 predictors, and the outcome variable is the median house price (MEDV). Table 6.11 describes each of the predictors and the outcome variable.

![TABLE_6.11](img/TABLE_6.11.PNG)


pd.options.display.max_rows=506
housing_df = pd.read_csv('BostonHousing-1.csv')
housing_df


housing_df.isnull().sum()

# load the data into a dataframe called housing data
#MISSING 1 line of code
housing_df = pd.read_csv('BostonHousing-1.csv')

# display column/variable names
#Create a list called columns with all of the housing_df columns names in it
#MISSING 1 line of code
columns=housing_df.columns.values


print("Variables in the data are: ")
print(columns)

# review first 5 records in the data
print("\nFirst 5 records in the data are:")
#MISSING 1 line of code
housing_df.head(5)


__6.1.a.__ (Double click and add answers here)

Why should the data be partitioned into training and validation sets? 
<br><b>Answer:</b><br>?????

What will the training set be used for? 
<br><b>Answer:</b><br>?????

What will the validation set be used for?
<br><b>Answer:</b><br> ?????



__6.1.b.__ Fit a multiple linear regression model to the median house price (MEDV) as a function of CRIM, CHAS, and RM. Write the equation for predicting the median house price from the predictors in the model.

__Answer:__


# select columns for regression analysis
outcome = housing_df[['MEDV']]
predictors = housing_df[['CRIM', 'CHAS', 'RM']]

#Create a dataframe called x containing the predictor columns
#MISSING 1 line of code
x=pd.DataFrame(predictors)



#Create a dataframe (technically a series) containing the outcome variable. Call it y
#MISSING 1 line of code
y=pd.DataFrame(outcome)




x.ndim

y.ndim

#Create a model called housing_lm and set it to be a LinearRegression() model
#MISSING 1 line of code
model=LinearRegression(fit_intercept=True,
    normalize=True,
    copy_X=True,
    n_jobs=5,
    positive=False)

# fit the regression model y on x
#MISSING 1 line of code
model.fit(x,y)

# print the intercept
#MISSING 1 line of code
print(model.predict(x))

#print the list of predictor columns and the coefficients
#MISSING 1 line of code
print(model.coef_)



__6.1.c.__ Using the estimated regression model, what median house price is predicted for a tract in the Boston area that does not bound the Charles River, has a crime rate of 0.1, and where the average number of rooms per house is 6?

__Answer:__ 

new_df = pd.DataFrame(
                        [[0.1, 0, 6]],
                        columns=['CRIM', 'CHAS', 'RM']
                        )
new_df

#Run the prediction model that you created using the above created dataframe containing
# the new predictor values the the results housing_lm_pred
#MISSING 1 line of code


print('Predicted value for median house price based on the model built using  dataset is:', housing_lm_pred)

The predicted value of median house price based on the model built using training dataset is $20832.32.

__6.1.d.i.__ Reduce the number of predictors:
Which predictors are likely to be measuring the same thing among the 13 predictors? Discuss the relationships among INDUS, NOX, and TAX.

__Answer:__ 
?????


__6.1.d.ii.__ Compute the correlation table for the 11 <b><font color="red">numerical</font></b> predictors and search for highly correlated pairs. These have potential redundancy and can cause multicollinearity. Choose which ones to remove based on this table.

__Answer:__
?????

# variables in the data
housing_df.columns

# Create a new dataframe called predictors_df with only numerical predictors
#MISSING 1-5 lines of code (many different ways we have done this before)
a=housing_df[["CRIM","ZN","INDUS","NOX","RM","AGE","DIS","RAD","TAX","PTRATIO","LSTAT"]]
predictors_df=pd.DataFrame(a)


predictors_df.columns

#Create a correlation table called corr using the predictors_df dataframe
#MISSING 1 line of code
corr=predictors_df.corr()

# Use seaborn to create a correlation/heatmap (color and such does not matter)
#MISSING 1 line of code
sns.heatmap(corr.corr())



# correlation table
corr

If the correlation between variables if greater than 0.7 we can say that the two variables are highly correlated. From the above table, the pairs of highly correlated variables are:

1) ?????

2) ?????

3) ?????

4) ?????

5) ?????

6) ?????

According to the correlation table, we might be able to remove some variables that do not add much information to others that we keep. We might remove ?????

__6.1.d.iii.__ Use three subset selection algorithms: backward, forward, and stepwise to reduce the remaining predictors as follows: Run stepwise on the training set. Choose the top model from each stepwise run. Then use each of these models separately to predict the validation set. Compare RMSE, MAPE, and mean error, as well as lift charts. Finally, describe the best model.

__Answer:__

# partition the data into training (60%) and validation (40%) sets
predictors = housing_df[['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'LSTAT']]
outcome = housing_df[['MEDV']]

# partition the data
#Create a dataframe called X with the columns in the predictors[] list above
# Make sure to turn text columns (categorical) values into dummy variable columns
#MISSING 1 line of code
x=pd.DataFrame(predictors)

#Create a dataframe (technically a Series) called y containing the outcome column
#MISSING 1 line of code
pd.DataFrame(outcome)

#Split the data into 40/60 validation and training datasets with a random state of 1
#MISSING 1 line of code
train_X,valid_X,train_y,valid_y=train_test_split(x,y,test_size=0.4,random_state=1)

print('Training set:', train_X.shape, 'Validation set:', valid_X.shape)


predictors = housing_df[['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'LSTAT']]
outcome = housing_df[['MEDV']]


x=pd.DataFrame(predictors)


pd.DataFrame(outcome)


train_X,valid_X,train_y,valid_y=train_test_split(x,y,test_size=0.4,random_state=1)

print('Training set:', train_X.shape, 'Validation set:', valid_X.shape)

NOTE: There is currently no support in scikit-learn or statsmodels for stepwise regression. It is however straightforward to 
implement such an approach in a few lines of code using the `dmba` implementation.

# backward elimination

def train_model(variables):
    model = LinearRegression()
    model.fit(train_X[variables], train_y)
    return model

def score_model(model, variables):
    return AIC_score(train_X, model.predict(train_X[variables]), model)

#Run the backward_elimination function created by the author and shown in Ch06
#MISSING 1 line of code
best_variables, best_model = backward_elimination(x.columns, train_model, score_model, verbose=True)


print("Best Subset:", best_variables)

# backward elimination

def train_model(variables):
    model = LinearRegression()
    model.fit(train_X[variables], train_y)
    return model

def score_model(model, variables):
    return AIC_score(train_X, model.predict(train_X[variables]), model)

#Run the backward_elimination function created by the author and shown in Ch06
#MISSING 1 line of code
#best_variables, best_model = backward_elimination(x.columns, train_model, score_model, verbose=True)
best_score, best_variables = backward_elimination(train_X.columns, train_y)



print("Best Subset:", best_variables)


def train_model(variables):
    model = LinearRegression()
    model.fit(train_X[variables], train_y)
    return model

def score_model(model, variables):
    return AIC_score(train_y, model.predict(train_X[variables]), model)

# Run the backward_elimination function to find the best variables to include in the model
best_variables, best_model = backward_elimination(x.columns, train_model, score_model, verbose=True)
print("Best Subset", best_variables)

Variables: CRIM, ZN, INDUS, CHAS, NOX, RM, AGE, DIS, RAD, TAX, PTRATIO, LSTAT
Start: score=1807.23
Step: score=1805.30, remove AGE
Step: score=1803.57, remove INDUS
Step: score=1803.57, remove None
Best Subset: ['CRIM', 'ZN', 'CHAS', 'NOX', 'RM', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'LSTAT']

# forward selection
# The initial model is the constant model - this requires special handling in train_model and score_model

# define train_model function
def train_model(variables):
    if len(variables) == 0:
        return None
model = LinearRegression()
model.fit(train_X[variables], train_y)
return model

# define score_model function
def score_model(model, variables):
    if len(variables) == 0:
        return

# forward selection
# The initial model is the constant model - this requires special handling in train_model and score_model

#Write the train_model function (starting with "def") shown in the examples for Ch 06
#MISSING 6 lines of code
def .......

#Write the score_model function (starting with "def") shown in the examples for Ch 06
#MISSING 4 lines of code
def .....

#Run the forward_selection function created by the author and shown in Ch06
#MISSING 1 line of code

print("Best Subset:", best_variables)

Variables: CRIM, ZN, INDUS, CHAS, NOX, RM, AGE, DIS, RAD, TAX, PTRATIO, LSTAT
Start: score=4542.98
Step: score=4542.86, remove LSTAT
Step: score=4542.63, remove RM
Step: score=4542.50, remove TAX
Step: score=4542.28, remove PTRATIO
Step: score=4542.02, remove INDUS
Step: score=4541.83, remove NOX
Step: score=4541.57, remove AGE
Step: score=4541.28, remove RAD
Step: score=4540.32, remove CRIM
Step: score=4540.07, remove DIS
Step: score=4539.04, remove ZN
Best Subset: LinearRegression()

# stepwise (both) method

#Run the stepwise_selection function created by the author and shown in Ch06
#MISSING 1 line of code

print("Best Subset:", best_variables)

Above we see that all three subset selection methods gave the same best subset/model:
Each model is designed using the following 10 variables: LSTAT, RM, PTRATIO, CHAS, NOX, DIS, ZN, CRIM, RAD, TAX. 

# Re-run the Regression but this time fit the model with best subset variables from the 
#  subset reductions from above

#Define the outcome and predictor variables
outcome = 'MEDV'
predictors = ['LSTAT', 'RM', 'PTRATIO', 'CHAS', 'NOX', 'DIS', 'ZN', 'CRIM', 'RAD', 'TAX']

#Create a dataframe called X containing the new predictor columns
#MISSING 1 line of code

#Create a dataframe (Series) called y containing the outcome column.
#MISSING 1 line of code

# fit the regression model y on X
#MISSING 2 lines of code


# print the intercept
#MISSING 1 line of code

#print the predictor column names and the coefficients
#MISSING 1 line of code


# print performance measures (training set)
print("\nModel performance on training data:")
#MISSING 1 line of code


# predict prices in validation set, print first few predicted/actual values and residuals
#MISSING 1 line of code


result = pd.DataFrame({'Predicted': house_lm_pred, 'Actual': valid_y, 'Residual': valid_y - house_lm_pred})

# print performance measures (validation set)
print("\nModel performance on validation data:")
#MISSING 1 line of code


