Linear regression notation for models with multiple predictors¶

This page will introduce the notation for linear regression models with multiple predictor variables, expanding on the notation which we have already seen for linear regression models with a single predictor.

We will add the same caveat as before: this page will most likely be challenging if you are new to this, so please do not worry if you do not understand some concepts the first time you read them.

But do keep reading, and re-reading. You'll get there.

We will again use the Duncan (1961) occupational prestige dataset. The dataset combines information from the 1950 U.S. Census with data collected by the National Opinion Research Centre (NORC). The Census data contained information about different occupations, such as the percentage of people working in that occupation who earned over a certain amount per year. The NORC data was from a survey which asked participants to rate how prestigious they considered each occupation. Here are the descriptions of the variables in the dataset, which covers 45 occupations (adapted from here):

• name - the name of the occupation, from the 1950 US Census
• type- type of occupation, with the following categories prof, professional and managerial; wc, white-collar; bc, blue-collar. (E.g. how the occupation was classified in the 1950 US Census)
• income - percentage of census respondents within the occupation who earned 3,500 dollars or more per year (about 36,000 US dollars in 2017)
• education - percentage of census respondents within the occupation who were high school graduates
• prestige - percentage of respondents in the NORC survey who rated the occupation as “good” or better in prestige

See the dataset page for more detail.

Again, we will restrict our attention to the first 15 observations in the dataset (to make it easier to visualise/think about the concepts):

The observational units in this dataset are professions, and you'll remember from the previous page, that the linear regression model we fit was as follows:

prestige = $b * $ education + $ \text{c} + \vec{\varepsilon} $

The model equation reads as "we are modeling prestige as a linear function of education, plus error".

The error is the randomness that remains after the straight line relationship between education and prestige is accounted for e.g. the y-axis distance between each data point and the line of best fit.

Remember that prestige and education are vectors of 15 values.

Just for ease of typing, let's store the vectors again as separate python variables.

We learned the mathematical notation for the model we fit previously:

$ \vec{y} = b \vec{x} + \text{c} + \vec{\varepsilon} $

$\begin{bmatrix}{} \text{$y_{1}$} \\ \text{$y_{2}$} \\ \text{$y_{3}$} \\ \text{$y_{4}$} \\ \text{$y_{5}$} \\ \text{$y_{6}$} \\ \text{$y_{7}$} \\ \text{$y_{8}$} \\ \text{$y_{9}$} \\ \text{$y_{10}$} \\ \text{$y_{11}$} \\ \text{$y_{12}$} \\ \text{$y_{13}$} \\ \text{$y_{14}$} \\ \text{$y_{15}$} \\ \end{bmatrix} = b * \begin{bmatrix}{} \text{$x_{1}$} \\ \text{$x_{2}$} \\ \text{$x_{3}$} \\ \text{$x_{4}$} \\ \text{$x_{5}$} \\ \text{$x_{6}$} \\ \text{$x_{7}$} \\ \text{$x_{8}$} \\ \text{$x_{9}$} \\ \text{$x_{10}$} \\ \text{$x_{11}$} \\ \text{$x_{12}$} \\ \text{$x_{13}$} \\ \text{$x_{14}$} \\ \text{$x_{15}$} \\ \end{bmatrix} + c + \begin{bmatrix}{} \text{$\varepsilon_{1}$} \\ \text{$\varepsilon_{2}$} \\ \text{$\varepsilon_{3}$} \\ \text{$\varepsilon_{4}$} \\ \text{$\varepsilon_{5}$} \\ \text{$\varepsilon_{6}$} \\ \text{$\varepsilon_{7}$} \\ \text{$\varepsilon_{8}$} \\ \text{$\varepsilon_{9}$} \\ \text{$\varepsilon_{10}$} \\ \text{$\varepsilon_{11}$} \\ \text{$\varepsilon_{12}$} \\ \text{$\varepsilon_{13}$} \\ \text{$\varepsilon_{14}$} \\ \text{$\varepsilon_{15}$} \\ \end{bmatrix}$

Run the cell below which will remind us how the model looks with the data from our actual prestige and education vectors:

Just to remind ourselves about the parameter estimates $b$ and $c$ of the best-fitting line, let's fit the model again.

As we will be fitting several different models on this page, we will store the slope and intercept using python variable names which indicate which model they are from. In this case we will use b_education and c_education to indicate that they are the slope and intercept from a linear regression using education as the predictor variable:

Note: prestige will always be the outcome variable, for all models fit on this page.

Here are the data, and the fitted values from the model, shown graphically:

Linear regression using scipy.optimize¶

In order to make sense of the notation for having multiple predictors in the model, we need to pay attention to the cost function used in linear regression. The method of fitting regression models we used above, and on the previous page, involves a pre-written function (scipy.stats.linregress) which uses its own cost function.

Because the scipy library is open source, we can look at the cost function it uses. But, to deepen our understanding of the notation for linear regression with multiple predictors, let's re-visit another method of fitting linear regression models, which you've already seen on a previous page.

Here we write our own cost function. The cost function for a line is a function that returns a quality-of-fit value for a particular line. That is, it accepts parameters of the slope and intercept defining the line, and returns a value that will be low (low cost) for a well-fitting line, and high (high cost) for a poorly fitting line.

We have seen different cost functions for linear regression before on this page. You'll recall that we can use the sum of the squared errors or the root mean squared error - they will identify the same slope and intercept as giving the lowest cost. As we are working with a small number of observations, let's use the sum of the squared errors.

The cell below defines a function which takes a sequence (b_and_c) containing a slope and an intercept value. It then calculates fitted values for that slope and intercept pair, using the formula we saw on the linear regression notation page

