18.1 – Multiple Linear Regression
Introduction
Last time we introduced simple linear regression:
- one independent X variable
- one dependent Y variable.
The linear relationship between Y and X was estimated by the method of Ordinary Least Squares (OLS). OLS minimizes the sum of squared distances between the observed responses,
and responses predicted
by the line. Simple linear regression is analogous to our one-way ANOVA — one outcome or response variable and one factor or predictor variable (Chapter 12.2).
But the world is complicated and so, our one-way ANOVA was extended to the more general case of two or more predictor (factor) variables (Chapter 14). As you might have guessed by now, we can extend simple regression to include more than one predictor variable. In fact, combining ANOVA and regression gives you the general linear model! And, you should not be surprised that statistics has extended this logic to include not only multiple predictor variables, but also multiple response variables. Multiple response variables falls into a category of statistics called multivariate statistics.
Like multi-way ANOVA, multiple regression is the extension of simple linear regression from one independent predictor variable to include two or more predictors. The benefit of this extension is obvious — our models gain realism. All else being equal, the more predictors, the better the model will be at describing and/or predicting the response. Things are not all equal, of course, and we’ll consider two complications of this basic premise, that more predictors are best; in some cases they are not.
However, before discussing the exceptions or even the complications of a multiple linear regression model, we begin by obtaining estimates of the full model, then introduce aspects of how to evaluate the model. We also introduce and whether a reduced model may be the preferred model.
R code
Multiple regression is easy to do in Rcmdr — recall that we used the general linear model function, lm(), to analyze one-way ANOVA and simple linear regression. In R Commander, we access lm() by
Rcmdr: Statistics → Fit model → Linear model
You may, however, access linear regression through R Commander
We use the same general linear model function for cases of multi-way ANOVA and for multiple regression problems. Simply enter more than one ratio-scale predictor variable and boom!
You now have yourself a multiple regression. You would then proceed to generate the ANOVA table for hypothesis testing
Rcmdr: Models → Hypothesis testing → ANOVA tables
From the output of the regression command, estimates of the coefficients along with standard errors for the estimate and results of t-tests for each coefficient against the respective null hypotheses for each coefficient are also provided. In our discussion of simple linear regression we introduced the components: the intercept, the slope, as well as the concept of model fit, as evidenced by R2, the coefficient of determination. These components exist for the multiple regression problem, too, but now we call the slopes partial regression slopes because there are more than one.
Our full multiple regression model becomes
![]()
where the coefficients β1, β2, … βn are the partial regression slopes and β0 is the Y-intercept for a model with 1 – n predictor variables. Each coefficient has a null hypothesis, each has a standard error, and therefore, each coefficient can be tested by t-test.
Now, regression, like ANOVA, is an enormous subject and we cannot do it justice in the few days we will devote to it. We can, however, walk you through a fairly typical example. I’ve posted a small data set diabetesCholStatin at the end of this page. Scroll down or click here. View the data set and complete your basic data exploration routine: make scatterplots and box plots. We think (predict) that body size and drug dose cause variation in serum cholesterol levels in adult men. But do both predict cholesterol levels?
Selecting the best model
We have two predictor variables, and we can start to see whether none, one, or both of the predictors contribute to differences in cholesterol levels. In this case, both contribute significantly. The power of multiple regression approaches is that it provides a simultaneous test of a model which may have many explanatory variables deemed appropriate to describe a particular response. More generally, it is sometimes advisable to think more philosophically about how to select a best model.
In model selection, some would invoke Occam’s razor — given a set of explanations, the simplest should be selected — to justify seeking simpler models. There are a number of approaches (forward selection, backward selection, or stepwise selection), and the whole effort of deciding among competing models is complicated with a number of different assumptions, strengths and weaknesses. I refer you to the discussion below, which of course is just a very brief introduction to a very large subject in (bio)statistics!
Let’s get the full regression model
The statistical model is
![]()
As written in R format, our model is ChLDL ~ BMI + Dose + Statin.
Note 1. BMI is ratio scale and Statin is categorical (two levels: Statin1, Statin2). Dose can be viewed as categorical, with five levels (5, 10, 20, 40, 80 mg), interval scale, or ratio scale. If we are make the assumption that the difference between 5, 10, up to 80 is meaningful, and that the effect of dose is at least proportional if not linear with respect to ChLDL, then we would treat Dose as ratio scale, not interval scale. That’s what we did here.
We can now proceed in R Commander to fit the model.
Rmdr: Statistics → Fit models → Linear model
How the model is inputted into linear model menu is shown in Figure 1.

