Stats with Python: Simple Linear Regression

We’ve seen several aspects of the correlation coefficient in the previous posts. The correlation coefficient treats two variables equally; they are symmetrical. When two variables are not symmetrical, that is, when you want to explain $y$ by $x$, correlation analysis alone is not sufficient. Instead, you might want to conduct a regression analysis.

The simplest approach, simple linear regression, considers a single explanatory variable (independent variable) $x$ for explaining the objective variable (dependent variable) $y$.

Least Square Estimates

How to determine the parameters $\beta_0$ and $\beta_1$ in the above equation? Given the paired data $\{(x_i,y_i)\}_{i=1}^n$, they are determined by the method of least squares. That is, they are chosen to minimize the sum of the squared error between the predicted $\hat{y}_i = b_0+b_1x_i$ and the actual $y_i$:

Therefore, the least squares estimates are:

, where

Proof

Considering $\sum_{i=1}^n(x_i-\bar{x})=0$ and $\sum_{i=1}^n(y_i-\bar{y})=0$,

From the above calculation, $\mathcal{L}(b_0,b_1)$ takes its minimum value when $b_1=S_{xy}/S_x^2$ and $b_0 = \bar{y}-b_1\bar{x}$.

Coefficient of Determination

Now we have the predicted values $\{\hat{y}_i\}_{i=1}^n$. How good are these predictions? To evaluate the goodness, coefficient of determination $R^2$ is frequently used.

The coefficient of determination is the ratio of ESS (explained sum of squares) to TSS (total sum of squares). As you may imagine from its notation, the coefficient of determination $R^2$ is the square of the Pearson correlation coefficient $r$.

Proof

Using the equation: $\hat{y}_i - \bar{y} = \beta_1(x_i - \bar{x})$,

Experiment

Lastly, let’s confirm that $R^2=r^2$, introducing how to use linear regression with Python. As done in the previous post, I generated 100 pairs of correlated random samples (x and y).

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

sns.set_style("darkgrid")

n = 100
x = np.random.rand(n)
y = x + 0.5*np.random.rand(n)

Scikit-learn implements linear regression as LinearRegression and coefficient of determination as r2_score. After fitting the model, we can plot the regression line like below:

from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score

model = LinearRegression()
model.fit(x.reshape(-1, 1), y.reshape(-1, 1))

sns.scatterplot(x, y)
plt.plot(x, model.predict(x.reshape(-1, 1)), color="k")

As expected, $R^2=r^2$ is confirmed.

r2_score(y, model.predict(x.reshape(-1, 1)))
# >> 0.7922606713476185

np.corrcoef(x, y)[0,1]**2
# >> 0.7922606713476184

References

[1] 倉田 博史, 星野 崇宏. “入門統計解析”（第3章）. 新世社. 2009.