In [1]:
from lec_utils import *

Lecture 15 Supplementary Notebook¶

Simple Linear Regression¶

EECS 398-003: Practical Data Science, Fall 2024¶

practicaldsc.org • github.com/practicaldsc/fa24

Note: This notebook is only a supplementary notebook to the main lecture slides, which are in PDF form.

The main lecture slides can be found at practicaldsc.org under Lecture 15. (After the live lecture, an annotated version of the slides will be made available as well.)

Understanding the Data¶

Let's load in the commute times dataset as a DataFrame.

In [2]:
df = pd.read_csv('data/commute-times.csv')
df.head()
Out[2]:
date day home_departure_time home_departure_mileage ... mileage_to_work minutes_to_home work_departure_time_hr mileage_to_home
0 5/15/2023 Mon 2023-05-15 10:49:00 15873.0 ... 53.0 72.0 17.17 53.0
1 5/16/2023 Tue 2023-05-16 07:45:00 15979.0 ... 54.0 NaN NaN NaN
2 5/22/2023 Mon 2023-05-22 08:27:00 50407.0 ... 54.0 82.0 15.90 54.0
3 5/23/2023 Tue 2023-05-23 07:08:00 50535.0 ... 54.0 NaN NaN NaN
4 5/30/2023 Tue 2023-05-30 09:09:00 50664.0 ... 54.0 76.0 17.12 54.0

5 rows × 18 columns

There are many columns in here, but the only ones we're interested in for now are 'departure_hour' and 'minutes'.

In [3]:
df[['departure_hour', 'minutes']]
Out[3]:
departure_hour minutes
0 10.82 68.0
1 7.75 94.0
2 8.45 63.0
... ... ...
62 7.58 68.0
63 7.45 90.0
64 7.60 83.0

65 rows × 2 columns

In [4]:
fig = px.scatter(df,
           x='departure_hour',
           y='minutes',
           size=np.ones(len(df)) * 50,
           size_max=8)
fig.update_xaxes(title='Home Departure Time (AM)')
fig.update_yaxes(title='Minutes')
fig.update_layout(title='Commuting Time vs. Home Departure Time')
fig.update_layout(width=700)

Implementing $w_0^*$ and $w_1^*$¶

Let's implement the formulas for the best slope, $w_1^*$, and intercept, $w_0^*$, we found in the lecture slides:

\begin{align*} w_1^* &= \frac{ \displaystyle \sum_{i=1}^n (x_i - \bar x)(y_i - \bar y) }{ \displaystyle \sum_{i=1}^n (x_i - \bar x)^2 } \qquad \qquad w_0^* = \bar y - w_1^* \bar x \end{align*}
In [5]:
def slope(x, y):
    # Assume x and y are two Series.
    numerator = ((x - np.mean(x)) * (y - np.mean(y))).sum()
    denominator = ((x - np.mean(x)) ** 2).sum()
    return numerator / denominator
def intercept(x, y):
    return y.mean() - slope(x, y) * x.mean()
In [6]:
w1_star = slope(df['departure_hour'], df['minutes'])
w1_star
Out[6]:
-8.186941724265552
In [7]:
w0_star = intercept(df['departure_hour'], df['minutes'])
w0_star
Out[7]:
142.4482415877287

The results above tell us that the linear hypothesis function with the lowest mean squared error on our dataset is:

$$\text{predicted commute time (minutes)} = 142.45 - 8.19 \cdot \text{departure hour}$$

We can use it to make predictions:

In [8]:
def predict_commute(x_new):
    return w0_star + w1_star * x_new

What if I leave at 8AM? 10:45AM?

In [9]:
predict_commute(8)
Out[9]:
76.95270779360428
In [10]:
predict_commute(10 + 45 / 60)
Out[10]:
54.438618051874016

What do all of our predictions look like?

In [11]:
hline = px.line(x=[5.5, 11.5], y=[predict_commute(5.5), predict_commute(11.5)]).update_traces(line={'color': 'red', 'width': 4})
fline1 = go.Figure(fig.data + hline.data)
fline1.update_xaxes(title='Home Departure Time (AM)')
fline1.update_yaxes(title='Minutes')
fline1.update_layout(title='<span style="color:red">Predicted Commute Time</span> = 142.25 - 8.19 * Departure Hour')
fline1.update_layout(width=700, margin={'t': 60})

Aside: What does $R_{\text{sq}}(w_0, w_1)$ look like?¶

Let's draw a plot of $R_{\text{sq}}(w_0, w_1)$, the empirical risk that we're trying to minimize.

  • When we only had a single parameter, $h$, $R(h)$ was in 2D.
    • One axis for $h$, one axis for $R(h)$.
  • Now that we have two parameters, $w_0$ and $w_1$, $R(w_0, w_1)$ will be in 3D!
    • One axis for $w_0$, one axis for $w_1$, one axis for $R(w_0, w_1)$.
    • The bottom plane consists of all possible combinations of slope and intercept.
    • The height of the function above any pair of points on the bottom plane represents the MSE for that combination of slope and intercept.
In [12]:
def mse(y_actual, y_pred):
    return np.mean((y_actual - y_pred)**2)
def mse_for_departure_model(w):
    w0, w1 = w
    return mse(df['minutes'], w0 + w1 * df['departure_hour'])
num_points = 50 # increase for better resolution, but it will run more slowly. 
# if (num_points <= 100):
uvalues = np.linspace(90, 190, num_points)
vvalues = np.linspace(-13, -3, num_points)
(u,v) = np.meshgrid(uvalues, vvalues)
thetas = np.vstack((u.flatten(),v.flatten()))
MSE = np.array([mse_for_departure_model(t) for t in thetas.T])
loss_surface = go.Surface(x=u, y=v, z=np.reshape(MSE, u.shape))
minimizer = go.Scatter3d(x=[w0_star], y=[w1_star], z=[mse_for_departure_model([w0_star, w1_star])], 
                         mode='markers', name='optimal parameters',
                         marker=dict(size=10, color='gold'))
fig = go.Figure(data=[loss_surface, minimizer])
# fig.add_trace(opt_point)
fig.update_layout(title='Loss Surface', scene = dict(
    xaxis_title = "w0",
    yaxis_title = "w1",
    zaxis_title = r"R(w0, w1)"))
fig.show()
# else:
#     print("Picking num points > 100 can be really slow. If you really want to try, edit the code above so that this if statement doesn't trigger.")

We used partial derivatives to minimize the graph above!

Correlation¶

$$\begin{align*} r &= \text{the average of the product of $x$ and $y$, when both are standardized} \\ &= \frac{1}{n} \sum_{i = 1}^n \left( \frac{x_i - \bar{x}}{\sigma_x} \right) \left( \frac{y_i - \bar{y}}{\sigma_y} \right) \end{align*}$$
In [13]:
def correlation(x, y): 
    x_su = (x - np.mean(x)) / np.std(x)
    y_su = (y - np.mean(y)) / np.std(y)
    return np.mean(x_su * y_su)
In [14]:
correlation(df['departure_hour'], df['minutes'])
Out[14]:
-0.6486426165832002
In [15]:
# Symmetric!
correlation(df['minutes'], df['departure_hour'])
Out[15]:
-0.6486426165832002
In [16]:
# Doesn't change if we multiply x or y by constants!
correlation(df['departure_hour'] * 1000, df['minutes'] * 545)
Out[16]:
-0.6486426165832
In [17]:
# DataFrames have a built-in correlation method.
df[['departure_hour', 'minutes']].corr()
Out[17]:
departure_hour minutes
departure_hour 1.00 -0.65
minutes -0.65 1.00
In [18]:
# numpy has a built-in corrcoef method.
np.corrcoef(df['departure_hour'], df['minutes'])
Out[18]:
array([[ 1.  , -0.65],
       [-0.65,  1.  ]])

Implementing $w_0^*$ and $w_1^*$, Again¶

Recall, the formulas for the optimal intercept and slope are:

$$w_1^* = r \frac{\sigma_y}{\sigma_x}$$$$w_0^* = \bar{y} - w_1^* \bar{x}$$

Let's define two new functions, slope_again and intercept_again, which use these slightly updated formulas. (Really, only the formula for $w_1^*$ has changed.)

In [19]:
def slope_again(x, y):
    return correlation(x, y) * np.std(y) / np.std(x)
In [20]:
def intercept_again(x, y):
    return y.mean() - slope_again(x, y) * x.mean()
In [21]:
w1_star_again = slope_again(df['departure_hour'], df['minutes'])
w1_star_again
Out[21]:
-8.186941724265553
In [22]:
w0_star_again = intercept_again(df['departure_hour'], df['minutes'])
w0_star_again
Out[22]:
142.44824158772872

We get the same optimal intercept and slope as before!

In [23]:
# From before:
(w1_star, w0_star)
Out[23]:
(-8.186941724265552, 142.4482415877287)
In [24]:
# Now:
(w1_star_again, w0_star_again)
Out[24]:
(-8.186941724265553, 142.44824158772872)

Implementing $w_0^*$ and $w_1^*$ using sklearn¶

In practice, you wouldn't manually implement formulas for $w_0^*$ and $w_1^*$. Instead, you'd use a pre-built implementation.

The Python package we'll use for machine learning is sklearn. We'll start seeing it more in lectures next week.

In [25]:
from sklearn.linear_model import LinearRegression

To build a linear regression model that we can use for prediction, we first need to instantiate a LinearRegression object.

In [26]:
model = LinearRegression()

Then, we need to fit the model by telling it what our $x$'s and $y$'s are.

In [27]:
model.fit(X=df[['departure_hour']], y=df['minutes'])
Out[27]:
LinearRegression()
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
LinearRegression()

Once the model is fit, we can look at its intercept_ and coef_ attributes to see $w_0^*$ and $w_1^*$, respectively.

In [28]:
model.intercept_
Out[28]:
142.4482415877287
In [29]:
model.coef_
Out[29]:
array([-8.19])

These are exactly the same as we found with our manual calculations! This means that sklearn is doing the same three step modeling process that we outlined in lecture.

Now that model is fit, we can use it for making predictions:

In [30]:
# We'll discuss this warning more in coming lectures.
model.predict([[8]])

/Users/surajrampure/miniforge3/envs/pds/lib/python3.10/site-packages/sklearn/base.py:493: UserWarning:

X does not have valid feature names, but LinearRegression was fitted with feature names
Out[30]:
array([76.95])
In [31]:
# Using our hand-build predict_commute function from earlier in the lecture:
predict_commute(8)
Out[31]:
76.95270779360428

Aside: Pitfalls of Correlation¶

In [32]:
anscombe = pd.read_csv('data/anscombe.csv')
In [33]:
plt.figure(figsize=(12, 10))
for i, n in enumerate(['I', 'II', 'III', 'IV']):
    rows = anscombe[anscombe.get('dataset') == n]
    x = rows['x']
    y = rows['y']
    plt.subplot(2, 2, i+1)
    plt.scatter(x, y, label=f'Dataset {n}', alpha=0.65, s=65)
    plt.title(f'Dataset {n}');
No description has been provided for this image

What do all four of these datasets have in common?

In [34]:
for i, n in enumerate(['I', 'II', 'III', 'IV']):
    rows = anscombe[anscombe.get('dataset') == n]
    x = rows['x']
    y = rows['y']
    r = correlation(x, y)
    outstr = f'''
    <b>Dataset {n}</b><br>
    $\\bar x$: {np.round(np.mean(x), 2)}<br>
    $\\bar y$: {np.round(np.mean(y), 2)}<br>
    $\\sigma_x$: {np.round(np.std(x), 2)}<br>
    $\\sigma_y$: {np.round(np.std(y), 2)}<br>
    $r$: {np.round(r, 2)}
    '''
    display(HTML(outstr))
Dataset I
$\bar x$: 9.0
$\bar y$: 7.5
$\sigma_x$: 3.16
$\sigma_y$: 1.94
$r$: 0.82
Dataset II
$\bar x$: 9.0
$\bar y$: 7.5
$\sigma_x$: 3.16
$\sigma_y$: 1.94
$r$: 0.82
Dataset III
$\bar x$: 9.0
$\bar y$: 7.5
$\sigma_x$: 3.16
$\sigma_y$: 1.94
$r$: 0.82
Dataset IV
$\bar x$: 9.0
$\bar y$: 7.5
$\sigma_x$: 3.16
$\sigma_y$: 1.94
$r$: 0.82

They all share the exact same mean and standard deviation of $x$ and $y$, and the same correlation coefficient $r$! This means they all have the same best linear hypothesis function, in the sense of minimizing squared loss.

However, that linear hypothesis function looks better for some datasets than it does for others:

In [35]:
plt.figure(figsize=(12, 10))
for i, n in enumerate(['I', 'II', 'III', 'IV']):
    rows = anscombe[anscombe.get('dataset') == n]
    x = rows['x']
    y = rows['y']
    w0_ans = intercept(x, y)
    w1_ans = slope(x, y)
    plt.subplot(2, 2, i+1)
    plt.scatter(x, y, label=f'Dataset {n}', alpha=0.65, s=65)
    plt.plot(x, w0_ans + w1_ans * x, color='red');
    plt.title(f'Dataset {n}');
No description has been provided for this image

Moral of the story – visualize your data before trying to fit a prediction rule!

Another example of this phenomenon is the Datasaurus Dozen 🦕.