Table of Contents

Inferential Statistics - Simple Linear Regression¶

This notebook explores various implementations of statistical inferences of a simple linear regression model.

As the derivations of the model and the associated statistics have been covered in numerous books, lectures and notebooks (see sources) we will focus instead on summary of the main formula and a case study and its Python implementation using three main libraries.

This notebook is divided into four main sections:

1. Introduction: notation, model description, formula, hypothesis
2. Implementation using Numpy and Pandas
3. Implementation using Statsmodel
4. Implementation using Scikit Learn

Introduction¶

Libraries¶

Dataset¶

The dataset illustrates the well known problem of trying to predict appartment prices given data on their surface.

• The explanatory (or independent, exogeneous) variable $x$ is the appartment surface area in $m^2$
• The dependent (or endogeneous) variable $y$ is the appartment price in French Francs

Notation¶

The following notation convention will be used throughout:

• $x_i$ is an observation of the explanatory variable
• $y_i$ is an observation of the dependent variable
• $\bar y$ and $\bar x$ are the mean of a set of observations
• $\hat y_i$ is the value predicted by the model, given the observation $x_i$

Decomposition of the sum of squares and notation¶

$SS_{Total} = SS_{Residual} + SS_{Explained}$

• $SS_{Total} = \large \sum_{i-1}^n(y_i - \bar y)^2$ also noted $\large S_{yy}$
• $SS_{Residual} = \large \sum_{i-1}^n(y_i - \hat y_i)^2$ also noted $\large \sum e_i^2$
• $SS_{Expl} = \large \sum_{i-1}^n(\hat y_i - \bar y)^2$ also noted $\large S_{\hat y \hat y}$

Simple linear regression¶

The model based on the true (unobservable) parameters is:

$y_i = \alpha + \beta x_i + \epsilon_i$

Given the observations, we estimate the model and the coefficients of regression with $a$ and $b$ such that:

$y_i = a + b x_i + \hat \epsilon_i$

$\hat y_i = a + bx_i$

$\hat \epsilon_i = y_i - \hat y_i$

Note that $\epsilon_i$ is the unobservable statistical error between the value of $y_i$ and the true population mean. The value $\hat \epsilon_i$ is the observable residual and is an estimate of the statistical error.

Coefficient of determination $R^2$¶

$R^2 = \large \frac{SS_{exp}}{SS_{tot}} = \frac{S_{\hat y \hat y}}{S_{yy}} = \Large \frac{\sum_{i-1}^n(\hat y_i - \bar y)^2}{\sum_{i-1}^n(y_i - \bar y)^2}$

Leverage¶

The leverage of an observcation $x_i$ can be viewed intuitively as a measure of the "influence" of the point on the regression model. It is defined as:

$h_i = \Large \frac{1}{n} + \frac{(x_i - \bar x)^2}{\sum(x_i - \bar x)^2}$

Leverage values are bounded between $1/n$ and $1$. The average is $2/n$ and we can consider that the leverage is significant if it is greater than $4/n$

Maximum Likelihood Model - Normally distributed errors¶

Model description¶

Given a sample of size $n$: $\{(Y_i, x_i)\}$ for $i = 1..n$ taken from a population described by the following model:

$Y_i = Y \ | \ x_i = \alpha + \beta x_i + \epsilon_i$

Where:

• $\alpha$ is an unknown constant describing the y intercept of the regression line (called constant term)
• $\beta$ is an unknown constant describing the slope of the regression line (called coefficient of regression)
• $x_i$ is the value of the explanatory variable
• $\epsilon_i$ is a random variable which follows a Normal distribution $\epsilon_i \ i.i.d. \ \sim N(0, \sigma^2)$. It is the statistical error i.e. the difference between the random variable $Y_i$ and the true population mean
• $Y_i$ is a random variable which also follows a Normal distribution $Y_i \ i.i.d. \ \sim N(\alpha + \beta x_i, \sigma^2)$

Assumptions¶

• $E(\epsilon_i) = 0 \ \forall \ i$ Expectation of the random variable $\epsilon_i$ is zero
• $var(\epsilon_i) = \sigma^2 \ \forall \ i$ Variance of the random variable $\epsilon_i$ is finite and defined
• $\epsilon_i$ are i.i.d. (independent, identically distributed) $\Rightarrow \forall \ i \neq j \ cov(\epsilon_i,\epsilon_j) = 0$
• $\epsilon_i$ are normally distributed: $\epsilon_i \ i.i.d. \ \sim N(0, \sigma^2)$
• $x_i$ are observed with no error

