Linear Regression: Boston House Price Prediction

---

Problem Statement

---

The problem on hand is to predict the housing prices of a town or a suburb based on the features of the locality provided to us. In the process, we need to identify the most important features in the dataset. We need to employ techniques of data preprocessing and build a linear regression model that predicts the prices for us.

---

Data Information

---

Each record in the database describes a Boston suburb or town. The data was drawn from the Boston Standard Metropolitan Statistical Area (SMSA) in 1970. Detailed attribute information can be found below-

Attribute Information (in order):

• CRIM: per capita crime rate by town
• ZN: proportion of residential land zoned for lots over 25,000 sq.ft.
• INDUS: proportion of non-retail business acres per town
• CHAS: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
• NOX: nitric oxides concentration (parts per 10 million)
• RM: average number of rooms per dwelling
• AGE: proportion of owner-occupied units built before 1940
• DIS: weighted distances to five Boston employment centers
• RAD: index of accessibility to radial highways
• TAX: full-value property-tax rate per 10,000 dollars
• PTRATIO: pupil-teacher ratio by town
• LSTAT: %lower status of the population
• MEDV: Median value of owner-occupied homes in 1000 dollars

Let us start by importing the required libraries

# import libraries for data manipulation
import pandas as pd
import numpy as np

# import libraries for data visualization
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.graphics.gofplots import ProbPlot

# import libraries for building linear regression model
from statsmodels.formula.api import ols
import statsmodels.api as sm
from sklearn.linear_model import LinearRegression

# import library for preparing data
from sklearn.model_selection import train_test_split

# import library for data preprocessing
from sklearn.preprocessing import MinMaxScaler

import warnings
warnings.filterwarnings("ignore")

Read the dataset

df = pd.read_csv("Boston.csv")
df.head()

	CRIM	ZN	INDUS	CHAS	NOX	RM	AGE	DIS	RAD	TAX	PTRATIO	LSTAT	MEDV
0	0.00632	18.0	2.31	0	0.538	6.575	65.2	4.0900	1	296	15.3	4.98	24.0
1	0.02731	0.0	7.07	0	0.469	6.421	78.9	4.9671	2	242	17.8	9.14	21.6
2	0.02729	0.0	7.07	0	0.469	7.185	61.1	4.9671	2	242	17.8	4.03	34.7
3	0.03237	0.0	2.18	0	0.458	6.998	45.8	6.0622	3	222	18.7	2.94	33.4
4	0.06905	0.0	2.18	0	0.458	7.147	54.2	6.0622	3	222	18.7	5.33	36.2

Observations

• The price of the house indicated by the variable MEDV is the target variable and the rest are the independent variables based on which we will predict house price.

Get information about the dataset using the info() method

df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 506 entries, 0 to 505
Data columns (total 13 columns):
 #   Column   Non-Null Count  Dtype  
---  ------   --------------  -----  
 0   CRIM     506 non-null    float64
 1   ZN       506 non-null    float64
 2   INDUS    506 non-null    float64
 3   CHAS     506 non-null    int64  
 4   NOX      506 non-null    float64
 5   RM       506 non-null    float64
 6   AGE      506 non-null    float64
 7   DIS      506 non-null    float64
 8   RAD      506 non-null    int64  
 9   TAX      506 non-null    int64  
 10  PTRATIO  506 non-null    float64
 11  LSTAT    506 non-null    float64
 12  MEDV     506 non-null    float64
dtypes: float64(10), int64(3)
memory usage: 51.5 KB

Observations

• There are a total of 506 non-null observations in each of the columns. This indicates that there are no missing values in the data.

• Every column in this dataset is numeric in nature.

Let's now check the summary statistics of this dataset

Step 1: Find the summary statistics and write observations

df.describe()

	CRIM	ZN	INDUS	CHAS	NOX	RM	AGE	DIS	RAD	TAX	PTRATIO	LSTAT	MEDV