Let's see what the sum of squared error is for an arbitrarily chosen slope and intercept value.

The graph below shows the line that results from this slope and intercept pair, along with the original data (blue dots) and the errors (black dashed lines).

Remember the mathematical form of a linear regression model with a single predictor:

$ \vec{y} = b \vec{x} + \text{c} + \vec{\varepsilon} $

$\begin{bmatrix}{} \text{$y_{1}$} \\ \text{$y_{2}$} \\ \text{$y_{3}$} \\ \text{$y_{4}$} \\ \text{$y_{5}$} \\ \text{$y_{6}$} \\ \text{$y_{7}$} \\ \text{$y_{8}$} \\ \text{$y_{9}$} \\ \text{$y_{10}$} \\ \text{$y_{11}$} \\ \text{$y_{12}$} \\ \text{$y_{13}$} \\ \text{$y_{14}$} \\ \text{$y_{15}$} \\ \end{bmatrix} = b * \begin{bmatrix}{} \text{$x_{1}$} \\ \text{$x_{2}$} \\ \text{$x_{3}$} \\ \text{$x_{4}$} \\ \text{$x_{5}$} \\ \text{$x_{6}$} \\ \text{$x_{7}$} \\ \text{$x_{8}$} \\ \text{$x_{9}$} \\ \text{$x_{10}$} \\ \text{$x_{11}$} \\ \text{$x_{12}$} \\ \text{$x_{13}$} \\ \text{$x_{14}$} \\ \text{$x_{15}$} \\ \end{bmatrix} + c + \begin{bmatrix}{} \text{$\varepsilon_{1}$} \\ \text{$\varepsilon_{2}$} \\ \text{$\varepsilon_{3}$} \\ \text{$\varepsilon_{4}$} \\ \text{$\varepsilon_{5}$} \\ \text{$\varepsilon_{6}$} \\ \text{$\varepsilon_{7}$} \\ \text{$\varepsilon_{8}$} \\ \text{$\varepsilon_{9}$} \\ \text{$\varepsilon_{10}$} \\ \text{$\varepsilon_{11}$} \\ \text{$\varepsilon_{12}$} \\ \text{$\varepsilon_{13}$} \\ \text{$\varepsilon_{14}$} \\ \text{$\varepsilon_{15}$} \\ \end{bmatrix}$

Let's make sure we can relate all aspects of the notation to the graph below (the error vector for this slope/intercept pair, as well as the sum of the squared error is shown below the graph):

Note: if you are using this notebook interactively, try and use a few different slope and intercept pairs by changing the value of b and c in the input to the make_scatter function in the cell below:

Let's compare this to the graph, the error vector and the sum of the squared error value from our arbitrarily guessed pairing to those from the parameter estimates for the best-fitting line:

We can see that when we compare the badly fitting line with the best-fitting line the value of the cost function (the sum of squared error, calculated via our ss_any_line function) is smaller for the best fitting line.

Just now, and on the regression notation page we used the scipy.stats.linregress function to quickly fit our regression. This essentially fits the model using a cost function that is internal to scipy.stats.linregress — we do not get to control the cost function.

Let's refresh our memories on how to find the parameter estimates using a cost function which we have written. We do this by passing our cost function (ss_any_line) to the minimize function from the scipy.optimize. scipy.optimize uses an efficient method to try a variety of slope/intercept pairings until it finds the pair which give the minimum value of the cost function:

Let's pull out our parameter estimates, from the output of minimize, so that we can compare them to the values we got above from scipy.stats.linregress:

Both of these approaches (using scipy.stats.linregress or using minimize with our cost function ss_any_line) return the same (near as dammit) parameter estimates. This is because linregress uses the mathematical shortcut to the calculation of the least-squares slope and intercept. That is, the mathematics behind linregress is using the sum-of-squares (or mean square or root mean square) cost function.

The benefit of writing our own cost function is firstly that it gives us a better understanding of how we are getting our parameter estimates, and secondly that we can customize the cost function. Minimizing cost functions is a key concept in machine learning, so you can think of this approach as a foundational machine learning concept, applying to newer machine learning methods, as well as to classical statistical methods like linear regression.

Let's just recap how the process of fitting a linear regression model (with one predictor variable) works, with reference to the mathematical notation:

$\begin{bmatrix}{} \text{$y_{1}$} \\ \text{$y_{2}$} \\ \text{$y_{3}$} \\ \text{$y_{4}$} \\ \text{$y_{5}$} \\ \text{$y_{6}$} \\ \text{$y_{7}$} \\ \text{$y_{8}$} \\ \text{$y_{9}$} \\ \text{$y_{10}$} \\ \text{$y_{11}$} \\ \text{$y_{12}$} \\ \text{$y_{13}$} \\ \text{$y_{14}$} \\ \text{$y_{15}$} \\ \end{bmatrix} = b * \begin{bmatrix}{} \text{$x_{1}$} \\ \text{$x_{2}$} \\ \text{$x_{3}$} \\ \text{$x_{4}$} \\ \text{$x_{5}$} \\ \text{$x_{6}$} \\ \text{$x_{7}$} \\ \text{$x_{8}$} \\ \text{$x_{9}$} \\ \text{$x_{10}$} \\ \text{$x_{11}$} \\ \text{$x_{12}$} \\ \text{$x_{13}$} \\ \text{$x_{14}$} \\ \text{$x_{15}$} \\ \end{bmatrix} + c + \begin{bmatrix}{} \text{$\varepsilon_{1}$} \\ \text{$\varepsilon_{2}$} \\ \text{$\varepsilon_{3}$} \\ \text{$\varepsilon_{4}$} \\ \text{$\varepsilon_{5}$} \\ \text{$\varepsilon_{6}$} \\ \text{$\varepsilon_{7}$} \\ \text{$\varepsilon_{8}$} \\ \text{$\varepsilon_{9}$} \\ \text{$\varepsilon_{10}$} \\ \text{$\varepsilon_{11}$} \\ \text{$\varepsilon_{12}$} \\ \text{$\varepsilon_{13}$} \\ \text{$\varepsilon_{14}$} \\ \text{$\varepsilon_{15}$} \\ \end{bmatrix}$

Let's think about what the minimize function is doing here, with reference to this notation:

• Our $x$ and $y$ vectors are fixed, they are given by our dataset (in our case education is our $x$ vector and prestige is our $y$ vector).
• We write a cost function (ss_any_line) which takes a value of $b$ and $c$ and then returns the "cost" for those parameters, given our our $x$ and $y$ vectors.
  • It does this by generating fitted values using $b$ and $c$,
  subtracting these fitted values from the actual prestige scores to generate the error vector, and then squaring the values in the error vector and adding them up (to calculate the sum of squared error for this specific $b$ and $c$ pairing).
• We pass our cost function to minimize and it tries multiple $b$ and $c$ pairs until it finds the pair that gives the lowest value of the cost function, out of all possible $b$ and $c$ pairs.

Exploring other models¶

So far we've refreshed our knowledge of how we can use minimize on the sum of squared errors cost function (ss_any_line) to get our parameter estimates ($b$ and $c$) for linear regression.

We've done this using education as our predictor, and prestige as our outcome.

However, there are other potential predictor variables in the Duncan dataset, which might be relevant to predicting occupational prestige.

Let's investigate the relationship between income and prestige. Reasonably, we might predict a priori (before analysing the data) that higher paid jobs might be more prestigious.

Let's store our income scores as a separate python variable.

Before we continue: we would like you to do an exercise, for some practice.

We would like you to fit the following model, using the scipy.optimize

method first, then checking your results with the scipy.stats.linregress method:

prestige = $b * $ income + $ \text{c} + \vec{\varepsilon} $

The cell below shows the model, with the actual data in the vectors:

Let's use minimize and the ss_any_line cost function to find values of the slope $b$ and the intercept $c$ which give us the best-fitting line, when we use linear regression to predict prestige from income:

Let's compare these, once more, to the parameter estimates we get from scipy.stats.linregress, when we use linear regression to predict prestige from income

The graph below shows the best fitting line for the model prestige = $b *$ income + $c$ + $\vec{\varepsilon}$

Expanding the model - multiple predictors¶

So far we have fit two separate linear regression models.

Each separate model is shown below, in vector form, with the actual data and error values in the relevant vectors, as well as the parameter estimates obtained from each linear regression:

The fitted values from each separate linear regression are shown in the output of the cell below. You'll notice that the two models give different fitted values. You can think of each of these regressions as modeling the linear relationship between each predictor and the outcome variable if we ignore the other predictor.

Generally, when exploring datasets and doing data science, we will be interested in more than one predictor. We may want to take into account another predictor when we get our parameter estimates.

In most cases we will want a model where the fitted values take into account information from both predictor variables.

You've already seen that it is possible to include multiple predictors in a linear regression model on the previous multiple regression page.

In fact, the original model Duncan fit included both education and income as predictor variables - (and as mentioned previously, the parameter estimates were calculated with a mechanical calculator, a labor which we will not subject you to!).

But how might we express the inclusion of multiple predictors in the mathematical notation that we've already seen for linear regression models?

You can think of it this way, we want to include a second predictor variable in our model. One way of saying this is that we are including a second predictor vector in our model. So we need to include a second predictor vector in the mathematical notation.

Remember the meaning of the slope? It relates the predictor variable to the outcome variable by describing how the outcome variable changes as the predictor variable scores change (where "change" refers to a comparison of the predictor variable and outcome variable scores of different observational units).

Therefore, if we include a second predictor vector, we will also want a second slope in the model. So if our original single-predictor model was:

prestige = $b * $ education + $ \text{c} + \vec{\varepsilon} $

If we want to include income as a second predictor in the model, we could write this as:

prestige = $b_1 * $ education + $b_2 * $ income + $ \text{c} + \vec{\varepsilon} $

You'll notice there are now two slopes $b_1$ and $b_2$.

• $b_1$ relates changes in education to changes in prestige.
• $b_2$ relates changes in income to changes in prestige.

We keep the form of the relationship between each slope and each predictor the same (e.g. the predictor vector gets multiplied by the slope for that predictor) in order that the model is still a linear regression - if we changed the form of the slope/predictor relationship, say, by raising the slope to the power of the predictor values, then the model would be nonlinear.

We have only one intercept term for reasons we will explain later (for now you can just think of the intercept as allowing us to be more flexible in the model we fit, allowing us to fit it better to a wider variety of datasets).

Remember that our predictor vector education can also be referred to by its mathematical name $\vec{x}$? We might want to use different symbols to refer to our second slope and second predictor, so we could write something like:

$ \vec{y} = b_1 \vec{x} + b_2 \vec{w} + c + \vec{\varepsilon} $

Where in our case:

• $\vec{y}$ would be our prestige vector (e.g. our outcome variable).
• $\vec{x}$ would be our education vector (e.g. our first predictor variable - the slope for this predictor is $b_1$).
• $\vec{w}$ would be our income vector (e.g. our second predictor variable — the slope for this predictor is $b_2$).

