Lab-7 (Simple Linear Regression)¶

(Intro. to Data Science: IIT Delhi - Abu Dhabi)¶

In this lab, we will experiment with the concepts we discussed for single-variable linear regression. Recall that this is the statistical method to model the linear relationship between the response variable (typically denoted $y$) and a single predictor variable (typically denoted $x$). The relationship is modeled as a line $y = \beta_0 + \beta_1 x + \varepsilon$, where $\varepsilon$ is the residual error. We will use real datasets for our experiments in this lab. A classical dataset used for linear regression consists of records on various aspects of a house (e.g., Living Area, Front Area, etc.), along with the sale price. We will use the sale price as the response variable and one of the recorded attributes as the predictor. This dataset is available in CSV format at the following location. Please download this dataset to your working directory for further analysis. We start by loading this dataset into a pandas dataframe.

In [ ]:
import pandas as pd
import numpy as np
from matplotlib.pyplot import subplots

housing = pd.read_csv('housing.csv')
housing.head(3)

For the first part of the analysis, we shall work with 'Gr Liv Area' as the prector and 'SalePrice' as the reponse.

In [ ]:
x = housing['Gr Liv Area']
y = housing['SalePrice']

Noting the the Pearson correlation should be the very first step in analysing the linear relationship withween the response and the predictor.

In [ ]:
print('The Pearson correlation beween Liv Area and Sale Price is: ', y.corr(x))

Task 7.1 (15 points)¶

The Pearson correlation $r$ closer to 1 indicates a linear relationship between our variables. The next step is to visualise the relationship using a scatter plot. You are asked to write a Python function that takes as input the predictor and response variables and produces a scatter plot in standard units, and also draws the line with slope $r$ at the origin (i.e., $(0, 0)$).

In [ ]:
# Write a function that outputs a scatter plot y versus x 
# after conversion to standard units. It also plots a line 
# through the origin (i.e., x=0, y=0) with slope equal to the 
# Pearson Correlation r. Some of the other inputs are described below:
#   - xlabel and ylabel: String labels of y and x axis 
#   - plot_title: String title of the plot
def plot_scatter_in_standard(x, y, xlabel, ylabel, plot_title):
    # Calculate the Pearson correlation r
    r = y.corr(x)
    # Fill in the other details below
    #
    #

Once you have completed the above function, we will use it to see a scatter plot using the command below:

In [ ]:
#plot_scatter_in_standard(x, y, 'Gr Liv Area', 'Sale Price', 'price-vs-area')

Task 7.2 (10 points)¶

In the class, we saw how to plot the regression line within the scatter plot in the original units. You are asked to use the derivations for the slope and intercept for the regression line to write a function that produces the scatter plot in original units along with the regression line, where the slope and intercept has been calculated using the formulas derived in class. Please do not use polyfit.

In [ ]:
# This function should output the scatter plot in original units along with the regression line
def plot_scatter_in_original(x, y, xlabel, ylabel, plot_title):
    # Calculate the Pearson correlation r
    r = y.corr(x)
    # Fill in the other details below
    #
    #

Once you have completed the above function, we will use it to see a scatter plot using the command below:

In [ ]:
#plot_scatter_in_original(x, y, 'Gr Liv Area', 'Sale Price', 'price-vs-area')

Task 7.3 (20 points)¶

You will observe that certain data points are far from the regression line and can be regarded as 'outliers' that impact how well the regression line fits the remaining points. Even though some data points stand out as outliers, the important question is whether there is a way to quantify the notion of outliers. A more formal notion can be used to automate the process of finding the outliers in the data. For this, we define the Cook's distance ($D_i$), which is a diagnostic measure in least-squares regression that quantifies the influence of a single data point on the model's predicted values. The key idea is to define a measure of the impact of removing the $i^{th}$ data point. The formula for the Cook's distance is given below: $$D(i) = \frac{\sum_{j=1}^{n} (\hat{y}_j - \hat{y}_{(i)j})^2}{2 \cdot MSE}$$ where $\hat{y}_j$'s are the fitted values when using the entire sample, $\hat{y}_{(i)j}$'s are the fitted values when using the sample without the $i^{th}$ data point, and $MSE$ is the mean squared error, i.e., $MSE = \frac{1}{n} \cdot \sum_{j=1}^{n} (y_i - \hat{y}_i)^2$. The rule of thumb is to consider a data point $i$ to be an outlier if $D(i) > \frac{4}{n}$. However, this is only a high-level guideline. We will use a fixed $threshold = 0.05$ to mark outlier points. Check $D(i) > threshold$ to isolate the outlier points. You are asked to mark these points in a different color in the scatter plot of 'Sale Price' versus 'Gr Liv Area'. Follow the instructions given as comments in the code block below:

In [ ]:
# You are allowed to use the norm function of linalg
from numpy import linalg as LA

# Define a function that when given x and y values returns a vector of
# Cook's distances. For this, you are allowed to use numpy's polyfit function
def calculate_cook(x, y):
    # This function must return a numpy array with Cook's distances
    n = x.shape[0]
    D = np.zeros(n)
    # Fill in the rest of the function below
    #
    #
    return(D)


# This function gives a list of indices of data points that are outliers
# based on the criteria defined, i.e., D[i] > threshold
def return_outlier_indices(D, threshold):
    # Write the function below
    #
    #
    return([])# This is dummy code. Please replace this


# Once the above functions are defined, we will use the following command to find the 
# list of outlier indices. The threshold used is 0.05
outlier_indices = return_outlier_indices(calculate_cook(x, y), 0.05)

# We then provide a scatter plot with the outlier points colored in red.
# The command is given below.
fig, ax = subplots(figsize = (5, 5))
ax.scatter(x.drop(outlier_indices), y.drop(outlier_indices), c='blue', alpha = 0.5);

Let us now compare the Pearson coefficients before and after dropping the outliers

In [ ]:
print('Correlation before dropping outliers', x.corr(y))
print('Correlation after dropping outliers', x.drop(outlier_indices).corr(y.drop(outlier_indices)))

Henceforth, we shall be working with the dataset without the outlier points.

In [ ]:
x_wo_outliers = x.drop(outlier_indices)
y_wo_outliers = y.drop(outlier_indices)

# Let us give a scatter plot of the points without outliers. 
# Fortunately, we now have a function to do this. We just need to make the function call
plot_scatter_in_original(x_wo_outliers, y_wo_outliers, 'Liv Area', 'Sale Price',  'price-vs-area')

Task 7.4 (10 points)¶

Let's now use visual diagnostics to evaluate the linear prediction of the response in terms of the predictor. For this we first do a scatter plot of the residuals.

In [ ]:
# Write the function below that outputs the residuals versus the predictor scatter plot
def plot_residuals(x, y, xlabel, ylabel, plot_title):
    # Write the function below
    #
    #

# We will use the above-defined function to plot the residuals
#plot_residuals(x_wo_outliers, y_wo_outliers, 'Liv Area', 'Sale Price', 'price-vs-area')

We notice some Heteroscedasticity in the residual scatter plot. There is an uneven variation around the horizontal line. The variance of the residuals seems to increase as we go towards the right. Let us try using the log transformation to remove the Heteroscedasticity. So, instead of plotting the residuals with respect to $x$ and $y$, we will plot the residuals with respect to $x$ and $\log{y}$.

In [ ]:
#plot_residuals(x_wo_outliers, np.log(y_wo_outliers), 'Liv Area', 'log(Sale Price)',  'price-vs-area')

Task 7.5 (15 points)¶

We now observe that the log transformation makes it more appropriate to apply the regression model, where the errors have a normal distribution with a constant variance across values of the predictor variable. This means we can apply confidence intervals for prediction in the log(price)-versus-area model. Suppose we want a 95% CI for the Sale Price prediction based on this linear regression model with $\log{y}$ as the response and $x$ as the predictor. Use the bootstrap method to give a 95% confidence interval on the sale price of a house with an area of $2300$. Use $5000$ repetitions in your bootstrap method.

In [ ]:
# Write the code to implement the bootstrap method and use it to build a 95% confidence interval on the price of a house with an area of 2300.

