The title of this project is R Workshop - Linear Models v1 28-05-2024. The slides for the presentation accompanying this workflow can be found here.
In this workflow we focus on practical data analysis by presenting statistical workflows applicable in any software for four of the most common univariate analyses: linear regression, ANOVA, ANCOVA, and repeated measures (a simple mixed model) – all assuming a normal (gaussian) residual. These workflows can be easily extended to more complex models. The R code used to create output is also included.
In order to re-run this .qmd file, either open it with RStudio and press “render” or, if you want both the .pdf and the .html file to be produced, run the following code with this .qmd file being present in your working directory: quarto::quarto_render("linear_modelsI_workflow.qmd", output_format = "all")
If you used this workflow, please don’t forget to acknowledge us. You can use the following sentence:
The authors acknowledge the Statistical Workshops and Workflows provided by the Sydney Informatics Hub, a Core Research Facility of the University of Sydney.
Before we delve into the content, it is necessary to load certain libraries. If you want to use them, please install them beforehand with install.packages("library_name").
suppressPackageStartupMessages({
library("tibble")
library("magrittr")
library("ggplot2")
library("emmeans")
library("lme4")
library("lmerTest")
library("gglm")
library("patchwork")
library("writexl")
})
# GGplot theme
theme_set(theme_bw())
# No warnings
options(warn = -1)For this part, we assume that we have a dataset called dataset with two variables:
response which is the response or target variable, i.e. the ‘variable of interest’.predictor1 which is the predictor variable, i.e. the variable that we are using to explain the variance in response.predictor1 <- rnorm(100, 6, 0.25)
variance <- rnorm(100, 0, 0.05)
response <- 1 + 0.5 * predictor1 + variance
dataset <- data.frame(response, predictor1)In this section, we first check for validity of assumptions prior to the formal model fitting procedure using statistical diagnostics.
A simple linear model assumes a linear relationship between our predictor(s) and the target variable. In order to verify this, we create a scatterplot between predictor1 and response:
ggplot(data = dataset, aes(x = predictor1, y = response)) +
geom_point() +
labs(x = "predictor1 values", y = "response values")
As we can see, the response values rise nicely with increased values of predictor1, indicating a linear relationship. Have a look at a graph, in which linearity is clearly not present:
tibble(x = rnorm(100), y = 1 + x^2 + rnorm(100, sd = 0.5)) %>%
ggplot(data = ., aes(x = x, y = y)) +
geom_point() +
labs(x = "predictor1 values", y = "response values")
In this above graph, the response is actually dependent on \(x^2\), not \(x\).
Another assumption of linear models is that the observations are not dependent on each other. One way to check this is to do a serial plot, which lines up all observations in order of appearance. If any pattern can be detected, they are likely to not be independent.
ggplot(data = dataset, aes(
x = seq_along(response),
y = response
)) +
geom_point() +
labs(x = "index", y = "response values") +
ggplot(data = dataset, aes(
x = seq_along(response),
y = predictor1
)) +
geom_point() +
labs(x = "index", y = "predictor.linear1 values")
There are no visible patterns in the above graph. A problematic pattern could occur in time-dependent observations, like below:
y0 <- 10
y_1 <- stats::filter(c(y0, runif(99, -0.5, 0.5)), 0.75, method = "recursive")
qplot(seq_along(y_1), y_1, geom = "point") +
labs(x = "index values", y = "response values")Don't know how to automatically pick scale for object of type <ts>. Defaulting
to continuous.
In the graph above, we can clearly see how each observations depend on the ones that came before.
The assumption of normality is technically only relevant for residuals, which can only be obtained post-model fit. But it’s still a good idea to check normality of the response variable , since they are related.
First, we create a histogram of all response observations, and have a look at the predictor1 values as well:
ggplot(data = dataset, aes(x = response, y = after_stat(density))) +
geom_histogram(fill = "cornflowerblue", bins = 15, col = "black") +
ggplot(data = dataset, aes(x = predictor1, y = after_stat(density))) +
geom_histogram(fill = "cornflowerblue", bins = 15, col = "black")
What we’re looking for here is a roughly symmetric distributions with a single peak.
This is a also very poorly understood assumption. We want a model represent the bulk of the data. We don’t want it biased towards 1 or 2 outlying influential points. Just like checking the normality assumption we can only test this for sure once we have fit a model. However, it is always worth looking at all our data to see if there are any outliers we might need to deal with. The best way to do this is via histograms, re-using the one created above.
After a sucessful Exploratory Data Analysis, we use lm() to fit a linear model.
model <- lm(response ~ predictor1, data = dataset)In order to gauge the model fit, we first visually display the residuals. Residuals are variability in our response that cannot be explained by the model. Visually, they appear as the distance between the fitted line and the observed datapoints. Have a look:
ggplot(data = dataset, aes(x = predictor1, y = response)) +
geom_point() +
geom_line(data = NULL, aes(x = predictor1, y = fitted(model)), col = "firebrick") +
labs(x = "predictor1 values", y = "response values")
We are looking for a few things here: Firstly, that the deviation from the line is roughly equal across the spread of predictor1. Secondly, that the residuals are normally distributed. We can also display the deviations from the line in a histogram including a kernel density graph, like so:
model_res <- residuals(model)
ggplot(data = NULL, aes(x = model_res, y = after_stat(density))) +
geom_histogram(fill = "cornflowerblue", bins = 15, col = "black") +
labs(x = "predictor1 values", y = "density") +
ggplot(data = NULL, aes(x = model_res, y = after_stat(density))) +
geom_density(fill = "firebrick") +
geom_rug(aes(y = NULL)) +
labs(x = "predictor1 values", y = "density")
Using a histogram it is often easier to spot irregularities. In this case, using the graph above and below we can confirm that the residuals are normally and equally distributed across the line.
Other diagnostic plots can be produced with plot(model_object_name). The gglm package produces the same plots but using the graphing package ggplot2, the plots of which look a little more up-to-date:
gglm(model, theme = theme_bw()) +
ggplot(data = model) +
stat_cooks_obs() +
geom_hline(yintercept = 0.5, linetype = "dashed", col = "red") +
plot_layout(nrow = 2, heights = c(2, 1))
The following plots are produced above:
All of these plots look ideal.
In this part, we want to check the goodness-of-fit of the model. Above, we already displayed the fitted values vs the residuals of the model. It also makes sense to compare the values of each predictor (explanatory variable), as well as the response with the residuals, to make sure the linearity assumption holds for each variable that we’re predicting with. In our case, we only have one predictor variable.
# Plot 1
ggplot(dataset, aes(x = predictor1, y = residuals(model))) +
geom_point() +
geom_hline(yintercept = 0, col = "red", linetype = "dashed") +
labs(x = "predictor1", y = "model residuals", title = "predictor1 vs. residuals") +
# Plot 2
ggplot(dataset, aes(x = response, y = residuals(model))) +
geom_point() +
geom_hline(yintercept = 0, col = "red", linetype = "dashed") +
labs(x = "response", y = "model residuals", title = "response vs. residuals") +
# Plot 3
ggplot(dataset, aes(x = predict(model), y = residuals(model))) +
geom_point() +
geom_hline(yintercept = 0, col = "red", linetype = "dashed") +
labs(x = "fitted values", y = "model residuals", title = "fitted values vs residuals") +
plot_layout(ncol = 2)
As we can see, the linearity assumption holds nicely. Let’s also check the response vs predicted plot:
ggplot(dataset, aes(x = response, y = predict(model))) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, col = "red", formula = y ~ x) +
labs(y = "predicted")
This plot is a good visual representation of model fit. If the response is being exactly predicted than we expect it to fall along the 1:1 line.
Now, we look at the summary output of our linear model:
summary(model)
Call:
lm(formula = response ~ predictor1, data = dataset)
Residuals:
Min 1Q Median 3Q Max
-0.115978 -0.030085 -0.001874 0.028932 0.161342
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.94554 0.11584 8.162 1.14e-12 ***
predictor1 0.50825 0.01927 26.370 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.04691 on 98 degrees of freedom
Multiple R-squared: 0.8765, Adjusted R-squared: 0.8752
F-statistic: 695.4 on 1 and 98 DF, p-value: < 2.2e-16A model summary output has three parts:
We already looked at the residuals a lot, so we’re starting with the second section here. Each predictor and the intercept has one line in the table, which holds information on both the coefficient estimate and a significance test. A p value smaller than 0.05 indicates that we see a significant partial correlation between that predictor and the response variable. In this case, we have one predictor, the coefficient of which is significantly different from 0. It’s also a good idea to look at the confidence interval of coefficients:
confint(model) 2.5 % 97.5 %
(Intercept) 0.7156560 1.1754214
predictor1 0.4700031 0.5465012In the third section of the summary output, we can see the \(R^2\), which ranges from 0 to 1 and gives an indication on much variance of the response variable is explained by the predictor variables. The last line is important, as it tests for significance of the entire model.
In this case, a p value of less than 0.05 indicates that the model is significantly different from 0 (no model). We can therefore conclude that the model is significant. Refer to the presentation here for more detail.
To interpret the paramters, we use the fitted regression coefficient. In our case, using the model table above, we can say (in this case we use “fake” predictor and response meanings):
There is strong evidence to show that predictor1 influences response (p<2e-16), with each 1 unit of increase in predictor1 adding between 0.47-0.54 units of response (95% CI). This effect on weight has been estimated very accurately [as 95% CI is quite narrow].
The model is a good fit to the data with an \(R^2=88\%\). There were no outliers or unexplained structure. The error was normal.
In this section, we have a data.frame object with two variables:
treatment: The treatment variable (binary)response: the response variable (metric)set.seed(171975)
variance <- rnorm(100, 0, 0.05)
b0.control <- 3
b1.treatment <- 0.5
data1 <- tibble(
treatment = factor(c(rep("Control", 50), rep("Treatment", 50))),
response = b0.control + b1.treatment*ifelse(treatment=="Treatment", 1, 0) + variance
)We assume that the response variable is dependent on the treatment. Due to the binary nature of treatment, we are fitting an Analysis of Variance (a subtype of linear model).
First, let’s visually display both variables using boxplots (as we see in the powerpoint):
ggplot(data1, aes(x = treatment, y = response)) +
# geom_violin(alpha=0.4, position = position_dodge(width = .75),size=1,color="black") +
geom_point(
position = position_jitter(),
color = "black"
) +
labs(x = "Treatment", y = "Response") +
ggplot(data1, aes(x = treatment, y = response)) +
# geom_violin(alpha=0.4, position = position_dodge(width = .75),size=1,color="black") +
geom_boxplot(
notch = FALSE,
outlier.size = -1,
color = "black",
lwd = 1.2,
alpha = 0.7
) +
labs(x = "Treatment", y = "Response")
As we can see above, all observations are clearly separated, indicating a strong treatment effect.
Note: The “jitter” creates a nice visual scatter, but is purely for exploration. It should be removed in the final publication.
First, we check for independence using a serial plot:
ggplot(data = data1, aes(
x = seq_along(response),
y = response
)) +
geom_point() +
labs(x = "index", y = "response values")
We see the separation of treatment values here, but within those there doesn’t appear to be any (auto) correlation. Let’s check normality next:
ggplot(data1, aes(
x = response,
y = after_stat(density)
)) +
geom_histogram(fill = "cornflowerblue", col = "black", bins = 30) +
geom_density(aes(fill = NULL)) +
labs(title = "Full dataset") +
ggplot(data1, aes(
x = response,
y = after_stat(density)
)) +
facet_wrap(~ treatment) +
geom_histogram(fill = "cornflowerblue", col = "black", bins = 30) +
geom_density(aes(fill = NULL)) +
labs(title = "Divided by treatment") +
plot_layout(ncol = 1)
We can see clear separation here as well, otherwise normality. If we only looked in the above plot, we would have assumed non-normality. This is a great example of how we shouldn’t worry about violations of normality in the dependent variable too much, since the residuals have to be normally distributed, not the outcome.
The combined data (as seen above) exhibits a clear bimodal distribution and deviates significantly from normality. Thus, we need to address whether this poses a problem. However, the error distribution should be normal, not the response, and upon examination, the errors around the mean of each treatment appear to be approximately normal.
There are two ways to fit an ANOVA in R. You can either use the lm() function as above, or the aov() function. In this case, I’ll use the aov() function, mainly because I prefer the summary output.
anova_model <- aov(response ~ treatment, data = data1)gglm(anova_model, theme = theme_bw()) +
ggplot(data = anova_model) +
stat_cooks_obs() +
geom_hline(yintercept = 0.5, linetype = "dashed", col = "red") +
plot_layout(nrow = 2, heights = c(2, 1))
The model diagnostics plot is less informative than with the linear model, due to the binary nature of the treatment variable, but we can still observe interesting trends, like with the QQ plot for example.
Let’s investigate the residuals a bit more:
ggplot(data1, aes(
x = residuals(anova_model),
y = after_stat(density)
)) +
labs(x = "Model residuals") +
geom_histogram(fill = "cornflowerblue", col = "black", bins = 30) +
ggplot(data1, aes(
x = residuals(anova_model),
y = after_stat(density)
)) +
geom_density(fill = "firebrick") +
geom_rug(aes(y = NULL)) +
labs(x = "Model residuals")
This time, we cannot observe large differences betwen the residuals of the both treatment effects, which indicates that the differences are nicely taken into account.
Let’s have another look at the Residuals vs Fitted plot:
ggplot(data1, aes(x = treatment, y = residuals(anova_model))) +
# geom_violin(alpha=0.4, position = position_dodge(width = .75),size=1,color="black") +
geom_boxplot(
notch = FALSE,
outlier.size = -1,
color = "black",
lwd = 1.2,
alpha = 0.7
) +
geom_point(
shape = 21,
size = 2,
position = position_jitter(),
color = "black",
alpha = 1
) +
labs(x = "Treatment", y = "Residuals")
This time we can again see that the residuals are nicely spread around 0 for both treatment and control.
Another great plot is Fitted vs Residuals as as well as Response vs Residuals :
ggplot(data1, aes(
x = predict(anova_model),
y = residuals(anova_model),
col = treatment)) +
geom_point() +
labs(x = "Predictions", y = "Residuals") +
geom_hline(yintercept = 0, col = "red") +
theme(legend.position = "none") +
# Plot 2
ggplot(data1, aes(
x = response,
y = residuals(anova_model),
col = treatment)) +
geom_point() +
labs(x = "Predictions", y = "Residuals") +
geom_hline(yintercept = 0, col = "red")
We expect the ‘lines’ of data rather than a random ‘cloud’ of data which we saw in the regression. This is because rather than a range of predictions for each different value of the predictor we only get 1 prediction for control and another for treatment, hence 2 vertical lines in the upper chart.
Let’s also check the response vs. predictions:
ggplot(data1, aes(x = response, y = predict(anova_model))) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, col = "red", formula = y ~ x) +
labs(y = "predicted")
To check the goodness-of-fit, we have a look at the model summary output:
summary(anova_model) Df Sum Sq Mean Sq F value Pr(>F)
treatment 1 6.481 6.481 2400 <2e-16 ***
Residuals 98 0.265 0.003
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In this case, we purely observe the significance of our treatment effect, which is seen in the first line of the Anova-Output. A p value of less than 0.05 indicates a significant difference between the categories of the treatment variable (of which there are only two).
In order to find out how large the difference is, we need to convert the anova object back to a linear model like so:
coef(lm(anova_model)) (Intercept) treatmentTreatment
3.0017546 0.5091524 confint(lm(anova_model)) 2.5 % 97.5 %
(Intercept) 2.9871714 3.0163378
treatmentTreatment 0.4885287 0.5297761We can now see that the average difference of response between the treatment and control effect is 0.509. Since our effect is significantly different from 0, we can now conclude that our treatment is successful.
Refer to the presentation here for more detail.
You can use the following sentence (we again assume a “fake” meaning):
There is strong evidence to show that the treatment influences response (p<2e-16). It increases the response by between 0.49-0.53 units (95% CI), from an average of approximately 3 (95% CI=2.98-3.01). This effect on the response has been estimated very accurately [as 95% CI is quite narrow].
The model is a good fit to the data with an R2=97%. There were no outliers or unexplained structure. The error was normal”
In ANCOVA model situations, we are interested in a treatment effect like in ANOVA, but we also want to account for other covariates. This could, for example, be the age of clinical trial candidates.
Let’s create some data first (I collapsed this element because it’s a lot of code):
set.seed(171974)
predictor.linear1 <- rnorm(100, 6, 0.25)
variance <- rnorm(100, 0, 0.05)
treatment <- factor(c(rep("Control", 50), rep("Treatment", 50)))
# Model 1) if no difference between groups i.e. lines are parrallel
b0.control <- 1
b0.treatment <- 0.5
b1.control <- 1
b1.treatment <- 0
data2 <- data.frame()[1:100, ]
data2$treatment <- treatment
data2$predictor.linear1 <- predictor.linear1
data2$response <- b0.control +
b1.control*predictor.linear1 +
b0.treatment*ifelse(data2$treatment=="Treatment", 1, 0) +
b1.treatment*ifelse(data2$treatment=="Treatment", 1, 0)*predictor.linear1 +
variance
row.names(data2) <- NULL
# Model 2) if difference between groups i.e. lines aren't parallel
b0.control <- 1
b0.treatment <- -8.5
b1.control <- 1
b1.treatment <- 1.5
data3 <- data.frame()[1:100, ]
data3$treatment <- treatment
data3$predictor.linear1 <- predictor.linear1
data3$response <- b0.control +
b1.control*predictor.linear1 +
b0.treatment*ifelse(data3$treatment=="Treatment", 1, 0) +
b1.treatment*ifelse(data3$treatment=="Treatment", 1, 0)*predictor.linear1 +
variance
row.names(data3) <- NULLWe now have three covariates in a data.frame object called data2, and the same in a different object called data3:
treatment: A categorical variable depicting the treatment/control groupspredictor.linear1: Our numeric predictorresponse: The numeric responseLet’s look at our datasets graphically:
d2graph <- ggplot(data2, aes(x = predictor.linear1, y = response, col = treatment)) +
geom_point() +
guides(col = "none") +
labs(x = "predictor.linear1 values", title = "Scatterplot data2")
d3graph <- ggplot(data3, aes(x = predictor.linear1, y = response, col = treatment)) +
geom_point() +
labs(x = "predictor.linear1 values", title = "Scatterplot data3")
d2graph + d3graph
We can now see the relationship between predictor.linear1 and response in two different scenarios (data2 vs data3). The difference is that on the left graph, the differently coloured points have a parallel increase, whereas on the right graph the slopes of the increase between the points is also different. We will be further using data3, because it has a more complex (and interesting) model structure.
From the (right) graph above, we can see that there is a binary treatment effect as well as a linear predictor called predictor.linear1. This is a classic ANCOVA scenario.
First, let’s check whether the response observations seem to be (auto) correlated using a serial plot:
ggplot(data = data3, aes(
x = seq_along(response),
y = response,
col = treatment
)) +
geom_point() +
labs(x = "index", y = "response values") +
ggplot(data = data3, aes(
x = seq_along(response),
y = predictor.linear1,
col = treatment
)) +
geom_point() +
labs(x = "index", y = "predictor.linear1 values") +
plot_layout(guides = "collect")