However, if we include extra predictors (three, or four, or five...) then we will have to use new symbols for each one, which will get cumbersome.

Instead we add a subscript to the mathematical notation for each predictor variable, like this:

$ \vec{x_1} = \vec{x} = $ education

$ \vec{x_2} = \vec{w} = $ income

We can then refer to each predictor vector as "the $x$ one vector" for $\vec{x_1}$ and the "the $x$ two vector" for $\vec{x_2}$.

If we had additional predictor variables (e.g. more than two) we can extend this notation to as many variables as we like, $\vec{x_3}$ for the third predictor variable, $\vec{x_4}$ for the fourth predictor variable and so on...

We can refer to each of the values within each vector using this notation:

education = $ \vec{x_1} = \begin{bmatrix}{} x_{1_1}, \\ x_{1_2} \\ ... \\ x_{1_n}\end{bmatrix} $

income = $\vec{x_2} = \begin{bmatrix}{} x_{2_1}, \\ x_{2_2} \\ ... \\ x_{2_n}\end{bmatrix}$

So each $x$ value now has a number indicating which predictor vector it is part of (e.g. $x_{1}$ or $x_{2}$), and then has a second subscript indicating it's index location within each vector. So:

• $x_{1_1}$ is the first value from $\vec{x_1}$
• $x_{1_2}$ is the second value from $\vec{x_1}$
• $x_{1_3}$ is the third value from $\vec{x_1}$
• ...
• $x_{1_{13}}$ is the thirteenth value from $\vec{x_1}$
• $x_{1_{14}}$ is the fourteenth value from $\vec{x_1}$
• $x_{1_{15}}$ is the fifteenth value from $\vec{x_1}$

And:

• $x_{2_1}$ is the first value from $\vec{x_2}$
• $x_{2_2}$ is the second value from $\vec{x_2}$
• $x_{2_3}$ is the third value from $\vec{x_2}$
• ...
• $x_{2_{13}}$ is the thirteenth value from $\vec{x_2}$
• $x_{2_{14}}$ is the fourteenth value from $\vec{x_2}$
• $x_{2_{15}}$ is the fifteenth value from $\vec{x_2}$

Remember the advantage of mathematical notation like this is that it's general. It can be used to refer to datasets with more than 15 observations. You'll remember that n represents the total number of values within each vector, in our case $n = 15$, as there are 15 observations. So we can also write the final value in each vector as:

• $\Large x_{1_n}$ is the final value from $ \Large \vec{x_1}$, which in this case has 15 values, so can also be written as $ \Large x_{1_{15}}$

Now we can refer to our education as our $x_1$ vector and income as our $x_2$ vector, the whole model, with both predictors, in mathematical notations is:

$ \vec{y} = b_1 \vec{x_1} + b_2 \vec{x_2} + \text{c} + \vec{\varepsilon} $

We can write this using the Python names for the vectors we are using:

prestige = $b_1 * $ education + $b_2 * $ income + $ \text{c} + \vec{\varepsilon} $

And here is the mathematical form, showing the mathematical notation for the values in each vector:

$\begin{bmatrix}{} \text{$y_{1}$} \\ \text{$y_{2}$} \\ \text{$y_{3}$} \\ \text{$y_{4}$} \\ \text{$y_{5}$} \\ \text{$y_{6}$} \\ \text{$y_{7}$} \\ \text{$y_{8}$} \\ \text{$y_{9}$} \\ \text{$y_{10}$} \\ \text{$y_{11}$} \\ \text{$y_{12}$} \\ \text{$y_{13}$} \\ \text{$y_{14}$} \\ \text{$y_{15}$} \\ \end{bmatrix} = b_1 * \begin{bmatrix}{} \text{$x_{1_1}$} \\ \text{$x_{1_2}$} \\ \text{$x_{1_3}$} \\ \text{$x_{1_4}$} \\ \text{$x_{1_5}$} \\ \text{$x_{1_6}$} \\ \text{$x_{1_7}$} \\ \text{$x_{1_8}$} \\ \text{$x_{1_9}$} \\ \text{$x_{1_{10}}$} \\ \text{$x_{1_{11}}$} \\ \text{$x_{1_{12}}$} \\ \text{$x_{1_{13}}$} \\ \text{$x_{1_{14}}$} \\ \text{$x_{1_{15}}$} \\ \end{bmatrix} + b_2 * \begin{bmatrix}{} \text{$x_{2_{1}}$} \\ \text{$x_{2_{2}}$} \\ \text{$x_{2_{3}}$} \\ \text{$x_{2_{4}}$} \\ \text{$x_{2_{5}}$} \\ \text{$x_{2_{6}}$} \\ \text{$x_{2_{7}}$} \\ \text{$x_{2_{8}}$} \\ \text{$x_{2_{9}}$} \\ \text{$x_{2_{10}}$} \\ \text{$x_{2_{11}}$} \\ \text{$x_{2_{12}}$} \\ \text{$x_{2_{13}}$} \\ \text{$x_{2_{14}}$} \\ \text{$x_{2_{15}}$} \\ \end{bmatrix} +c + \begin{bmatrix}{} \text{$\varepsilon_{1}$} \\ \text{$\varepsilon_{2}$} \\ \text{$\varepsilon_{3}$} \\ \text{$\varepsilon_{4}$} \\ \text{$\varepsilon_{5}$} \\ \text{$\varepsilon_{6}$} \\ \text{$\varepsilon_{7}$} \\ \text{$\varepsilon_{8}$} \\ \text{$\varepsilon_{9}$} \\ \text{$\varepsilon_{10}$} \\ \text{$\varepsilon_{11}$} \\ \text{$\varepsilon_{12}$} \\ \text{$\varepsilon_{13}$} \\ \text{$\varepsilon_{14}$} \\ \text{$\varepsilon_{15}$} \\ \end{bmatrix}$

The code cell below "fills in" the notation with the actual values from our vectors:

Let's look at our vectors individually, to compare them to the notation:

We will now pause for an exercise...¶

Multiple linear regression using arbitrary guesses¶

Now we are armed with the notation for linear regression models which include multiples predictors, let's look at the mechanics of fitting such models.

Essentially, we now need to modify our cost function to accept two slopes and an extra predictor vector, as per the new model specification. Let's quickly recap the new model, with multiple predictors:

prestige = $b_1 * $ education + $b_2 * $ income + $ \text{c} + \vec{\varepsilon} $

Which we can write in mathematical notation as:

$ \vec{y} = b_1 \vec{x_1} + b_2 \vec{x_2} + \text{c} + \vec{\varepsilon} $

• $\vec{x_1}$ is the education vector.
• $\vec{x_2}$ is the income vector.
• $b_1$ is the slope for education.
• $b_2$ is the slope for income.
• $c$ is the intercept (again, including this lets us fit our model more flexibly, just as when we use a single predictor).
• $\vec{\varepsilon}$ is the error vector.

You'll notice that there is still only one error vector despite there now being two slopes. This is a key concept to appreciate, when we include extra predictors in the model.

Recall that the procedure for fitting a linear regression with a single predictor variable is:

• we write a cost function which takes a value of $b$ and $c$, and returns a "cost" - e.g. a single value summary of how well the fitted values from that value of $b$ and $c$ fit the data we have.
• Specifically, in our case, the cost function generates fitted values ($\hat{y}$), using $b$ and $c$, and then subtracts these from our actual $y$ values, to give us the error vector $\vec{\varepsilon}$. The values $\vec{\varepsilon}$ and then squared and added together, to give us the sum of the squared error. Our cost function is shown below:

• We then pass our cost function to minimize and it tries multiple $b$ and $c$ pairs until it finds the pair that gives the lowest value of the cost function (e.g. the lowest sum of the squared error), out of all possible $b$ and $c$ pairs.

How do we include an extra predictor variable in this process? - we include it in the computation of the fitted values. When we have two predictors, we calculate our fitted values $\hat{y}$ as follows:

$\vec{\hat{y}} = b_1 \vec{x_1} + b_2 \vec{x_2} + \text{c} $

Or, showing the mathematical notation for the values within the vectors:

$\begin{bmatrix}{} \text{$\hat{y_{1}}$} \\ \text{$\hat{y_{2}}$} \\ \text{$\hat{y_{3}}$} \\ \text{$\hat{y_{4}}$} \\ \text{$\hat{y_{5}}$} \\ \text{$\hat{y_{6}}$} \\ \text{$\hat{y_{7}}$} \\ \text{$\hat{y_{8}}$} \\ \text{$\hat{y_{9}}$} \\ \text{$\hat{y_{10}}$} \\ \text{$\hat{y_{11}}$} \\ \text{$\hat{y_{12}}$} \\ \text{$\hat{y_{13}}$} \\ \text{$\hat{y_{14}}$} \\ \text{$\hat{y_{15}}$} \\ \end{bmatrix} = b_1 * \begin{bmatrix}{} \text{$x_{1_1}$} \\ \text{$x_{1_2}$} \\ \text{$x_{1_3}$} \\ \text{$x_{1_4}$} \\ \text{$x_{1_5}$} \\ \text{$x_{1_6}$} \\ \text{$x_{1_7}$} \\ \text{$x_{1_8}$} \\ \text{$x_{1_9}$} \\ \text{$x_{1_{10}}$} \\ \text{$x_{1_{11}}$} \\ \text{$x_{1_{12}}$} \\ \text{$x_{1_{13}}$} \\ \text{$x_{1_{14}}$} \\ \text{$x_{1_{15}}$} \\ \end{bmatrix} + b_2 * \begin{bmatrix}{} \text{$x_{2_{1}}$} \\ \text{$x_{2_{2}}$} \\ \text{$x_{2_{3}}$} \\ \text{$x_{2_{4}}$} \\ \text{$x_{2_{5}}$} \\ \text{$x_{2_{6}}$} \\ \text{$x_{2_{7}}$} \\ \text{$x_{2_{8}}$} \\ \text{$x_{2_{9}}$} \\ \text{$x_{2_{10}}$} \\ \text{$x_{2_{11}}$} \\ \text{$x_{2_{12}}$} \\ \text{$x_{2_{13}}$} \\ \text{$x_{2_{14}}$} \\ \text{$x_{2_{15}}$} \\ \end{bmatrix} +c $

You'll note that, despite there now being two predictors and two slopes, we still only get one set of fitted values. That is what makes this a multiple regression model rather than two separate linear regressions. We performed two separate linear regressions near the start of this page, and they produced two separate vectors of fitted values, and two separate error vectors.

Because the multiple predictor model takes into account both predictors when finding its parameter estimates, for any set of parameter estimates ($b_1, b_1, c $) we only get one vector of fitted values, and therefore we also have only one vector of error values (in the error vector $\vec{\varepsilon}$).

We square the values in $\vec{\varepsilon}$ and add them up to get our sum of squared error, for a specific combination of parameter estimates ($b_1, b_1, c $).

We then use minimize to search for the combination of parameter estimates ($b_1, b_1, c $) which give the lowest sum of the squared error, out of all possible combinations.

