7.2. Multiple Linear Regression

Incorporating Multiple Features¶

Motivation¶

The commute times dataset that we were reintroduced to in Chapter 7.1 has several columns in it, but we’ve only used one – departure_hour – as an input variable.

Let’s take a look at a subset of them.

For example, the first row above tells us that the first recorded commute was on a Monday, which was the 15th of that month.

Using day seems like it would be a good idea, since traffic patterns likely differ based on the day of the week. But, it’s not clear how to use strings, like 'Tue', in a linear model. We’ll address this in a moment. For now, suppose we’d like to fit a linear model that predicts commute time in minutes using both departure_hour and day_of_month (“dom” for short).

Such a hypothesis function is of the form

which is a plane in $\mathbb{R}^3$. Once we figure out how to find $w_0^*$, $w_1^*$, and $w_2^*$, we’ll visualize the resulting plane, not to worry.

(To be crystal clear, the $\cdot$’s in the above equations refer to “regular” scalar multiplication, not the dot product. There are no vectors in the equation above.)

But, how do we find these three optimal model parameters? If we consult the modeling recipe, and choose squared loss, we need to find the values of $w_0^*$, $w_1^*$, and $w_2^*$ that minimize

The solution is to use what we know about spans, projections, and the design matrix. The design matrix for this problem will have 3 columns: a column of 1’s for the intercept, a column of departure_hour values, and a column of dom values.

To find $\vec w^*$, all we need to do is define the observation vector, $\color{orange} \vec y = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}$, and find the $\vec w^*$ that minimizes

which is a solved problem at this point: it’s the $\vec w^*$ that satisfies the normal equation, ${\color{#3d81f6} X^TX} \vec w^* = {\color{#3d81f6} X^T} \color{orange} \vec y$.

Generalized Notation¶

Let me try and state this problem in slightly more general terms, where we have $d$ features rather than just 1 or 2.

As before, subscripts distinguish between individuals (rows) in our dataset, and we will try and use superscripts (in parentheses) to distinguish between features (columns). Specifically, we’ll use the notation $x_i^{(j)}$ to represent the $j$-th feature for the $i$-th individual. In the example above, $x_5^{(1)}$ is the departure hour for the 5th row in the dataset. Think of $x^{(1)}$, $x^{(2)}$, ... as new variable names, like new letters.

Consider the following example dataset.

There are $n = 3$ rows, each of which has $d = 2$ features (minutes is not a feature, it’s what we’re trying to predict).

We can represent each row (day) with a feature vector, $\vec x_i = \begin{bmatrix} x_i^{(1)} \\ x_i^{(2)} \end{bmatrix}$. A feature vector contains all of the information we’re using to make predictions for the $i$-th individual.

In general,

Note that if our model has $d$ features, then there are $d+1$ parameters, $w_0$, $w_1$, $w_2$, ..., $w_d$. There is one parameter for each feature, plus one for the intercept. So, the generalized multiple linear regression model has a hypothesis function of the form

It’d be great if we could express our hypothesis function $h$ as a dot product between a parameter vector $\vec w$ and a feature vector $\vec x_i$, but we can’t: $\vec w$ has $d+1$ components while $\vec x_i$ has $d$ components.

To address this issue, we’ll define the augmented feature vector, $\text{Aug}({\vec x_i})$, which is the vector obtained by adding a 1 to the front of feature vector $\vec x_i$ (much like the design matrix $\color{#3d81f6} X$ has a column of 1’s).

Then, the multiple linear regression model can be written as

The General Problem¶

Now we can finally state the problem of multiple linear regression in its most general form. Suppose we have $n$ data points, $\left({ \vec x_1}, {y_1}\right), \left({ \vec x_2}, {y_2}\right), \ldots, \left({ \vec x_n}, {y_n}\right)$, where each $\vec x_i$ is a feature vector of $d$ features:

We’d like to find a good linear model,

By good, we mean that we’d like to find optimal parameters $w_0^*$, $w_1^*$, ..., $w_d^*$ that minimize mean squared error:

The solution is to define the $n \times (d + 1)$ design matrix $\color{#3d81f6} X$ and $n$-dimensional observation vector $\color{orange} \vec y$:

and to solve the normal equation to find the optimal parameter vector, $\vec{w}^*$:

The $\vec w^*$ that satisfies the equations above minimizes mean squared error, $R_\text{sq}(\vec w)$, so use this $\vec w^*$ to make predictions in $h(\vec x_i) = \vec w^* \cdot \text{Aug}(\vec x_i)$. Once you find that $\vec w^*$,

• ${\color{#004d40} \vec p} = {\color{#3d81f6} X} \vec w^*$ is a prediction vector of predicted $y$-values for the entire dataset. It is also the projection of $\color{orange} \vec y$ onto the span of the columns of $\color{#3d81f6} X$.

• $h(\vec x_i) = \vec w^* \cdot \text{Aug}(\vec x_i)$ is the predicted $y$-value for just the input $\vec x_i$.

Using sklearn¶

Earlier in the course, we saw how to use sklearn to fit a simple linear regression model. Performing multiple linear regression is just as easy.

After calling fit, we can access $\vec w^*$.

The outputs above tell us that the “best way” (according to squared loss) to make commute time predictions using a linear model is using:

This is the plane of best fit given historical data; it is not a causal statement.

Fit LinearRegression objects have a predict method, which can be used to predict commute times given a departure_hour and day_of_month.

Comparing Models¶

While it’s not difficult to implement ourselves, sklearn has a built-in mean_squared_error function.

So, the mean squared error of our plane of best fit is 96.78 $\text{minutes}^2$.

Throughout this section, we’ll fit various different models to the commute times dataset, and it’ll be useful to keep track of their mean squared errors in one place. I’ll do so in a Python dictionary. (I’m intentionally showing you more code than I typically have, so you have some sense of how to do this all yourself – I hope you don’t mind!)

As you can see, adding day_of_month barely reduced our model’s mean squared error, which matches our intuition that knowing whether it’s the 15th or 19th or 3rd of the month doesn’t seem like it would be helpful in predicting commute time. The activity below has you think about why the mean squared error still did decrease, and whether this is always the case.

Activity 1¶

Feature Engineering¶

So far, we’ve limited ourselves to using departure_hour and day_of_month, two features that already came to us in the dataset. In practice, we’ll often need to engineer features, which means creating new features from existing ones, or “transforming” raw data into good features. We might do this to improve model performance, for example, in case the relationship between the features and output variable isn’t truly linear, or maybe even to enable the use of categorical features.

Example: One Hot Encoding¶

Suppose we’d like to use the day column from commutes_full as a feature.

One naive approach would be to convert each day value to a number, e.g. "Mon" is 1, "Tue" is 2, "Wed" is 3, etc. The reason this is a bad idea is that this approach implies that Monday is somehow “less than” Tuesday and should impact our model’s predictions less than Tuesday (remember that this column will be multiplied by a parameter in order to make predictions). This is what’s called an ordinal encoding, and these only make sense when the features have some inherent order that is meaningful to the prediction problem.

Instead, we’ll perform one hot encoding, which turns a categorical feature into several binary features: one per unique value of the categorical feature.

Since day has 5 unique values, we’ll need to create 5 new binary features.

Suppose we want to build a model that predicts minutes using departure_hour, day_of_month, and one hot encoded day. That model would have $1 + 1 + 5 = 7$ features, so its design matrix would have $7 + 1 = 8$ columns!

The actual process of one hot encoding can be done efficiently using sklearn.preprocessing.OneHotEncoder. You’ll see how to use this in a homework assignment.

Note that while it looks like the day == Wed, day == Thu, and day == Fri columns are all 0’s, this is just because there were no Wednesday, Thursday, or Friday values in the first five rows of the dataset. Somewhere later on in the dataset, there were rows with Wednesday, Thursday, or Friday values, and so those columns have non-zero values later on. If there were never any Wednesday, Thursday, or Friday values, then we wouldn’t make those columns in the first place (just like we never made a Saturday column, or a Zebra column, etc.).

Question: What is the rank of the design matrix, whose first five rows are shown above?

You might notice that among the 5 binary columns, in a particular row, exactly one of them is 1, and the other 4 are 0. This means that if we sum together the 5 binary columns, we’ll get a column of 1’s – which is already the first column in our design matrix!

This means that, by default, when we perform one hot encoding, the design matrix doesn’t have full rank. This doesn’t stop us from projecting $\color{orange} \vec y$ onto $\text{colsp}({\color{#3d81f6} X})$, but it does mean that there are infinitely many possible optimal parameters $\vec w^*$. Again, all of these would give us the same predictions and same mean squared error, but it’s annoying to have to deal with this setting.

So, a common solution when performing one hot encoding is to drop one of the binary columns. This doesn’t change the information conveyed by the other features, and doesn’t change $\text{colsp}({\color{#3d81f6} X})$, which is what matters for the quality of the model’s predictions. And intuitively, even if we don’t have the value of day == Fri, we know it’s Friday if day == Mon, day == Tue, day == Wed, and day == Thu are all 0.

Now that we’ve converted day to a numerical variable, we can use it as input in a regression model. Here’s the model we’ll try to fit:

This model has 6 features (1 + 1 + 4) and hence 7 parameters. Its design matrix is $n \times 7$.

Notice that below, I’ve set fit_intercept=False, because I manually added a column of 1’s to the design matrix, so I don’t need sklearn to add another one.

So, our linear model to predict commute time given departure_hour, day_of_month, and day (Mon, Tue, Wed, or Thu) is:

While this model has 6 features, and thus requires 7 dimensions to graph, it’s really a collection of five parallel planes in 3D, all with slightly different $z$-intercepts! Remember that $\text{day}_i == \text{Mon}$, $\text{day}_i == \text{Tue}$, etc. aren’t all free variables: if one of them is 1, the others must be 0. There are 5 cases to consider, and each one corresponds to one of these 5 parallel planes. You’re exploring this idea in Homework 7, Question 3.

If we do want to visualze the model in $\mathbb{R}^2$, we need to pick a single feature to place on the $x$-axis.

Despite being a linear model, this model doesn’t look like a straight line, since it’s linear in terms of all 6 features, but when you project it into 2D, it no longer appears linear.

How does this new model compare to our previous models?

Note that adding the day of the week decreased our model’s mean squared error significantly: this one hot encoded model is our best yet.

Activity 2¶

model = LinearRegression()
model.fit(X, y)

# Used in the answer choices below.
ws = np.append(model.intercept_, model.coef_)
preds = model.predict(X)
squares = X.shape[0] * mean_squared_error(y, preds)

Example: Numerical Transformations¶

Another common step in feature engineering is to produce new numerical features out of existing ones.

As an example, I could fit a degree 2 polynomial model

by producing the design matrix

Below, we’ll apply this thinking to the departure_hour feature. We’ll create these features manually, though in homework assignments you’ll see how to do this easily and more efficiently using sklearn.preprocessing.PolynomialFeatures.

This tells us that the best-fitting quadratic model is

which looks like:

This is a linear model, even though the features themselves are made up of non-linear functions.

If we treat $x_i$ and $x_i^2$ as two separate features, then each data point lives in $\mathbb{R}^3$ with coordinates

In that feature space, the model is

so the predictions lie on a plane. It only looks curved when we project back down to the usual 2D plot of departure_hour versus minutes, because the points we draw must satisfy the relationship $x_i^2 = (x_i)^2$.

We could follow the same idea to fit a degree 3 polynomial model, or a degree 4 polynomial model, or a degree 5 polynomial model, and so on. Just to try this one more time, let’s fit a degree 3 model.

Overfitting¶

Is this model better than the simple linear regression model or quadratic model from earlier?

75.27, the mean squared error of the cubic model, is lower than the simple linear regression model’s mean squared error of 97.04 (but worse than we got when one hot encoding day) and lower than the quadratic model’s mean squared error of 75.32.

Keep in mind that adding complex features doesn’t always equate to a “better model”, even if it lowers the mean squared error on the training data. To illustrate what I mean, what if we fit a degree 10 polynomial to the data?

Without computing it, we can tell this model’s mean squared error is lower than that of the simple linear, quadratic, or cubic models, since it passes closer to the training data. But, this model is overfit to the training data. Remember, our eventual goal is to build models that generalize well to new data, and it seems unlikely that this degree 10 polynomial reflects the true relationship between departure hour and commute time in nature.

The question, then, is how to choose the right model, e.g. how to choose the right polynomial degree, if we’re committed to making polynomial features. The solution is to intentionally split our data into training data and test data, and use only the training data to pick the features we’d like to include (e.g. polynomial degree). Then, we pick the combination of features that performed best on the test data, since that combination of features seems most likely to generalize well on unseen data. You’ll explore this paradigm in homework assignments.

In the above examples, I created polynomial features out of departure_hour only, but there’s nothing stopping me from creating features out of day_of_month too, or out of both of them simultaneously. A perfectly valid model is

which, using our more generic syntax, would look like

You may wonder, how do I know which features to create? Through a variety of techniques:

• Domain knowledge

• Trial and error

• By visualizing the data

In the polynomial examples above, the data didn’t particularly look like it was quadratic or cubic, so a linear model sufficed. But if we’re given a dataset that clearly has some non-linear relationship, visualization can help us identify what features might improve model performance.

As the dimensionality of our data increases, though, visualization becomes less possible. (How do we visualize a dataset with 100 features?) Strategies for reducing the dimensionality of our data will be explored in Chapter 10.

Linear in the Parameters¶

What is the extent to which I can use linear regression to fit a model? As long as a model can be expressed in the form

for some parameter vector $\vec w$ and some feature vector $\vec x_i$, then it can be fit using linear regression (i.e. by creating a design matrix $X$ and finding $\vec w^*$ by solving the normal equations). We say such a model is linear in the parameters. The choice of features to include in $\vec x_i$ is up to us.

For example,

• If $x_i$ is a scalar, then

  $$h(x_i) = w_0 + w_1 x_i + w_2 \cos(x_i) + w_3 e^{x_i} = \vec w \cdot \begin{bmatrix} 1 \\ x_i \\ \cos(x_i) \\ e^{x_i} \end{bmatrix}$$

  is linear in the parameters. It would be linear in the parameters even if the $w_0$ term wasn’t there; then, we wouldn’t need to augment $\vec x_i$ with a 1.

• If $x_i^{(1)}$ and $x_i^{(2)}$ are scalars (e.g. height and weight), then

  $$h(\vec x_i) = w_0 + w_1 x_i^{(1)} + w_2 x_i^{(2)} + w_3 x_i^{(1)} x_i^{(2)} + w_4 \frac{e^{x_i^{(1)}}}{x_i^{(2)}} = \vec w \cdot \begin{bmatrix} 1 \\ x_i^{(1)} \\ x_i^{(2)} \\ x_i^{(1)} x_i^{(2)} \\ \frac{e^{x_i^{(1)}}}{x_i^{(2)}} \end{bmatrix}$$

  is linear in the parameters.

• If $x_i$ is a scalar, then

  $$h(x_i) = w_0 + w_1 \cos(x_i) + \sin(w_2 x_i)$$

  is not linear in the parameters, because $w_2$ is within the $\sin$ function. This model can’t be written as $\vec w \cdot \text{Aug}(\vec x_i)$ for any choice of $\vec x_i$.