Example: Simple Bivariate Linear Regression
Linear regression modeling is one of the more common tasks in data analysis and the following example will cover the very basic topic of bivariate linear regression. The storyline follows the one from Zuur et al. (2007) to a certain degree.
While one could use actual data sets, we keep it controlled by using an artificial data set originally compiled by Francis Anscombe. The anscombe dataset comes as part of base R. For now, we will use x1 as independent variable and y1 as dependent variable.
ind <- anscombe$x1
dep <- anscombe$y1
par_org <- par(no.readonly = TRUE)
par(mfrow = c(1,2))
hist(ind, main = "Independent variable")
hist(dep, main = "Dependent variable")
A look at the relationship between the variables by using a scatterplot justifies a linear modeling attempt. Fitting a bivariate linear model to the data set would result in the function shown in the plot and visualized by the red line.
par(par_org)
plot(ind, dep)
lmod <- lm(dep ~ ind) # Compute linear regression model
abline(lmod, col = "red") # Add the fitted line using base R
text(10, 6, bquote(paste("y"["i"], " = ", .(round(lmod$coefficients[1], 3)),
" + ", .(round(lmod$coefficients[2], 3)), "* x"["i"])))
The fitted line describes the observed relationship. Under suitable model assumptions, an ANOVA or t-test can assess evidence against a zero slope. Neither test checks all those assumptions or establishes a causal relationship. We calculate the summaries first, then use the separate assumptions walkthrough to judge how cautiously they should be interpreted.
Testing a linear regression relationship by an analysis of variance
Let’s start with the ANOVA. Here we divide the variation in the response into variation described by the fitted line and variation left in the residuals.
Given a (linear) model, each actual data value can be calculated by adding the fitted value and the corresponding residual value:
data value = fitted value + residual value (or y = y’ + res)
The corresponding differences are:
- Observed differences: each observed y minus the mean response, describing total variation.
- Model differences: each fitted y’ minus the mean response, describing variation explained by the model.
- Residual differences: each observed y minus its fitted y’, describing the remaining variation.
Square the differences and sum them to obtain the total, model and residual sums of squares. For this ordinary least-squares model with an intercept, the model and residual sums of squares add to the total sum of squares.
In order to calculate the variances, one can use the lm class of the model since - among others - it contains the original independent and dependent values as well as the predicted ones:
ss_obsrv <- sum((lmod$model$dep - mean(lmod$model$dep))**2)
ss_model <- sum((lmod$fitted.values - mean(lmod$model$dep))**2)
ss_resid <- sum((lmod$model$dep - lmod$fitted.values)**2)
Since the sum increases with increasing numbers of observations, the resulting sums of squares are normalized by the respective degrees of freedom. This gives us:
- the mean observed sum of squares
- the mean model sum of squares
- the mean residual sum of squares, i.e. the mean squared error if the model is a simple linear regression model.
mss_obsrv <- ss_obsrv / (length(lmod$model$dep) - 1)
mss_model <- ss_model / 1
mss_resid <- ss_resid / (length(lmod$model$dep) - 2)
Under the zero-slope null hypothesis and the usual model assumptions, both the model and residual mean squares estimate the error variance. Their ratio fluctuates around a reference distribution; it is not exactly one in every sample. A systematic relationship can make the model mean square larger relative to the residual mean square.
This provides us with the test statistic for the null-hypothesis that the true slope is not different from 0.
f_value <- mss_model / mss_resid
The F-statistic quantifies how much larger the variation explained by the model is compared to the variation left unexplained (residuals).
- A large F value indicates that the model explains a lot more variation than is left unexplained, suggesting that the slope is likely not 0.
- A small F value indicates that the model doesn’t explain much more variation than the residuals, suggesting that the slope might be 0.
By comparing the computed F value to an F distribution with 1 numerator degree of freedom and (number of observations - 2) denominator degrees of freedom, the p value is determined. If the p value is very small, the null-hypothesis can be rejected.
p_value <- pf(f_value, 1, length(dep) - 2, lower.tail = FALSE)
In summary, this gives us the following:
data.frame(Name = c("ind", "residuals"),
Deg_of_freedom = c(1, length(lmod$model$dep) - 2),
Sum_of_squares = c(ss_model, ss_resid),
Mean_sum_of_squares = c(mss_model, mss_resid),
F_value = f_value,
p_value = p_value)
## Name Deg_of_freedom Sum_of_squares Mean_sum_of_squares F_value
## 1 ind 1 27.51000 27.510001 17.98994
## 2 residuals 9 13.76269 1.529188 17.98994
## p_value
## 1 0.002169629
## 2 0.002169629
Of course, one does not have to compute this every time. A simple call to the anova function will do it:
anova(lmod)
## Analysis of Variance Table
##
## Response: dep
## Df Sum Sq Mean Sq F value Pr(>F)
## ind 1 27.510 27.5100 17.99 0.00217 **
## Residuals 9 13.763 1.5292
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Exercise: Recreate the regression plot
Now you can check your knowledge and strengthen your plotting skills by writing code to rebuild the linear regression scatterplot from the previous page. When coding, think about the meaning of slope, intercept, observed variance, model variance, residual variance, and the different types of sums of squares.
R squared
The variance explained by the model is one of the most often used variables to explain the relationship in simple linear models. It is computed by normalizing either the model sum of squares by the observed sum of squares or by subtracting the normalization of the residual sum of squares from 1.
ss_model / ss_obsrv
## [1] 0.6665425
1 - ss_resid / ss_obsrv
## [1] 0.6665425
Convention. p and r-squared often denote results of bivariate linear regressions while P and R-squared (capital letters) often denote results of multiple linear regressions. The letters p (P) and r2 (R2) are often written in italics.
Finished?
Well, the above looks like a really good example of linear regression analysis, right? And the r-squared of about 0.67 is also quite OK not to mention the significance of the independent variable.
Before we clap our hands, let’s just have a look at the other variable combinations of the Anscombe data set.
par(mfrow = c(2,2))
for(i in seq(4)){
plot(anscombe[, i], anscombe[, i+4],
xlab = names(anscombe[i]), ylab = names(anscombe[i+4]))
}
Reading the panels from top left to bottom right: the first dataset is approximately linear; the second follows a curve; the third is approximately linear apart from one unusually high response; and the fourth has nearly all predictor values at the same position, with one point far to the right. These different patterns explain why similar regression summaries do not replace inspecting the observations.
While x3/y3 might still justify a linear regression if we remove the outlier, the two plots on the right side do not. So what to do? Unfortunately, all of these data combinations result in almost the same regression statistics:
lmodels <- lapply(seq(4), function(i){
lm(anscombe[, i+4] ~ anscombe[, i])
})
lapply(lmodels, summary)
## [[1]]
##
## Call:
## lm(formula = anscombe[, i + 4] ~ anscombe[, i])
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.92127 -0.45577 -0.04136 0.70941 1.83882
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.0001 1.1247 2.667 0.02573 *
## anscombe[, i] 0.5001 0.1179 4.241 0.00217 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.237 on 9 degrees of freedom
## Multiple R-squared: 0.6665, Adjusted R-squared: 0.6295
## F-statistic: 17.99 on 1 and 9 DF, p-value: 0.00217
##
##
## [[2]]
##
## Call:
## lm(formula = anscombe[, i + 4] ~ anscombe[, i])
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9009 -0.7609 0.1291 0.9491 1.2691
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.001 1.125 2.667 0.02576 *
## anscombe[, i] 0.500 0.118 4.239 0.00218 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.237 on 9 degrees of freedom
## Multiple R-squared: 0.6662, Adjusted R-squared: 0.6292
## F-statistic: 17.97 on 1 and 9 DF, p-value: 0.002179
##
##
## [[3]]
##
## Call:
## lm(formula = anscombe[, i + 4] ~ anscombe[, i])
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.1586 -0.6146 -0.2303 0.1540 3.2411
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.0025 1.1245 2.670 0.02562 *
## anscombe[, i] 0.4997 0.1179 4.239 0.00218 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.236 on 9 degrees of freedom
## Multiple R-squared: 0.6663, Adjusted R-squared: 0.6292
## F-statistic: 17.97 on 1 and 9 DF, p-value: 0.002176
##
##
## [[4]]
##
## Call:
## lm(formula = anscombe[, i + 4] ~ anscombe[, i])
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.751 -0.831 0.000 0.809 1.839
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.0017 1.1239 2.671 0.02559 *
## anscombe[, i] 0.4999 0.1178 4.243 0.00216 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.236 on 9 degrees of freedom
## Multiple R-squared: 0.6667, Adjusted R-squared: 0.6297
## F-statistic: 18 on 1 and 9 DF, p-value: 0.002165
So be careful when interpreting p and R-squared values! We will learn how to treat non-linear relationships in the subsequent units.
par(par_org)
Continue with model checking
The Anscombe examples show why similar regression summaries are not enough. Next, work through Checking Regression Assumptions for residual plots, the graphical plausibility check, normality-test limitations and possible responses to concerns. That page prepares you for Assignment 05.