count	506.000000	506.000000	506.000000	506.000000	506.000000	506.000000	506.000000	506.000000	506.000000	506.000000	506.000000	506.000000	506.000000
mean	3.613524	11.363636	11.136779	0.069170	0.554695	6.284634	68.574901	3.795043	9.549407	408.237154	18.455534	12.653063	22.532806
std	8.601545	23.322453	6.860353	0.253994	0.115878	0.702617	28.148861	2.105710	8.707259	168.537116	2.164946	7.141062	9.197104
min	0.006320	0.000000	0.460000	0.000000	0.385000	3.561000	2.900000	1.129600	1.000000	187.000000	12.600000	1.730000	5.000000
25%	0.082045	0.000000	5.190000	0.000000	0.449000	5.885500	45.025000	2.100175	4.000000	279.000000	17.400000	6.950000	17.025000
50%	0.256510	0.000000	9.690000	0.000000	0.538000	6.208500	77.500000	3.207450	5.000000	330.000000	19.050000	11.360000	21.200000
75%	3.677083	12.500000	18.100000	0.000000	0.624000	6.623500	94.075000	5.188425	24.000000	666.000000	20.200000	16.955000	25.000000
max	88.976200	100.000000	27.740000	1.000000	0.871000	8.780000	100.000000	12.126500	24.000000	711.000000	22.000000	37.970000	50.000000

Observations:

• CRIM is left skewed - most of the towns have low crime rate.
• ZN is left skewed
• Most of the towns do not bound river

Before performing the modeling, it is important to check the univariate distribution of the variables.

Univariate Analysis

Check the distribution of the variables

# let's plot all the columns to look at their distributions
for i in df.columns:
    plt.figure(figsize=(7, 4))
    sns.histplot(data=df, x=i, kde = True)
    plt.show()

Observations

• The variables CRIM and ZN are positively skewed. This suggests that most of the areas have lower crime rates and most residential plots are under the area of 25,000 sq. ft.
• The variable CHAS, with only 2 possible values 0 and 1, follows a binomial distribution, and the majority of the houses are away from Charles river (CHAS = 0).
• The distribution of the variable AGE suggests that many of the owner-occupied houses were built before 1940.
• The variable DIS (average distances to five Boston employment centers) has a nearly exponential distribution, which indicates that most of the houses are closer to these employment centers.
• The variables TAX and RAD have a bimodal distribution., indicating that the tax rate is possibly higher for some properties which have a high index of accessibility to radial highways.
• The dependent variable MEDV seems to be slightly right skewed.

As the dependent variable is sightly skewed, we will apply a log transformation on the 'MEDV' column and check the distribution of the transformed column.

df['MEDV_log'] = np.log(df['MEDV'])

sns.histplot(data=df, x='MEDV_log', kde = True)

<AxesSubplot:xlabel='MEDV_log', ylabel='Count'>

Observations

• The log-transformed variable (MEDV_log) appears to have a nearly normal distribution without skew, and hence we can proceed.

Before creating the linear regression model, it is important to check the bivariate relationship between the variables. Let's check the same using the heatmap and scatterplot.

Bivariate Analysis

Let's check the correlation using the heatmap

Step 2:

• Plot the correlation heatmap between the variables

plt.figure(figsize=(12,8))
cmap = sns.diverging_palette(230, 20, as_cmap=True)
sns.heatmap(df.corr(),annot=True,fmt='.2f',cmap=cmap ) #write your code here
plt.show()

Observations:

• Tax(TAX) is highly correlated with distance from radial highways(RAD).
• NOX is highly corealted with older houses(AGE)
• Old houses(AGE) are highly correalted with INDUS. Old houses are near industrial places.
• RM is highly corelated with Median Cost. More the rooms in house, more costly it is.

Now, we will visualize the relationship between the pairs of features having significant correlations.

Visualizing the relationship between the features having significant correlations (> 0.7)

Step 3 :

• Create a scatter plot to visualize the relationship between the features having significant correlations (>0.7)

# scatterplot to visualize the relationship between NOX and INDUS
plt.figure(figsize=(6, 6))
sns.scatterplot(x='INDUS',y='NOX',data=df) #write you code here
plt.show()

Observations:

• NOX is postively corelated with INDUS. There seem to be few outlierrs in the data though around INDUS=18.
• More industrial the area, more the NO2 content.
• We need to explore, values of NOX at around INDUS=18.

df[df.INDUS>18.00]['INDUS'].value_counts()

18.10    132
19.58     30
21.89     15
25.65      7
27.74      5
Name: INDUS, dtype: int64

# Explore data at INDUS==18.10
df[df.INDUS==18.10]['NOX'].value_counts()

0.713    18
0.693    14
0.740    13
0.700    11
0.584     8
0.770     8
0.679     8
0.671     7
0.614     7
0.597     6
0.718     6
0.631     5
0.532     5
0.583     4
0.580     4
0.655     3
0.668     3
0.659     2
Name: NOX, dtype: int64

df[df.INDUS==18.10][df.NOX==.713]

	CRIM	ZN	INDUS	CHAS	NOX	RM	AGE	DIS	RAD	TAX	PTRATIO	LSTAT	MEDV	MEDV_log
433	5.58107	0.0	18.1	0	0.713	6.436	87.9	2.3158	24	666	20.2	16.22	14.3	2.660260
434	13.91340	0.0	18.1	0	0.713	6.208	95.0	2.2222	24	666	20.2	15.17	11.7	2.459589
448	9.32909	0.0	18.1	0	0.713	6.185	98.7	2.2616	24	666	20.2	18.13	14.1	2.646175
449	7.52601	0.0	18.1	0	0.713	6.417	98.3	2.1850	24	666	20.2	19.31	13.0	2.564949
450	6.71772	0.0	18.1	0	0.713	6.749	92.6	2.3236	24	666	20.2	17.44	13.4	2.595255
451	5.44114	0.0	18.1	0	0.713	6.655	98.2	2.3552	24	666	20.2	17.73	15.2	2.721295
452	5.09017	0.0	18.1	0	0.713	6.297	91.8	2.3682	24	666	20.2	17.27	16.1	2.778819
453	8.24809	0.0	18.1	0	0.713	7.393	99.3	2.4527	24	666	20.2	16.74	17.8	2.879198
454	9.51363	0.0	18.1	0	0.713	6.728	94.1	2.4961	24	666	20.2	18.71	14.9	2.701361
455	4.75237	0.0	18.1	0	0.713	6.525	86.5	2.4358	24	666	20.2	18.13	14.1	2.646175
456	4.66883	0.0	18.1	0	0.713	5.976	87.9	2.5806	24	666	20.2	19.01	12.7	2.541602
457	8.20058	0.0	18.1	0	0.713	5.936	80.3	2.7792	24	666	20.2	16.94	13.5	2.602690
458	7.75223	0.0	18.1	0	0.713	6.301	83.7	2.7831	24	666	20.2	16.23	14.9	2.701361
459	6.80117	0.0	18.1	0	0.713	6.081	84.4	2.7175	24	666	20.2	14.70	20.0	2.995732
460	4.81213	0.0	18.1	0	0.713	6.701	90.0	2.5975	24	666	20.2	16.42	16.4	2.797281
461	3.69311	0.0	18.1	0	0.713	6.376	88.4	2.5671	24	666	20.2	14.65	17.7	2.873565
462	6.65492	0.0	18.1	0	0.713	6.317	83.0	2.7344	24	666	20.2	13.99	19.5	2.970414
463	5.82115	0.0	18.1	0	0.713	6.513	89.9	2.8016	24	666	20.2	10.29	20.2	3.005683

# Lets check distribution of NOX in boxplot
sns.boxplot(df.NOX)

<AxesSubplot:xlabel='NOX'>

• These look like lower income suburbs - mostly very old dwellings
• But there are no outliers in data.
• I think we can drop NOX

# scatterplot to visualize the relationship between AGE and NOX
plt.figure(figsize=(6, 6))
sns.scatterplot(x='AGE',y='NOX',data=df) #Write your code here
plt.show()

Observations:

• AGE is positively corelated with NOX
• We can drop NOX

# scatterplot to visualize the relationship between DIS and NOX
plt.figure(figsize=(6, 6))
sns.scatterplot(x='DIS',y='NOX',data=df) #Write your code here
plt.show()

Observations:

• NOX is negatively correlated with DIS
• Need to check if NOX can be dropped

# scatterplot to visualize the relationship between AGE and DIS
plt.figure(figsize=(6, 6))
sns.scatterplot(x = 'AGE', y = 'DIS', data = df)
plt.show()

Observations:

• The distance of the houses to the Boston employment centers appears to decrease moderately as the the proportion of the old houses increase in the town. It is possible that the Boston employment centers are located in the established towns where proportion of owner-occupied units built prior to 1940 is comparatively high.