This is exactly the same process as for a single predictor model - only now the fitted values are generated using information from two predictor variables.

For instance, let's say we take a completely arbitrary guess at the combination of parameter estimates, to see how the fitted values and error vector are generated, using two predictor variables:

With our guessed parameter estimates, the model with two predictors could be written as follows:

Remember the formula for the fitted values with a single predictor?:

$\vec{\hat{y}} = b \vec{x} + \text{c}$

As mentioned above, just as the formula for the model, this can be modified for two predictors:

$\vec{\hat{y}} = b_1 \vec{x_1} + b_2\vec{x_2} + \text{c}$

In our case this is:

$\vec{\hat{y}} = b_1 * $ education + $b_2 * $ income + $\text{c} $

Let's calculate the fitted values for the slopes and intercepts we just guessed:

The cell below prints out the process of computing the fitted values, using two predictors, in vector form.

It first shows the model, and then shows the values within the vectors in the model.

Remember that we get our error vector ($\vec{\varepsilon}$) by subtracting our fitted values ($\hat{y}$) from the values of our actual outcome variable:

$$ \vec{\varepsilon} = \begin{bmatrix} y_1 - \hat{y}_1 \\ y_2 - \hat{y}_2 \\ \vdots \\ y_n - \hat{y}_n \end{bmatrix} $$

This is performed in the cell below, using the fitted models based on our guessed parameter estimates. We reiterate again that these fitted values are now generated using two predictor vectors and two slopes.

The output of the cell below shows the vectors in the model, and now also shows the error vector, along with our guessed parameter estimates.

We can then square the values in this error vector, and add them up. This gives us the sum of the squared error for our particular combination of slopes and intercepts ($b_1, b_1, c $) which we arbitrarily guessed.

The code cell below lets us make different guesses about the parameter estimates (the slopes and intercepts $b_1, b_1, c $) and shows us the error vector and the sum of the squared error, for that particular combination.

If you are using this page interactively, try using the sliders below to change the parameter estimates, and see how this changes the error vector, and the resultant sum of the squared error:

• The first slider controls the $x_1$ (education) slope ($b_1$).
• The second slider controls the $x_2$ (income) slope ($b_1$).
• The third slider controls the intercept.

Note: it's best to change these parameters "slowly" e.g. take a deep note of how the vectors look with a certain parameter combination, and what the sum of the squared error is for that combination. Then change one parameter and see how it affects the error vector.

To recap:

• we are now using two predictor variables and two slopes to generate our fitted values.
• Everything else in the fitting process is the same as for a single predictor model.

The value of the sum of the squared error depends on our particular combination of parameter estimates ($b_1, b_2, c $). We can use minimize to try lots of combinations (many triplets) in order to identify the combination which gives the smallest sum of the squared error.

In order to do that, we need to modify our cost function so that it accepts two predictors, and generates fitted values which result from those:

• two predictors,
• two corresponding slopes and
• an intercept.

Let's test our function out, with the last set of guessed parameters we used:

Multiple linear regression using minimize¶

Once we have modified our cost function, we can pass it to minimize in exactly the same way as with a single predictor model.

