13: Linear Regression

Scatterplots in Standard Units

In Chapter 12 we learned about standard units or the z-score. A value in standard units tells us how many standard deviations it is above or below the mean of the dataset.

We wrote the standard_units function to convert every value in an array to standard units:

def standard_units(array):
    return (array - np.mean(array))/np.std(array)

Output?

We will now apply this function to both variables in a scatterplot. Although the numbers on the axes will change, the overall pattern of the scatterplot will remain the same.

We work with the weight and horsepower of the cars.

We convert the weight and horsepower data into standard units and create a new DataFrame to store them.

su = pd.DataFrame({'weight (z)': standard_units(cars['weight']),
                   'horsepower (z)': standard_units(cars['horsepower'])})
display(su)

Output?

Note that the data for both weight (z) and horsepower (z) are now centered around 0. The positive z values come from actual weight or horsepower values that are greater than the mean, and the negative z values come from actual weight or horsepower values that are less than the mean. The numbers have changed, but the relative position of the data values remains the same.

We now display the weight (z) vs. horsepower (z) in a scatterplot to show that the plot in standard units still has the same shape as the plot with actual units in pounds and horsepower shown above.

plt.figure(figsize=(4,4))
plt.scatter(su["weight (z)"], su["horsepower (z)"])
plt.xlabel("Weight (z)")
plt.ylabel("Horsepower (z)")
plt.grid()
plt.show()

Output?

Although the axes are now measured in standard units instead of in actual units of pounds and horsepower, the scatterplot has the same overall shape as before. This shows that converting data to standard units does not change the relationship between the two variables.

The Correlation Coefficient

A scatterplot visually shows us the relationship between two variables. To measure the strength of that relationship, we calculate the correlation coefficient.

The correlation coefficient is denoted by r and has the following characteristics:

• r is a value between 1 and -1.
• When r = 1, all the points lie on a straight line with a positive correlation.
• When r = 0, there is no linear correlation.
• When r = -1, all the points lie on a straight line with a negative correlation.

We observe the shape of the scatterplot and its relationship to the r value below.

# You don't need to study the following code.
# It's used to demo the different r values and its direct
# relationship with the clustering of the data.

plt.figure(figsize=(6,6))

# r = 0.95
plt.subplot(2, 2, 1)
x1 = np.linspace(0, 10, 100)
y1 = x1 * 0.95 + np.random.randn(100)
plt.scatter(x1, y1)
plt.title('subplot 1: r = 0.95')
plt.grid()

# r = 0.5
plt.subplot(2, 2, 2)
x2 = np.random.rand(100) * 10
y2 = x2 * 0.5 + np.random.randn(100)
plt.scatter(x2, y2)
plt.title('subplot 2: r = 0.5')
plt.grid()

# r = 0
plt.subplot(2, 2, 3)
x3 = np.random.rand(100) * 10
y3 = np.random.rand(100) * 10
plt.scatter(x3, y3)
plt.title('subplot 3: r = 0')
plt.grid()

# r = -0.8
plt.subplot(2, 2, 4)
x4 = np.linspace(0, 10, 100)
y4 = -x4 * 0.8 + np.random.randn(100)
plt.scatter(x4, y4)
plt.title('subplot 4: r = -0.8')
plt.grid()

plt.tight_layout()
plt.show()

Output?

The correlation coefficient r is the average of the products of the two variables, after both variables have been converted to standard units, or have been standardized.

Multiplying the two standardized variables captures two aspects of their relationship:

• Direction. If both variables are above their means or both are below their means, their product is positive. If one variable is above its mean and the other is below its mean, their product is negative. When we average these products, a positive average indicates a positive correlation, while a negative average indicates a negative correlation.
• Strength. A large product means that both variables are far from their means at the same time. A small product means that one or both variables are close to their means. Therefore, for the correlation coefficient r to have a large magnitude, many data points must move far from their means together in a consistent direction, which means that they have a strong linear correlation.

We define a function correlation that calculates the r value for 2 sequences of data:

def correlation(array1, array2):
    return np.mean(standard_units(array1) * standard_units(array2))

Output?

We test the correlation function, with the x1 and y1 data from subplot 1 above, and with the x3 and y3 data from subplot 3 above.

print("r for subplot 1:", round(correlation(x1, y1), 2))
print("r for subplot 3:", round(correlation(x3, y3), 2))

Output?

r for subplot 1: 0.94
r for subplot 3: -0.02

The calculated r results correspond with the scatterplots' shapes.

Generally, the closer the value of r is to 1 or −1, the stronger the linear relationship between the variables. The closer r is to 0, the weaker the linear relationship.

There is no universal cutoff for deciding what is considered "strong" or "weak." The scatterplot should always be examined together with the value of r.

Calculate the Regression Line

In the previous example, we estimated the regression line visually. However, different people might draw slightly different lines. Linear regression gives us a mathematical method for calculating the line that best fits the data.

In this case weight is the predictor variable because it is used to make the prediction, and mpg is the response variable because it is the value being predicted. The predictor variable is placed on the x-axis, and the response variable is placed on the y-axis.

• weight = predictor, on x-axis
• mpg = response, on y-axis

x = standard_units(cars["weight"])
y = standard_units(cars["mpg"])

Output?

r = correlation(x, y)
print("r:", r)

Output?

r: -0.8322442148315758

plt.figure(figsize=(4,4))

plt.scatter(x, y)
regression_line = r * x
plt.plot(x, regression_line, color='red')
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
plt.show()

Output?

Convert to Original Units

The regression line: ŷ = r × x predicts values only in standard units. To make predictions with the original data, we must convert the regression line back to the original units.

Recall the definition of the standard units: $$ \text{standard units} = \frac{ data - mean } { \text{standard deviation}} $$

Substitute the definition of standard_units for ŷ and x in the equation for the regression line: $$ \frac{\hat{y}- mean_{\hat{y}}}{SD_{y}} = r \times \frac{x- mean_{x}}{SD_{x}} $$

Solve for ŷ: $$ \hat{y} = r \times \frac{SD_{y}}{SD_{x}} \times (x-mean_{x}) + mean_{\hat{y}} $$

Distribute (x − mean_x) $$ \hat{y} = (r \times \frac{SD_{y}}{SD_{x}} \times x) - (r \times \frac{SD_{y}}{SD_{x}} \times mean_{x}) + mean_{\hat{y}} $$

Rearrange into the standard form of a line: y = mx + b $$ \hat{y} = (r \times \frac{SD_{y}}{SD_{x}}) \times x + (mean_{\hat{y}} - (r \times \frac{SD_{y}}{SD_{x}} \times mean_{x})) $$

In the equation of the line: y = mx + b
m is the slope and b is the intercept.

Therefore the slope of the regression line is: $$ slope = r \times \frac{SD_{y}}{SD_{x}} $$

The intercept is: $$ intercept = mean_{\hat{y}} - (r \times \frac{SD_{y}}{SD_{x}} \times mean_{x}) $$ Or: intercept = mean_ŷ − slope × mean_x

We can now write a function to find the slope and the intercept:

def slope(x, y):
    return correlation(x, y) * np.std(y) / np.std(x)

Output?

def intercept(x, y):
    return np.mean(y) - slope(x, y) * np.mean(x)

Output?

We now apply the functions to the cars weight and mpg data. First we review the dataset.

weight_mpg_slope = slope(cars["weight"], cars["mpg"])
weight_mpg_intercept = intercept(cars["weight"], cars["mpg"])

Output?

Displaying the scatterplot for weight vs mpg, along with the regression line.

plt.figure(figsize=(4,4))
plt.scatter(cars["weight"], cars["mpg"], alpha=0.9)
# regression line y = slope * x + intercept
plt.plot(cars["weight"], weight_mpg_slope * cars["weight"] + weight_mpg_intercept, color='red')
plt.xlabel("Weight")
plt.ylabel("MPG")
plt.grid()
plt.show()

Output?

With the linear regression line, we now have a consistent way to predict the MPG or ŷ, given a car's weight, which is x.

weight = 4000
predicted_mpg = weight_mpg_slope * weight + weight_mpg_intercept
print("Predicted MPG for a car weighing", weight, "lbs:", round(predicted_mpg, 1))

Output?

Predicted MPG for a car weighing 4000 lbs: 15.6

Meaning of the Slope

With the weight vs mpg data, the slope is:

print(round(weight_mpg_slope, 2))

Output?

-0.01

The slope of the regression line is −0.01. This means that, on average, for every additional pound of weight, the predicted MPG decreases by about 0.01 MPG. Equivalently, a car that is 100 pounds heavier is predicted to get about 1 MPG less.

In general, the slope of a regression line represents the average change in the predicted value of y for each one-unit increase in x. A positive slope means the predicted value of y increases as x increases, while a negative slope means the predicted value of y decreases as x increases.