# scatterplot to visualize the relationship between AGE and INDUS
plt.figure(figsize=(6, 6))
sns.scatterplot(x = 'AGE', y = 'INDUS', data = df)
plt.show()

Observations:

• No trend between the two variables is visible in the above plot.

# scatterplot to visulaize the relationship between RAD and TAX
plt.figure(figsize=(6, 6))
sns.scatterplot(x = 'RAD', y = 'TAX', data = df)
plt.show()

Observations:

• The correlation between RAD and TAX is very high. But, no trend is visible between the two variables. This might be due to outliers.

Let's check the correlation after removing the outliers.

# remove the data corresponding to high tax rate
df1 = df[df['TAX'] < 600]
# import the required function
from scipy.stats import pearsonr
# calculate the correlation
print('The correlation between TAX and RAD is', pearsonr(df1['TAX'], df1['RAD'])[0])

The correlation between TAX and RAD is 0.24975731331429196

So the high correlation between TAX and RAD is due to the outliers. The tax rate for some properties might be higher due to some other reason.

# scatterplot to visualize the relationship between INDUS and TAX
plt.figure(figsize=(6, 6))
sns.scatterplot(x = 'INDUS', y = 'TAX', data = df)
plt.show()

Observations:

• The tax rate appears to increase with an increase in the proportion of non-retail business acres per town. This might be due to the reason that the variables TAX and INDUS are related with a third variable.

# scatterplot to visulaize the relationship between RM and MEDV
plt.figure(figsize=(6, 6))
sns.scatterplot(x = 'RM', y = 'MEDV', data = df)
plt.show()

Observations:

• The price of the house seems to increase as the value of RM increases. This is expected as the price is generally higher for more rooms.

• There are a few outliers in a horizontal line as the MEDV value seems to be capped at 50.

# scatterplot to visulaize the relationship between LSTAT and MEDV
plt.figure(figsize=(6, 6))
sns.scatterplot(x = 'LSTAT', y = 'MEDV', data = df)
plt.show()

Observations:

• The price of the house tends to decrease with an increase in LSTAT. This is also possible as the house price is lower in areas where lower status people live.
• There are few outliers and the data seems to be capped at 50.

We have seen that the variables LSTAT and RM have a linear relationship with the dependent variable MEDV. Also, there are significant relationships among a few independent variables, which is not desirable for a linear regression model. Let's first split the dataset.

Split the dataset

Let's split the data into the dependent and independent variables and further split it into train and test set in a ratio of 70:30 for train and test set.

# separate the dependent and independent variable
Y = df['MEDV_log']
X = df.drop(columns = {'MEDV', 'MEDV_log'})

# add the intercept term
X = sm.add_constant(X)

# splitting the data in 70:30 ratio of train to test data
X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=0.30 , random_state=1)

Next, we will check the multicollinearity in the train dataset.

Check for Multicollinearity

We will use the Variance Inflation Factor (VIF), to check if there is multicollinearity in the data.

Features having a VIF score > 5 will be dropped/treated till all the features have a VIF score < 5

from statsmodels.stats.outliers_influence import variance_inflation_factor

# function to check VIF
def checking_vif(train):
    vif = pd.DataFrame()
    vif["feature"] = train.columns

    # calculating VIF for each feature
    vif["VIF"] = [
        variance_inflation_factor(train.values, i) for i in range(len(train.columns))
    ]
    return vif


print(checking_vif(X_train))

    feature         VIF
0     const  535.372593
1      CRIM    1.924114
2        ZN    2.743574
3     INDUS    3.999538
4      CHAS    1.076564
5       NOX    4.396157
6        RM    1.860950
7       AGE    3.150170
8       DIS    4.355469
9       RAD    8.345247
10      TAX   10.191941
11  PTRATIO    1.943409
12    LSTAT    2.861881

Observations:

• There are two variables with a high VIF - RAD and TAX. Let's remove TAX as it has the highest VIF values and check the multicollinearity again.

Step 4: Drop the column 'TAX' from the training data and check if multicollinearity is removed?

# create the model after dropping TAX
X_train = X_train.drop(columns={'TAX'})

# check for VIF
print(checking_vif(X_train))

    feature         VIF
0     const  532.025529
1      CRIM    1.923159
2        ZN    2.483399
3     INDUS    3.270983
4      CHAS    1.050708
5       NOX    4.361847
6        RM    1.857918
7       AGE    3.149005
8       DIS    4.333734
9       RAD    2.942862
10  PTRATIO    1.909750
11    LSTAT    2.860251

Now, we will create the linear regression model as the VIF is less than 5 for all the independent variables, and we can assume that multicollinearity has been removed between the variables.

Step 5: Write the code to create the linear regression model and print the model summary.

# create the model
model1 = sm.OLS(y_train, X_train)

# fitting the model
model1 = model1.fit()

# get the model summary
model1.summary()

OLS Regression Results

Dep. Variable:	MEDV_log	R-squared:	0.769
Model:	OLS	Adj. R-squared:	0.761
Method:	Least Squares	F-statistic:	103.3
Date:	Mon, 22 Nov 2021	Prob (F-statistic):	1.40e-101
Time:	18:12:35	Log-Likelihood:	76.596
No. Observations:	354	AIC:	-129.2
Df Residuals:	342	BIC:	-82.76
Df Model:	11
Covariance Type:	nonrobust

	coef	std err	t	P>|t|	[0.025	0.975]
const	4.6324	0.243	19.057	0.000	4.154	5.111
CRIM	-0.0128	0.002	-7.445	0.000	-0.016	-0.009
ZN	0.0010	0.001	1.425	0.155	-0.000	0.002
INDUS	-0.0004	0.003	-0.148	0.883	-0.006	0.005
CHAS	0.1196	0.039	3.082	0.002	0.043	0.196
NOX	-1.0598	0.187	-5.675	0.000	-1.427	-0.692
RM	0.0532	0.021	2.560	0.011	0.012	0.094
AGE	0.0003	0.001	0.461	0.645	-0.001	0.002
DIS	-0.0503	0.010	-4.894	0.000	-0.071	-0.030
RAD	0.0076	0.002	3.699	0.000	0.004	0.012
PTRATIO	-0.0452	0.007	-6.659	0.000	-0.059	-0.032
LSTAT	-0.0298	0.002	-12.134	0.000	-0.035	-0.025

Omnibus:	30.699	Durbin-Watson:	1.923
Prob(Omnibus):	0.000	Jarque-Bera (JB):	83.718
Skew:	0.372	Prob(JB):	6.62e-19
Kurtosis:	5.263	Cond. No.	2.09e+03

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.09e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

Observations:

• Our R-Squared value is .769. This seem to be a good fit.
• Considering p-values, we can drop insignificant(p-value > .05).
• Look like only ZN , INDUS and AGE are the insignificant features. We will drop these.

Step 6: Drop insignificant variables from the above model and create the regression model again.

Examining the significance of the model

It is not enough to fit a multiple regression model to the data, it is necessary to check whether all the regression coefficients are significant or not. Significance here means whether the population regression parameters are significantly different from zero.

From the above it may be noted that the regression coefficients corresponding to ZN, AGE, and INDUS are not statistically significant at level α = 0.05. In other words, the regression coefficients corresponding to these three are not significantly different from 0 in the population. Hence, we will eliminate the three features and create a new model.

# create the model after dropping columns 'MEDV', 'MEDV_log', 'TAX', 'ZN', 'AGE', 'INDUS' from df dataframe
Y = df['MEDV_log']
X = df.drop(columns=['AGE','ZN','INDUS','MEDV_log','MEDV','TAX']) #write your code here
X = sm.add_constant(X)

#splitting the data in 70:30 ratio of train to test data
X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=0.30 , random_state=1)

# create the model
model2 = sm.OLS(y_train, X_train)

# fitting the model
model2 = model2.fit()

# get the model summary
model2.summary()

OLS Regression Results

Dep. Variable:	MEDV_log	R-squared:	0.767
Model:	OLS	Adj. R-squared:	0.762
Method:	Least Squares	F-statistic:	142.1
Date:	Tue, 23 Nov 2021	Prob (F-statistic):	2.61e-104
Time:	05:55:07	Log-Likelihood:	75.486
No. Observations:	354	AIC:	-133.0
Df Residuals:	345	BIC:	-98.15
Df Model:	8
Covariance Type:	nonrobust

	coef	std err	t	P>|t|	[0.025	0.975]