Figure 1. Screenshot of Rcmdr linear model menu with our model elements in place.
The output
summary(LinearModel.1)
Call:
lm(formula = ChLDL ~ BMI + Dose + Statin, data = cholStatins)
Residuals:
Min 1Q Median 3Q Max
-3.7756 -0.5147 -0.0449 0.5038 4.3821
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.016715 1.178430 0.863 0.39041
BMI 0.058078 0.047012 1.235 0.21970
Dose -0.014197 0.004829 -2.940 0.00411 **
Statin[Statin2] 0.514526 0.262127 1.963 0.05255 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.31 on 96 degrees of freedom
Multiple R-squared: 0.1231, Adjusted R-squared: 0.09565
F-statistic: 4.49 on 3 and 96 DF, p-value: 0.005407
Question. What are the estimates of the model coefficients (rounded)?
= intercept = 1.017
= slope for variable BMI = 0.058
= slope for variable Dose = -0.014
= slope for variable Statin = -0.515
Question. Which of the three coefficients were statistically different from their null hypothesis (p-value < 5%)?
Answer: Only
coefficient was judged statistically significant at the Type I error level of 5% (p = 0.0041). Of the four null hypotheses we have for the coefficients:
= 0
= 0
= 0; we only reject the null hypothesis for Dosecoefficient.
= 0)
The important take-home concept is about the lack of a direct relationship between the magnitude of the estimate of the coefficient and the likelihood that it will be statistically significant! In absolute value terms
, but
was not close to statistical significance (p = 0.220).
We generate a 2D scatterplot and include the regression lines (by group=Statin) to convey the relationship between at least one of the predictors (Fig 2).

Figure 2. Scatter plot of predicted LDL against dose of a statin drug. Regression lines represent the different statin drugs (Statin1, Statin2).
Question. Based on the graph, can you explain why there will be no statistical differences between levels of the statin drug type, Statin1 (shown open circles) vs. Statin2 (shown closed red circles).
Because we have two predictors (BMI and Statin Dose), you may also elect to use a 3D-scatterplot. Here’s one possible result (Fig 3).

Figure 3. 3D plot of BMI and dose of Statin drugs on change in LDL levels (green Statin2, blue Statin1).
R code for Figure 3.
Graph made in Rcmdr: Graphs → 3D Graph → 3D scatterplot …
scatter3d(ChLDL~BMI+Dose|Statin, data=diabetesCholStatin, fit="linear", residuals=TRUE, parallel=FALSE, bg="white", axis.scales=TRUE, grid=TRUE, ellipsoid=FALSE)
Note 2: Figure 3 is a challenging graphic to interpret. I wouldn’t use it because it doesn’t convey a strong message. With some effort we can see the two planes representing mean differences between the two statin drugs across all predictors, but it’s a stretch. No doubt the graph can be improved by changing colors, for example, but I think the 2d plot Figure 2 works better). Alternatively, if the platform allows, you can use animation options to help your reader see the graph elements. Interactive graphics are very promising and, again, unsurprisingly, there are several R packages available. For this example, plot3d() of the package rgl can be used. Figure 4 is one possible version; I saved images and made animated gif.

Figure 4. An example of a possible interactive 3D plot; the file embedded in this page is not interactive, just an animation.
Diagnostic plots
While visualization concerns are important, let’s return to the statistics. All evaluations of regression equations should involve an inspection of the residuals. Inspection of the residuals allows you to decide if the regression fits the data; if the fit is adequate, you then proceed to evaluate the statistical significance of the coefficients.
The default diagnostic plots (Fig 5) R provides are available from
Rcmdr: Models → Graphs → Basic diagnostics plots
Four plots are returned