Only now we need to supply it with an extra slope in the list of initial parameter estimates, and an extra predictor vector in the args = argument (e.g. we now use args=(education, income, prestige) rather than just args=(education, prestige).

We can view our parameter estimates just as before, only now you will see that there are three of them (because of the extra slope $b_2$).

We can save them as separate python variables, just as before:

Let's compare this to the sum of squared error from the last set of guessed parameter estimates that we used:

We can see that minimize found a much better fitting model!

So, just to recap:

We added an extra predictor into our model:

prestige = $b_1 * $ education + $b_2 * $ income + $ \text{c} + \vec{\varepsilon} $

The mathematical notation for this is:

$ \vec{y} = b_1 \vec{x_1} + b_2 \vec{x_2} + \text{c} + \vec{\varepsilon} $

$\begin{bmatrix}{} \text{$y_{1}$} \\ \text{$y_{2}$} \\ \text{$y_{3}$} \\ \text{$y_{4}$} \\ \text{$y_{5}$} \\ \text{$y_{6}$} \\ \text{$y_{7}$} \\ \text{$y_{8}$} \\ \text{$y_{9}$} \\ \text{$y_{10}$} \\ \text{$y_{11}$} \\ \text{$y_{12}$} \\ \text{$y_{13}$} \\ \text{$y_{14}$} \\ \text{$y_{15}$} \\ \end{bmatrix} = b_1 * \begin{bmatrix}{} \text{$x_{1_1}$} \\ \text{$x_{1_2}$} \\ \text{$x_{1_3}$} \\ \text{$x_{1_4}$} \\ \text{$x_{1_5}$} \\ \text{$x_{1_6}$} \\ \text{$x_{1_7}$} \\ \text{$x_{1_8}$} \\ \text{$x_{1_9}$} \\ \text{$x_{1_{10}}$} \\ \text{$x_{1_{11}}$} \\ \text{$x_{1_{12}}$} \\ \text{$x_{1_{13}}$} \\ \text{$x_{1_{14}}$} \\ \text{$x_{1_{15}}$} \\ \end{bmatrix} + b_2 * \begin{bmatrix}{} \text{$x_{2_{1}}$} \\ \text{$x_{2_{2}}$} \\ \text{$x_{2_{3}}$} \\ \text{$x_{2_{4}}$} \\ \text{$x_{2_{5}}$} \\ \text{$x_{2_{6}}$} \\ \text{$x_{2_{7}}$} \\ \text{$x_{2_{8}}$} \\ \text{$x_{2_{9}}$} \\ \text{$x_{2_{10}}$} \\ \text{$x_{2_{11}}$} \\ \text{$x_{2_{12}}$} \\ \text{$x_{2_{13}}$} \\ \text{$x_{2_{14}}$} \\ \text{$x_{2_{15}}$} \\ \end{bmatrix} +c + \begin{bmatrix}{} \text{$\varepsilon_{1}$} \\ \text{$\varepsilon_{2}$} \\ \text{$\varepsilon_{3}$} \\ \text{$\varepsilon_{4}$} \\ \text{$\varepsilon_{5}$} \\ \text{$\varepsilon_{6}$} \\ \text{$\varepsilon_{7}$} \\ \text{$\varepsilon_{8}$} \\ \text{$\varepsilon_{9}$} \\ \text{$\varepsilon_{10}$} \\ \text{$\varepsilon_{11}$} \\ \text{$\varepsilon_{12}$} \\ \text{$\varepsilon_{13}$} \\ \text{$\varepsilon_{14}$} \\ \text{$\varepsilon_{15}$} \\ \end{bmatrix}$

We modified our cost function so that it uses fitted values ($\hat{y}$) which are generated using two slopes and two predictor variables:

$\vec{\hat{y}} = b_1 * $ education + $b_2 * $ income + $\text{c} $

Which can be shown with the mathematical notation for the individual values within the vectors:

$\begin{bmatrix}{} \text{$\hat{y_{1}}$} \\ \text{$\hat{y_{2}}$} \\ \text{$\hat{y_{3}}$} \\ \text{$\hat{y_{4}}$} \\ \text{$\hat{y_{5}}$} \\ \text{$\hat{y_{6}}$} \\ \text{$\hat{y_{7}}$} \\ \text{$\hat{y_{8}}$} \\ \text{$\hat{y_{9}}$} \\ \text{$\hat{y_{10}}$} \\ \text{$\hat{y_{11}}$} \\ \text{$\hat{y_{12}}$} \\ \text{$\hat{y_{13}}$} \\ \text{$\hat{y_{14}}$} \\ \text{$\hat{y_{15}}$} \\ \end{bmatrix} = b_1 * \begin{bmatrix}{} \text{$x_{1_1}$} \\ \text{$x_{1_2}$} \\ \text{$x_{1_3}$} \\ \text{$x_{1_4}$} \\ \text{$x_{1_5}$} \\ \text{$x_{1_6}$} \\ \text{$x_{1_7}$} \\ \text{$x_{1_8}$} \\ \text{$x_{1_9}$} \\ \text{$x_{1_{10}}$} \\ \text{$x_{1_{11}}$} \\ \text{$x_{1_{12}}$} \\ \text{$x_{1_{13}}$} \\ \text{$x_{1_{14}}$} \\ \text{$x_{1_{15}}$} \\ \end{bmatrix} + b_2 * \begin{bmatrix}{} \text{$x_{2_{1}}$} \\ \text{$x_{2_{2}}$} \\ \text{$x_{2_{3}}$} \\ \text{$x_{2_{4}}$} \\ \text{$x_{2_{5}}$} \\ \text{$x_{2_{6}}$} \\ \text{$x_{2_{7}}$} \\ \text{$x_{2_{8}}$} \\ \text{$x_{2_{9}}$} \\ \text{$x_{2_{10}}$} \\ \text{$x_{2_{11}}$} \\ \text{$x_{2_{12}}$} \\ \text{$x_{2_{13}}$} \\ \text{$x_{2_{14}}$} \\ \text{$x_{2_{15}}$} \\ \end{bmatrix} +c $

We then used minimize to find the combination of slopes and intercepts ($b_1, b_2, c$) which gives the lowest sum of the squared error, out of any combination.

The model is shown below, with the best-fitting parameter estimates ($b_1, b_2, c$), as well as the prestige, education and income vectors, and the error vector $\vec{\varepsilon}$.

Lines in 3D? WTF?¶

Note: (WTF? = What To Fit?)

When we fit a linear regression with one predictor variable, we used minimize to find the line of best fit.

Now we have two predictors, are we still fitting a line?

When using one predictor, we used plt.scatter to create a scatterplot to show the data.

Now we have two predictors, we will need a 3D scatterplot (do not worry about the specifics of how to generate this plot, it is just to visualise what is happening when we fit our model):

Each datapoint (in blue) on the scatterplot is one observational unit (in this case an occupation). The location on the education axis of each point is given by that occupation's education score. The location on the income axis of each point is given by that occupation's income score. The location on the prestige axis is given by that occupation's prestige score.

Let's look at the formula for the fitted values, with two predictors:

$\vec{\hat{y}} = b_1 * $ education + $b_2 * $ income + $\text{c} $

This equation let's us generate a fitted value for each observational unit. The fitted value for a specific observational unit is generated by the education score and income score of that observational unit, as well as the parameter estimates $b_1, b_2$ and $c$.

But remember that the foundational purposes of regression models are description and prediction? If we were to use our model to make predictions about new, unseen data from the same underlying population, then we might get an observational unit with a new pairing of education and income scores (for instance, scoring 100 on each).

As you look at the 3D scatterplot, imagine it as a physical object. As if you're literally viewing the points floating in sapce. A specific pairing of an education score and income score can occur at any point on the "floor" of the scatterplot e.g at any point in the square "floor" of the graph.

Fortunately, once we have our parameter estimates $b_1, b_2$ and $c$ we can generate fitted values for any education and income score pairing.

This means that we are fitting a plane of best fit rather than a line of best fit.

This is easier to appreciate graphically.

The output of the cell below shows the plane of best fit, obtained using minimize:

Let's try some arbitrary slope/intercept combinations, so we can see what it looks like when the plane doesn't fit the data well.

If you are using this page interactively, use the sliders below to try some combinations for the slopes and intercept, so you can see how they affect the plane and how it fits the data.

The vector notation for the model, along with your custom parameter estimates and the resultant sum of the squared errors, are shown beneath the scatter plot.

• The first slider controls the size of the education slope $b_1$
• The second slider controls the size of the income slope $b_2$
• The third slider controls the size of the intercept

Let's recap how we fit this plane to the data, in a model with two predictors (make sure you can relate each step of this process to the graph and the mathematical notation shown above this cell):

• our $x_1$, $x_2$ and $y$ vectors are fixed, they are given by our dataset (in our case education is our $x_1$ vector, income is our $x_1$ vector, and prestige is our $y$ vector)
• we write a cost function (ss_any_line) which takes a value of $b_1$, $b_2$ and $c$ and then returns the "cost" for those parameters, given our our $x_1$, $x_2$ and $y$ vectors.
  • Specifically, our cost function does this by generating fitted values using $b_1$, $b_2$ and $c$, subtracting these fitted values from the actual prestige scores to generate the error vector. The values in the error vector are then squared and added up (to calculate the sum of squared error for this specific $b_1$, $b_2$ and $c$ combination)
• we pass our cost function to minimize and it tries multiple $b_1$, $b_2$ and $c$ combinations until it finds the combination that gives the lowest value of the cost function, out of all possible $b_1$, $b_2$ and $c$ combinations

How does this compare to linear regression with one predictor?¶

As mentioned above, using multiple regression takes into account multiple predictors when finding it's parameter estimates (e.g. the $b_1$, $b_2$ and $c$ combination which best fit the data).

In contrast, we mentioned earlier that when we conduct separate single predictor linear regressions, each regression ignores other predictors.

The plot below shows the 3D regression plane if we use the slope from a single predictor regression, and ignore the second predictor.

It does this by using the slope b_education and intercept c_education which we got from the very first regression model we fit on this page - e.g. a single predictor regressing, predicting prestige from education.

The slope for the second predictor (income) is set to 0, so that it is ignored - only information about the education scores is used to fit this plane:

Again, as you look at the 3D scatterplot, imagine it as a physical object.

Imagine you are standing in front of the education axis, where the axis label is, looking towards the "back wall" of the graph.

The edge of the plane would appear to you as a line, and assuming you were the right height, you would see something like this:

Where the red line would be the edge of the regression plane.

This is the line of best fit from the single predictor linear regression - it ignores the other potential predictor variable income. It does not factor it into it's parameter estimates.

Likewise, the 3D scatterplot below is generated from the slope (b_income) and the intercept (c_income) which were obtained from the second single-predictor linear regression model that we fit on this page.

Now the slope for education is set to 0 (ignored) when generating the plane:

Now imagine you are standing in front of the income axis, where the axis label is, looking towards the "left wall" of the graph.

Again, the edge of the plane would appear to you as a line, and assuming you were the right height, you would see something like this:

When we fit a multiple regression model both of these slopes are adjusted in light of the other.

We get a regression plane that is fit to the datapoints using information from both predictor variables.

You can see that this changes the plane considerably, relatively to when we only consider one predictor variable. The graph below shows the best-fitting plane, fit by considering both predictor variables:

What is the meaning of the slopes?¶

Remember that in linear regression with a single predictor, the interpretation of the slope is:

"The slope is the expected change in the outcome variable if we compared two observational units which differed only by 1 point on the predictor variable".

Again, if you imagine the graph as a physical object/space, if you imagine walking in a straight line up the regression plane in one direction, say you were walking in the direction of the education axis, from 0 to 100, but climbing the plane as you went. The slope of the plane in the other direction (along the income axis) would remain constant under your feet.

Similarly, if you imagine walking in the direction of the income axis, from 0 to 100, but climbing the plane as you went. The slope of the plane in the other direction (along the education axis) would remain constant under your feet.

Thus, the interpretation of each slope (now we have two slopes) in a multi-predictor linear regression is:

"The slope of a predictor variable (e.g. $b_1$) is the expected change in the outcome variable if we compared two observational units which differed only by 1 point on that predictor variable, holding the other predictor constant".

Probably the best (and simplest) way of thinking about this is that the estimate of the slope of one predictor ($x_1$) - which describes the relationship of that predictor and the outcome variable ($y$) - has been adjusted in light of the relationship between the other predictor ($x_2$) and the outcome ($y$).

This gives us an estimate of the expected change in the outcome variable if we compared two observational units which differed only by 1 point on that predictor variable, but had the same score on the other predictor variable.

How good this estimate is depends on how well our regression plane fits the underlying population of data that the sample was drawn from...a topic for a later class!

A quick note on the intercept¶

The intercept is the point where the plane crosses the vertical axis (in our case the prestige axis) at the point where both the other axes (education and income) both equal 0. This is why the intercept is a single term, as there is only one point where this crossing occurs! If we didn't include an intercept then our plane would have to pass through the "origin" (e.g. 0 on all three axes) which would limit the number of planes we can fit. Including an intercept term let's us move the plane up and down the vertical axis when fitting it to the data points...

The plot below let's you use some sliders to manipulate a plane by changing the values of the $x_1$ slope, the $x_2$ slope and the intercept. It also shows the origin (where all three axes == 0) as a black dot.

Note: the axes are shown with values relevant to the dataset used on this page - the plane itself is technically infinite, it runs on infinitely outside of all the borders of the graph - a pretty cool concept!