const	4.6494	0.242	19.242	0.000	4.174	5.125
CRIM	-0.0125	0.002	-7.349	0.000	-0.016	-0.009
CHAS	0.1198	0.039	3.093	0.002	0.044	0.196
NOX	-1.0562	0.168	-6.296	0.000	-1.386	-0.726
RM	0.0589	0.020	2.928	0.004	0.019	0.098
DIS	-0.0441	0.008	-5.561	0.000	-0.060	-0.028
RAD	0.0078	0.002	3.890	0.000	0.004	0.012
PTRATIO	-0.0485	0.006	-7.832	0.000	-0.061	-0.036
LSTAT	-0.0293	0.002	-12.949	0.000	-0.034	-0.025

Omnibus:	32.514	Durbin-Watson:	1.925
Prob(Omnibus):	0.000	Jarque-Bera (JB):	87.354
Skew:	0.408	Prob(JB):	1.07e-19
Kurtosis:	5.293	Cond. No.	690.

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Observations:

• We can see that the R-squared value has decreased by 0.002, since we have removed variables from the model, whereas the adjusted R-squared value has increased by 0.001, since we removed statistically insignificant variables only.

Now, we will check the linear regression assumptions.

Check the below linear regression assumptions

1. Mean of residuals should be 0
2. No Heteroscedasticity
3. Linearity of variables
4. Normality of error terms

Step 7: Check the above linear regression assumptions and provide insights.

Check for mean residuals

residuals = model2.resid

residuals.mean()

-2.7385501274087196e-15

Observations:

Mean of residulas is very close to zero.

Check for homoscedasticity

• Homoscedasticity - If the residuals are symmetrically distributed across the regression line, then the data is said to homoscedastic.

• Heteroscedasticity - If the residuals are not symmetrically distributed across the regression line, then the data is said to be heteroscedastic. In this case, the residuals can form a funnel shape or any other non-symmetrical shape.

• We'll use Goldfeldquandt Test to test the following hypothesis with alpha = 0.05:

  • Null hypothesis: Residuals are homoscedastic
  • Alternate hypothesis: Residuals have heteroscedastic

from statsmodels.stats.diagnostic import het_white
from statsmodels.compat import lzip
import statsmodels.stats.api as sms

name = ["F statistic", "p-value"]
test = sms.het_goldfeldquandt(y_train, X_train)
lzip(name, test)

[('F statistic', 1.083508292342528), ('p-value', 0.3019012006766869)]

Observations:

p-value is > .05, so we fail to reject NULL hypothesis and assume that residuals are homescedastic.

Linearity of variables

It states that the predictor variables must have a linear relation with the dependent variable.

To test the assumption, we'll plot residuals and fitted values on a plot and ensure that residuals do not form a strong pattern. They should be randomly and uniformly scattered on the x-axis.

# predicted values
fitted = model2.fittedvalues

# sns.set_style("whitegrid")
sns.residplot(x = fitted, y = residuals, color="lightblue", lowess=True) #write your code here
plt.xlabel("Fitted Values")
plt.ylabel("Residual")
plt.title("Residual PLOT")
plt.show()

Observations:

• We see that the residuals are evenly scattered and linear assumption is satisfied

Normality of error terms

The residuals should be normally distributed.

# Plot histogram of residuals
sns.histplot(residuals,kde=True)

<AxesSubplot:ylabel='Count'>

# Plot q-q plot of residuals
import pylab
import scipy.stats as stats

stats.probplot(residuals, dist="norm", plot=pylab)
plt.show()

Observations:

• Residuals are pretty much making straight line plot and sataisfy normality condition.

Check the performance of the model on the train and test data set

Step 8: Compare model performance of train and test dataset

# RMSE
def rmse(predictions, targets):
    return np.sqrt(((targets - predictions) ** 2).mean())


# MAPE
def mape(predictions, targets):
    return np.mean(np.abs((targets - predictions)) / targets) * 100


# MAE
def mae(predictions, targets):
    return np.mean(np.abs((targets - predictions)))


# Model Performance on test and train data
def model_pref(olsmodel, x_train, x_test):

    # In-sample Prediction
    y_pred_train = olsmodel.predict(x_train)
    y_observed_train = y_train

    # Prediction on test data
    y_pred_test = olsmodel.predict(x_test)
    y_observed_test = y_test

    print(
        pd.DataFrame(
            {
                "Data": ["Train", "Test"],
                "RMSE": [
                    rmse(y_pred_train, y_observed_train),
                    rmse(y_pred_test, y_observed_test),
                ],
                "MAE": [
                    mae(y_pred_train, y_observed_train),
                    mae(y_pred_test, y_observed_test),
                ],
                "MAPE": [
                    mape(y_pred_train, y_observed_train),
                    mape(y_pred_test, y_observed_test),
                ],
            }
        )
    )


# Checking model performance
model_pref(model2, X_train, X_test)  

    Data      RMSE       MAE      MAPE
0  Train  0.195504  0.143686  4.981813
1   Test  0.198045  0.151284  5.257965

Observations:

• We see that model perform pretty well across test data.
• RMSE, MAE and MAPE are quite consistent across train and test data.

Apply cross validation to improve the model and evaluate it using different evaluation metrics

Step 9: Apply the cross validation technique to improve the model and evaluate it using different evaluation metrics.

# import the required function

from sklearn.model_selection import cross_val_score

# build the regression model and cross-validate
linearregression = LinearRegression()                                    

cv_Score11 = cross_val_score(linearregression, X_train, y_train, cv = 10)
cv_Score12 = cross_val_score(linearregression, X_train, y_train, cv = 10,scoring = 'neg_mean_squared_error')                               


print("RSquared: %0.3f (+/- %0.3f)" % (cv_Score11.mean(), cv_Score11.std() * 2))
print("Mean Squared Error: %0.3f (+/- %0.3f)" % (-1*cv_Score12.mean(), cv_Score12.std() * 2))

RSquared: 0.729 (+/- 0.232)
Mean Squared Error: 0.041 (+/- 0.023)

Observations

• The R-squared on the cross validation is 0.729, whereas on the training dataset it was 0.769
• And the MSE on cross validation is 0.041, whereas on the training dataset it was 0.038

We may want to reiterate the model building process again with new features or better feature engineering to increase the R-squared and decrease the MSE on cross validation.

Step 10: Get model Coefficients in a pandas dataframe with column 'Feature' having all the features and column 'Coefs' with all the corresponding Coefs.

coef = model2.params
coef.shape

(9,)

# Let us write the equation of the fit
Equation = "log (Price) ="
print(Equation, end='\t')
for i in range(len(coef)):
    print('(', coef[i], ') * ', coef.index[i], '+', end = ' ')

log (Price) =	( 4.649385823266648 ) *  const + ( -0.012500455079104059 ) *  CRIM + ( 0.11977319077019681 ) *  CHAS + ( -1.0562253516683295 ) *  NOX + ( 0.058906575109278964 ) *  RM + ( -0.04406889079940562 ) *  DIS + ( 0.00784847460624374 ) *  RAD + ( -0.04850362079499892 ) *  PTRATIO + ( -0.029277040479796852 ) *  LSTAT + 

Step 11: Write the conclusions and business recommendations derived from the model.

Write Conclusions here

• We performed Univariate and BiVariate analysis on the data and check if independent and dependent variables satisfy linearity requirements.
• We performed VIF(Variable Inflation Faactor) test to check for multicolinearity and eliminated TAX as a feature sicne its VIF > 5.
• We then split data 70-30 between train and test and fitted the model on train data (model1).
• We reviewed model summary and observed some insignificant variables(p-value > .05). We eliminated AGE, ZN and INDUS features based on these stats.
• We then fitted the model again(model2) on train data.
• We then checked if the model satisfy 4 assumptions of linearity
• We then evaluated model performance based on train and test data and found MAPE,MAE and RMSE across train and test dataset are pretty consistent.
• We finally applied cross-validation on the model

Write Recommendations here

• We notice that house values are positively corelated with location near the river. More the house near the river, more the value.
• As exepcted , more the rooms in house and more the proximity to radial highways more the value of house.
• Look like house prices are heavily impacted by Nitric Oxide concentration. This can be understood by presence of old industrial units and there might be prior evidence of contamination.