Notebook

Lab 7: Regression¶

Welcome to Lab 7!

Today we will get some hands-on practice with linear regression. You can find more information about this topic in section 15.2.

In [ ]:

# Run this cell, but please don't change it.

# These lines import the Numpy and Datascience modules.
import numpy as np
from datascience import *

# These lines do some fancy plotting magic.
import matplotlib
%matplotlib inline
import matplotlib.pyplot as plots
plots.style.use('fivethirtyeight')
import warnings
warnings.simplefilter('ignore', FutureWarning)

from scipy.stats import pearsonr as cor

import scipy.stats as stats

import scipy.stats

1. How Faithful is Old Faithful?¶

Old Faithful is a geyser in Yellowstone National Park that is famous for eruption on a fairly regular schedule. Run the cell below to see Old Faithful in action!

In [ ]:

# For the curious: this is how to display a YouTube video in a
# Jupyter notebook.  The argument to YouTubeVideo is the part
# of the URL (called a "query parameter") that identifies the
# video.  For example, the full URL for this video is:
#   https://www.youtube.com/watch?v=wE8NDuzt8eg
from IPython.display import YouTubeVideo
YouTubeVideo("wE8NDuzt8eg")

Some of Old Faithful's eruptions last longer than others. Whenever there is a long eruption, it is usually followed by an even longer wait before the next eruption. If you visit Yellowstone, you might want to predict when the next eruption will happen, so that you can see the rest of the park instead of waiting by the geyser.

Today, we will use a dataset on eruption durations and waiting times to see if we can make such predictions accurately with linear regression.

The dataset has one row for each observed eruption. It includes the following columns:

duration: Eruption duration, in minutes
wait: Time between this eruption and the next, also in minutes

Run the next cell to load the dataset.

In [ ]:

faithful = Table.read_table("faithful.csv")
faithful

We would like to use linear regression to make predictions, but that won't work well if the data aren't roughly linearly related. To check that, we should look at the data.

Question 1.1. Make a scatter plot of the data. It's conventional to put the column we want to predict on the vertical axis and the other column on the horizontal axis. Add fit_line= True to make the regression line appear on the plot.

In [ ]:

...

Question 1.2. Are eruption duration and waiting time roughly linearly related based on the scatter plot above? Is this relationship positive?

Write your answer here, replacing this text.

We're going to continue with the assumption that they are linearly related, so it's reasonable to use linear regression to analyze this data.

We'd next like to plot the data in standard units. If you don't remember the definition of standard units, textbook section 14.2 might help!

Question 1.3. Use the means and standard deviations provided and convert both the eruption durations and waiting times into standard units. Then create a table called faithful_standard containing the eruption durations and waiting times in standard units. The columns should be named duration (standard units) and wait (standard units).

$$SU = \frac{x - \overline{x}}{S_x}$$

In [ ]:

import scipy.stats as stats

duration_mean = np.mean(faithful.column("duration"))
duration_std = stats.tstd(faithful.column("duration"))
wait_mean = np.mean(faithful.column("wait"))
wait_std = stats.tstd(faithful.column("wait"))

duration_in_su = ...
wait_in_su = ...

faithful_standard = Table().with_columns("duration (standard units)",
                                        ...,
                                        "wait (standard units)",
                                        ...)

faithful_standard

Question 1.4. Plot the data again, but this time in standard units.

In [ ]:

...

You'll notice that this plot looks the same as the last one! However, the data and axes are scaled differently. So it's important to read the ticks on the axes.

Recall that $$r = \frac{1}{n-1} \sum_{i=1}^n\left(\frac{x_i-\overline{x}}{S_x}\right)\left(\frac{y_i-\overline{y}}{S_y}\right)$$

So if x and y are in standard unit, then $$r = \frac{1}{n-1} \sum_{i=1}^nx_iy_i$$

Question 1.5. Use this fact to compute the correlation r from the faithful_standard table.

Hint: Use faithful_standard. Section 15.1 explains how to do this.

In [ ]:

correlation = ...
correlation

Question 1.6. Use the function cor (which was imported from the scipy.stats module) to compute the correlation between wait and duration.

In [ ]:

r = ...
r

Notice that the pearsonr (or cor) returns two values. The second number is the p-value for a hypothesis test. What hypothesis test? This is an aside, but here's the details.

$H_o: \rho = 0$ The correlation is 0.

$H_a: \rho \not= 0$ The correlation is not 0.

If the p-value is small (smaller than $\alpha$), then we reject the null hypothesis, meaning we conclude that the correlation is NOT 0. This means that there is some relationship between the two variables.

2. The regression line¶

Recall that the Root Mean Squared Error is

