?airquality2 Linear Regression
NoteIn this chapter
In this chapter, you will learn
- how to fit a linear model, read the regression table, visualize the result, and check the residuals and adjust the functional form if the model does not fit the mean (Section 2.1, Section 2.3, Section 2.4, Section 2.5);
- how categorical predictors are coded (contrasts), and how to follow up with ANOVA, post-hoc tests and compact letter displays (Section 2.6);
- how multiple regression, collinearity, centering, scaling and interactions change what the coefficients mean (Section 2.7, Section 2.8, Section 2.9, Section 2.10);
- how to make predictions and how to report an effect that is not constant, e.g. a curved effect (Section 2.11);
- how to get uncertainty from the bootstrap when the analytical results are not available (Section 2.12).
Workflow steps (Chapter 12): step 3 (a-priori model), step 5 (residual checks), step 6 (inference and predictions).
Prerequisites: basic R and the idea of a regression (Section 1.2).
2.1 Fitting a simple linear regression
Let’s start with the basics. For this chapter, we will rely a lot on the airquality data, which is one of the built-in datasets in R.
TipMore info on the airquality data
Most datasets in R have a help file. Thus, if you don’t know the data set, have a look at the description via pressing F1 with the cursor on the word, or via typing
As you can read, the datasets comprises daily air quality measurements in New York, May to September 1973. You can see the structure of the data via
str(airquality)Variables are described in the help. Note that Month / Day are currently coded as integers, would be better coded as date or (ordered) factors (something you could do later, if you use this dataset, as we do for fMonth in Section 2.7).
Let’s say we want to examine the relationship between Ozone and Wind in these data.
plot(Ozone ~ Wind, data = airquality)
OK, I would say there is some dependency here. To quantify this numerically, we want to fit a linear regression model through the data with the lm() function of R. We can do this in R by typing
fit = lm(Ozone ~ Wind, data = airquality)This command creates a linear regression model that fits a straight line, adjusting slope and intercept to get the best fit through the data. We can see the fitted coefficients if we print the model object
fit
Call:
lm(formula = Ozone ~ Wind, data = airquality)
Coefficients:
(Intercept) Wind
96.873 -5.551
Note
The function name lm is short for “linear model”. Remember from your basic stats course: This model is not called linear because we necessarily fit a straight line. It’s called linear because the response (in our case Ozone) is expressed as a linear combination of the parameters (we call the parameters “effects”). Being linear in the parameters is what ensures that they all affect predictions linearly, so that the estimation can be solved as a system of linear equations. Note that the predictors may still be transformed (e.g. \(x^2\), \(\log x\), or a spline basis) - as long as the effects enter linearly, the model remains a linear model.
So, an lm is any regression of the form \(y = \operatorname{f}(x) + \varepsilon\) where \(\varepsilon \sim \mathcal{N}(0, \sigma^2)\), with \(\operatorname{f}\) linear in the parameters, e.g. \({a}_{0} + {a}_{1} \cdot x + {a}_{2} \cdot {x}^{2}\) (the effects \({a}_{0}, {a}_{1}, {a}_{2}\) enter linearly, even though \(x^2\) is a nonlinear transform of \(x\)). Here, \(\varepsilon \sim \mathcal{N}(0, \sigma^2)\) means that the residuals (deviations from \(\operatorname{f}(x)\)) follow a normal (Gaussian) distribution with unknown variance \(\sigma^2\).
TipUnderstanding how the lm works
To understand how the lm works, we will estimate an lm using a custom likelihood function:
- Create some toy data:
Xobs = rnorm(100, sd = 1.0)
Y = 0.0+1.0*Xobs + rnorm(100,sd = 0.2)- Define the likelihood function and optimize the parameters (slope+intercept):
likelihood = function(par) { # three parameters now
lls = dnorm(Y, mean = Xobs*par[2]+par[1], sd = par[3], log = TRUE)
# calculate for each observation the probability to observe the data given our model
# we use the logLikilihood because of numerical reasons
return(sum(lls))
}
likelihood(c(0, 0, 0.2))[1] -1328.254opt = optim(c(0.0, 0.0, 1.0), fn = function(par) -likelihood(par), hessian = TRUE )
opt$par[1] -0.01770037 1.00544470 0.17984788- Test the estimated parameters against 0 (standard errors can be derived from the hessian of the MLE): Our true parameters are 0.0 for the intercept, 1.0 for the slope, and 0.2 for the sd of the observational error:
- calculate standard error
- calculate t-statistic
- calculate p-value
st_errors = sqrt(diag(solve(opt$hessian)))
st_errors[2][1] 0.01735754t_statistic = opt$par[2] / st_errors[2]
2 * pt(abs(t_statistic), df = 100-2, lower.tail = FALSE)[1] 1.280913e-77
Tip
For a model \(y = \beta_0 + \beta_1 x + \varepsilon\), there are three unknowns: \(\beta_0\), \(\beta_1\), and \(\sigma^2\). However, the degrees of freedom for the t-test of the regression coefficients are \(n\) minus the number of parameters in the mean function (i.e. not counting \(\sigma^2\), which is estimated separately from the residuals). For this model that is \(n - 2 = 100 - 2 = 98\). More generally, for a model with \(p\) predictors plus an intercept, the residual df is \(n - (p + 1)\).
The p-value is smaller than \(\alpha\), so the effect is significant! Let’s compare it against the output of the lm function:
model = lm(Y~Xobs) # 1. Get estimates, MLE
model
Call:
lm(formula = Y ~ Xobs)
Coefficients:
(Intercept) Xobs
-0.01767 1.00543 summary(model) # 2. Calculate standard errors, CI, and p-values
Call:
lm(formula = Y ~ Xobs)
Residuals:
Min 1Q Median 3Q Max
-0.37768 -0.10133 0.00245 0.10821 0.57248
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.01767 0.01818 -0.972 0.333
Xobs 1.00543 0.01753 57.345 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1817 on 98 degrees of freedom
Multiple R-squared: 0.9711, Adjusted R-squared: 0.9708
F-statistic: 3288 on 1 and 98 DF, p-value: < 2.2e-162.2 Interpreting the fitted model
We can get a more detailed summary of the statistical estimates with the summary() function.
summary(fit)
Call:
lm(formula = Ozone ~ Wind, data = airquality)
Residuals:
Min 1Q Median 3Q Max
-51.572 -18.854 -4.868 15.234 90.000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 96.8729 7.2387 13.38 < 2e-16 ***
Wind -5.5509 0.6904 -8.04 9.27e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 26.47 on 114 degrees of freedom
(37 observations deleted due to missingness)
Multiple R-squared: 0.3619, Adjusted R-squared: 0.3563
F-statistic: 64.64 on 1 and 114 DF, p-value: 9.272e-13In interpreting this, recall that:
- “Call” repeats the regression formula.
- “Residuals” gives you an indication about how far the observed data scatters around the fitted regression line / function.
- The regression table (starting with “Coefficients”) provides the estimated parameters, one row for each fitted parameter. The first column is the estimate, the second is the standard error (approximate 95% CI: \(\hat{\beta} \pm 1.96 \cdot \text{SE}\); for small samples use the exact \(t_{n-p}\) quantile instead of 1.96), and the fourth column is the p-value for a two-sided test with \({H}_{0}\): “Estimate is zero”. The t-value is used for calculation of the p-value and can usually be ignored. Note that each coefficient answers one specific question, e.g. for a factor the difference to the reference level, and with interactions or curved effects the effect at one particular point; often this is not the number you want to report (Section 2.11.2).
- The last section of the summary provides information about the model fit.
- Residual error = Standard deviation of the residuals,
- Residual degrees of freedom (here 114) = number of observations used in the fit minus number of fitted mean parameters.
- R-squared \(\left({R}^{2}\right)\) = How much of the signal, respective variance is explained by the model, calculated by \(\displaystyle 1 - \frac{\text{residual variance}}{\text{total variance}}\).
- Adjusted R-squared = Adjusted for model complexity.
- F-test = Test against intercept only model, i.e. is the fitted model significantly better than the intercept only model (most relevant for models with > 1 predictor).
2.3 Visualizing the results
We visualize fitted models with the effects package. This works in the same way for simple and multiple regression (Section 2.7) and later for generalized linear models (GLM, Chapter 9). Effects plots are for looking at a model; for estimates, comparisons and tests of what you see, we will use the emmeans package (Section 2.6.3). The base command to visualize the fitted model is:
library(effects)
plot(allEffects(fit, partial.residuals = T))
The blue line is the fitted model (with confidence interval in light blue). When using the argument partial.residuals = T, purple circles are the data, and the purple line is a nonparametric fit to the data. If the two lines deviate strongly, this indicates that we may have a misspecified model (see the section on residual checks, Section 2.4).
TipSolution
fit = lm(Ozone ~ Temp, data = airquality)
plot(allEffects(fit, partial.residuals = T))
summary(fit)
Call:
lm(formula = Ozone ~ Temp, data = airquality)
Residuals:
Min 1Q Median 3Q Max
-40.729 -17.409 -0.587 11.306 118.271
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -146.9955 18.2872 -8.038 9.37e-13 ***
Temp 2.4287 0.2331 10.418 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 23.71 on 114 degrees of freedom
(37 observations deleted due to missingness)
Multiple R-squared: 0.4877, Adjusted R-squared: 0.4832
F-statistic: 108.5 on 1 and 114 DF, p-value: < 2.2e-16Temperature has a significant positive effect of Ozone. However, effect seems nonlinear, could consider adding a quadratic term of a GAM to improve the fit. The model explains nearly 49% of the variance of the given data (R2). 37 observations have missing data and are omitted. The model explains the data significantly better compared a model with only an intercept (F-statistic).
2.4 Residual checks
The regression line of an lm, including p-values and all that, is estimated under the assumptions that the residuals scatter iid normal; p-values and CIs are calculated under this assumption, and violations of the functional form, independence and constant variance usually matter more than non-normality. If these assumptions are not (approximately) met, the results may be nonsense. Looking at the results of the effect plots in Section 2.3, we can already see a first indication that there may be a problem.
What we see highlighted here is that the data seems to follow a completely different curve than the fitted model. The conclusion here would be: The model we are fitting does not fit to the data, we should not interpret its outputs, but rather say that we reject it, it’s the wrong model, we have to search for a more appropriate description of the data (see the next section on adjusting the functional form, Section 2.5).
2.4.1 LM residual plots
To better analyze these residuals (and potential problems), R offers a function for residual plots. It produces 4 plots. I think it’s most convenient plotting them all into one figure, via the command par(mfrow = c(2, 2)), which produces a figure with 2 x 2 = 4 panels. For the residual checks, we continue with the model Ozone ~ Temp from the exercise in the visualization section (Section 2.3); the same checks work for any lm.
fit = lm(Ozone ~ Temp, data = airquality)
par(mfrow = c(2, 2))
plot(fit)
Interpretation of the panels:
- Residuals vs Fitted: Shows misfits and wrong functional form. Scattering should be uniformly distributed.
- Normal Q-Q: Checks if residuals follow an overall normal distribution. Bullets should lie on the line in the middle of the plot and may scatter a little bit at the ends.
- Scale - Location: Checks for heteroskedasticity. Does the variance change with predictions/changing values? Scattering should be uniformly distributed.
- Residuals vs Leverage: How much impact do outliers have on the regression? Data points with high leverage should not have high residuals and vice versa. Bad points lie in the upper right or in the lower right corner. This is measured via the Cook’s distance. Distances higher than 0.5 indicate candidates for relevant outliers or strange effects.
Note
The plot(model) function should ONLY be used for linear models, because it explicitly tests for the assumptions of normally distributed residuals. Unfortunately, it also works for the results of the glm() function, but when used on a GLM beyond the linear model, it is not interpretable. Use the DHARMa residuals instead (Section 2.4.2, Section 9.3.1).
2.4.2 DHARMa residual plots
Because we will relax the assumption of normal residuals once we go to more complex models, we need a more general check. The DHARMa package provides a general way to perform residual checks for practically any model from the generalized linear mixed-effect model (GLMM) family.
library(DHARMa)
res <- simulateResiduals(fit, plot = T)
The simulateResiduals() function uses simulations to generate standardized residuals for any GLMM. Standardized means that a uniform distribution of residuals signifies perfect fit. Here, the residuals are saved in the variable res.
The argument plot = T creates a plot with two panels (alternatively, you could also type plot(res)). The left panel is a uniform qq plot (calling plotQQunif), and the right panel shows residuals against predicted values (calling plotResiduals), with outliers highlighted in red.
Very briefly, we would expect that a correctly specified model shows:
a straight 1-1 line, as well as n.s. of the displayed tests in the qq-plot (left) -> evidence for the correct overall residual distribution (for more details on the interpretation of this plot, see plotQQunif)
visual homogeneity of residuals in both vertical and horizontal direction, as well as n.s. of quantile tests in the res ~ predictor plot (for more details on the interpretation of this plot, see
?plotResiduals)
Deviations from these expectations can be interpreted similar to a linear regression. See the DHARMa vignette for detailed examples.
2.5 Adjusting the functional form
2.5.1 Mean vs. variance problems
Looking at the residual plots of our regression model, there are two types of problems, and there is a clear order in which you should solve them:
Systematic misfit of the mean: the first thing we should worry about is if the model describes the mean of the data well, i.e. if our function form (we assumed a straight line) fits to the data. You can see misfit for the model
Ozone ~ Tempin a) the effects plot (Section 2.3), b) the res ~ fitted plot (Section 2.4.1), and c) the DHARMa res ~ fitted plot (Section 2.4.2).Distributional problems: only if your model fits the mean of the data well should you look at distributional problems. Distributional problems means that the residuals do not scatter as assumed by the regression (for an lm iid normal). The reason why we look only after solving 1) at the distribution is that a misfit can easily create distributional problems.
Here, and in this part of the book in general, we will first be concerned with adjusting the model to describe the mean correctly. In the third part of the book (starting with Chapter 9; variance modeling in Chapter 10), we will then also turn to advanced options to model the variance.
Important
Residuals are always getting better for more complex models. They should therefore NOT solely be used for automatic model selection, i.e. you shouldn’t assume that the model with the best residuals is necessarily preferable (on how to select models, see the chapter on model selection, Chapter 8). The way you should view residual checks is as a rejection test: the question you are asking is if the residuals are so bad that the model needs to be rejected!
2.5.2 R regression syntax
If we see a misfit of the mean, we should adjust the functional form of the model. Here a few options (see Table B.2 for more details):
fit = lm(Ozone ~ Wind, data = airquality) # Intercept + slope.
fit = lm(Ozone ~ 1, data = airquality) # Only intercept.
fit = lm(Ozone ~ Wind - 1 , data = airquality) # Only slope.
fit = lm(Ozone ~ log(Wind), data = airquality) # Predictor variables can be transformed.
fit = lm(Ozone^0.5 ~ Wind, data = airquality) # Output variables can also be transformed.
fit = lm(Ozone ~ Wind + I(Wind^2), data = airquality) # Mathematical functions with I() command.
fit = lm(Ozone ~ poly(Wind, 2), data = airquality) # Polynomials: the same curve, with poly().Basically, you can use all these tools to search for a better fit through the data. I() works for any transformation of the predictor inside a formula; for polynomials, poly() is usually the better choice (see the box on coding polynomials, Caution 2.1).
Caution 2.1: Coding polynomials:
poly() vs I() syntax
Rule: for a polynomial effect, use poly(x, 2) as the default. Use poly(x, 2, raw = TRUE) or x + I(x^2) only if you need the equation or the turning point of the curve. For a curved effect there is no single slope; how to report it is explained in Section 2.11.2.
dat = airquality[!is.na(airquality$Ozone), ] # rows without missing values (see below)
fitI = lm(Ozone^0.35 ~ Wind + I(Wind^2), data = dat)
fitP = lm(Ozone^0.35 ~ poly(Wind, 2), data = dat)
summary(fitI)$coefficients Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.63092132 0.451187270 14.696605 5.520190e-28
Wind -0.49607080 0.086992904 -5.702428 9.588689e-08
I(Wind^2) 0.01560796 0.003955471 3.945916 1.384441e-04summary(fitP)$coefficients Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.454417 0.0704794 49.013139 5.623049e-78
poly(Wind, 2)1 -6.203325 0.7590864 -8.172094 4.885679e-13
poly(Wind, 2)2 2.995291 0.7590864 3.945916 1.384441e-04all.equal(fitted(fitI), fitted(fitP))[1] TRUE- Same model: both formulas give the same fitted curve, residuals, R², AIC and predictions, and the same test for the quadratic term (last row, t = 3.95). They differ in the linear term and in the meaning of the coefficients.
- The linear row: with
I(), theWindrow refers to the slope at Wind = 0, a value far outside the data, and its test changes if you add the quadratic term or center Wind. Withpoly(), the first row answers the question that you usually have: is there an overall increase or decrease over the observed range? The reason is thatpoly()builds two columns that are uncorrelated with each other, whereas Wind and Wind² are almost the same variable (r = 0.97). - Do not read
poly()coefficients as effects per unit of x. The columns are rescaled, so the numbers have no direct meaning in terms of Wind. A positive quadratic term means a U-shaped curve and a negative one a hump, whatever the linear term; the plot of the curve (plot(allEffects(fitP))) andpredict(fitP, newdata)show the size of the effect. raw = TRUE:poly(x, 2, raw = TRUE)returns the plain powers x and x² instead of orthogonal columns, so it is the same model with the same columns and coefficients asx + I(x^2). Use it if you need the equation or the turning point of the curve; for higher degrees it is shorter than writing every power withI(), e.g.poly(x, 4, raw = TRUE). The coefficients are then directly interpretable, but the terms are strongly correlated again (see the figure below).- Missing values:
poly()builds its columns from all rows of the data beforelm()removes rows with NA. We therefore removed them first; otherwise the columns are only approximately uncorrelated in the rows used (r = 0.06 here), which does not change the model, only the numbers in the linear row slightly. - What to report: the test of the quadratic term, the direction of the curve (U or hump), and the effect plot.
What do the coefficients mean? Both models above describe the same curve (black line in the lower panels), but they build it from different columns. The upper panels show the columns that enter the model; the lower panels show how the fitted curve is the sum of the intercept (dotted) and the two columns, each multiplied with its coefficient (dashed).
dat = airquality[!is.na(airquality$Ozone), ]
fitI = lm(Ozone^0.35 ~ Wind + I(Wind^2), data = dat)
fitP = lm(Ozone^0.35 ~ poly(Wind, 2), data = dat)
w = seq(min(dat$Wind), max(dat$Wind), length.out = 200)
XI = cbind(w, w^2) # columns used by I()
XP = predict(poly(dat$Wind, 2), w) # columns used by poly()
bI = coef(fitI); bP = coef(fitP)
yl = range(c(dat$Ozone^0.35, bI[1], bI[2] * XI[, 1], bI[3] * XI[, 2]))
par(mfrow = c(2, 2), mar = c(4, 4, 3.5, 1), cex.main = 0.95)
# Row 1: the two columns that enter the model
matplot(w, XI, type = "l", lty = 1, lwd = 2, col = c(2, 4), xlab = "Wind", ylab = "value of the column",
main = "I(): columns Wind and Wind^2")
legend("topleft", c("Wind", "Wind^2"), col = c(2, 4), lty = 1, lwd = 2, bty = "n")
matplot(w, XP, type = "l", lty = 1, lwd = 2, col = c(2, 4), xlab = "Wind", ylab = "value of the column",
main = "poly(): columns 1 and 2")
legend("topleft", c("column 1 (linear trend)", "column 2 (curvature)"), col = c(2, 4), lty = 1, lwd = 2, bty = "n")
# Row 2: fitted curve = intercept + coefficient 1 x column 1 + coefficient 2 x column 2
parts = function(X, b, ttl) {
plot(dat$Wind, dat$Ozone^0.35, col = "grey60", pch = 16, cex = 0.6, ylim = yl,
xlab = "Wind", ylab = "Ozone^0.35", main = ttl)
abline(h = b[1], col = "grey40", lty = 3)
lines(w, X[, 1] * b[2], col = 2, lwd = 2, lty = 2)
lines(w, X[, 2] * b[3], col = 4, lwd = 2, lty = 2)
lines(w, b[1] + X %*% b[-1], lwd = 3)
}
parts(XI, bI, sprintf("y = %.2f %+.3f * Wind %+.4f * Wind^2", bI[1], bI[2], bI[3]))
parts(XP, bP, sprintf("y = %.2f %+.2f * column 1 %+.2f * column 2", bP[1], bP[2], bP[3]))
legend("bottomright", c("fitted curve (sum)", "intercept", "coefficient 1 x column 1", "coefficient 2 x column 2"),
col = c(1, "grey40", 2, 4), lty = c(1, 3, 2, 2), lwd = c(3, 1, 2, 2), bty = "n", cex = 0.8)
- With
I()(left), the two columns are almost the same shape, so large positive and negative contributions have to cancel each other to give the curve. This is why the individual coefficients and their tests depend on each other, and why high-degree polynomials built withI()become numerically unstable (for a polynomial of degree 8 on this data, the condition number of the model matrix is about 2e12 withI()and 11 withpoly(); a large condition number means that the estimates are sensitive to rounding errors). - With
poly()(right), column 1 is a straight line centered on the mean of Wind and column 2 is the curvature that remains after the straight line has been taken out. Because they are uncorrelated, each contribution is small and can be judged on its own, and the intercept is simply the mean of the response. - The curve and its fit are identical. Only the way of writing the same curve differs, and therefore the meaning of the coefficients.
Note that if you transform the response variable, you change both the functional form and the spread of the residuals. Thus, this option can change (and often improve) both the mean and distributional problems, and is a common strategy to deal with non-symmetric residuals etc. There is a small convenience function in the package MASS that automatically searches for the best-fitting power transformation of the response.
library(MASS)
fitLin = lm(Ozone ~ Wind, data = airquality)
boxcox(fitLin)
In the result, you can see, that the best fit assuming a power transformation is Ozone^k is delivered by k is about 0.3 (95% interval roughly 0.2-0.4); we use 0.35.
TipSolution
Possible solution, adding a quadratic predictor (with poly(), see Caution 2.1) and choosing a power of 0.35 transformation based on the boxcox function:
fit1 = lm(Ozone^0.35 ~ poly(Wind, 2), data = airquality)
plot(allEffects(fit1, partial.residuals = T), selection = 1)
Note that if you plot(fit1), there seems to be still a bit of a pattern in the scale-location plot (heteroskedasticity). I would say that this is not particularly concerning, but if you are concerned and don’t manage to address it via changing the functional form, you could fit a regression with variable variance. We will discuss this in the section on modeling the variance (Section 10.1.2).
You could get an even better fit by adding more and more predictors, for example:
fit2 = lm(Ozone^0.35 ~ poly(Wind, 8), data = airquality)
plot(allEffects(fit2, partial.residuals = T), selection = 1)
However, as noted in Section 2.5.1, and as we will further discuss in the chapter on model selection (Chapter 8), this is not a good idea, because while a more complicated model always improves residuals, it has other disadvantages. This model is likely too complex for the data (aka it overfits). We can see this by looking at common model selection indicators (again, more in the section on AIC model selection, Section 8.2.2).
AIC comparison (lower = better)
AIC(fit1)[1] 270.2059AIC(fit2)[1] 274.7512Nested model comparison with anova() (for an lm an F-test, which is the likelihood ratio test for normal errors; see Tip 3.1): is there evidence for the more complex model?
anova(fit1, fit2)Analysis of Variance Table
Model 1: Ozone^0.35 ~ poly(Wind, 2)
Model 2: Ozone^0.35 ~ poly(Wind, 8)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 113 65.112
2 107 61.059 6 4.0528 1.1837 0.3205Both AIC (lower for the simpler model) and the test (n.s.) give no support for the more complex model.
2.5.3 Generalized additive models (GAMs)
Another option to get the fit right is Generalized Additive Models (GAM). The idea is to fit a smooth function to data, to automatically find the “right” functional form. The smoothness of the function is automatically optimized.
library(mgcv)
fit = gam(Ozone ~ s(Wind) , data = airquality)
# allEffects doesn't work here.
plot(fit, pages = 0, residuals = T, pch = 20, lwd = 1.8, cex = 0.7,
col = c("black", rep("red", length(fit$residuals))))
In the summary(), you still get significance for the smoothing term (i.e. a p-value, testing the smooth against a flat line), but given that you fit a line that could go up and down, you don’t get an effect direction (how to report such a curved effect is explained in Section 2.11.2). The edf (effective degrees of freedom) tell you how wiggly the curve is: an edf near 1 is a straight line, larger values a wigglier curve.
summary(fit)
Family: gaussian
Link function: identity
Formula:
Ozone ~ s(Wind)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 42.129 2.194 19.2 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(Wind) 2.995 3.763 28.76 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.487 Deviance explained = 50%
GCV = 578.45 Scale est. = 558.53 n = 1162.6 Categorical predictors
The lm() function can handle both numerical and categorical variables. To understand what happens if the predictor is categorical, we’ll use another data set here, PlantGrowth (type ?PlantGrowth or F1 help if you want details). We visualize the data via:
boxplot(weight ~ group, data = PlantGrowth)
2.6.1 Understanding contrasts
Let’s fit an lm() now with the categorical explanatory variable group. The syntax is the same as before
fit = lm(weight ~ group, data = PlantGrowth)
summary(fit)
Call:
lm(formula = weight ~ group, data = PlantGrowth)
Residuals:
Min 1Q Median 3Q Max
-1.0710 -0.4180 -0.0060 0.2627 1.3690
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.0320 0.1971 25.527 <2e-16 ***
grouptrt1 -0.3710 0.2788 -1.331 0.1944
grouptrt2 0.4940 0.2788 1.772 0.0877 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6234 on 27 degrees of freedom
Multiple R-squared: 0.2641, Adjusted R-squared: 0.2096
F-statistic: 4.846 on 2 and 27 DF, p-value: 0.01591but the interpretation of the results often leads to confusion. We see effects for trt1 and trt2, but where is ctrl? The answer is: R is calculating so-called treatment contrasts for categorical predictors. In treatment contrasts, we have a reference level (in our case ctrl), which is the intercept of the model, and then we estimate the difference to the reference for all other factor levels (in our case trt1 and trt2).
Tip
By default, R orders factor levels alphabetically (whether upper or lower case comes first can depend on your computer settings), and the first level is used as reference level in lm(). It is therefore always worth checking the order with levels() before fitting. In our example, ctrl is before trt in the alphabet and thus used as reference, which also makes sense, because it is scientifically logical to compare treatments to the control. In other cases, it may be that the logical reference level does not correspond to the alphabetical order of the factors.
- To change only the reference level, use
relevel(x, ref = "ctrl"). - To set the order of all levels, e.g. for a treatment none / low / high (which R would sort as high, low, none), use
factor(x, levels = c("none", "low", "high")).
When you change the reference level, estimates and p-values will change, because different comparisons (contrasts) are being made, but model predictions will stay the same! Also, avoid choosing a reference level with very few observations: all other estimates are differences to this level, so they will all be imprecise.
This also becomes clear if we plot the results.
plot(allEffects(fit))
Note that in this case, the effects package depicts the CI by purple error bars.
So, in our example, the intercept estimates the mean weight in the ctrl group (with a p-value tested against zero, usually not interesting), and the effects of trt1 and trt2 are the differences of weight compared to the control, with a p-value for a test against an effect size of zero.
Tip 2.1: More details on contrasts
Treatment contrasts are not the only option to deal with categorical variables. A simple alternative is to remove the intercept:
fit = lm(weight ~ 0 + group, data = PlantGrowth)
summary(fit)
Call:
lm(formula = weight ~ 0 + group, data = PlantGrowth)
Residuals:
Min 1Q Median 3Q Max
-1.0710 -0.4180 -0.0060 0.2627 1.3690
Coefficients:
Estimate Std. Error t value Pr(>|t|)
groupctrl 5.0320 0.1971 25.53 <2e-16 ***
grouptrt1 4.6610 0.1971 23.64 <2e-16 ***
grouptrt2 5.5260 0.1971 28.03 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6234 on 27 degrees of freedom
Multiple R-squared: 0.9867, Adjusted R-squared: 0.9852
F-statistic: 665.5 on 3 and 27 DF, p-value: < 2.2e-16This results in fitting the mean for each factor level (cell means), with p-values for a test of each mean against zero. Unfortunately, this syntax option doesn’t generalize to multiple regressions.
More generally, the coding of a factor is controlled by its contrasts. The coding never changes the fitted model (predictions, residuals, R² stay the same); it only changes what the intercept and the coefficients mean. For PlantGrowth (group means ctrl = 5.03, trt1 = 4.66, trt2 = 5.53; average of the means \(\bar\mu\) = 5.07), the common codings give:
| Coding (R) | Intercept | Coefficients | Typical use |
|---|---|---|---|
Treatment (contr.treatment, default) |
\(\mu_1\) = 5.03 | \(\mu_2-\mu_1\), \(\mu_3-\mu_1\) | regression tables, comparison to a control |
Sum-to-zero (contr.sum) |
\(\bar\mu\) = 5.07 | \(\mu_1-\bar\mu\), \(\mu_2-\bar\mu\) | type III ANOVA, effects averaged over levels |
Helmert (contr.helmert) |
\(\bar\mu\) = 5.07 | \((\mu_2-\mu_1)/2\), \((\mu_3-(\mu_1+\mu_2)/2)/3\) | sequential comparisons (note the scaling) |
Polynomial (contr.poly) |
\(\bar\mu\) = 5.07 | linear and quadratic trend (scaled) | ordered levels (ordered() uses this) |
You can set the contrasts for a factor, or locally in the model call:
fit = lm(weight ~ group, data = PlantGrowth,
contrasts = list(group = contr.helmert))
summary(fit)
Call:
lm(formula = weight ~ group, data = PlantGrowth, contrasts = list(group = contr.helmert))
Residuals:
Min 1Q Median 3Q Max
-1.0710 -0.4180 -0.0060 0.2627 1.3690
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.07300 0.11381 44.573 < 2e-16 ***
group1 -0.18550 0.13939 -1.331 0.19439
group2 0.22650 0.08048 2.814 0.00901 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6234 on 27 degrees of freedom
Multiple R-squared: 0.2641, Adjusted R-squared: 0.2096
F-statistic: 4.846 on 2 and 27 DF, p-value: 0.01591Note that the Helmert coefficients are scaled: -0.19 is half the difference between trt1 and ctrl. Polynomial contrasts assume equally spaced levels; if your levels are e.g. doses of 1, 2 and 5, use contr.poly(3, scores = c(1, 2, 5)).
Sum-to-zero vs. orthogonal contrasts. Two properties of contrasts are often confused. Sum-to-zero means that each contrast column sums to zero over the factor levels. Then the intercept is the unweighted average of the group means, and in models with interactions, main effects are averages over the levels of the other factor; this is the property a type III ANOVA requires (see Section 3.2.4). Orthogonal contrasts are sum-to-zero contrasts whose columns are, in addition, uncorrelated with each other. You can see both properties in the cross-product of the coding matrix, including a column for the intercept: zeros in the first row mean sum-to-zero, zeros between the other columns mean orthogonal:
crossprod(cbind(intercept = 1, contr.sum(3))) intercept
intercept 3 0 0
0 2 1
0 1 2crossprod(cbind(intercept = 1, contr.helmert(3))) intercept
intercept 3 0 0
0 2 0
0 0 6Thus, contr.sum sums to zero but is not orthogonal, while Helmert and polynomial contrasts are both (which is why all three have the same intercept in the table above). For a type III ANOVA, sum-to-zero is enough: all sum-to-zero codings give the same results. Orthogonality matters if you want to split the sum of squares of a factor into independent parts (see planned comparisons below), and only translates into independent estimates if the group sizes are equal. Treatment contrasts do not sum to zero, and are therefore no contrasts in this strict sense.
Planned comparisons. If you have specific hypotheses, e.g. “do the treatments differ from the control?” and “do the two treatments differ from each other?”, you can test them as custom contrasts. Don’t write the weights of your hypotheses directly into contrasts(): this only gives the intended tests if the contrasts are mutually orthogonal and the design is balanced. The safer route is the emmeans package, where you specify the weights of each hypothesis on the group means directly:
library(emmeans)
emm = emmeans(lm(weight ~ group, data = PlantGrowth), ~ group)
contrast(emm, list(ctrl_vs_trt = c(-1, 0.5, 0.5),
trt1_vs_trt2 = c(0, -1, 1))) contrast estimate SE df t.ratio p.value
ctrl_vs_trt 0.0615 0.241 27 0.255 0.8009
trt1_vs_trt2 0.8650 0.279 27 3.103 0.0045The two contrasts here are orthogonal (the products of their weights sum to zero) and the design is balanced. In this case, their sums of squares add up exactly to the sum of squares of the factor in the ANOVA, which you can see with the split argument of summary.aov():
PlantGrowth$group2 = PlantGrowth$group
contrasts(PlantGrowth$group2) = cbind(ctrl_vs_trt = c(-2, 1, 1),
trt1_vs_trt2 = c(0, -1, 1))
summary(aov(weight ~ group2, data = PlantGrowth),
split = list(group2 = list(ctrl_vs_trt = 1, trt1_vs_trt2 = 2))) Df Sum Sq Mean Sq F value Pr(>F)
group2 2 3.766 1.883 4.846 0.01591 *
group2: ctrl_vs_trt 1 0.025 0.025 0.065 0.80086
group2: trt1_vs_trt2 1 3.741 3.741 9.627 0.00446 **
Residuals 27 10.492 0.389
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1With unequal group sizes, the parts still add up (the split is sequential), but the result now depends on the order of the contrasts. Note also that orthogonal contrasts do not protect against multiple testing. Planned comparisons should be defined before seeing the data, and you should decide in advance if and how you correct for multiple testing.
For more details, see Schad et al. (2020), chapter 10 in the book by Berg (2025) and also Schielzeth (2010).
Berg, S.M. van den (2025) Analysing data using linear models: With applications in R. bookdown.org. Available at: https://bookdown.org/pingapang9/linear_models_bookdown/ (Accessed: October 3, 2026).
Schad, D.J. et al. (2020) “How to capitalize on a priori contrasts in linear (mixed) models: A tutorial,” Journal of Memory and Language, 110, p. 104038.
How the coding affects the ANOVA, in particular the type III ANOVA in models with interactions, is discussed in the section on contrasts and ANOVA results (Section 3.2.4).
2.6.2 ANOVA
For categorical predictors, one nearly always performs an ANOVA on top of the linear regression estimates. The underlying ideas behind this will be explained in more detail in the section on the idea of ANOVA (Section 3.1), but just very quickly:
An ANOVA starts with a base model (in this case intercept only) and adds the variable group. It then measures:
- How much the model improves in terms of \({R}^{2}\) (this is in the column Sum Sq).
- If this increase of model fit is significant.
We can run an ANOVA on the fitted model object using aov(). For a linear model, summary(aov(fit)) and anova(fit) give the same (sequential, type I) ANOVA table. In general, we will use anova(), R’s general function for comparing models, which also works for GLMs and mixed models; we use aov() here only because the post-hoc test TukeyHSD() used below requires an aov object as input:
fit = lm(weight ~ group, data = PlantGrowth)
anov = aov(fit)
summary(anov) Df Sum Sq Mean Sq F value Pr(>F)
group 2 3.766 1.8832 4.846 0.0159 *
Residuals 27 10.492 0.3886
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In this case, we would conclude that the variable group (3 levels) significantly improves model fit, i.e. the group overall has a significant effect, even though the individual contrasts in the original model were not significant.
Note
The situation that the ANOVA is significant, but the linear regression contrasts are not, is often perceived as confusing, but this is perfectly normal and not a contradiction.
The ANOVA tests a different hypothesis (at least one group mean differs) than the individual coefficients (each treatment vs the reference). In the summary of a lm(), for >2 levels, not all possible contrasts are tested: in our example, we estimate p-values for ctrl-trt1 and ctrl-trt2, which are both n.s., but we do not test for differences of trt1-trt2 (which happens to be significant). The largest difference is between trt1 and trt2; therefore the omnibus test can be significant although both treatment-vs-control contrasts are not.
Thus, the general idea is: 1. ANOVA tells us if there is an effect of the variable at all, 2. summary() tests for specific group differences.
2.6.3 Post-hoc tests
Another common test performed on top of an ANOVA or a categorical predictor of an lm() is a so-called post-hoc test. Basically, a post-hoc test computes p-values for all possible combinations of factor levels, usually corrected for multiple testing:
TukeyHSD(anov) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = fit)
$group
diff lwr upr p adj
trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
trt2-ctrl 0.494 -0.1972161 1.1852161 0.1979960
trt2-trt1 0.865 0.1737839 1.5562161 0.0120064In the result, we see, as hinted at in Section 2.6.2, that there is a significant difference between trt1 and trt2. Note that the p-values of the post-hoc tests are larger than those of the lm summary(). This is because post-hoc tests nearly always correct p-values for multiple testing, while regression estimates are usually not corrected for multiple testing (you could of course if you wanted).
Note
Remember: a significance test with a significance threshold of alpha = 0.05 has a type I error rate of 5%. Thus, if you make 20 independent tests, you would expect about one significant result by chance (probability of at least one: 1 - 0.95^20, about 64%), even if there are no effects whatsoever. Corrections for multiple testing are used to compensate for this, usually aiming at controlling the family-wise error rate (FWER), the probability of making one or more type I errors when performing multiple hypotheses tests. For any set of p-values, p.adjust(p, method = "holm") (FWER) or p.adjust(p, method = "BH") (false discovery rate) applies such a correction; the family of tests is defined by the question the analysis answers.
TukeyHSD() only works for simple ANOVA models fitted with aov(). A more general alternative is the emmeans package, which works for lm, glm and mixed models, and also for models with several predictors. emmeans() calculates the estimated marginal mean for each factor level: the model’s prediction for that level, with all other predictors in the model held fixed (numeric predictors at their mean, other factors averaged with equal weight; “marginal” here only means averaged over the other predictors). It is a model-based mean, which equals the raw group mean only in simple balanced designs like this one. pairs() compares all levels, with Tukey adjustment as default:
library(emmeans)
fit = lm(weight ~ group, data = PlantGrowth)
emm = emmeans(fit, ~ group)
emm group emmean SE df lower.CL upper.CL
ctrl 5.03 0.197 27 4.63 5.44
trt1 4.66 0.197 27 4.26 5.07
trt2 5.53 0.197 27 5.12 5.93
Confidence level used: 0.95 pairs(emm) contrast estimate SE df t.ratio p.value
ctrl - trt1 0.371 0.279 27 1.331 0.3909
ctrl - trt2 -0.494 0.279 27 -1.772 0.1980
trt1 - trt2 -0.865 0.279 27 -3.103 0.0120
P value adjustment: tukey method for comparing a family of 3 estimates If you are only interested in comparing each treatment against a control, you make fewer comparisons, and thus have more power. This is the so-called Dunnett test:
contrast(emm, "trt.vs.ctrl", ref = "ctrl") contrast estimate SE df t.ratio p.value
trt1 - ctrl -0.371 0.279 27 -1.331 0.3296
trt2 - ctrl 0.494 0.279 27 1.772 0.1582
P value adjustment: dunnettx method for 2 tests Make sure to name the control with ref; otherwise, the first factor level is used as control. The default adjustment for this test (dunnettx) is an approximation; for exact p-values, add adjust = "mvt". We will use emmeans more extensively in the chapter on ANOVA (Chapter 3; see in particular Note 3.1).
2.6.4 Compact letter display
It is common to visualize the results of the post-hoc tests with the so-called Compact Letter Display (cld). This doesn’t work with the base TukeyHSD function, so we will use the multcomp package:
library(multcomp)
fit = lm(weight ~ group, data = PlantGrowth)
tuk = glht(fit, linfct = mcp(group = "Tukey"))
summary(tuk) # Standard display.
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lm(formula = weight ~ group, data = PlantGrowth)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
trt1 - ctrl == 0 -0.3710 0.2788 -1.331 0.3909
trt2 - ctrl == 0 0.4940 0.2788 1.772 0.1979
trt2 - trt1 == 0 0.8650 0.2788 3.103 0.0119 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)tuk.cld = cld(tuk) # Letter-based display.
plot(tuk.cld)The cld gives a new letter for each group of factor levels that are statistically undistinguishable. You can see the output via tuk.cld, here I only show the plot:
2.7 Specifying a multiple regression
Multiple (linear) regression means that we consider more than 1 predictor in the same model. The syntax is very easy: Just add your predictors (numerical or categorical) to your regression formula, as in the following example for the airquality dataset. To be able to also add a factor, I created a new variable fMonth to have month as a factor (categorical):
airquality$fMonth = factor(airquality$Month)
fit = lm(Ozone ~ Temp + Wind + Solar.R + fMonth, data = airquality)In principle, everything is interpreted as before: First of all, residual plots:
par(mfrow =c(2,2))
plot(fit)
Tip
Maybe you have noted that the residual plots used fitted values as the x axis. This is in general the most parsimonious choice, as there are many predictors in a multiple regression, and each influences the response. However, if running a multiple regression, you should additionally plot residuals against each predictor separately - if doing so, you often see a misfit that doesn’t occur in the res ~ predicted plot.
The resulting regression table already looks a bit intimidating,
summary(fit)
Call:
lm(formula = Ozone ~ Temp + Wind + Solar.R + fMonth, data = airquality)
Residuals:
Min 1Q Median 3Q Max
-40.344 -13.495 -3.165 10.399 92.689
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -74.23481 26.10184 -2.844 0.00537 **
Temp 1.87511 0.34073 5.503 2.74e-07 ***
Wind -3.10872 0.66009 -4.710 7.78e-06 ***
Solar.R 0.05222 0.02367 2.206 0.02957 *
fMonth6 -14.75895 9.12269 -1.618 0.10876
fMonth7 -8.74861 7.82906 -1.117 0.26640
fMonth8 -4.19654 8.14693 -0.515 0.60758
fMonth9 -15.96728 6.65561 -2.399 0.01823 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.72 on 103 degrees of freedom
(42 observations deleted due to missingness)
Multiple R-squared: 0.6369, Adjusted R-squared: 0.6122
F-statistic: 25.81 on 7 and 103 DF, p-value: < 2.2e-16Luckily, we also have the effect plots to make sense of this:
library(effects)
plot(allEffects(fit, partial.residuals = T) )
Tip
Note that for the multiple regression, that scatter around the regression line displayed by the effects package is not the raw data, but the so-called partial residuals. Essentially, what this means is that we plot the fitted line for the respective predictor, and then we calculate the residual for the full model (including the other predictors) and add them on the regression line. The advantage is that this subtracts the effect of other predictors, so the scatter shows the part of the response not explained by the other predictors, i.e. the conditional effect of this predictor (see Section 2.8 for collinearity).
| Formula | Meaning | Details |
|---|---|---|
y~x_1 |
\(y=a_0 +a_1*x_1\) | Slope+Intercept |
y~x_1 - 1 |
\(y=a_1*x_1\) | Slope, no intercept |
y~x_1+x_2 |
\(y=a_0+a_1*x_1+a_2*x_2\) | Multiple linear regression (two variables) |
y~x_1:x_2 |
\(y=a_0+a_1*(x_1*x_2)\) | Interaction between x1 and x2 |
y~x_1*x_2 |
\(y=a_0+a_1*x_1+a_2*x_2+a_3*(x_1*x_2)\) | Interaction and main effects |
2.8 The effect of collinearity
A common misunderstanding is that a multiple regression is just a convenient way to run several independent univariate regressions. This, however, is wrong. A multiple regression is doing a fundamentally different thing. Let’s compare the results of the multiple regression in Section 2.7 to the simple regression.
fit = lm(Ozone ~ Wind , data = airquality)
summary(fit)
Call:
lm(formula = Ozone ~ Wind, data = airquality)
Residuals:
Min 1Q Median 3Q Max
-51.572 -18.854 -4.868 15.234 90.000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 96.8729 7.2387 13.38 < 2e-16 ***
Wind -5.5509 0.6904 -8.04 9.27e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 26.47 on 114 degrees of freedom
(37 observations deleted due to missingness)
Multiple R-squared: 0.3619, Adjusted R-squared: 0.3563
F-statistic: 64.64 on 1 and 114 DF, p-value: 9.272e-13The simple regression estimates an effect of Wind that is -5.55, while in the multiple regression, we had -3.1.
The reason is that a multiple regression estimates the effect of Wind, corrected for the effect of the other variables, while a univariate regression estimates the raw correlation of Ozone and Wind. And this can make a big difference if Wind has correlations with other variables(the technical term is collinear). We can see, e.g., that Wind and Temp are indeed correlated by running
plot(Wind ~ Temp, data = airquality)
How to quantify collinearity? For two predictors, their correlation (cor()) is sufficient (|r| > 0.7 is a common warning sign). With more predictors, variance inflation factors (car::vif(fit) on the fitted model) are more informative: they also detect a predictor that is well explained by a combination of the others, and √VIF is the factor by which its standard error is inflated. As a rough guide, values above 5–10 are often considered problematic.
In such a case, if Temp is not included in the regression, the effect of Temp will be absorbed in the Wind effect to the extent to which both variables are collinear. In other words, if one variable is missing, there will be a spillover of the effect of this variable to all collinear variables. This spillover is omitted-variable bias, the core of confounding; whether adding a variable to the regression is the right thing to do depends on the causal structure (colliders, mediators), see Chapter 7.
Tip
To give a simple example of this spillover, let’s say
- Wind has a positive effect of 1
- Temp has a negative effect of -1
- Wind and Temp are correlated 50%
A multiple regression would estimate 1,-1. In a single regression of the response on Wind alone (with Wind and Temp standardized and correlated at 0.5), the collinear part of the Temp effect would be absorbed by Wind, so we would estimate a Wind effect of 1 + (-1)(0.5) = 0.5.
Note that the same happens in univariate plots of the correlation, i.e. this is not a problem of the linear regression, but rather of looking at the univariate correlation of Ozone ~ Wind without correcting for the effect of Temp.
Also, note that collinearity is not fundamentally a problem for the regression - regression estimates stay unbiased and p-values and CIs calibrated under arbitrary high collinearity. However, the more collinearity we have, the more uncertain regression estimates will become (with associated large p-values). Here an example, where we simulate data under the assumption that x1 and x2 have each an effect of 1
set.seed(123)
x1 = runif(100)
x2 = 0.95 *x1 + 0.05 * runif(100)
y = x1 + x2 + rnorm(100)
summary(lm(y ~ x1 + x2))
Call:
lm(formula = y ~ x1 + x2)
Residuals:
Min 1Q Median 3Q Max
-1.8994 -0.6821 -0.1086 0.5749 3.3663
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.1408 0.2862 -0.492 0.624
x1 1.6605 7.0742 0.235 0.815
x2 0.4072 7.4695 0.055 0.957
Residual standard error: 0.9765 on 97 degrees of freedom
Multiple R-squared: 0.2668, Adjusted R-squared: 0.2517
F-statistic: 17.65 on 2 and 97 DF, p-value: 2.91e-07We see that the multiple regression is not terribly wrong, but very uncertain about the true values. Some authors therefore advise to remove collinear predictors (for an overview of methods to deal with collinearity, see Dormann et al., 2013). This can sometimes be useful, in particular for predictions (see the bias-variance trade-off, Section 8.1); however, note that the effect of the removed predictor will then be absorbed by the remaining predictor(s). In our case, the univariate regression would be very sure that there is an effect, note that the CI does not cover the true effect of x1 (which was 1).
Dormann, C.F. et al. (2013) “Collinearity: A review of methods to deal with it and a simulation study evaluating their performance,” Ecography, 36(1), pp. 27–46.
summary(lm(y ~ x1))
Call:
lm(formula = y ~ x1)
Residuals:
Min 1Q Median 3Q Max
-1.9000 -0.6757 -0.1067 0.5799 3.3630
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.1295 0.1965 -0.659 0.511
x1 2.0456 0.3426 5.971 3.78e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9715 on 98 degrees of freedom
Multiple R-squared: 0.2668, Adjusted R-squared: 0.2593
F-statistic: 35.65 on 1 and 98 DF, p-value: 3.785e-08
TipUnderstanding collinearity spillover in more detail
The simulation above showed that a univariate regression absorbs the effect of a collinear predictor that is left out. We can look at this spillover in more detail, and visually, with two collinear predictors that both have an effect of 1 on a normal response \(y\) (intercept 0). First, the case where the predictors are positively correlated, i.e. their effects have equal signs:
x1 = runif(100, -5, 5)
x2 = x1 + 0.2*runif(100, -5, 5)
cor(x1, x2)[1] 0.9810774y = 0 + 1*x1 + 1*x2 + rnorm(100)If you look at the formula, you can nearly visually see what the problem of the univariate regression is: as x2 is nearly identical to x1, the apparent slope between y and either x1 or x2 will be around 2 (although the true effect is 1). Consequently, the univariate models will estimate too high effect sizes, while the multiple regression gets the right estimates (with a bit of error):
coef(lm(y ~ x1))(Intercept) x1
0.002362246 2.003283129 coef(lm(y ~ x2))(Intercept) x2
-0.08366353 1.91355593 coef(lm(y ~ x1 + x2))(Intercept) x1 x2
-0.03627352 1.05770683 0.92146207 You can also see this visually. As we need the same plot again below, we write it as a small function:
plotSpillover = function(x1, x2, y) {
par(mfrow = c(1, 2))
plot(x1, y, main = "x1 effect", ylim = c(-12, 12))
abline(lm(y ~ x1))
# Draw a line with intercept 0 and slope 1,
# just like we simulated the true dependency of y on x1:
abline(0, 1, col = "red")
legend("topleft", c("fitted", "true"), lwd = 1, col = c("black", "red"))
plot(x2, y, main = "x2 effect", ylim = c(-12, 12))
abline(lm(y ~ x2))
abline(0, 1, col = "red")
legend("topleft", c("fitted", "true"), lwd = 1, col = c("black", "red"))
}
plotSpillover(x1, x2, y)
Check what happens if the two predictors are negatively correlated (both effects are still 1).
x1 = runif(100, -5, 5)
x2 = -x1 + 0.2*runif(100, -5, 5)
y = 0 + 1*x1 + 1*x2 + rnorm(100)
cor(x1, x2)[1] -0.9804515coef(lm(y ~ x1))(Intercept) x1
-0.03467159 -0.01170653 coef(lm(y ~ x2))(Intercept) x2
-0.04036562 0.04795055 coef(lm(y ~ x1 + x2))(Intercept) x1 x2
-0.04622019 0.91103840 0.94186191 plotSpillover(x1, x2, y)
In the univariate regressions, the two effects cancel out (slopes close to 0), while the multiple regression recovers both effects.
2.9 Centering and scaling of predictor variables
2.9.1 Centering
In our multiple regression (Section 2.7), we saw an intercept of -74. Per definition, the intercept is the predicted value for \(y\) (Ozone) at a value of 0 for all other variables. It’s fine to report this, as long as we are interested in this value, but often the value at predictor = 0 is not particularly interesting. In this specific case, a value of -74 clearly doesn’t make sense, as Ozone concentrations can only be positive.
To explain what’s going on, let’s look at the univariate regression against Temp:
fit = lm(Ozone ~ Temp, data = airquality)
summary(fit)
Call:
lm(formula = Ozone ~ Temp, data = airquality)
Residuals:
Min 1Q Median 3Q Max
-40.729 -17.409 -0.587 11.306 118.271
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -146.9955 18.2872 -8.038 9.37e-13 ***
Temp 2.4287 0.2331 10.418 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 23.71 on 114 degrees of freedom
(37 observations deleted due to missingness)
Multiple R-squared: 0.4877, Adjusted R-squared: 0.4832
F-statistic: 108.5 on 1 and 114 DF, p-value: < 2.2e-16Now, the intercept is -147. We can see the reason more clearly when we plot the results:
plot(Ozone ~ Temp, data = airquality, xlim = c(-10, 110), ylim = c(-200, 170))
abline(fit)
abline(h = 0, lty = 2)
abline(v = 0, lty = 2)
That shows us that the value 0 is far outside of the set of our observed values for Temp, which is measured in Fahrenheit. Thus, we are extrapolating the Ozone far beyond the observed data. To avoid this, we can simply re-define the x-Axis, by subtracting the mean, which is called centering:
airquality$cTemp = airquality$Temp - mean(airquality$Temp)Alternatively, you can center with the built-in R command scale
airquality$cTemp = scale(airquality$Temp, center = T, scale = F)
Tip
scale() returns a matrix that remembers the transformation
scale() returns a one-column matrix that stores the values used for the transformation as attributes: attr(airquality$cTemp, "scaled:center") gives the mean that was subtracted (77.9), and if you also scale (see Section 2.9.2), "scaled:scale" gives the standard deviation that was divided by. This is useful to translate results back to the original units, e.g. the estimate for a scaled predictor divided by attr(x, "scaled:scale") is the effect per original unit. The downside is that some functions expect plain numeric columns. For example, predict() with newdata stops with the error ‘variable … was fitted with type “nmatrix.1” but type “numeric” was supplied’ if you supply plain numbers for such a variable. In this case, transform the new values in the same way, e.g. scale(newdata$Temp, center = attr(airquality$cTemp, "scaled:center"), scale = FALSE) (for a scaled variable, also pass scale = attr(..., "scaled:scale")). Some other functions (e.g. from the effects package) convert such columns to plain vectors and print a warning about it, which is harmless.
Fitting the model with the centered variable moves the intercept line in the middle of the predictor range
fit = lm(Ozone ~ cTemp, data = airquality)
plot(Ozone ~ cTemp, data = airquality)
abline(fit)
abline(v = 0, lty = 2)
which produces a more interpretable value for the intercept: when we center, the intercept of the centered variable can be interpreted as the Ozone concentration at the mean temperature. For a univariate regression, this value will also typically be very similar to the grand mean mean(airquality$Ozone, na.rm = TRUE) (42.1).
summary(fit)
Call:
lm(formula = Ozone ~ cTemp, data = airquality)
Residuals:
Min 1Q Median 3Q Max
-40.729 -17.409 -0.587 11.306 118.271
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 42.1576 2.2018 19.15 <2e-16 ***
cTemp 2.4287 0.2331 10.42 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 23.71 on 114 degrees of freedom
(37 observations deleted due to missingness)
Multiple R-squared: 0.4877, Adjusted R-squared: 0.4832
F-statistic: 108.5 on 1 and 114 DF, p-value: < 2.2e-16
Caution
For categorical predictors, the analog of centering is sum-to-zero coding (contr.sum, see Section 2.6.1): the intercept is then the average over the levels. With the default treatment contrasts, the intercept is predicted for the reference group of all categorical variables, evaluated at the mean of all numeric variables (if centered). Random effects also estimate a grand mean over all groups (see Section 4.3.1).
2.9.2 Scaling
Another very common transformation is to divide the x axis by a certain value to bring all variables on a similar scale. This is called scaling.
A frequent motivation for scaling (which is optional, see the summary In short at the end of this section) is to make effect sizes in a multiple regression more comparable. As an example, look again at the multiple regression with 4 predictors (Section 2.7):
airquality$fMonth = factor(airquality$Month)
fit = lm(Ozone ~ Temp + Wind + Solar.R + fMonth, data = airquality)
summary(fit)
Call:
lm(formula = Ozone ~ Temp + Wind + Solar.R + fMonth, data = airquality)
Residuals:
Min 1Q Median 3Q Max
-40.344 -13.495 -3.165 10.399 92.689
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -74.23481 26.10184 -2.844 0.00537 **
Temp 1.87511 0.34073 5.503 2.74e-07 ***
Wind -3.10872 0.66009 -4.710 7.78e-06 ***
Solar.R 0.05222 0.02367 2.206 0.02957 *
fMonth6 -14.75895 9.12269 -1.618 0.10876
fMonth7 -8.74861 7.82906 -1.117 0.26640
fMonth8 -4.19654 8.14693 -0.515 0.60758
fMonth9 -15.96728 6.65561 -2.399 0.01823 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.72 on 103 degrees of freedom
(42 observations deleted due to missingness)
Multiple R-squared: 0.6369, Adjusted R-squared: 0.6122
F-statistic: 25.81 on 7 and 103 DF, p-value: < 2.2e-16So, which of the predictors is the strongest (= largest effect on the response)? Superficially, it looks as if Month has the highest values. But we have to remember that the effect on the response is \(y = \text{regression estimate} * \text{predictor}\), i.e if we have a predictor with a large range (difference between min/max values), it may have a strong effect even though the estimate is small. So, if the predictors have a different range, we cannot compare effect sizes directly regarding their “effective influence” on y.
To make the effect sizes comparable, we can scale all numeric predictors by dividing them by their standard deviation, which will bring them all roughly on a range between -2, 2. You can do this by hand, or use the scale() function in R:
By default, scale() centers and scales; scale(x, scale = FALSE) only centers.
airquality$sTemp = scale(airquality$Temp)
airquality$sWind = scale(airquality$Wind)
airquality$sSolar.R = scale(airquality$Solar.R)These are one-column matrices that keep their mean and standard deviation as attributes (see the tip box in Section 2.9.1).
Tip
Here, I create a new variable for each scaled predictor. It is possible to use the scale() command inside the regression formula as well, but that sometimes creates problems with downstream plotting functions, so I recommend to create a new variable in your data for scaling.
Running the regression:
fit = lm(Ozone ~ sTemp + sWind + sSolar.R + fMonth, data = airquality)
summary(fit)
Call:
lm(formula = Ozone ~ sTemp + sWind + sSolar.R + fMonth, data = airquality)
Residuals:
Min 1Q Median 3Q Max
-40.344 -13.495 -3.165 10.399 92.689
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 50.558 5.526 9.150 5.74e-15 ***
sTemp 17.748 3.225 5.503 2.74e-07 ***
sWind -10.952 2.325 -4.710 7.78e-06 ***
sSolar.R 4.703 2.131 2.206 0.0296 *
fMonth6 -14.759 9.123 -1.618 0.1088
fMonth7 -8.749 7.829 -1.117 0.2664
fMonth8 -4.197 8.147 -0.515 0.6076
fMonth9 -15.967 6.656 -2.399 0.0182 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.72 on 103 degrees of freedom
(42 observations deleted due to missingness)
Multiple R-squared: 0.6369, Adjusted R-squared: 0.6122
F-statistic: 25.81 on 7 and 103 DF, p-value: < 2.2e-16As you see, effect sizes have now changed considerably. The difference in interpretation is the following: for the unscaled variables, we estimate the effect of 1 unit change (e.g. 1 degree for temp). For the scaled variable, we estimate the effect of a change of 1 standard deviation. If we look at the sd of temperature
sd(airquality$Temp)[1] 9.46527we can conclude that the scaled regression is estimating an effect of a 9.46527 degree change in temperature.
While the latter is numerically a bit confusing, the advantage is that so we can interpret the Ozone effect scaled to typical temperature differences in the data, and the same for all other variables. Consequently, it makes sense to compare the scaled effect sizes between variables, which suggests, judged by the point estimates, that Temp is the most important predictor (check the CIs before ranking). “Most important” here means the largest effect of a change by one SD, which depends on the variation of the predictor in your sample.
What about categorical variables? They are not scaled: their coefficient is the difference between two levels, which is easy to interpret as it is, but it is not an effect per SD. An effect per SD of a numeric predictor and a difference between two levels are therefore not directly comparable. To make them more comparable, Gelman (2008) suggests dividing numeric predictors by 2 SD instead of 1 SD.
Gelman, A. (2008) “Scaling regression inputs by dividing by two standard deviations,” Statistics in Medicine, 27(15), pp. 2865–2873.
In short:
- Centering: centering numeric predictors at their mean is almost always useful: the intercept then refers to an average observation, and in models with interactions, the main effects refer to the other predictors at their center (Section 2.10.3). The exception is when you want the intercept (and the main effects) to refer to a particular reference value: then center at this value instead of the mean (e.g. a reference temperature), or do not center at all if x = 0 is the meaningful reference (e.g. dose 0 = untreated control). For strongly skewed predictors, the median is often a better center than the mean.
- Scaling is optional. It often makes the results easier to interpret, because the coefficients become effects per SD of the predictor, which can be compared between numeric predictors, and for more complex models such as GLMMs, it often helps (and is sometimes needed) for the model to converge (Section 4.4.1). The downside is that the effects are no longer in natural units (per degree, per mm of rain); they can be transformed back (see the tip in Section 2.9.1), or you report both versions.
- In practice: create the centered or scaled predictors as new variables (e.g.
sTempfor scaled,cTempfor centered only; see the tip above on why not inside the formula), and say in the methods which version you report.
2.10 Interactions
2.10.1 Syntax
When we have multiple variables, we can have the situation that the value of one variable influences the effect of the other(s). Technically, this is called an interaction. In situations where the causal direction is known, this is also called a moderator. An example: Imagine we observe that the effect of aspirin differs depending on the weight of the subject. Technically, we have an interaction between aspirin and weight. Physiologically, we know the causal direction is “weight -> effect of aspirin”, so we can say weight is a moderator for the effect of aspirin.
fit = lm(Ozone ~ Temp * Wind, data = airquality)
plot(allEffects(fit))
We will have a look at the summary later, but for the moment, let’s just look at the output visually. In the effect plots, we see the effect of Temperature on Ozone for different values of Wind. We also see that the slope changes. For low Wind, we have a strong effect of Temperature. For high Wind, the effect is basically gone.
Let’s look at the interaction syntax in more detail. The “*” operator in an lm() is a shorthand for main effects + interactions. You can write equivalently:
fit = lm(Ozone ~ Wind + Temp + Wind:Temp, data = airquality)What is fit here is literally a third predictor that is specified as Wind * Temp (normal multiplication). The above syntax would allow you to also have interactions without main effects, e.g.:
fit = lm(Ozone ~ Wind + Wind:Temp, data = airquality)Although this is generally never advisable, as the main effect influences the interaction, unless you are sure that the main effect must be zero.
There is another important syntax in R:
fit = lm(Ozone ~ (Wind + Temp + Solar.R)^2 , data = airquality)
summary(fit)
Call:
lm(formula = Ozone ~ (Wind + Temp + Solar.R)^2, data = airquality)
Residuals:
Min 1Q Median 3Q Max
-38.685 -11.727 -2.169 7.360 91.244
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.408e+02 6.419e+01 -2.193 0.03056 *
Wind 1.055e+01 4.290e+00 2.460 0.01555 *
Temp 2.322e+00 8.330e-01 2.788 0.00631 **
Solar.R -2.260e-01 2.107e-01 -1.073 0.28591
Wind:Temp -1.613e-01 5.896e-02 -2.735 0.00733 **
Wind:Solar.R -7.231e-03 6.688e-03 -1.081 0.28212
Temp:Solar.R 5.061e-03 2.445e-03 2.070 0.04089 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 19.17 on 104 degrees of freedom
(42 observations deleted due to missingness)
Multiple R-squared: 0.6863, Adjusted R-squared: 0.6682
F-statistic: 37.93 on 6 and 104 DF, p-value: < 2.2e-16plot(allEffects(fit), selection = 1)
plot(allEffects(fit), selection = 2)
plot(allEffects(fit), selection = 3)
This creates all main effect and second order (aka two-way) interactions between variables. You can also use ^3 to create all possible 2-way and 3-way interactions between the variables in the parentheses. By the way: the ()^2 syntax for interactions is the reason why powers of a single predictor need I(x^2) or poly(x, 2) in the formula (see Caution 2.1).
2.10.2 Interactions with categorical variables
When you include an interaction with a categorical variable, that means a separate effect will be fit for each level of the categorical variable, as in
fit = lm(Ozone ~ Wind * fMonth, data = airquality)
summary(fit)
Call:
lm(formula = Ozone ~ Wind * fMonth, data = airquality)
Residuals:
Min 1Q Median 3Q Max
-54.528 -12.562 -2.246 10.691 77.750
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 50.748 15.748 3.223 0.00169 **
Wind -2.368 1.316 -1.799 0.07484 .
fMonth6 -41.793 31.148 -1.342 0.18253
fMonth7 68.296 20.995 3.253 0.00153 **
fMonth8 82.211 20.314 4.047 9.88e-05 ***
fMonth9 23.439 20.663 1.134 0.25919
Wind:fMonth6 4.051 2.490 1.627 0.10680
Wind:fMonth7 -4.663 2.026 -2.302 0.02329 *
Wind:fMonth8 -6.154 1.923 -3.201 0.00181 **
Wind:fMonth9 -1.874 1.820 -1.029 0.30569
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 23.12 on 106 degrees of freedom
(37 observations deleted due to missingness)
Multiple R-squared: 0.5473, Adjusted R-squared: 0.5089
F-statistic: 14.24 on 9 and 106 DF, p-value: 7.879e-15The interpretation is like for a single categorical predictor, i.e. we see the effect of Wind as the effect for the first Month 5, and the Wind:fMonth6 effect, for example, tests for a difference in the Wind effect between month 5 (reference) and month 6. As before, you could change this behavior by changing the contrasts (see the section on contrasts, Section 2.6.1). For example, removing the intercept gives a separate intercept for each month, while the Wind effect is still estimated for month 5, with differences for the other months:
fit = lm(Ozone ~ Wind * fMonth - 1, data = airquality)and removing also the Wind main effect gives a separate intercept and a separate Wind slope for each month:
fit = lm(Ozone ~ Wind * fMonth - 1 - Wind, data = airquality)2.10.3 Interactions and centering
A super important topic when working with numeric interactions is centering. To get this started, let’s make a small experiment:
What you should have seen in the models above is that the main effects for Wind / Temp change significantly when the interaction is introduced.
Maybe you know the answer already. If not, consider the following simulation, where we create data with known effect sizes:
# Create predictor variables.
x1 = runif(100, -1, 1)
x2 = runif(100, -1, 1)
# Create response for lm, all effects are 1.
y = x1 + x2 + x1*x2 + rnorm(100, sd = 0.3)
# Fit model, but shift the mean of the predictor.
fit = lm(y ~ x1 * I(x2 + 5))
summary(fit)
Call:
lm(formula = y ~ x1 * I(x2 + 5))
Residuals:
Min 1Q Median 3Q Max
-0.78590 -0.24028 -0.01261 0.23085 0.72614
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.52248 0.32056 -17.227 < 2e-16 ***
x1 -3.56776 0.60688 -5.879 5.98e-08 ***
I(x2 + 5) 1.09592 0.06277 17.460 < 2e-16 ***
x1:I(x2 + 5) 0.89799 0.11935 7.524 2.85e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3248 on 96 degrees of freedom
Multiple R-squared: 0.8632, Adjusted R-squared: 0.8589
F-statistic: 201.9 on 3 and 96 DF, p-value: < 2.2e-16plot(allEffects(fit))
Play around with the shift in x2, and observe how the effects change. Try how the estimates change when centering the variables via the scale() command. If you understand what’s going on, you will realize that you should almost always center your predictors (at the mean or at a meaningful reference value, Section 2.9.2) whenever you use interactions; scale(x, scale = FALSE) centers without scaling.
Excellent explanations of the issues also in Schielzeth (2010).
Schielzeth, H. (2010) “Simple means to improve the interpretability of regression coefficients,” Methods in Ecology and Evolution, 1(2), pp. 103–113.
Describing interactions between numeric predictors. With an interaction, the slope of one predictor changes with the value of the other. Instead of interpreting the main effects, you can describe this directly by the slope of one predictor at a few typical values of the other, e.g. its 10%, 50% and 90% quantiles, with confidence intervals. This is easy with the emmeans package (see Section 2.6.3), whose functions all work with model predictions at fixed values of the other predictors: emtrends() gives the slopes, and emmip() plots the predictions:
library(emmeans)
aq = airquality[complete.cases(airquality[, c("Ozone", "Wind", "Solar.R")]), ]
aq$cWind = aq$Wind - mean(aq$Wind)
aq$cSolar = aq$Solar.R - mean(aq$Solar.R)
fitNN = lm(Ozone ~ cWind * cSolar, data = aq)
q = round(quantile(aq$cSolar, c(0.1, 0.5, 0.9)))
emtrends(fitNN, ~ cSolar, var = "cWind", at = list(cSolar = q)) cSolar cWind.trend SE df lower.CL upper.CL
-148 -2.34 1.270 107 -4.87 0.184
22 -5.79 0.667 107 -7.11 -4.465
100 -7.37 0.960 107 -9.27 -5.468
Confidence level used: 0.95 emmip(fitNN, cSolar ~ cWind, at = list(cSolar = q, cWind = seq(-7, 10, 1)), CIs = TRUE)
A possible results sentence: “Ozone decreased with wind speed, more strongly at high solar radiation: the slope was -2.3 ppb per mph (95% CI -4.9 to 0.2) at low radiation (10% quantile), -5.8 (-7.1 to -4.5) at median and -7.4 (-9.3 to -5.5) at high radiation (90% quantile).”
Which number is “the effect” of wind? In fitNN, the slope of cWind on a given day is b1 + b3 * cSolar. With centered predictors, the coefficient of cWind (-5.3) is the slope at the mean solar radiation. Because the model is linear, it is also the average of the day-specific slopes (we call this the data-averaged effect, Section 2.11.2) (check: mean(coef(fitNN)["cWind"] + coef(fitNN)["cWind:cSolar"] * aq$cSolar)). Uncentered, it would be the slope at Solar.R = 0 (-1.6). The same holds for a quadratic term with a centered predictor: mean-centering makes the coefficient of the linear term the data-averaged effect, exactly for quadratic terms and bilinear interactions, as a good approximation otherwise.
2.10.4 Interactions in GAMs
The equivalent of an interaction in a GAM is a tensor spline. It describes how the response varies as a smooth function of 2 predictors. You can specify a tensor spline with the te() command, as in
library(mgcv)
fit = gam(Ozone ~ te(Wind, Temp) , data = airquality)
plot(fit, pages = 0, residuals = T, pch = 20, lwd = 1.9, cex = 0.4)
2.11 Predictions and effects
This section shows how to predict from a fitted model, and how to use predictions to report an effect that is not a single number.
2.11.1 Predictions with confidence intervals
To predict with a fitted model
fit = lm(Ozone ~ Wind, data = airquality)we can use the predict() function.
newData = data.frame(Wind = c(5, 10))
predict(fit) # predicts on current data
predict(fit, newdata = newData) # predicts on new data
predict(fit, newdata = newData, se.fit = T)With se.fit = T, predict() also returns the standard error of the predicted mean. An approximate 95% confidence interval for the mean prediction is fit ± 1.96 * se (large-sample approximation; interval = "confidence" in predict() uses the exact t quantile). For new individual observations use interval = "prediction", which also includes the residual variation and is therefore much wider.
So, let’s predict Ozone with a 95% CI for wind between 2.5 and 10 (within the observed range of 2.3 to 20.7, so that we do not extrapolate):
newData = data.frame(Wind = seq(2.5, 10, 0.1))
pred = predict(fit, newdata = newData, se.fit = T)
lower = pred$fit - 1.96 * pred$se.fit
upper = pred$fit + 1.96 * pred$se.fit
plot(newData$Wind, pred$fit, type = "l", ylim = range(lower, upper), xlab = "Wind", ylab = "Ozone")
lines(newData$Wind, lower, lty = 2)
lines(newData$Wind, upper, lty = 2)
2.11.2 Reporting effects that are not constant
The regression table gives one coefficient per term, but this coefficient is often not the number that we want to report. We have met three reasons in this chapter: for a factor, a coefficient is the difference to the reference level, while we may want other comparisons, e.g. between all pairs of levels or against a control; in a model with an interaction, a coefficient is the effect at one value of the other predictor, usually its mean after centering, while we may want the effect at other values, or an average over them (Section 2.10.3); and for a curved effect (a polynomial or a GAM smooth), the slope changes along the predictor (on a curve, the slope at 60 °F is not the slope at 90 °F), so that no single coefficient is the effect at all. The same question returns with GLMs, where the link function makes every effect vary on the response scale (Section 9.4).
In all these cases, we have to decide which number describes the effect. Beyond the comparisons between factor levels, which we have covered, there are three ways to report an effect that is not constant:
- The curve: an effect plot with its confidence band (Section 2.3). Always show this first.
- Slopes at chosen values: the slope at, e.g., the 10%, 50% and 90% quantiles of the predictor, or at biologically meaningful values, each with a confidence interval (
emtrends(), as for interactions in Section 2.10.3; it works for lm and gam objects, and holds the other predictors at their means). - The data-averaged effect: the slope averaged over all observations (called the average marginal effect, AME, in econometrics). This is one number that summarizes the effect, e.g. for comparisons between predictors or between studies. For a factor, the corresponding quantity is the data-averaged difference between its levels; in a linear model, if the covariates are balanced between the groups (or there are none), this is simply the
emmeanscontrast.
As an example, we model the curved effect of temperature on ozone with a polynomial and with a GAM, and compute the slopes at the quantiles:
library(emmeans)
library(mgcv)
aq = airquality[complete.cases(airquality[, c("Ozone", "Temp", "Wind")]), ]
fitQ = lm(Ozone ~ poly(Temp, 2) + Wind, data = aq)
fitG = gam(Ozone ~ s(Temp) + Wind, data = aq)
q = quantile(aq$Temp, c(0.1, 0.5, 0.9))
emtrends(fitQ, ~ Temp, var = "Temp", at = list(Temp = q)) Temp Temp.trend SE df lower.CL upper.CL
64.5 0.146 0.519 112 -0.882 1.17
79.0 2.202 0.257 112 1.693 2.71
90.0 3.762 0.575 112 2.623 4.90
Confidence level used: 0.95 emtrends(fitG, ~ Temp, var = "Temp", at = list(Temp = q)) Temp Temp.trend SE df lower.CL upper.CL
64.5 -0.0078 0.858 111 -1.71 1.69
79.0 2.8183 0.650 111 1.53 4.11
90.0 2.7968 0.898 111 1.02 4.58
Confidence level used: 0.95 The slope rises from about zero at low temperatures to 2-4 ppb per °F at high temperatures; at 90 °F, where there are few data, the two models disagree (the polynomial slope keeps rising to 3.8, the GAM levels off at about 2.8), which is a good reason to show the curve.
The data-averaged effect is easy to compute from predictions: imagine nudging the predictor up a little for every observation, recording how much each prediction changes, and averaging. Formally, the slope at an observation is the change of the prediction for a small change h of the predictor, and we average these slopes over all observations (transform() returns a copy of the data with the shifted column; any small h works, e.g. 1% of the standard deviation of the predictor):
h = 0.01
d = (predict(fitQ, transform(aq, Temp = Temp + h)) - predict(fitQ, transform(aq, Temp = Temp - h))) / (2 * h)
mean(d)[1] 2.039344
TipConfidence interval for the data-averaged slope
To get a confidence interval for the data-averaged slope, let emtrends() compute the slopes at all observed temperatures and average them with the frequencies of these values in the data:
ut = sort(unique(aq$Temp))
w = as.numeric(table(factor(aq$Temp, levels = ut))) / nrow(aq) # frequency of each value
confint(contrast(emtrends(fitQ, ~ Temp, var = "Temp", at = list(Temp = ut)), list(average = w))) contrast estimate SE df lower.CL upper.CL
average 2.04 0.244 112 1.56 2.52
Confidence level used: 0.95 This works as long as temperature is the only predictor with a curved or interacting effect; with several such predictors, use the finite differences above and a bootstrap for the confidence interval (Section 2.12.1).
Caution
A data-averaged slope near zero does not mean that there is no effect: a hump-shaped response (rising, then falling) averages to about zero. If the slope changes sign, report the curve and the position of the optimum, not the average. The average also depends on where the observations lie: it describes the fitted curve over the covariate values of this sample, and its confidence interval assumes that the model is right.
A possible results sentence: “Ozone increased on average by 2.0 ppb per °F (95% CI 1.6 to 2.5), with the slope rising from 0.1 ppb per °F at 64 °F (10% quantile) to 3.8 at 90 °F (90% quantile).”
2.12 Nonparametrics and the bootstrap
The linear model (LM) has the beautiful property that it delivers unbiased estimates, calibrated p-values and all that based on analytical expressions that remain valid for any sample size, provided residuals are normally distributed. Unfortunately, there are very few extensions of the LM that completely maintain all of these desirable properties for limited sample sizes. The underlying problem is that for more complicated models, there are usually no closed-form solutions for the expected null or sampling distributions, or those distributions depend on the observed data. For that reason, non-parametric methods to generate confidence intervals and p-values are increasingly used when going to more complicated models, which motivates this first small excursion into non-parametric methods, and in particular cross-validation and the bootstrap.
There are two non-parametric methods that you should absolutely know and that form the backbone of basically all non-parametric for regression analyses with statistical and machine learning models. Those two are:
- the bootstrap: generate variance of estimates, allows us to calculate non-parametric CIs and p-values
- cross-validation: generate non-parametric estimates of the predictive error of the model (connects to ANOVA and AIC, see the chapters on ANOVA, Chapter 3, and model selection, Chapter 8)
Here, we will look at the bootstrap. Cross-validation is introduced in the chapter on model selection (Section 8.4).
TipExcursion: randomization null models
Another non-parametric method that you may encounter are randomization null models. Let’s look at this with an example: here, we create data from a normal and a lognormal distribution. We are interested in the difference of the mean of the two (test statistic).
set.seed(1337)
groupA = rnorm(50)
groupB = rlnorm(50)
dat = data.frame(value = c(groupA, groupB), group = factor(rep(c("A", "B"), each = 50)))
plot(value ~ group, data = dat)
# test statistic: difference of the means
reference = mean(groupA) - mean(groupB)We can show mathematically that the mean of normal draws, standardized by its estimated standard error, follows a Student-t distribution, which is why we can use this distribution in the t-test to calculate the probability of a certain deviation from a null expectation. If we work with different distributions, however, we may not be so sure what the null distribution should be. Either we have to do the (possibly very complicated) math to get a closed-form expression for this problem, or we try to generate a non-parametric expectation for the null distribution via a randomization null model.
The idea is that we generate a null expectation for the test statistic by re-shuffling the data in a way that conforms to the null assumption. In this case, we can just randomize the group labels, which means that we assume that both distributions (A,B) are the same (note that this is a bit more than asking if they have the same mean, but for the sake of the example, we are happy with this simplification).
nSim = 5000
nullDistribution = rep(NA, nSim)
for(i in 1:nSim){
sel = dat$value[sample.int(100, size = 100)]
nullDistribution[i] = mean(sel[1:50]) - mean(sel[51:100])
}
hist(nullDistribution, xlim = c(-2,2))
abline(v = reference, col = "red")
(1 + sum(abs(nullDistribution) >= abs(reference))) / (nSim + 1) # two-sided p-value[1] 0.00019996Voila, we have a p-value. None of the randomizations was as extreme as the observed difference, so p < 1/nSim; the smallest p-value that we can report is 1/(nSim + 1) (never p = 0), and it is significant. Randomization null models are used in many R packages where analytical p-values are not available, e.g., in the packages vegan (ordinations, community analysis) and bipartite (networks).
2.12.1 Non-parametric bootstrap
The key problem to calculate parametric p-values or CIs is that we need to know the distribution of the statistics of interest that we would obtain if we (hypothetically) repeated the experiment many times. The bootstrap generates new data so that we can approximate this distribution. There are two variants of the bootstrap:
- Non-parametric bootstrap: we sample new data from the old data with replacement
- Parametric bootstrap: we sample new data from the fitted model
Let’s first look at the non-parametric bootstrap, using a simple regression of Ozone ~ Wind
fit = lm(Ozone ~ Wind, data = airquality)First, for convenience, we generate a function that extracts the estimate that we want to bootstrap from a fitted model. I decided to extract the slope estimate for Wind, but you could make difference choices
getEstimate <- function(model){
coef(model)[2]
}
getEstimate(fit) Wind
-5.550923 Next, we perform the non-parametric bootstrap
performBootstrap = function(){
sel = sample.int(n = nrow(airquality), replace = T)
fitNew = lm(Ozone ~ Wind, data = airquality[sel,])
return(getEstimate(fitNew))
}
bootstrappedEstimate = replicate(1000, performBootstrap())The bootstrapped estimates of the slope can basically be interpreted as the uncertainty of the slope estimate.
hist(bootstrappedEstimate, breaks = 50)
abline(v = mean(bootstrappedEstimate), col = "red", lwd = 2)
abline(v = getEstimate(fit), col = "red", lwd = 2, lty = 2)
sd(bootstrappedEstimate)[1] 0.8314342We estimate a standard error of 0.83, slightly higher than the parametric error of 0.69 returned by summary.lm()
summary(fit)
Call:
lm(formula = Ozone ~ Wind, data = airquality)
Residuals:
Min 1Q Median 3Q Max
-51.572 -18.854 -4.868 15.234 90.000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 96.8729 7.2387 13.38 < 2e-16 ***
Wind -5.5509 0.6904 -8.04 9.27e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 26.47 on 114 degrees of freedom
(37 observations deleted due to missingness)
Multiple R-squared: 0.3619, Adjusted R-squared: 0.3563
F-statistic: 64.64 on 1 and 114 DF, p-value: 9.272e-13In R, you don’t have to code the bootstrap by hand, you can also use the boot package. The code will not necessarily be shorter, because you still need to code the statistics, but it offers a few convenient downstream functions for CIs etc.
TipCode for same example using the boot package
# see help, needs arguments d = data, k = index
getBoot = function(d, k){
fitNew = lm(Ozone ~ Wind, data = d[k,])
return(getEstimate(fitNew))
}
library(boot)
b1 = boot(airquality, getBoot, R = 1000)
b1
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = airquality, statistic = getBoot, R = 1000)
Bootstrap Statistics :
original bias std. error
t1* -5.550923 -0.04177248 0.81135222.12.2 Parametric bootstrap
The parametric bootstrap does the same thing, but the data is not re-sampled, but generated from the fitted model. This can be done conveniently in R as long as the regression package implements a simulate() function that works on the fitted model.
performBootstrapParametric = function(){
newData = model.frame(fit)
newData$Ozone = unlist(simulate(fit))
fitNew = lm(Ozone ~ Wind, data = newData)
return(getEstimate(fitNew))
}
bootstrappedEstimateParametric = replicate(1000,
performBootstrapParametric())In the results, we see that the parametric bootstrap is much closer to the analytical results
hist(bootstrappedEstimateParametric, breaks = 50)
abline(v = mean(bootstrappedEstimateParametric), col = "red", lwd = 2)
abline(v = getEstimate(fit), col = "red", lwd = 2, lty = 2)
sd(bootstrappedEstimateParametric)[1] 0.6835152
TipCode for same example using the boot package
For the parametric bootstrap, we just need a function of the data to return the estimate
getBoot = function(d){
fitNew = lm(Ozone ~ Wind, data = d)
return(getEstimate(fitNew))
}Plus, we need a function that generates the simulations.
rgen = function(dat, mle){
newData = model.frame(mle)
newData$Ozone = unlist(simulate(mle))
return(newData)
}The call to boot is then
b2 = boot(airquality, getBoot, R = 1000, sim = "parametric", ran.gen = rgen, mle = fit)
b2
PARAMETRIC BOOTSTRAP
Call:
boot(data = airquality, statistic = getBoot, R = 1000, sim = "parametric",
ran.gen = rgen, mle = fit)
Bootstrap Statistics :
original bias std. error
t1* -5.550923 0.02205086 0.67357572.13 Take-home messages
Tip
- Check the residuals before you interpret any p-value or estimate. If the checks reject the model, do not interpret it, change it: first fix a misfit of the mean (functional form), and only then look at distributional problems.
- A coefficient in a multiple regression is the effect adjusted for the other predictors in the model. If a collinear predictor is missing, part of its effect is absorbed by the remaining predictors. Collinearity itself does not bias the estimates, but it makes them more uncertain.
- For a categorical predictor, the coefficients are differences to the reference level (alphabetically first, change it with
relevel()). Test the factor as a whole with an ANOVA, and compare levels with post-hoc tests (e.g.emmeans) that correct for multiple testing. - Center numeric predictors, in particular if the model contains interactions, so that the intercept and main effects refer to a meaningful point (usually the mean); with centered predictors, a main effect of a linear model is also the data-averaged effect (exactly for quadratic terms and bilinear interactions). Scaling (by 1 SD) is optional: it makes effects comparable and often helps complex models to converge, but the effects are then no longer in natural units; say which version you report.
- Do not conclude “no effect” from a non-significant p-value: say that there is no significant evidence for an effect, and look at the confidence interval to see which effect sizes remain plausible.
- For a curved effect, report the curve, slopes at chosen values and, if one number is needed, the data-averaged slope with its confidence interval, not the polynomial coefficients alone.
2.14 Case studies
2.14.1 Univariate regression
Caution 2.2: Exercise — Plant height
Look at the plantHeight dataset in EcoData. Let’s assume we want to analyze whether height of plant species from around the world depends on temperature at the location of occurrence. Note that “loght” = log10(height).
library(EcoData)
plot(loght ~ temp, data = plantHeight)
Fit a linear model for this relationship. Perform residual checks (Section 2.4) and modify the model if you think it is necessary (Section 2.5). Interpret the fitted model - how would you present the results in a scientific publication?
The data set also includes a categorical variable “growthform”. Test if growthform has an effect on the plant height. First check how many observations each growth form has, and remove the growth forms that have only one observation. Also, consider the reference level of the treatment contrasts (Section 2.6.1) in this case.
TipSolution
1. Temperature dependence
model = lm(loght ~ temp, data = plantHeight)par(mfrow = c(2, 2))
plot(model)
We’re seeing a bit of heteroskedasticity, maybe a tiny bit of misfit, otherwise it looks OK! If you want to do something:
- Heteroskedasticity - see the chapter on dispersion models (Chapter 10)
- Misfit: could change functional form (Section 2.5) or use a GAM (Section 2.5.3)
But that doesn’t really seem to be a problem, so we can now interpret this model!
model = lm(loght ~ temp, data = plantHeight)
summary(model)
Call:
lm(formula = loght ~ temp, data = plantHeight)
Residuals:
Min 1Q Median 3Q Max
-1.97903 -0.42804 -0.00918 0.43200 1.79893
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.225665 0.103776 -2.175 0.031 *
temp 0.042414 0.005593 7.583 1.87e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6848 on 176 degrees of freedom
Multiple R-squared: 0.2463, Adjusted R-squared: 0.242
F-statistic: 57.5 on 1 and 176 DF, p-value: 1.868e-122. Including “growthform”
Before fitting the model, we check how many observations each growth form has:
table(plantHeight$growthform)
Fern Herb Herb/Shrub Shrub Shrub/Tree Tree
1 44 1 49 12 61 Fern and Herb/Shrub have only one observation each. The mean of such a group rests on a single point, and interactions with temperature cannot be estimated for it at all, so we remove these two growth forms:
plantHeight2 = plantHeight[!plantHeight$growthform %in% c("Fern", "Herb/Shrub"), ]
plantHeight2$growthform = as.factor(plantHeight2$growthform) # drops the unused levels; a factor is needed for the post-hoc tests below
table(plantHeight2$growthform)
Herb Shrub Shrub/Tree Tree
44 49 12 61 boxplot(loght ~ growthform, data = plantHeight2)
Fit the model
model2 = lm(loght ~ growthform, data = plantHeight2)
summary(model2)
Call:
lm(formula = loght ~ growthform, data = plantHeight2)
Residuals:
Min 1Q Median 3Q Max
-1.3229 -0.2609 0.0381 0.2515 1.5076
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.41419 0.06281 -6.594 5.76e-10 ***
growthformShrub 0.58222 0.08653 6.728 2.81e-10 ***
growthformShrub/Tree 1.23945 0.13569 9.135 2.57e-16 ***
growthformTree 1.67510 0.08241 20.327 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4166 on 162 degrees of freedom
(10 observations deleted due to missingness)
Multiple R-squared: 0.7363, Adjusted R-squared: 0.7314
F-statistic: 150.8 on 3 and 162 DF, p-value: < 2.2e-16The reference level (intercept) is the first level in alphabetical order, here Herb. It is well sampled, which is a good choice: the comparisons of the other growth forms with it are estimated precisely. Had we kept Fern, which comes first alphabetically, the reference level would have been a single observation. With another reference level that is better for your question, use relevel(): it is the same model, only the coefficients change.
The overall test of whether growth form has an effect is the ANOVA (Section 2.6.2), possible post-hoc tests (Section 2.6.3) follow up on it.
anova(model2)Analysis of Variance Table
Response: loght
Df Sum Sq Mean Sq F value Pr(>F)
growthform 3 78.530 26.1768 150.8 < 2.2e-16 ***
Residuals 162 28.121 0.1736
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(multcomp)
tuk = glht(model2, linfct = mcp(growthform = "Tukey"))
summary(tuk) # Standard display.
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lm(formula = loght ~ growthform, data = plantHeight2)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
Shrub - Herb == 0 0.58222 0.08653 6.728 < 0.001 ***
Shrub/Tree - Herb == 0 1.23945 0.13569 9.135 < 0.001 ***
Tree - Herb == 0 1.67510 0.08241 20.327 < 0.001 ***
Shrub/Tree - Shrub == 0 0.65723 0.13419 4.898 < 0.001 ***
Tree - Shrub == 0 1.09287 0.07993 13.673 < 0.001 ***
Tree - Shrub/Tree == 0 0.43565 0.13157 3.311 0.00584 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)tuk.cld = cld(tuk) # Letter-based display.
par(mar = c(5,3,10,3))
plot(tuk.cld)
TipSolution
1. Fitting the species-area curve
model = lm(Species ~ log(Area), data = solomonislands)
summary(model)
Call:
lm(formula = Species ~ log(Area), data = solomonislands)
Residuals:
Min 1Q Median 3Q Max
-31.422 -3.705 3.460 6.309 13.014
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 30.105 1.613 18.67 <2e-16 ***
log(Area) 4.935 0.359 13.75 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.6 on 50 degrees of freedom
Multiple R-squared: 0.7907, Adjusted R-squared: 0.7866
F-statistic: 188.9 on 1 and 50 DF, p-value: < 2.2e-16plot(allEffects(model))
Area has a small positive significant effect on the number of species.
2. Residual checks
par(mfrow = c(2, 2))
plot(model)
Residuals don’t look particularly normal … would get better fitting ~ log(Area+1), but that changes the model of the exercise. In this case, it’s a trade-off between fitting the model that was asked for and getting better residuals.
With advanced methods, we could also think about GLMs (Chapter 9) or quantile regressions (Section 10.3.2).
3. Interpretation
Significant increase of species richness with area that follows a log(Area) relationship.