Figure 5. R’s default regression diagnostic plots.
Each of these diagnostic plots in Figure 4 gives you clues about the model fit.
- Plot of residuals vs. fitted helps you identify patterns in the residuals
- Normal Q-Q plot helps you to see if the residuals are approximately normally distributed
- Scale-location plot provides a view of the spread of the residuals
- The residuals vs. leverage plot allows you to identify influential data points.
We introduced these plots in Chapter 17.8 when we discussed fit of simple linear model to data. My conclusion? No obvious trend in residuals, so linear regression is a fit to the data; data not normally distributed Q-Q plot shows S-shape.
Interpreting the diagnostic plots for this problem
The “Normal Q-Q” plot allows us to view our residuals against a normal distribution (the dotted line). Our residuals do no show an ideal distribution: low for the first quartile, about on the line for intermediate values, then high for the 3rd and 4th quartile residuals. If the data were bivariate normal we would see the data fall along a straight line. The “S-shape” suggests log-transformation of the response and or one or more of the predictor variables.
There also seems to be a pattern in residuals vs the predicted (fitted) values. There is a trend of increasing residuals as cholesterol levels increase, which is particularly evident in the “scale-location” plot. Residuals tended to be positive at low and high doses, but negative at intermediate doses. This suggests that the relationship between predictors and cholesterol levels may not be linear, and it demonstrates what statisticians refer to as a monotonic spread of residuals.
The last diagnostic plot looks for individual points that influence, change, or “leverage” the regression — in other words, if a point is removed, does the general pattern change? If so, then the point had “leverage” and thus we need to decide whether or not to include the datum. diagnostic plots Cook’s distance is a measure of the influence of a point in regression. Points with large Cook’s distance values warrant additional checking.
The multicollinearity problem
Statistical model building is a balancing act by the statistician. While simpler models may be easier to interpret and, perhaps, to use, it is a basic truism that the more predictor variables the model includes, the more realistic the statistical model. However, each additional parameter that is added to the statistical model must be independent of all other parameters already in the model. To the extent that this assumption is violated, the problem is termed multicollinearity. If predictor variables are highly correlated, then they are essentially just linear combinations and do not provide independent evidence. For example, one would naturally not include two core body temperature variables in a statistical model on basal metabolic rate, one in degrees Fahrenheit and the other in degrees Celsius, because it is a simple linear conversion between the two units. This would be an example of structural collinearity: the collinearity is because of misspecification of the model variables. In contrast, collinearity among predictor variables may because the data are themselves correlated. For example, if multiple measures of body size are included: weight, height, length of arm, etc., then we would expect these to be correlated, ie, data multicollinearity.
Collinearity in statistical models may have a number of undesirable effects on a multiple regression model. These include
- estimates of coefficients not stable: with collinearity, values of coefficients depend on other variables in the model; if collinear predictors, then the assumption of independent predictor variables is violated.
- precision of the estimates decreases (standard error of estimates increase).
- statistical power decreases.
- p-values for individual coefficients not trust worthy.
Tolerance and Variance Inflation Factor
Absence of multicollinearity is important assumption of multiple regression. A partial test is to calculate product moment correlations among predictor variables. For example, when we calculate the correlation between BMI and Dose for our model, we get r = 0.101 (P = 0.3186), and therefore would tentatively conclude that there was little correlation between our predictor variables.
A number of diagnostic statistics have been developed to test for multicollinearity. Tolerance for a particular independent variable (Xi) is defined as 1 minus the proportion of variance it shares with the other independent variables in the regression analysis (1 − R2i) (O’Brien 2007). Tolerance reports the proportion of total variance explained by adding the Xi–th predictor variable that is unrelated to the other variables in the model. A small value for tolerance indicates multicollinearity — and that the predictor variable is nearly a perfect combination (linear) of the variables already in the model and therefore should be omitted from the model. Because tolerance is defined in relation to the coefficient of determination you can interpret a tolerance score as the unique variance accounted for by a predictor variable.
A second, related diagnostic of multicollinearity is called the Variance Inflation Factor, VIF. VIF is the inverse of tolerance.
![]()
VIF shows how much of the variance of a regression coefficient is increased because of collinearity with the other predictor variables in the model. VIF is easy to interpret: a tolerance of 0.01 has a VIF of 100; a tolerance of 0.1 has a VIF of 10; a tolerance of 0.5 has a VIF of 2, and so on. Thus, small values of tolerance and large values of VIF are taken as evidence of multicollinearity.
Rcmdr: Models → Numerical diagnostics → Variation-inflation factors
vif(RegModel.2) BMI Dose 1.010256 1.010256
A rule of thumb is that if VIF is greater than 5 then there is multicollinearity; with VIF values close to one we would conclude, like our results from the partial correlation estimate above, that there is little evidence for a problem of collinearity between the two predictor variables. They can therefore remain in the model.
Solutions for multicollinearity
If there is substantial multicollinearity then you cannot simply trust the estimates of the coefficients. Assuming that there hasn’t been some kind of coding error on your part, then you may need to find a solution. One solution is to drop one of the predictor variables and redo the regression model. Another option is to run what is called a Principal Components Regression (PCR). One takes the predictor variables and runs a Principal Component Analysis to reduce the number of variables, then the regression is run on the PCA components. By definition, the PCA components are independent of each other. Another option is to use ridge regression approach.
Note 3. Principal Component Analysis (PCA) is a dimensionality reduction technique that finds new uncorrelated variables (principal components) by maximizing variance. Principal Component Regression (PCR) uses PCA as a preprocessing step to create a regression model. The key difference is that PCA is an unsupervised method for data transformation, whereas PCR is a supervised regression technique that uses the components from PCA to predict an outcome. It’s probably not a great idea to use the acronym PCR in our biology reports; easily confused with polymerase chain reaction
Like any diagnostic rule, however, one should not blindly apply a rule of thumb. A VIF of 10 or more may indicate multicollinearity, but it does not necessarily lead to the conclusion that the linear regression model requires that the researcher reduce the number of predictor variables or analyze the problem using a different statistical method to address multicollinearity as the sole criteria of a poor statistical model. Rather, the researcher needs to address all of the other issues about model and parameter estimate stability, including sample size. Unless the collinearity is extreme (like a correlation of 1.0 between predictor variables!), larger sample sizes alone will work in favor of better model stability (by lowering the sample error) (O’Brien 2007).
Questions
- Write up three learning outcomes for this page. Hint: Point your favorite generative AI to this page and ask for help.
- Please explain why the magnitude of the slope is not the key to statistical significance of a slope? Hint: look at the equation of the t-test for statistical significance of the slope.
- Consider the following scenario. A researcher measures his subjects multiple times for blood pressure over several weeks, then plots all of the values over time. In all, the data set consists of thousands of readings. He then proceeds to develop a model to explain blood pressure changes over time. What kind of collinearity is present in his data set? Explain your choice.
- We noted that
Dosecould be viewed as categorical variable. ConvertDoseto factor variable (fDose) and redo the linear model. Compare the summary output and discuss the additional coefficients.- Use Rcmdr: Data → Manage variables in active data set → Convert numeric Variables to Factors to create a new factor variable
fDose. It’s ok to use the numbers as factor levels.
- Use Rcmdr: Data → Manage variables in active data set → Convert numeric Variables to Factors to create a new factor variable
- We flagged the change in LDL as likely to be not normally distributed. Create a log10-transformed variable for
ChLDLand perform the multiple regression again.- Write the new statistical model
- Obtain the regression coefficients — are they statistically significant?
- Run basic diagnostic plots and evaluate for
Quiz Chapter 18.1
Multiple Linear Regression
Data set
| ID | Statin | Dose | BMI | LDL | ChLDL |
|---|---|---|---|---|---|
| 1 | Statin2 | 5 | 19.5 | 3.497 | 2.714778 |
| 2 | Statin1 | 20 | 20.2 | 4.268 | 1.276483 |
| 3 | Statin2 | 40 | 20.3 | 3.989 | 2.677377 |
| 4 | Statin2 | 20 | 20.3 | 3.502 | 2.430618 |
| 5 | Statin2 | 80 | 20.4 | 3.766 | 1.79463 |
| 6 | Statin2 | 20 | 20.6 | 3.44 | 2.234295 |
| 7 | Statin1 | 20 | 20.7 | 3.414 | 2.635305 |
| 8 | Statin1 | 10 | 20.8 | 3.222 | 0.809181 |
| 9 | Statin1 | 10 | 21.1 | 4.04 | 3.259599 |
| 10 | Statin1 | 40 | 21.2 | 4.429 | 1.763997 |
| 11 | Statin1 | 5 | 21.2 | 3.528 | 3.369377 |
| 12 | Statin1 | 40 | 21.5 | 3.01 | -0.827154 |
| 13 | Statin2 | 20 | 21.6 | 3.393 | 2.11172 |
| 14 | Statin1 | 10 | 21.7 | 4.512 | 3.166238 |
| 15 | Statin1 | 80 | 22 | 5.449 | 3.00833 |
| 16 | Statin2 | 10 | 22.2 | 4.03 | 3.05013 |
| 17 | Statin2 | 40 | 22.2 | 3.911 | 2.646034 |
| 18 | Statin2 | 10 | 22.2 | 3.724 | 2.945656 |
| 19 | Statin1 | 5 | 22.2 | 3.238 | 3.209584 |
| 20 | Statin2 | 10 | 22.5 | 4.123 | 3.088763 |
| 21 | Statin1 | 20 | 22.6 | 3.859 | 5.152548 |
| 22 | Statin1 | 10 | 23 | 4.926 | 2.58483 |
| 23 | Statin2 | 20 | 23 | 3.512 | 2.291975 |
| 24 | Statin1 | 5 | 23 | 3.838 | 1.469 |
| 25 | Statin2 | 20 | 23.1 | 3.548 | 2.34079 |
| 26 | Statin1 | 5 | 23.1 | 3.424 | 1.204346 |
| 27 | Statin1 | 40 | 23.2 | 3.709 | 3.238179 |
| 28 | Statin1 | 80 | 23.2 | 4.786 | 2.748643 |
| 29 | Statin1 | 20 | 23.3 | 4.103 | 1.250082 |
| 30 | Statin1 | 40 | 23.4 | 3.341 | 1.432292 |
| 31 | Statin1 | 10 | 23.5 | 3.828 | 1.381755 |
| 32 | Statin2 | 10 | 23.8 | 4.02 | 3.039187 |
| 33 | Statin1 | 20 | 23.8 | 3.942 | 0.848328 |
| 34 | Statin2 | 20 | 23.8 | 2.89 | 1.721163 |
| 35 | Statin1 | 80 | 23.9 | 3.326 | 1.939346 |
| 36 | Statin1 | 10 | 24.1 | 4.071 | 3.090741 |
| 37 | Statin1 | 40 | 24.1 | 4.222 | 1.322305 |
| 38 | Statin2 | 10 | 24.1 | 3.44 | 2.472223 |
| 39 | Statin1 | 5 | 24.2 | 3.507 | 0.076817 |
| 40 | Statin2 | 20 | 24.2 | 3.647 | 2.425758 |
| 41 | Statin2 | 80 | 24.3 | 3.812 | 1.710575 |
| 42 | Statin2 | 40 | 24.3 | 3.305 | 1.940572 |
| 43 | Statin2 | 5 | 24.3 | 3.455 | 2.502214 |
| 44 | Statin2 | 5 | 24.4 | 4.258 | 3.228008 |
| 45 | Statin1 | 5 | 24.4 | 4.16 | 3.477747 |
| 46 | Statin2 | 80 | 24.4 | 4.128 | 2.063247 |
| 47 | Statin1 | 80 | 24.5 | 4.507 | 3.784422 |
| 48 | Statin1 | 5 | 24.5 | 3.553 | 0.695709 |
| 49 | Statin2 | 10 | 24.5 | 3.616 | 2.69987 |
| 50 | Statin2 | 80 | 24.6 | 3.372 | 1.300401 |
| 51 | Statin2 | 80 | 24.6 | 3.667 | 1.418109 |
| 52 | Statin2 | 5 | 24.7 | 3.854 | 3.126671 |
| 53 | Statin1 | 80 | 24.7 | 3.32 | -1.286439 |
| 54 | Statin2 | 5 | 24.7 | 3.756 | 2.423664 |
| 55 | Statin1 | 40 | 24.8 | 4.398 | 2.907473 |
| 56 | Statin2 | 40 | 24.9 | 3.621 | 2.362429 |
| 57 | Statin1 | 10 | 25 | 3.17 | 1.264656 |
| 58 | Statin1 | 80 | 25.1 | 3.424 | -2.436908 |
| 59 | Statin2 | 10 | 25.1 | 3.196 | 2.001465 |
| 60 | Statin2 | 80 | 25.2 | 3.367 | 1.100704 |
| 61 | Statin1 | 80 | 25.2 | 3.067 | -0.23154 |
| 62 | Statin1 | 20 | 25.3 | 3.678 | 4.662866 |
| 63 | Statin2 | 5 | 25.5 | 4.077 | 2.611705 |
| 64 | Statin1 | 20 | 25.5 | 3.678 | 2.633053 |
| 65 | Statin2 | 5 | 25.6 | 4.994 | 4.180082 |
| 66 | Statin1 | 20 | 25.8 | 3.699 | 1.899031 |
| 67 | Statin1 | 10 | 25.9 | 3.507 | 4.063757 |
| 68 | Statin2 | 20 | 25.9 | 3.445 | 2.303761 |
| 69 | Statin1 | 5 | 26 | 4.025 | 2.501427 |
| 70 | Statin1 | 5 | 26.3 | 3.616 | 0.740863 |
| 71 | Statin2 | 40 | 26.4 | 3.937 | 2.573321 |
| 72 | Statin2 | 40 | 26.4 | 3.823 | 2.363839 |
| 73 | Statin1 | 10 | 26.7 | 4.46 | 2.174198 |
| 74 | Statin2 | 5 | 26.7 | 5.03 | 3.845271 |
| 75 | Statin2 | 10 | 26.7 | 3.73 | 2.708896 |
| 76 | Statin2 | 10 | 26.7 | 3.232 | 2.272627 |
| 77 | Statin1 | 80 | 26.8 | 3.693 | 1.751169 |
| 78 | Statin2 | 80 | 27 | 4.108 | 1.86131 |
| 79 | Statin2 | 40 | 27.2 | 5.398 | 4.028977 |
| 80 | Statin2 | 80 | 27.2 | 4.517 | 2.348903 |
| 81 | Statin2 | 20 | 27.3 | 3.901 | 2.790047 |
| 82 | Statin1 | 80 | 27.3 | 5.247 | 5.848545 |
| 83 | Statin2 | 80 | 27.4 | 3.507 | 1.247863 |
| 84 | Statin1 | 20 | 27.4 | 3.807 | -1.079928 |
| 85 | Statin2 | 80 | 27.6 | 3.574 | 1.486789 |
| 86 | Statin1 | 40 | 27.8 | 4.16 | 2.427753 |
| 87 | Statin2 | 20 | 28 | 4.501 | 3.284648 |
| 88 | Statin2 | 5 | 28.1 | 3.621 | 2.699007 |
| 89 | Statin1 | 40 | 28.2 | 3.652 | -1.091256 |
| 90 | Statin2 | 40 | 28.2 | 4.191 | 2.874231 |
| 91 | Statin2 | 40 | 28.4 | 5.791 | 4.445454 |
| 92 | Statin1 | 40 | 28.6 | 4.698 | 3.202874 |
| 93 | Statin1 | 5 | 29 | 4.32 | 4.070753 |
| 94 | Statin2 | 10 | 29.1 | 3.776 | 2.751281 |
| 95 | Statin2 | 5 | 29.2 | 4.703 | 3.64949 |
| 96 | Statin2 | 40 | 29.9 | 4.128 | 2.864691 |
| 97 | Statin1 | 40 | 30.4 | 4.693 | 4.983704 |
| 98 | Statin1 | 20 | 30.4 | 4.123 | 2.273898 |
| 99 | Statin1 | 80 | 30.5 | 3.921 | -0.903438 |
| 100 | Statin1 | 10 | 36.5 | 4.175 | 3.311437 |
Chapter 18 contents
12.7 – Many tests one model
Introduction
In our introduction to parametric tests we so far have covered one and two sample t-tests and now the multiple sample or one-way analysis of variance (ANOVA). In subsequent sections we will cover additional tests, each with their own name. It is time to let you in on a little secret. All of these tests, t-tests, ANOVA, and linear and multiple regression we will work on later in the book belong to one family of statistical models. That model is called the general Linear Model (LM), not to be confused with the Generalized Linear Model (GLM) (Burton et al 1998; Guisan et al 2002). This greatly simplifies our approach to learning how to implement statistical tests in R (or other statistical programs) — you only need to learn one approach: the general Linear Model (LM) function lm().
Brief overview of linear models
With the inventions of correlation, linear regression, t-tests, and analysis of variance in the period between 1890 and 1920, subsequent work led to the realization that these tests (and many others!) were special cases of a general model, the general linear model, or LM. The LM itself is a special case of the generalized linear model, or GLM; among the differences between LM and GLM, in LM, the dependent variable is ratio scale and errors in the response (dependent) variable(s) are assumed to come from a Gaussian (normal) distribution. In contrast, for GLM, the response variable may be categorical or continuous, and error distributions other than normal (Gaussian), may be applied. The GLM user must specify both the error distribution family (eg, Gaussian) and the link function, which specifies the relationship among the response and predictor variables. While we will use the GLM functions when we attempt to model growth functions and calculate EC50 in dose-response problems, we will not cover GLM this semester.
The general Linear Model, LM
In matrix form, the LM can be written as
![]()
where
is a matrix of response variables predicted by independent variables contained in matrix
and weighted by linear coefficients in the vector b. Basically, all of the predictor variables are combined to produce a single linear predictor
. By adding an error component we have the complete linear model.
![]()
In the linear model, the error distribution is assumed to be normally distributed, or “Gaussian.”
R code
The bad news is that LM in R (and in any statistical package, actually) is a fairly involved set of commands; the good news is that once you understand how to use this command, and can work with the Options, you will be able to conduct virtually all of the tests we will use this semester, from two-sample t-tests to multiple linear regression. In the end, all you need is the one Rcmdr command to perform all of these tests.
We begin with a data set, ohia.ch12. Scroll down this page or click here to get the R code.
I found a nice report on a common garden experiment with o`hia (Corn and Hiesey 1973). O`hia (Metrosideros polymorpha) is an endemic, but wide-spread tree in the Hawaiian islands (Fig 1). O`hia exhibits pronounced intraspecific variation: individuals differ from each other. O`hia grows over wide range of environments, from low elevations along the ocean right up the sides of the volcanoes, and takes on many different growth forms, from shrubs to trees. Substantial areas of o`hia trees on the Big Island are dying, attributed to two exotic fungal species of the genus Ceratocystis (Asner et al., 2018).

Figure 1. O’hia, Metrosideros polymorpha. Public domain image from Wikipedia.
The Biology. Individuals from distinct populations may differ because the populations differ genetically, or because the environments differ, or, and this is more realistic, both. Phenotypic plasticity is the ability of one genotype to produce more than one phenotype when exposed to different environments. Environmental differences are inevitable when populations are from different geographic areas. Thus, in population comparisons, genetic and environmental differences are confounded. A common garden experiment is a crucial genetic experiment to separate variation in phenotypes, P, among populations into causal genetic or environmental components.
If you recall from your genetics course, ![]()
where G stands for genetic (alleles) differences among individuals and E stands for environmental differences among individuals. In brief, the common garden experiment begins with individuals from the different populations are brought to the same location to control environmental differences. If the individuals sampled from the populations continue to differ despite the common environment, then the original differences between the populations must have a genetic basis, although the actual genetic scenario may be complicated (the short answer is that if “genotype by environment interaction exists, then results from a common garden experiment cannot be generalized back to the natural populations/locations — this will make more sense when we talk about two-way ANOVA). For more about common garden experiments, see de Villemereuil et al (2016). Nuismer and Gandon (2008) discuss statistical aspects of the common garden approach to studying local adaptation of populations and the more powerful “reciprocal translocation” experimental design.
Managing data for linear models
First, your data must be stacked in the worksheet. That means one column is for group labels (independent variable), the other column is for the response (dependent) variable.
If you have not already downloaded the data set, ohia.ch12, do so now. Scroll down this page or click here to get the R code.
Confirm that the worksheet is stacked. If it is not, then you would rearrange your data set using Rcmdr: Data → Active data set → Stack variables in data set…
The data set contains one factor, “Site” with three levels (M-1, 2, 3). M stands for Maui, and collection sites were noted in Figure 2 of Corn and Hiesey (1973). Once the dataset is in Rcmdr, click on View to see the data (Fig 2). There are two response variables, Height (shown in red below) and Width (shown in blue below).

Figure 2. The o`hia dataset as viewed in R Commander.
The data are from Table 5 of Corn and Hiesey 1973. (I simulated data based on their mean/SD reported in Table 5). This was a very cool experiment: they collected o`hia seeds from three elevations on Maui, then grew the seeds in a common garden in Honolulu. Thus, the researchers controlled the environment; what varied, then were the genotypes.
As always, you should look at the data. Box plots are good to compare central tendency (Fig 3).

Figure 3. Box plots of growth responses of o`hia seedlings collected from three Maui sites, M-1 (elevation 750 ft), M-2 (elevation 1100 ft), and M-3 (elevation 6600 ft). Data adapted from Table 5 of Corn and Hiersey 1973.
R code to make Figure 4 plots
par(mfrow=c(1,2)) Boxplot(Height ~ Site, data = ohia, id = list(method = "y")) Boxplot(Width ~ Site, data = ohia, id = list(method = "y"))
This dataset would typically be described as a one-way ANOVA problem. There was one treatment variable (population source) with three levels (M-1, M-2, M-3). From Rcmdr we select the one-way ANOVA: Statistics → Means → One-way ANOVA… and after selecting the Groups (from the Site variable) and the Response variable (eg, Height), we have
AnovaModel.1 <- aov(Height ~ Site, data = ohia.ch12)
summary(AnovaModel.1)
Df Sum Sq Mean Sq F value Pr(>F)
Site 2 4070 2034.8 22.63 0.000000131 ***
Residuals 47 4227 89.9
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Let us proceed to test the null hypothesis (what was it???) using instead the lm() function. Four steps in all.
Step 1. Rcmdr: Statistics → Fit models → Linear model … (Fig 4).

Figure 4. R Commander, select to fit a Linear model.
Step 2. The pop up menu below (Fig 5) follows.
First, What is our response (dependent) variable? What is our predictor (independent) variable? We then input our model. In this case, with only the one predictor variable, Sites, our model formula is simple to enter (Fig 5): Height ~ Site.

Figure 5. Input linear model formula, Height ~ Site.
Step 3. Click OK to carry out the command.
Here is the R output and the statistical results from the application of the linear model.
LinearModel.1 <- lm(Height ~ Site, data=ohia.ch12)
summary(LinearModel.1)
Call:
lm(formula = Height ~ Site, data = ohia.ch12)
Residuals:
Min 1Q Median 3Q Max
-18.808 -4.761 -1.755 4.758 29.257
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.314 2.121 7.222 0.00000000377 ***
Site[T.M-2] 19.261 2.999 6.423 0.00000006153 ***
Site[T.M-3] 2.924 3.673 0.796 0.43
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.483 on 47 degrees of freedom
Multiple R-squared: 0.4905, Adjusted R-squared: 0.4688
F-statistic: 22.63 on 2 and 47 DF, p-value: 0.0000001311
End R output
The linear model has produced a series of estimates of coefficients for the linear model, statistical tests of the significance of each component of the model, and the coefficient of determination, R2, which is a descriptive statistic of how well model fits the data. Instead of our single factor variable for Source Population like in ANOVA we have a series of what are called dummy variables or contrasts between the populations. Thus, there is a coefficient for the difference between M-1 and M-2. “Site[T.M-2]” in the output, between M-1 and M-3, and between M-2 and M-3.
Note 1: Brief description of linear model output; these will be discussed fuller in Chapter 17 and Chapter 18. The residual standard error is a measure of how well a model fits the data. The Adjusted R-squared is calculated by dividing the residual mean square error by the total mean square error. The result is then subtracted from 1.
It also produced our first statistic that assesses how well the model fits the data called the coefficient of determination, R2. An of of 1.0 would indicate that all variation in the data set can be explained by the predictor variable(s) in the model with no residual error remaining. Our value of 49% indicates that nearly 50% of the variation in height of the seedlings grown under common environments are due to the source population (= genetics).
Step 4. But we are not quite there — we want the traditional ANOVA results (recall the ANOVA table).
To get the ANOVA Table we have to ask Rcmdr (and therefore R) to give us this. Select
Rcmdr: Models → Hypothesis tests → ANOVA table … (Fig 6).

Figure 6. To retrieve an ANOVA table, select Models, Hypothesis tests, then ANOVA table…
Here’s the type of tests for the ANOVA table, select the default (Fig 7).

Figure 7. Options for types of tests.
Now, in the future when we work with more complicated experimental designs, we will also need to tell R how to conduct the test. For now, we will accept the default Type II type of test and ignore sandwich estimators. You should confirm that for a one-way ANOVA, Type I and Type II choices give you the same results.
The reason they do is because there is only one factor — when there are more than one factors, and if one or both of the factors are random effects, then Type I, II, and III will give you different answers. We will discuss this more as needed, but see the note below about default choices.
Note 2: Marginal or partial effects are slopes (or first derivatives): they quantify the change in one variable given change in one or more independent variables. Type I tests are sequential: sums of squares are calculated in the order the predictor variables are entered into the model. Type II tests the sums of squares are calculated after adjusting for some of the variables in the model. Type III, every sum of square calculation is adjusted for all other variables in the model. Sandwich estimator refers to algorithms for calculating the structure of errors or residuals remaining after the predictor variables are fitted to the data. The assumption for ordinary least square estimation (See Chapter 17) is that errors across the predictors are equal, ie, equal variances assumption. HC refers to “heteroscedasticity consistent” (Hayes and Chai 2007), which infers fitting of a model does not assume equal error variance across all observations.
By default, Rcmdr makes Type II. In most of the situations we will find ourselves this semester, this is the correct choice.
Below is the output from the ANOVA table request. Confirm that the information is identical to the output from the call to aov() function.
Anova(LinearModel.1, type = "II")
Anova Table (Type II tests)
Response: Height
Sum Sq Df F value Pr(>F)
Site 4069.7 2 22.626 0.0000001311 ***
Residuals 4226.9 47
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
And the other stuff we got from the linear model command? Ignore for now but make note that this is a hint that regression and ANOVA are special cases of the same model, the linear model.
We do have some more work to do with ANOVA, but this is a good start.
Why use the linear model approach?
Chief among the reasons to use the lm() approach is to emphasize that a model approach is in use. One purpose of developing a model is to provide a formula to predict new values. Prediction from linear models is more fully developed in Chapter 17 and Chapter 18, but for now, we introduce the predict() function with our O`hia example.
myModel <- predict(LinearModel.1, interval = "confidence") head(myModel, 3) #print out first 3 rows #Add the output to the data set ohiaPred <- data.frame(ohia,myModel) with(ohiaPred, tapply(fit, list(Site), mean, na.rm = TRUE)) #print out predicted values for the three sites
Output from R
myModel <- predict(LinearModel.1, interval = "confidence")
head(myModel, 3)
fit lwr upr
1 15.31374 11.04775 19.57974
2 15.31374 11.04775 19.57974
3 15.31374 11.04775 19.57974
with(myModel, tapply(fit, list(Site), mean, na.rm = TRUE))
M-1 M-2 M-3
15.31374 34.57474 18.23796
Questions
1. Revisit ANOVA problems in homework and questions from early parts of this chapter and apply lm() followed by Hypothesis testing (Rcmdr: Models → Hypothesis tests → ANOVA table) approach instead of one-way ANOVA command. Compare results using lm() to results from One-way ANOVA and other ANOVA problems.
Quiz Chapter 12.7
Many tests one model
Data set and R code used in this page
Corn and Hiesey (1973)
| Site | Height | Width |
|---|---|---|
| M-1 | 12.5567 | 19.1264 |
| M-1 | 13.2019 | 13.1547 |
| M-1 | 8.0699 | 16.032 |
| M-1 | 6.0952 | 22.8586 |
| M-1 | 11.3879 | 11.0105 |
| M-1 | 12.2242 | 21.8102 |
| M-1 | 16.0147 | 11.0488 |
| M-1 | 19.7403 | 25.9756 |
| M-1 | 36.4824 | 25.2867 |
| M-1 | 13.1233 | 20.0487 |
| M-1 | 21.7725 | 24.8511 |
| M-1 | 14.2013 | 43.7679 |
| M-1 | 37.7629 | 37.3438 |
| M-1 | 2.8652 | 2.5549 |
| M-1 | 0.6456 | 22.8013 |
| M-1 | 29.623 | 20.0194 |
| M-1 | 10.5812 | 29.0328 |
| M-1 | 18.3046 | 22.2867 |
| M-1 | 19.0528 | 24.684 |
| M-1 | 2.5693 | 35.74 |
| M-2 | 45.0162 | 14.3878 |
| M-2 | 40.8404 | 18.8396 |
| M-2 | 27.1032 | 21.0547 |
| M-2 | 29.8036 | 16.9327 |
| M-2 | 63.8316 | 30.7037 |
| M-2 | 42.107 | 3.2491 |
| M-2 | 30.0322 | 47.4412 |
| M-2 | 34.0516 | 42.239 |
| M-2 | 15.7664 | 32.8354 |
| M-2 | 35.1262 | 50.9698 |
| M-2 | 43.6988 | 19.3897 |
| M-2 | 26.7585 | 13.8168 |
| M-2 | 36.7895 | 0.5817 |
| M-2 | 30.9458 | 53.7757 |
| M-2 | 26.8465 | 15.4137 |
| M-2 | 40.3883 | 9.2161 |
| M-2 | 30.6555 | 56.8456 |
| M-2 | 19.9736 | 44.9411 |
| M-2 | 27.676 | 36.8543 |
| M-2 | 44.084 | 24.3396 |
| M-3 | 15.2646 | 11.4999 |
| M-3 | 19.6745 | 9.7757 |
| M-3 | 23.275 | 12.7825 |
| M-3 | 16.1161 | 2.4065 |
| M-3 | 16.8393 | 1.1253 |
| M-3 | 23.107 | 3.7349 |
| M-3 | 21.5322 | 6.9725 |
| M-3 | 13.4191 | 12.2867 |
| M-3 | 14.7273 | 11.4841 |
| M-3 | 18.4245 | 11.9078 |