$$RMSE = \sqrt{MSE} = \sqrt{\frac{1}{n-2} SSE} = \sqrt{\frac{1}{n-2} \sum_{i = 1}^n (y_i - \hat{y}_i)^2}$$

The linear regression line is the line that has the least possible value of the RMSE.

Also, recall that the correlation is the slope of the regression line when the data are put in standard units.

The next cell shows a scatterplot in standard units:

$$\text{waiting time in standard units} = r \times \text{eruption duration in standard units}.$$

Use the sliders to adjust the line shown. Find the slope/correlation by trying to find the smallest value of the RMSE, which is also displayed on the graph. (Adjust the slope and intercept slowly, there is a lag between when you move the sliders and when the RMSE is updated.)

In [ ]:

y1 = faithful.column('wait')
x1 = faithful.column('duration')

xx = (x1-np.average(x1))/stats.tstd(x1)
yy = (y1-np.average(y1))/stats.tstd(y1)

import ipywidgets as widgets

from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets

import matplotlib.pyplot as plt
%matplotlib inline

import numpy as np

def plot_func(correlation, intercept):
    x = np.linspace(-3.5, 3.5)
    y = x*correlation + intercept
    plt.scatter(xx,yy, s=50, color="green")
    plt.ylim(-2.5,2.5)
    plt.xlim(-2,2) 
    rmse = np.round((sum((yy - correlation*xx - intercept)**2)/(len(xx)-2))**0.5,2)
    plt.plot(x, y)
    plt.text(-1.5,1, f"RMSE = {rmse}", color="red", size='large')
    plt.ylabel("Duration")
    plt.xlabel("Wait")

interact(plot_func, correlation = widgets.FloatSlider(value=0.9001, min=0, max=1.0, step=0.01), 

         intercept=widgets.FloatSlider(value=0, min=-2, max=2, step=0.050));

Question 2.1 What is the least value of the RMSE?

In [ ]:

least_rmse_standard_units = ...
least_rmse_standard_units

How would you take a point in standard units and convert it back to original units? We'd have to "stretch" its horizontal position by duration_std and its vertical position by wait_std. That means the same thing would happen to the slope of the line.

Stretching a line horizontally makes it less steep, so we divide the slope by the stretching factor. Stretching a line vertically makes it more steep, so we multiply the slope by the stretching factor.

That is the slope of a regression line is: $\hat{\beta}_1 = r \frac{S_y}{S_x}$

Question 2.2. Calculate the slope of the regression line in original units, and assign it to slope.

(If the "stretching" explanation is unintuitive, consult section 15.2 in the textbook.)

In [ ]:

slope = ...
slope

We know that the regression line passes through the point (duration_mean, wait_mean). You might recall from high-school algebra that the equation for the line is therefore:

$$\text{waiting time} - \verb|wait_mean| = \texttt{slope} \times (\text{eruption duration} - \verb|duration_mean|)$$

The rearranged equation becomes:

$$\text{waiting time} = \texttt{slope} \times \text{eruption duration} + (- \texttt{slope} \times \verb|duration_mean| + \verb|wait_mean|)$$

Question 2.2. Calculate the intercept in original units and assign it to intercept.

In [ ]:

intercept = ...
intercept

3. Investigating the regression line¶

The slope and intercept tell you exactly what the regression line looks like. To predict the waiting time for an eruption, multiply the eruption's duration by slope and then add intercept.

Question 3.1. Compute the predicted waiting time for an eruption that lasts 2 minutes, and for an eruption that lasts 5 minutes.

In [ ]:

two_minute_predicted_waiting_time = ...
five_minute_predicted_waiting_time = ...

# Here is a helper function to print out your predictions.
# Don't modify the code below.
def print_prediction(duration, predicted_waiting_time):
    print("After an eruption lasting", duration,
          "minutes, we predict you'll wait", predicted_waiting_time,
          "minutes until the next eruption.")

print_prediction(2, two_minute_predicted_waiting_time)
print_prediction(5, five_minute_predicted_waiting_time)

The next cell plots the line that goes between those two points, which is (a segment of) the regression line.

In [ ]:

faithful.scatter( "duration", "wait") 
plots.plot([2, 5], [two_minute_predicted_waiting_time, five_minute_predicted_waiting_time]);

Question 3.2. Make predictions for the waiting time after each eruption in the faithful table. (Of course, we know exactly what the waiting times were! We are doing this so we can see how accurate our predictions are.) Put these numbers into a column in a new table called faithful_predictions. Its first row should look like this:

duration	wait	predicted wait
3.6	79	72.1011

Hint: Your answer can be just one line. There is no need for a for loop; use array arithmetic instead.

In [ ]:

faithful_predictions = 

faithful_predictions