The discussion on single-variable regression allowed us to highlight the key ideas. However, in applied settings, the goal is to predict the response variable from the predictors, and in many cases, using a linear combination of multiple predictors yields much better predictions than a single variable. Let us consider one such example.

Let us work with 'Gr Liv Area' and 'Lot Frontage' as the predictors for the 'SalePrice' as the response. We will try out various linear combinations of the predictors to create a single variable predictor. We then try to find out the linear combination that gives the highest correlation.

In [ ]:
x1 = housing['Gr Liv Area']
x2 = housing['Lot Frontage']
y = housing['SalePrice']

We note that the 'Lot Frontage' is not defined for a lot of houses. We remove those entries for the rest of the analysis.

In [ ]:
# Let us check how many houses have 'Lot Frontage' undefined
print('Number of undefined entries = ', np.sum(np.isnan(x2)))
In [ ]:
# Indices where the value is not defined
nan_indices = np.where(np.isnan(x2))[0]

# We drop the value at these indices
x1 = x1.drop(nan_indices).reset_index(drop=True)
x2 = x2.drop(nan_indices).reset_index(drop=True)
y = y.drop(nan_indices).reset_index(drop=True)

We'll use a linear combination of $x1$ and $x2$ as our predictor, i.e., we will use $x = p \cdot x1 + (1-p)\cdot x2$. We shall use the value of $p$ such that $x$ has a good correlation with $y$. By some trial-and-error, we arrive at the value of $p = 0.2$. So, we shall use $x = 0.2 \cdot x1 + 0.8 \cdot x2$ as our predictor variable.

In [ ]:
x = 0.2 * x1 + 0.8 * x2

Let us note down the Pearson correlation of this variable, $x = 0.2 \cdot x1 + 0.8 \cdot x2$, with the reponse, $y$

In [ ]:
print('Correlation of combined predictor with response: ', x.corr(y))

Note that the correlation is higher than with $x1$ alone.

Task 7.6 (25 points)¶

Repeat the analysis in Tasks 7.3, 7.4, and 7.5. After removing the outliers (using the threshold of 0.05) and using log transformation as in Task 7.5, you are asked to give a 95% CI on the prediction when the predictor variable $x = 600$.

In [ ]:
# Code cell for implementing Task 7.6. Please create more code/markup cells as required to complete this task.
#
#
#
#
#

Towards Multiple Linear Regression¶

Instead of trying to find an appropriate combined predictor variable for doing simple linear regression, we can perform mulitple linear regression, where we predict the response in terms of a linear combination of the predictors. For example, if there are two predictors $x_1$ and $x_2$, we predict $y$ using $\theta_0 + \theta_1 \cdot x_1 + \theta_2 \cdot x_3$. The coefficients $\theta_0, \theta_1, \theta_2$ are determined by the best fit plane through the data. We shall discuss how to go about finding these coefficients and other key aspects of multiple regression (such as which variables to choose for prediction) in the class soon. We close this lab by noting the tools that are available for determining the coefficients of multiple linear regression obtained by fitting the best plane through the data.

In [ ]:
from sklearn.linear_model import LinearRegression

# We first construct the predictor table where each column denotes the values of a predictor variable
# In our example, we will construct a table with x1 and x2 as columns
X = np.array([x1, x2])
X = X.T
print(X.shape)

reg = LinearRegression().fit(X, y)
print('Coefficients theta1 and theta 2:', reg.coef_)
print('Intercept theta0:', reg.intercept_)

So, we will use $-8794.098942603741 + 110.04890039 \cdot x1 + 361.61828017 \cdot x2$ as the predictor for the Sale Price. Let us quickly do a sanity check on our prediction. Let us see how good the prediction is for the first few data items.

In [ ]:
for i in range(10):
    print('response = ', y[i], end=', ')
    prediction = -8794.098942603741 + 110.04890039 * X[i][0] + 361.61828017 * X[i, 1]
    print('prediction = ', prediction, end = ', ')
    print('difference = ', (y[i] - prediction))

Task 7.7 (5 points)¶

Write code that calculates the Mean Squared Error for the above multiple linear predictor.

In [ ]:
# Write code for Task 7.7
#
#
In [ ]: