Let's jump in and take a look at some "real-life" examples in which a multiple linear regression model is used. Make sure you notice, in each case, that the model has more than one predictor. You might also try to pay attention to the similarities and differences among the examples and their resulting models. Most of all, don't worry about mastering all of the details now. In the upcoming lessons, we will revisit similar examples in greater detail. For now, my hope is that these examples leave you with an appreciation of the richness of multiple regression.

Are a person's brain size and body size predictive of his or her intelligence?

Interested in answering the above research question, some researchers (Willerman, et al, 1991) collected the following data (IQ Size data) on a sample of n = 38 college students:

• Response (y): Performance IQ scores (PIQ) from the revised Wechsler Adult Intelligence Scale. This variable served as the investigator's measure of the individual's intelligence.
• Potential predictor (\(x_{1}\)): Brain size based on the count obtained from MRI scans (given as count/10,000).
• Potential predictor (\(x_{2}\)): Height in inches.
• Potential predictor (\(x_{3}\)): Weight in pounds.

As always, the first thing we should want to do when presented with a set of data is to plot it. And, of course, plotting the data is a little more challenging in the multiple regression setting, as there is one scatter plot for each pair of variables. Not only do we have to consider the relationship between the response and each of the predictors, but we also have to consider how the predictors are related to each other.

A common way of investigating the relationships among all of the variables is by way of a "scatter plot matrix." Basically, a scatter plot matrix contains a scatter plot of each pair of variables arranged in an orderly array. Here's what one version of a scatter plot matrix looks like for our brain and body size example:

For each scatterplot in the matrix, the variable on the y-axis appears at the left end of the plot's row and the variable on the x-axis appears at the bottom of the plot's column. Try to identify the variables on the y-axis and x-axis in each of the six scatter plots appearing in the matrix. You can check your understanding by selecting the icons appearing in the above matrix.

Incidentally, in case you are wondering, the tick marks on each of the axes are located at 25%, and 75% of the data range from the minimum. That is:

• the first tick = ((maximum - minimum) * 0.25) + minimum
• the second tick = ((maximum - minimum) * 0.75) + minimum

Now, what does a scatterplot matrix tell us? Of course, one use of the plots is simple data checking. Are there any egregiously erroneous data errors? The scatter plots also illustrate the "marginal relationships" between each pair of variables without regard to the other variables. For example, it appears that brain size is the best single predictor of PIQ, but none of the relationships are particularly strong. In multiple linear regression, the challenge is to see how the response y relates to all three predictors simultaneously.

We always start a regression analysis by formulating a model for our data. One possible multiple linear regression model with three quantitative predictors for our brain and body size example is:

\(y_i=(\beta_0+\beta_1x_{i1}+\beta_2x_{i2}+\beta_3x_{i3})+\epsilon_i\) 

where:

• \(y_{i}\) is the intelligence (PIQ) of student i
• \(x_{i1}\) is the brain size (MRI) of student i
• \(x_{i2}\) is the height (Height) of student i
• \(x_{i3}\) is the weight (Weight) of student i

and the independent error terms \(\epsilon_{i}\) follow a normal distribution with mean 0 and equal variance \(\sigma_{2}\).

Of course, our interest in performing a regression analysis is almost always to answer some sort of research question. Can you think of some research questions that the researchers might want to answer here? How about the following set of questions? What procedure would you use to answer each research question? (Do the procedures that appear in parentheses seem reasonable?)

• Which, if any, predictors — brain size, height, or weight — explain some of the variations in intelligence scores? (Conduct hypothesis tests for individually testing whether each slope parameter could be 0.)
• What is the effect of brain size on PIQ, after taking into account height and weight? (Calculate and interpret a confidence interval for the brain size slope parameter.)
• What is the PIQ of an individual with given brain size, height, and weight? (Calculate and interpret a prediction interval for the response.)

Let's take a look at the output we obtain when we ask Minitab to estimate the multiple regression model we formulated above:

My hope is that you immediately observe that much of the output looks the same as before! The only substantial differences are:

• More predictors appear in the estimated regression equation and therefore also in the column labeled "Term" in the coefficients table.
• There is an additional row for each predictor term in the Analysis of Variance Table. By default in Minitab, these represent the reductions in the error sum of squares for each term relative to a model that contains all of the remaining terms (so-called Adjusted or Type III sums of squares). It is possible to change this using the Minitab Regression Options to instead use Sequential or Type I sums of squares, which represent the reductions in error sum of squares when a term is added to a model that contains only the terms before it.

We'll learn more about these differences later, but let's focus now on what you already know. The output tells us that:

• The \(R^{2}\) value is 29.49%. This tells us that 29.49% of the variation in intelligence, as quantified by PIQ, is reduced by taking into account brain size, height, and weight.
• The Adjusted \(R^{2}\) value — denoted "R-sq(adj)" — is 23.27%. When considering different multiple linear regression models for PIQ, we could use this value to help compare the models.
• The P-values for the t-tests appearing in the coefficients table suggest that the slope parameters for Brain (P = 0.001) and Height (P = 0.033) are significantly different from 0, while the slope parameter for Weight (P = 0.998) is not.
• The P-value for the analysis of variance F-test (P = 0.007) suggests that the model containing Brain, Height, and Weight is more useful in predicting intelligence than not taking into account the three predictors. (Note that this does not tell us that the model with the three predictors is the best model!)

So, we already have a pretty good start on this multiple linear regression stuff. Let's take a look at another example.

Some mammals burrow into the ground to live. Scientists have found that the quality of the air in these burrows is not as good as the air aboveground. In fact, some mammals change the way that they breathe in order to accommodate living in poor air quality conditions underground.

Some researchers (Colby, et al, 1987) wanted to find out if nestling bank swallows, which live in underground burrows, also alter how they breathe. The researchers conducted a randomized experiment on n = 120 nestling bank swallows. In an underground burrow, they varied the percentage of oxygen at four different levels (13%, 15%, 17%, and 19%) and the percentage of carbon dioxide at five different levels (0%, 3%, 4.5%, 6%, and 9%). Under each of the resulting 5 × 4 = 20 experimental conditions, the researchers observed the total volume of air breathed per minute for each of 6 nestling bank swallows. In this way, they obtained the following data (Baby birds) on the n = 120 nestling bank swallows:

• Response (y): percentage increase in "minute ventilation," (Vent), i.e., the total volume of air breathed per minute.
• Potential predictor (\(x_{1}\)): percentage of oxygen (O2) in the air the baby birds breathe.
• Potential predictor (\(x_{2}\)): percentage of carbon dioxide (CO2) in the air the baby birds breathe.

Here's a scatter plot matrix of the resulting data obtained by the researchers:

What does this particular scatter plot matrix tell us? Do you buy into the following statements?

• There doesn't appear to be a substantial relationship between minute ventilation (Vent) and the percentage of oxygen (O2).
• The relationship between minute ventilation (Vent) and the percentage of carbon dioxide (CO2) appears to be curved and with increasing error variance.
• The plot between the percentage of oxygen (O2) and the percentage of carbon dioxide (CO2) is the classical appearance of a scatter plot for the experimental conditions. The plot suggests that there is no correlation at all between the two variables. You should be able to observe from the plot the 4 levels of O2 and the 5 levels of CO2 that make up the 5×4 = 20 experimental conditions.

When we have one response variable and only two predictor variables, we have another sometimes useful plot at our disposal, namely a "three-dimensional scatter plot:"

If we added the estimated regression equation to the plot, what one word do you think describes what it would look like? Drag the slider on the bottom of the graph above to show the plot of the estimated regression equation for this data. Does it make sense that it looks like a "plane?" Incidentally, it is still important to remember that the plane depicted in the plot is just an estimate of the actual plane in the population that we are trying to study.

Here is a reasonable "first-order" model with two quantitative predictors that we could consider when trying to summarize the trend in the data:

\(y_i=(\beta_0+\beta_1x_{i1}+\beta_2x_{i2})+\epsilon_i\)

where:

• \(y_{i}\) is the percentage of minute ventilation of nestling bank swallow i
• \(x_{i1}\) is the percentage of oxygen exposed to nestling bank swallow i
• \(x_{i2}\) is the percentage of carbon dioxide exposed to nestling bank swallow i

and the independent error terms \(\epsilon_i\) follow a normal distribution with mean 0 and equal variance \(\sigma^{2}\).

The adjective "first-order" is used to characterize a model in which the highest power on all of the predictor terms is one. In this case, the power on \(x_{i1}\), although typically not shown, is one. And, the power on \(x_{i2}\) is also one, although not shown. Therefore, the model we formulated can be classified as a "first-order model." An example of a second-order model would be \(y=\beta_0+\beta_1x+\beta_2x^2+\epsilon\).

Do you have your research questions ready? How about the following set of questions? (Do the procedures that appear in parentheses seem appropriate in answering the research question?)

• Is oxygen related to minute ventilation, after taking into account carbon dioxide? (Conduct a hypothesis test for testing whether the O2 slope parameter is 0.)
• Is carbon dioxide related to minute ventilation, after taking into account oxygen? (Conduct a hypothesis test for testing whether the CO2 slope parameter is 0.)
• What is the mean minute ventilation of all nestling bank swallows whose breathing air is comprised of 15% oxygen and 5% carbon dioxide? (Calculate and interpret a confidence interval for the mean response.)

Here's the output we obtain when we ask Minitab to estimate the multiple regression model we formulated above:

What do we learn from the Minitab output?

• Only 26.82% of the variation in minute ventilation is reduced by taking into account the percentages of oxygen and carbon dioxide.
• The P-values for the t-tests appearing in the coefficients table suggest that the slope parameter for carbon dioxide level (P < 0.001) is significantly different from 0, while the slope parameter for oxygen level (P = 0.408) is not. Does this conclusion appear consistent with the above scatter plot matrix and the three-dimensional plot? Yes!
• The P-value for the analysis of variance F-test (P < 0.001) suggests that the model containing oxygen and carbon dioxide levels is more useful in predicting minute ventilation than not taking into account the two predictors. (Again, the F-test does not tell us that the model with the two predictors is the best model! For one thing, we have performed no model checking yet!)

Ugh! All of these definitions! Let's take a look at an example just to convince ourselves that, yes, indeed the least squares estimates are obtained by the following matrix formula:

\(b=\begin{bmatrix}
b_0\\
b_1\\
\vdots\\
b_{p-1}
\end{bmatrix}=(X^{'}X)^{-1}X^{'}Y\)

Let's consider the data in the Soap Suds dataset, in which the height of suds (y = suds) in a standard dishpan was recorded for various amounts of soap (x = soap, in grams) (Draper and Smith, 1998, p. 108). Using Minitab to fit the simple linear regression model to these data, we obtain:

Let's see if we can obtain the same answer using the above matrix formula. We previously showed that:

\(X^{'}X=\begin{bmatrix}
n & \sum_{i=1}^{n}x_i \\
\sum_{i=1}^{n}x_i  & \sum_{i=1}^{n}x_{i}^{2}
\end{bmatrix}\)

Using the calculator function in Minitab, we can easily calculate some parts of this formula:

That is, the 2 × 2 matrix X'X is:

\(X^{'}X=\begin{bmatrix}
7 & 38.5\\
38.5& 218.75
\end{bmatrix}\)

And, the 2 × 1 column vector X'Y is:

\(X^{'}Y=\begin{bmatrix}
\sum_{i=1}^{n}y_i\\
\sum_{i=1}^{n}x_iy_i
\end{bmatrix}=\begin{bmatrix}
347\\
1975
\end{bmatrix}\)

So, we've determined X'X and X'Y. Now, all we need to do is to find the inverse (X'X)^-1. As mentioned before, it is very messy to determine inverses by hand. Letting computer software do the dirty work for us, it can be shown that the inverse of X'X is:

\((X^{'}X)^{-1}=\begin{bmatrix}
4.4643 & -0.78571\\
-0.78571& 0.14286
\end{bmatrix}\)

And so, putting all of our work together, we obtain the least squares estimates:

\(b=(X^{'}X)^{-1}X^{'}Y=\begin{bmatrix}
4.4643 & -0.78571\\
-0.78571& 0.14286
\end{bmatrix}\begin{bmatrix}
347\\
1975
\end{bmatrix}=\begin{bmatrix}
-2.67\\
9.51
\end{bmatrix}\)

That is, the estimated intercept is \(b_{0}\) = -2.67 and the estimated slope is \(b_{1}\) = 9.51. Aha! Our estimates are the same as those reported by Minitab:

within rounding error!

A designed experiment is done to assess how the moisture content and sweetness of a pastry product affect a taster’s rating of the product (Pastry dataset). In a designed experiment, the eight possible combinations of four moisture levels and two sweetness levels are studied. Two pastries are prepared and rated for each of the eight combinations, so the total sample size is n = 16. The y-variable is the rating of the pastry. The two x-variables are moisture and sweetness. The values (and sample sizes) of the x-variables were designed so that the x-variables were not correlated.

A plot of moisture versus sweetness (the two x-variables) is as follows:

Notice that the points are on a rectangular grid so the correlation between the two variables is 0. (Please Note: we are not able to see that actually there are 2 observations at each location of the grid!)

The following figure shows how the two x-variables affect the pastry rating.

There is a linear relationship between rating and moisture and there is also a sweetness difference. The Minitab results given in the following output are for three different regressions - separate simple regressions for each x-variable and a multiple regression that incorporates both x-variables.

There are three important features to notice in the results:

1. The sample coefficient that multiplies Moisture is 4.425 in both the simple and the multiple regression. The sample coefficient that multiplies Sweetness is 4.375 in both the simple and the multiple regression. This result does not generally occur; the only reason that it does, in this case, is that Moisture and Sweetness are not correlated, so the estimated slopes are independent of each other. For most observational studies, predictors are typically correlated and estimated slopes in a multiple linear regression model do not match the corresponding slope estimates in simple linear regression models.

2. The \(R^{2}\) for the multiple regression, 95.21%, is the sum of the \(R^{2}\) values for the simple regressions (79.64% and 15.57%). Again, this will only happen when we have uncorrelated x-variables.

3. The variable Sweetness is not statistically significant in the simple regression (p = 0.130), but it is in the multiple regression. This is a benefit of doing a multiple regression. By putting both variables into the equation, we have greatly reduced the standard deviation of the residuals (notice the S values). This in turn reduces the standard errors of the coefficients, leading to greater t-values and smaller p-values.

(Data source: Applied Regression Models, (4th edition), Kutner, Neter, and Nachtsheim).

The data are from n = 214 females in statistics classes at the University of California at Davis (Stat Females dataset). The variables are y = student’s self-reported height, \(x_{1}\) = student’s guess at her mother’s height, and \(x_{2}\) = student’s guess at her father’s height. All heights are in inches. The scatterplots below are of each student’s height versus the mother’s height and the student’s height against the father’s height.

Both show a moderate positive association with a straight-line pattern and no notable outliers.

Interpretations

The first two lines of the Minitab output show that the sample multiple regression equations is predicted student height = 18.55 + 0.3035 × mother’s height + 0.3879 × father’s height:

To use this equation for prediction, we substitute specified values for the two parents’ heights.

We can interpret the “slopes” in the same way that we do for a simple linear regression model but we have to add the constraint that values of other variables remain constant. For example:

– When the father’s height is held constant, the average student height increases by 0.3035 inches for each one-inch increase in the mother’s height.

– When the mother’s height is held constant, the average student height increases by 0.3879 inches for each one-inch increase in the father’s height.

• The p-values given for the two x-variables tell us that student height is significantly related to each.
• The value of \(R^{2}\) = 43.35% means that the model (the two x-variables) explains 43.35% of the observed variation in student heights.
• The value S = 2.03115 is the estimated standard deviation of the regression errors. Roughly, it is the average absolute size of a residual.

Residual Plots

Just as in simple regression, we can use a plot of residuals versus fits to evaluate the validity of assumptions. The residual plot for these data is shown in the following figure:

It looks about as it should - a random horizontal band of points. Other residual analyses can be done exactly as we did for simple regression. For instance, we might wish to examine a normal probability plot of the residuals. Additional plots to consider are plots of residuals versus each x-variable separately. This might help us identify sources of curvature or non-constant variance.

Data from n = 113 hospitals in the United States are used to assess factors related to the likelihood that a hospital patient acquires an infection while hospitalized. The variables here are y = infection risk, \(x_{1}\) = average length of patient stay, \(x_{2}\) = average patient age, \(x_{3}\) = measure of how many x-rays are given in the hospital (Hospital Infection dataset). The Minitab output is as follows:

Interpretations for this example include:

• The p-value for testing the coefficient that multiplies Age is 0.330. Thus we cannot reject the null hypothesis \(H_0 \colon \beta_{2}\) = 0. The variable Age is not a useful predictor within this model that includes Stay and Xrays.
• For the variables Stay and X-rays, the p-values for testing their coefficients are at a statistically significant level so both are useful predictors of infection risk (within the context of this model!).
• We usually don’t worry about the p-value for Constant. It has to do with the “intercept” of the model and seldom has any practical meaning. It also doesn’t give information about how changing an x-variable might change y-values.

(Data source: Applied Regression Models, (4th edition), Kutner, Neter, and Nachtsheim).

For a sample of n = 20 individuals, we have measurements of y = body fat, \(x_{1}\) = triceps skinfold thickness, \(x_{2}\) = thigh circumference, and \(x_{3}\) = midarm circumference (Body Fat dataset). Minitab results for the sample coefficients, MSE (highlighted), and \(\left(X^{T} X \right)^{−1}\) are given below:

\(\left(X^{T} X \right)^{−1}\) - (calculated manually, see note below)

The variance-covariance matrix of the sample coefficients is found by multiplying each element in \(\left(X^{T} X \right)^{−1}\) by MSE. Common notation for the resulting matrix is either \(s^{2}\)(b) or \(se^{2}\)(b). Thus, the standard errors of the coefficients given in the Minitab output can be calculated as follows:

• Var(\(b_{0}\)) = (6.15031)(1618.87) = 9956.55, so se(\(b_{0}\)) = \(\sqrt{9956.55}\) = 99.782.
• Var(\(b_{1}\)) = (6.15031)(1.4785) = 9.0932, so se(\(b_{1}\)) = \(\sqrt{9.0932}\) = 3.016.
• Var(\(b_{2}\)) = (6.15031)(1.0840) = 6.6669, so se(\(b_{2}\)) = \(\sqrt{6.6669}\) = 2.582.
• Var(\(b_{3}\)) = (6.15031)(0.4139) = 2.54561, so se(\(b_{3}\)) = \(\sqrt{2.54561}\) = 1.595.

As an example of covariance and correlation between two coefficients, we consider \(b_{1 }\)and \(b_{2}\).

• Cov(\(b_{1}\), \(b_{2}\)) = (6.15031)(−1.2648) = −7.7789. The value -1.2648 is in the second row and third column of \(\left(X^{T} X \right)^{−1}\). (Keep in mind that the first row and first column give information about \(b_0\), so the second row has information about \(b_{1}\), and so on.)
• Corr(\(b_{1}\), \(b_{2}\)) = covariance divided by product of standard errors = −7.7789 / (3.016 × 2.582) = −0.999.

The extremely high correlation between these two sample coefficient estimates results from a high correlation between the Triceps and Thigh variables. The consequence is that it is difficult to separate the individual effects of these two variables.

If all x-variables are uncorrelated with each other, then all covariances between pairs of sample coefficients that multiply x-variables will equal 0. This means that the estimate of one beta is not affected by the presence of the other x-variables. Many experiments are designed to achieve this property. With observational data, however, we’ll most likely not have this happen.

(Data source: Applied Regression Models, (4th edition), Kutner, Neter, and Nachtsheim).