Estimators¶

Since $\alpha$, $\beta$ and $\sigma^2$ are unknown population parameters, we will estimate them using:

• B the estimator of $\beta$ with realization $b$ in a given sample
• A the estimator of $\alpha$ with realization $a$ in a given sample
• $S^2_{n-2}$ the estimator of $\sigma^2$ with realization $\hat \sigma^2$ in a given sample
• $\hat \epsilon_i$ or $e$ the residual, an estimator of the statistical error

Hence we have $Y_i = A + Bx_i + e_i$

Espectation and variance of the model¶

• $E(Y_i) = E(Y \ | \ x_i) = \alpha + \beta x_i$
• $Var(Y_i) = Var(Y \ | \ x_i) = \sigma^2$

Maximum Likelihood Estimators¶

Maximizing the likelihood function and solving the equations gives:

• $\large A = \bar y - B \bar x$
• $B = \large \frac{\sum (x_i - \bar x)(y_i - \bar y)}{\sum (x_i - \bar x)^2}=\frac{\sum y_i x_i - n \bar y \bar x}{\sum x_i^2 - n \bar x^2} = \frac{\text{Cov}(x,y)}{\text{Var}(x)}=\frac{S_{xy}}{S_{xx}}$

Properties of the estimators:¶

• Expectation: $E(A) = \alpha$ and $E(B) = \beta$
• Variance: $Var(A) = \sigma^2 \large \left(\frac{1}{n} + \frac{\bar x^2}{\sum (x_i -\bar x)^2} \right)$ and $Var(B) = \large \frac{\sigma^2}{\sum (x_i - \bar x)^2}$

Distribution of the estimators¶

Since A and B are linear combinations of i.i.d. normal random variables $Y_1, Y_2,...,Y_n$ they also follow a normal distribution:

• $A \sim N\left( \alpha, \sigma^2 \large \left(\frac{1}{n} + \frac{\bar x^2}{\sum (x_i -\bar x)^2} \right) \right)$ and $B \sim N \left(\beta, \large \frac{\sigma^2}{\sum (x_i - \bar x)^2} \right)$

It can be shown that the following estimator for the variance of the errors is unbiased and is calculated based on the observed residuals $e_i$:

• $S^2_{n-2} = \large \frac{1}{n-2} \sum(Y_i - A -B x_i)^2 = \frac{1}{n-2} \sum e_i^2$

Implications of the normal distribution hypothesis¶

Estimators can be Z standardized:

$\large \frac{A - E(A)}{\sqrt{Var(A)}} = \frac{A - \alpha}{\sqrt{\sigma^2 \large \left(\frac{1}{n} + \frac{\bar x^2}{\sum (x_i -\bar x)^2} \right)}} \small \sim N(0,1)$

$\large \frac{B - E(B)}{\sqrt{Var(B)}} = \frac{B - \beta}{\sqrt{\frac{\sigma^2}{\sum (x_i - \bar x)^2} }} \small \sim N(0,1)$

Replacing the true variance $\sigma^2$ with its estimator $S^2_{n-2}$ leads to a student $T$ distribution with $n-2$ degrees of freedom:

$\large \frac{A - \alpha}{S_A} \sim T_{n-2}$

$\large \frac{B - \beta}{S_B} \sim T_{n-2}$

Since A and B are independent of $S^2_{n-2}$

Under the assumption that $\beta = 0$ then the following statistics follow a Chi-squared distribution:

$\large \frac{\sum (\hat Y - \bar Y)^2}{\sigma^2} = \frac{SS_{exp}}{\sigma^2} \sim \chi_1^2$ and $\large \frac{\sum (Y_i - \bar Y)^2}{\sigma^2} =\frac{SS_{tot}}{\sigma^2} \sim \chi_{n-1}^2$

Also, the following follows a Chi-squared sitribution:

$\large \frac{\sum (Y_i - \bar Y_i)^2}{\sigma^2} = \frac{SS_{res}}{\sigma^2} \sim \chi^2_{n-2}$

And hence the ratio of the two follows a Fisher distribution

$\large F = \frac{SS_{exp}}{SS_{res} / (n-2)} \sim F_{1,n-2}$

Hypothesis testing for the parameters $\alpha$ and $\beta$¶

• $H_0 :\alpha = 0$
• $H_1: \alpha \neq 0$

Test statistic: $T_A = \frac{A}{S_A} \sim T_{n-2}$

• $H_0 :\beta = 0$
• $H_1: \beta \neq 0$

Test statistic: $T_B = \frac{B}{S_B} \sim T_{n-2}$

Hypothesis testing for the statistical significance of the model as a whole¶

• $H_0 : Y_i = \alpha + \epsilon_i$
• $H_1 : Y_i = \alpha + \beta x_i + \epsilon_i$

Test statistic: $\large F = \frac{SS_{exp}}{SS_{res} / (n-2)} \sim F_{1,n-2}$

Confidence intervals¶

Confidence interval around the mean of the regression line¶

For a given observation $x_0$ we want to construct a C.I. for $E(Y_0) = E(Y | x = x_0) = \alpha + \beta x_0$

• Expectation: $E[Y \ | \ x = x_0] = E(A) + E(B) x_0 = \alpha + \beta x_0$
• Variance: $\text{Var}[Y \ | \ x = x_0] = \sigma^2 \large \left( \frac{1}{n} + \frac{(x_0 - \bar x)^2 }{\sum_{i=1}^n (x_i - \bar x)^2} \right) = \small \sigma^2 h_0$ where $h_0$ is the leverage of point $x_0$.

Since $A + B x_0$ is a linear combination of Normal r.v. it also has a Normal distribution:

$A + B x_0 \sim N \left( \alpha + \beta x_0, \sigma^2 h_0 \right)$ and so $\large \frac{A + Bx_0 - (\alpha + \beta x_0)}{ \sigma \sqrt{h_0} } \small \sim N(0,1)$

Since $\sigma^2$ is unknown, we replace it by its estimator $S^2_{n-2}$ to obtain a pivotal quantity:

$\large \frac{A + Bx_0 - (\alpha + \beta x_0)}{ S_{n-2} \sqrt{h_0} } \small \sim T_{n-2}$

Re-arranging to build a confidence interval of level $1 - \gamma$ :

$IC_{1 - \gamma} [E(Y \ | \ x_0)] = [ A + B x_0 \pm t_{(n-2)(1 - \gamma /2)} S_{n-2} \sqrt{h_0} ]$

Prediction Confidence Interval around a new observation $x_p$¶

Here we are interested in the Expectation and Variance of the residual for a new observation (i.e. one which was not included in the model previously)

• Expectation of $E(Y_p - (A + B x_p)) = E(Y_p) - E(A + B x_p) = 0$
• Variance: $\text{Var}(Y_p - (A + B x_p)) = \text{Var}(Y_p) + \text{Var}(A + B x_p) = \sigma^2 \left(1 + h_p \right)$ Note here that the variance is greater than the variance in the previous C.I.

As previously, replacing the variance with $S_{n-2}^2$ and constructing the pivotal quantity:

$\large \frac{Y_p - (A + B x_p)}{S_{n-2} \sqrt{1 + h_p}} \sim T_{n-2}$

And the prediction confidence interval:

$IC_{1 - \gamma} [E(Y \ | \ x_p)] = [ A + B x_p \pm t_{(n-2)(1 - \gamma /2)} S_{n-2} \sqrt{1 + h_p} ]$

The uncertainty around the prediction interval is greater

Some comments on the residuals¶

If the hypotheses of the model are not violated (i.e. the errors are normally distributed etc..) then the following properties hold. Recall that residuals $e_i$ are the estimators of the true statistical error $\epsilon_i$

• Residuals are not independent: their sum is null
• The variance of the residuals is $\text{Var}(Y_i - (A + B x_i)) = \text{Var}(Y_i) + \text{Var}(A + Bx_i) - 2\text{Cov}(Y_i - (A + B x_i)) = \sigma^2 + \sigma^2 h_i - 2\sigma^2 h_i = \sigma^2(1 - h_i)$
• The variance can be estimated as $S^2_{n-2}(1 - h_i)$

Note: do not confuse the variance of the residuals of the sample $\sigma^2(1 - h_i)$ with the variance of the residuals of a new observation $x_p$ : $\sigma^2(1 + h_p)$

Cook's distance: a measure of influence¶

Defined as $D_i = \Large \frac{(e_i)^2}{2 \sigma^2} \left( \frac{h_i}{(1- h_i)^2} \right)$

Values larger than 1, 4 / n or 4 / (n - P -1) for multiple linear regression are considered to be large and point towards an influential observation

Implementation using Numpy and Pandas¶

Custom Python class¶

Note that Statsmodel is used for probability calculations

Initialize the model, perform the fitting and probability calculations and print the results in SAS format¶

Given the very small p-values for both the slope parameter $B$ and the model as a whole, we can confidently reject the null hypothesis and state that there is a statistically significant linear relationship between the explanatory variable $x$ and the dependent variable $y$

Print the resulting scatter plot and regression line with confidence intervals¶

Study of the residuals and influential points¶

Studying the residuals is fundamental as it allows to highlight observations which are

• Potential outliers
• Play an important role in the regression analysis

Moreover, it allows us to confirm the hypothesis of the model, in particular:

• Linearity
• Homoscedasticity
• Normal distribution of the errors $\epsilon_i$

Plotting the residual scatter plot and QQ plot¶

Comments: The residuals appear to be approximately symetrically distributed, with mean close to zero. The distribution is not perfectly normal, in particular the right tail is heavier than normal.

This points towards a violation of the hypothesis of heteroscedasticity, i.e. the residual increases as a function of x. A potential solution would be to perform linear regression on $log(Y)$ instead

Implementation using Statsmodel library¶

The python statsmodel library offers extensive functionalities and reports for the linear regression model. We will explore the main ones in this section

• As recommended by Statsmodel, import the statsmodel api to access the functions and models.
• Importing statsmodels.api will load most of the public parts of statsmodels. This makes most functions and classes conveniently available within one or two levels, without making the “sm” namespace too crowded.
• Use the add_constant utility function to add a column of 1s to the matrix of observation. This will allow the model to fit an intercept parameter
• Print results summary

http://www.statsmodels.org/dev/importpaths.html

Comments¶

As expected, the values for the parameters, test statistics, probabilities, R-squared values etc.. are the same as those calculated manually using numpy.

• The formatted table offers a more concise summary of the information
• Statsmodel also performs a number of tests by default: for example the Durbin Watson, Jarque Bera etc..

Study of influence¶

The Statsmodel result object has a built-in method which returns an Influence object with influence and outlier measures. This "influence" object contains many results, methods and methods which can be called as required. See below for a few interesting ones:

• .cooks_distance - Cooks distance
• .hat_matrix_diag - The leverage for observation points
• .summary_table() - Prints the available influence results
• .summary_frame() - Creates a DataFrame with all available influence results.

Documentation here:

http://www.statsmodels.org/dev/generated/statsmodels.regression.linear_model.OLSResults.get_influence.html

http://www.statsmodels.org/dev/generated/statsmodels.stats.outliers_influence.OLSInfluence.html#statsmodels.stats.outliers_influence.OLSInfluence

Influence plots¶

Influence plots show the externally studentized residuals vs the leverage of each observation as measured by the hat matrix (here this is only the leverage as there are no matrices in simple linear regression)

Externally studentized residuals are residuals that are scaled by their standard deviation where:

$var(\epsilon_i) = \sigma^2_i (1 - h_i)$

The influence of each point can be visualized by the criterion keyword argument. Options are Cook's distance and DFFITS, two measures of influence

Source: http://www.statsmodels.org/dev/examples/notebooks/generated/regression_plots.html#Influence-plots

Leverage-Resid2 Plot¶

Closely related to the influence_plot is the leverage-resid2 plot.

Implementation using Sklearn¶

There exists no R type regression summary report in sklearn. The main reason is that sklearn is used for predictive modelling / machine learning and the evaluation criteria are based on performance on previously unseen data (such as predictive r^2 for regression).

Nevertheless we can extract some key values and parameters using the built in methods and attributes, in particular the regression metrics:

Documentation: http://scikit-learn.org/stable/modules/model_evaluation.html#regression-metrics

Sources¶

Notebook Author: Xavier Bourret Sicotte

Date: May 2018

The dataset, formulas and insights were taken from the (excellent) course at CNAM in Paris:

• Course name: STA102 - Analyse des données : méthodes explicatives
• Prof: G. Russolillo
• Link1 , Link2

Other usefull resources are: Lecture 14 at Duke for detailed derivations of most formula