# A tibble: 928 × 2
parent child
<dbl> <dbl>
1 70.5 61.7
2 68.5 61.7
3 65.5 61.7
4 64.5 61.7
5 64 61.7
6 67.5 62.2
7 67.5 62.2
8 67.5 62.2
9 66.5 62.2
10 66.5 62.2
# ℹ 918 more rowsTopic 9 – Describing relationships
ENVX1002 Introduction to Statistical Methods
Apr 2024
About me
- Lecturer in Agricultural Data Science
- Spatial modelling and mapping, Soil science, Precision Agriculture
- This year back to being a student again, GradCert (Higher Education)
Brief history
Adrien-Marie Legendre, Carl Friedrich Gauss, Francis Galton, & Karl Pearson
Code
# Generate synthetic data
set.seed(123)
rainfall <- rnorm(100, mean = 50, sd = 10)
plant_growth <- rainfall + rnorm(100, mean = 0, sd = 5)
# Calculate correlation coefficient
correlation_coef <- cor(rainfall, plant_growth)
# Plot data using ggplot2
library(ggplot2)
p1 <- ggplot(data = data.frame(rainfall, plant_growth), aes(x = rainfall, y = plant_growth)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
labs(
title = "Positive correlation",
x = "Rainfall", y = "Plant growth",
subtitle = paste("Correlation coefficient: ", round(correlation_coef, 2))
)
# Generate synthetic data
set.seed(123)
phosphorous <- rnorm(100, mean = 50, sd = 10)
chlorophyll_a <- 100 - phosphorous + rnorm(100, mean = 0, sd = 5)
# Calculate correlation coefficient
correlation_coef <- cor(phosphorous, chlorophyll_a)
# Plot data using ggplot2
p2 <- ggplot(data = data.frame(phosphorous, chlorophyll_a), aes(x = phosphorous, y = chlorophyll_a)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
labs(
title = "Negative correlation",
x = "Phosphorous", y = "Chlorophyll-a",
subtitle = paste("Correlation coefficient: ", round(correlation_coef, 2))
)
library(patchwork)
p1 + p2
E.g. what is the correlation between all continuous variables in the iris dataset?
- We can essentially use correlation analysis to identify variables that are useful for predicting another variable, or if they will present issues for model fitting i.e. multicollinearity.
Pearson correlation coefficient
Formula:
r = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^n (x_i - \bar{x})^2 \sum_{i=1}^n (y_i - \bar{y})^2}}
covariance divided by the product of the standard deviations
Luckily we have Excel and R
- Excel:
=CORREL()formula, or use the analysis toolpak - R:
cor()function
We can also visually inspect the relationship between the two variables using a scatterplot:
… but this is not a very good way to assess the strength of the relationship between the two variables.
Interpretation: strong
A strong relationship is one where the correlation coefficient is close to 1 or -1.
Interpretation: weak
A weak or nonexistent relationship is one where the correlation coefficient is close to 0.
Anscombe’s quartet
All of these data have a correlation coefficient of about 0.8.
Datasaurus Dozen
All of these data have a correlation coefficient close to zero!
Correlation \neq causation
Spurious correlations: a relationship between two variables does not imply that one causes the other.
Visualising regression models
Visualising regression models
Code
library(ggplot2)
# Generate sample data
set.seed(136)
x <- 1:100
y <- x^2 + rnorm(100, sd = 100)
# Define regression functions
slr <- function(x, y) {
mod <- lm(y ~ x)
return(list(
data.frame(x = x, y = predict(mod), model_type = "Simple Linear Regression"),
paste("y =", round(coefficients(mod)[[2]], 2), "x +", round(coefficients(mod)[[1]], 2))
))
}
mlr <- function(x, y, z) {
mod <- lm(y ~ x + z)
return(list(
data.frame(x = x, z = z, y = predict(mod), model_type = "Multiple Linear Regression"),
paste("y =", round(coefficients(mod)[[3]], 2), "x +", round(coefficients(mod)[[2]], 2), "z +", round(coefficients(mod)[[1]], 2))
))
}
poly_reg <- function(x, y, degree) {
mod <- lm(y ~ poly(x, degree, raw = TRUE))
x_new <- seq(min(x), max(x), length.out = 100)
y_new <- predict(mod, newdata = data.frame(x = x_new))
return(list(
data.frame(x = x_new, y = y_new, model_type = paste("Polynomial Regression (", degree, ")", sep = "")),
paste(paste("x^", degree, sep = ""), ":", paste(round(coefficients(mod), 2), collapse = " + "))
))
}
exp_reg <- function(x, y) {
mod <- lm(log(y) ~ x)
x_new <- seq(min(x), max(x), length.out = 100)
y_new <- exp(predict(mod, newdata = data.frame(x = x_new)))
return(list(
data.frame(x = x_new, y = y_new, model_type = "Exponential Regression"),
paste("y =", round(exp(coefficients(mod)[[2]]), 2), "* e^(", round(coefficients(mod)[[1]], 2), "x", ")")
))
}
log_reg <- function(x, y) {
mod <- lm(y ~ log(x))
x_new <- seq(min(x), max(x), length.out = 100)
y_new <- predict(mod, newdata = data.frame(x = x_new))
return(list(
data.frame(x = x_new, y = y_new, model_type = "Logarithmic Regression"),
paste("y =", round(coefficients(mod)[[2]], 2), "* log(x) +", round(coefficients(mod)[[1]], 2))
))
}
# Create regression line dataframes and equations
reg_data <- list(slr(x, y), mlr(x, y, rnorm(100, sd = 10)), poly_reg(x, y, 3), exp_reg(x, y), log_reg(x, y))
reg_eqs <- sapply(reg_data, function(x) x[[2]])
# Plot the data and regression lines
ggplot(data.frame(x = x, y = y), aes(x, y)) +
lapply(seq_along(reg_data), function(i) geom_line(data = reg_data[[i]][[1]], aes(x, y, color = reg_data[[i]][[1]]$model_type), linewidth = 1.4)) +
ggtitle("Regression Models") +
geom_point(alpha = .2) +
xlab("X") +
ylab("Y") +
scale_color_discrete(name = "Model Type") +
ylim(-6000, 15000) +
theme(legend.position = "none")
Example: climate change modelling
Example: COVID-19 transmission modelling
Source: https://www.nature.com/articles/s41598-021-84893-4/figures/1
Defining a linear relationship
- Pearson correlation coefficient measures the linear correlation between two variables
- Does not distinguish different patterns of association, only the strength of the association
- Not quite usable for predictive modelling, or for inference about the relationship between variables
In pictures…
The variation in the response
Hypothesis testing
- How does our null (H_0: \beta_1=0) model compare to the linear (H_0: \beta_1 \neq 0) model?
- Null model thinks the data can be summarised by the mean \bar{y}, and the linear model thinks the data can be summarised by the estimate \hat{y}.
Code
null_model <- Galton %>%
lm(child ~ 1, data = .) %>%
broom::augment(Galton)
lin_model <- Galton %>%
lm(child ~ parent, data = .) %>%
broom::augment(Galton)
models <- bind_rows(null_model, lin_model) %>%
mutate(model = rep(c("Null model", "SLR model"), each = nrow(Galton)))
ggplot(data = models, aes(x = parent, y = child)) +
geom_smooth(
data = filter(models, model == "Null model"),
method = "lm", se = FALSE, formula = y ~ 1, size = 0.5
) +
geom_smooth(
data = filter(models, model == "SLR model"),
method = "lm", se = FALSE, formula = y ~ x, size = 0.5
) +
geom_segment(
aes(xend = parent, yend = .fitted),
arrow = arrow(length = unit(0.1, "cm")),
size = 0.3, color = "darkgray"
) +
geom_point(alpha = .2) +
facet_wrap(~model) +
xlab("Parent height (in)") +
ylab("Child height (in)")
Thanks!
This presentation is based on the SOLES Quarto reveal.js template and is licensed under a Creative Commons Attribution 4.0 International License.