Question 3.3. How close were we? Compute the residual for each eruption in the dataset. The residual is the actual waiting time minus the predicted waiting time. Add the residuals to faithful_predictions as a new column called residual and name the resulting table faithful_residuals.

Hint: Again, your code will be much simpler if you don't use a for loop.

In [ ]:

faithful_residuals_array = ...
faithful_residuals = ...
faithful_residuals

In [ ]:

Here is a plot of the residuals you computed. Each point corresponds to one eruption. It shows how much our prediction over- or under-estimated the waiting time.

In [ ]:

faithful_residuals.scatter("duration", "residual", color="r")

There isn't really a pattern in the residuals, which confirms that it was reasonable to try linear regression. It's true that there are two separate clouds; the eruption durations seemed to fall into two distinct clusters. But that's just a pattern in the eruption durations, not a pattern in the relationship between eruption durations and waiting times.

Question 3.4 Use the interact plot below to set the slope and intercept to the values that results in the least possible value for the RMSE. In the cell below that, record this least possible RMSE value.

In [ ]:

yy = faithful.column('wait')
xx = faithful.column('duration')

import ipywidgets as widgets

from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets

import matplotlib.pyplot as plt
%matplotlib inline

import numpy as np

def plot_func(slope, intercept):
    x = np.linspace(1, 6)
    y = x*slope + intercept
    plt.scatter(xx,yy, s=50, color="green")
    plt.ylim(40,100)
    plt.xlim(1,6) 
    rmse = np.round((sum((yy - slope*xx - intercept)**2)/(len(xx)-2))**0.5,2)
    plt.plot(x, y)
    plt.text(1.5,90, f"RMSE = {rmse}", color="red", size='large')
    plt.xlabel("Duration")
    plt.ylabel("Wait")

interact(plot_func, slope = widgets.FloatSlider(value=10.73, min=5, max=15, step=0.1), 

         intercept=widgets.FloatSlider(value=33.47, min=20, max=40, step=0.1));

In [ ]:

least_rmse_original_units = ...
least_rmse_original_units

4. How accurate are different predictions?¶

Earlier, you should have found that the correlation is fairly close to 1, so the line fits fairly well on the training data. That means the residuals are overall small (close to 0) in comparison to the waiting times.

We can see that visually by plotting the waiting times and residuals together:

In [ ]:

faithful_predictions.scatter("duration", "wait", label="actual waiting time", color="blue")
plots.scatter(faithful_predictions.column("duration"), faithful_residuals.column("residual"), label="residual", color="r")
plots.plot([2, 5], [two_minute_predicted_waiting_time, five_minute_predicted_waiting_time], label="regression line")
plots.legend(bbox_to_anchor=(1.7,.8));

However, unless you have a strong reason to believe that the linear regression model is true, you should be wary of applying your prediction model to data that are very different from the training data.

Question 4.1. In faithful, no eruption lasted exactly 0, 2.5, or 60 minutes. Using this line, what is the predicted waiting time for an eruption that lasts 0 minutes? 2.5 minutes? An hour?

In [ ]:

zero_minute_predicted_waiting_time = ...
two_point_five_minute_predicted_waiting_time = ...
hour_predicted_waiting_time = ...

print_prediction(0, zero_minute_predicted_waiting_time)
print_prediction(2.5, two_point_five_minute_predicted_waiting_time)
print_prediction(60, hour_predicted_waiting_time)

In [ ]:

Question 2. For each prediction, state whether you think it's reliable and explain your reasoning.

Write your answer here, replacing this text.

5. Divide and Conquer¶

It appears from the scatter diagram that there are two clusters of points: one for durations around 2 and another for durations between 3.5 and 5. A vertical line at 3 divides the two clusters.

In [ ]:

faithful.scatter("duration", "wait", label="actual waiting time", color="blue")
plots.plot([3, 3], [40, 100]);

Question 1. Separately compute the regression coefficients r for all the points with a duration below 3 and then for all the points with a duration above 3. To do so, first create two different tables of points, below_3 and above_3, then compute the correlations.

In [ ]:

below_3 = ...
above_3 = ...
below_3_r = ...
above_3_r = ...

print("For points below 3, r is", below_3_r, "; for points above 3, r is", above_3_r)

Question 5.2. Use the linregress function from scipy.stats to compute the slope and intercept for both the above_3 and below_3 datasets.

When you're done, record the slopes and intercepts in the cell below as below_3_slope, below_3_intercept, above_3_slope and above_3_intercept.

In [ ]:

...

In [ ]:

...

In [ ]:

below_3_slope = ...
below_3_intercept = ...
above_3_slope = ...
above_3_intercept = ...

The plot below shows two different regression lines, one for each cluster!

In [ ]:

faithful.scatter(0, 1)
plots.plot([1, 3], [below_3_slope*1+below_3_intercept, below_3_slope*3+below_3_intercept])
plots.plot([3, 6], [above_3_slope*3+above_3_intercept, above_3_slope*6+above_3_intercept]);

Question 3. In the cell below, we provide a function called predict_wait that takes a duration and returns the predicted wait time using the appropriate regression line, depending on whether the duration is below 3 or greater than (or equal to) 3. Use the .apply method and this new function to add the predicted wait times to the faithful table.

In [ ]:

def predict_wait(duration):
    if duration <=3:
        return below_3_slope*duration + below_3_intercept
    else:
        return above_3_slope*duration + above_3_intercept
        
predict_wait(2)   

The predicted wait times for each point appear below.

In [ ]:

faithful.with_column('predicted', faithful.apply(..., ...)).scatter(0)

Question 4. Do you think the predictions produced by predict_wait would be more or less accurate than the predictions from the regression line you created in section 2? How could you tell?

Write your answer here, replacing this text.

Question 5. Use the sliders to set the slopes and intercepts to the values from the linregress output.

In [ ]:

y1 = faithful.where('duration', are.below_or_equal_to(3)).column('wait')
x1 = faithful.where('duration', are.below_or_equal_to(3)).column('duration')

y2 = faithful.where('duration', are.above(3)).column('wait')
x2 = faithful.where('duration', are.above(3)).column('duration')

yy = faithful.column('wait')
xx = faithful.column('duration')

import ipywidgets as widgets

from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets

import matplotlib.pyplot as plt
%matplotlib inline

import numpy as np

def predict_wait_func(duration, lo_slope, hi_slope, lo_int, hi_int):
    if(duration <=3):
        return lo_slope * duration + lo_int
    else:
        return hi_slope * duration + hi_int

def plot_func(lo_slope, hi_slope, lo_int, hi_int):
    x3 = np.linspace(0, 3)
    y3 = make_array()
    x4 = np.linspace(3.01, 6)
    y4 = make_array()
    for i in np.arange(len(x3)):
        y3 =np.append(y3, predict_wait_func(x3.item(i),lo_slope, hi_slope, lo_int, hi_int))
    for i in np.arange(len(x4)):
        y4 =np.append(y4, predict_wait_func(x4.item(i),lo_slope, hi_slope, lo_int, hi_int))
    plt.scatter(x1,y1, s=50, color="green")
    plt.scatter(x2,y2, s=50, color="blue")
    plt.ylim(40, 100)
    plt.xlim(1, 5.5)
    new_residuals = make_array()
    for i in np.arange(len(xx)):
        new_residuals = np.append(new_residuals, yy.item(i) - predict_wait_func(xx.item(i),lo_slope, hi_slope, lo_int, hi_int) )
    rmse = np.round((sum(new_residuals**2)/(len(xx)-2))**0.5,2)
    plt.plot(x3, y3, color="green")
    plt.plot(x4, y4, color="blue")
    plt.text(1.5, 90, f"RMSE = {rmse}", color="red", size='large')
    plt.xlabel("Duration")
    plt.ylabel("Wait")

interact(plot_func, lo_slope = widgets.FloatSlider(value=6.35, min=5, max=15, step=0.1), 
        hi_slope = widgets.FloatSlider(value=5.44, min=0, max=15, step=0.1),
        lo_int=widgets.FloatSlider(value=41.6, min=20, max=50, step=0.1),
        hi_int=widgets.FloatSlider(value=56.6, min=20, max=70, step=0.1));

In [ ]:

split_rmse = ...

Which RMSE is lower? The one from the end of question 3.4 (least_rmse_original_units) or the one we just calculated (split_rmse)?

Summary¶

Run the cell below without changing anything in it. It should produce a nice summary of what we've learned during this lab.

In [ ]:

print("The mean duration of an eruption is ", np.round(duration_mean,1), "minutes with a standard deviation of", 
      np.round(duration_std,1),"." )
print("The mean wait time between eruptions is ", np.round(wait_mean,1), "minutes with a standard deviation of", 
      np.round(wait_std,1),"." )
print("These two variables have a correlation of", np.round(correlation,3), ".")

print("A single regression analysis has RMSE =", least_rmse_original_units,". \nThat model is predicted wait = ", 
      np.round(slope, 2), "duration + ", np.round(intercept,2), ".")
print("Splitting the data at durations of 3, we get an RMSE = ", split_rmse,". \nThis model uses two lines.  Before 3 we use ",
      below_3_slope,"duration + ", below_3_intercept, "and after 3 it is ", above_3_slope, "duration + ", above_3_intercept, "." )

In [ ]: