You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
```{r include = FALSE}
CWD <- read.table("./datasets/01-linear-regression/christ.csv", header = T, sep = ",", dec = ".")
model <- lm(CWD.BASA ~ RIP.DENS, data = CWD)
```
# Introduction
## Regression
Goal: describe the relationship between 2 series of observations $(X_i,Y_i)$, obtained for individual subjects $i=1,...,n$
Example:
Basal area of coarse woody debris (CWD) versus tree density along 16 North American lakes
- **Dependent variable, outcome, response** $Y$: CWD basal area
- **Independent variable, explanatory variable, predictor** $X$: tree density (in number per km)
{width="33%" fig-align="center"}
## CWD versus tree density
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 4.5
#| fig-align: center
plot(CWD$RIP.DENS, CWD$CWD.BASA, xlab = "Riparian tree density", ylab = "CWD basal area", pch = 19)
lines(lowess(CWD$RIP.DENS, CWD$CWD.BASA, f=3/4), lty = 2)
abline(model)
```
## Regression
For fixed $X$, $Y$ will be some function of $X$ plus random noise: **observation = signal + noise**
Mathematical modelling of observation:
$$
Y_i=f(X_i)+\epsilon_i
$$
where $f(x)$ is the expected outcome for subjects with $X_i=x$
$$
E(Y_i|X_i=x) = f(x)
$$
and $\epsilon_i$ is on average 0 for subjects with same $X_i$:
$$
E(\epsilon_i|X_i) = 0.
$$
## Linear regression
- To obtain accurate and interpretable results, $f(X)$ is often chosen as linear function of unknown parameters
- Use **linear regression model**
$$
E(Y|X=x)=\alpha + \beta x
$$
with unknown **intercept** $\alpha$ and **slope** $\beta$.
- Linear regression model makes assumption on distribution of $X$ and $Y$, so can be incorrect.
## Use of linear regression
- **Prediction**: when $Y$ unknown, but $X$ known, we can predict $Y$ based on $X$:
$$
E(Y|X=x)=\alpha + \beta x.
$$
- **Association**: describe biological relation between variable $X$ and continuous measurement $Y$
- Slope $\beta$: difference in mean outcome between subjects that differ 1 unit in the value of $X$:
\begin{align*}
E(Y|X=x+\delta) - E(Y|X=x) & = \alpha + \beta (x+\delta)
-\alpha-\beta x\\
& = \beta\delta.
\end{align*}
## Least squares estimates
- Least squares (regression) line: line that `best' fits data.
- Found by choosing values for $\alpha$ and $\beta$ that minimize sum of squares of **residuals**:
$$
\sum_{i=1}^n (\underbrace{Y_i-\alpha-\beta X_i}_{\text{Residual}})^2
$$
- Estimates for $\beta$ and $\alpha$:
$$
\hat{\beta} = \mathrm{Cor}(x,y) \frac{S_y}{S_x} \quad \text{and} \quad \hat{\alpha}=\bar Y - \hat{\beta} \bar X.
$$
with $\mathrm{Cor}(x, y)$ the sample correlation between $x$ and $y$ and $S_x$, $S_y$ the sample standard deviation.
## Residuals plot
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 4.5
#| fig-align: center
a <- -77.10
b <- 0.12
x <- CWD$RIP.DENS
y <- CWD$CWD.BASA
plot_residuals <- function(a, b) {
plot(y ~ x, xlab = "", ylab = "", pch = 19)
abline(a, b)
y_pred <- a + b*x
segments(x, y_pred, x, y, lwd=2)
}
plot_residuals(a, b)
text(1750, 100, paste("y =", a, "+", b, "x"))
```
See also: [residuals animation](https://yihui.org/animation/example/least-squares/).
## Output linear regression (coefficients only)
\footnotesize
```{r}
#| echo: true
model <- lm(CWD.BASA ~ RIP.DENS, data = CWD)
summary(model)$coefficients
```
\normalsize
Regression line:
$$
E(Y|X=x)=-77.10+0.12 x
$$
## Output linear regression (full)
\footnotesize
```{r}
#| echo: true
summary(model)
```
\normalsize
## Interpreting linear regression
- Model: $E(Y|X=x)=-77.10+0.12 x$
- Expected CWD basal area is 1.2 m$^2$ larger alongside lakes with 10 more trees per km
- Expected CWD basal area alongside lakes with 1,600 trees per km shoreline:
$$
-77.10 +0.12\times 1600=108 \ {\rm m}^2
$$
- Expected CWD basal area alongside lakes with 500 trees per km shoreline:
$$
-77.10 +0.12\times 500=-17 \ {\rm m}^2
$$
- **Be careful with extrapolation!** (linearity assumption can only be verified within range of data)
# Assumptions for linear regression
## Verifying linearity assumption
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 4.5
#| fig-align: center
plot(CWD$RIP.DENS, CWD$CWD.BASA, xlab = "Riparian tree density", ylab = "CWD basal area", pch = 19)
lines(lowess(CWD$RIP.DENS, CWD$CWD.BASA, f=3/4), lty = 2)
abline(model)
```
## Verifying linearity assumption
- An alternative (more convenient when there are multiple predictors) is a **residual plot**.
- Note: residuals are prediction errors:
$$
e_i = y_i-\hat{\alpha}-\hat{\beta}x_i
$$
- If linear model correct, then scatterplot of $e_i$ versus $x_i$ or predictions $\hat{\alpha}+\hat{\beta}x_i$ shows no pattern
```{r eval = FALSE}
model <- lm(CWD.BASA ~ RIP.DENS, data = CWD)
par(mfrow=c(2,2))
plot(model)
```
## Verifying linearity assumption
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 4.5
#| fig-align: center
model <- lm(CWD.BASA ~ RIP.DENS, data = CWD)
par(mfrow=c(2,2))
plot(model)
```
## Inference for simple linear regression
To be able to draw conclusions about the linear regression model
$$
E(Y|X)=\alpha+\beta X
$$
we need extra assumptions:
- **Homoscedasticity**: for fixed $X$, $Y$ has constant variance
$$
\text{Var}(Y|X)=\sigma^2,
$$
estimated by the residual mean square error:
$$
\text{MSE}=\sum_{i=1}^n e_i^2/(n-2)
$$
- **Normality**: for fixed $X$, $Y$ is normally distributed
$$
Y|X \sim N(\alpha+\beta X, \sigma^2)
$$
## Homoscedasticity versus heteroscedasticity
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 4.5
#| fig-align: center
set.seed(1234)
n <- 100
a <- -1
b <- 0.5
x <- seq(-3, 3, length.out = n)
y1 <- a + b*x + rnorm(n, sd = 0.5)
sd <- seq(1, 2, length.out = n)^3
y2 <- a + b*x + rnorm(n, sd = sd)
par(mfrow=c(1, 2))
plot(x, y1, xlab = "", ylab = "", main = "Homoscedasticity")
abline(a, b)
plot(x, y2, xlab = "", ylab = "", main = "Heteroscedasticity")
abline(a, b)
```
## Homoscedasticity?
Hard to check on regression plot directly!
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 4.5
#| fig-align: center
plot(CWD$RIP.DENS, CWD$CWD.BASA, xlab = "Riparian tree density", ylab = "CWD basal area", pch = 19)
lines(lowess(CWD$RIP.DENS, CWD$CWD.BASA, f=3/4), lty = 2)
abline(model)
```
## Assumption of homoscedasticity
- Squared residuals carry information on residual variability.
- If these are associated with explanatory variable, then indication of **heteroscedasticity**.
- Scatterplot of $e_i^2$ or $\sqrt{|e_i|}$ versus $x_i$ or predictions.
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 4
#| fig-align: center
plot(CWD$RIP.DENS, abs(residuals(model))^0.5, xlab = "Riparian tree density", ylab = "Square root of residuals", pch = 19)
lines(lowess(CWD$RIP.DENS, abs(residuals(model))^0.5, f=3/4))
```
## Assumption of homoscedasticity
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 5
#| fig-align: center
model <- lm(CWD.BASA ~ RIP.DENS, data = CWD)
par(mfrow=c(2,2))
plot(model)
```
## Normality assumption
- Assumption: outcomes normally distributed **for fixed values of explanatory variable**:
$$
Y|X \sim N(\alpha + \beta X, \sigma^2).
$$
- Can be checked using QQ-plot of the residuals.
## Normality assumption valid
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 4.5
#| fig-align: center
set.seed(456)
n <- 200
a <- 0
b <- 2
x <- seq(1, 10, length.out = n)[1:(n-1)]
y <- a + b*x + rnorm(n - 1)
bin <- floor(x)
par(mfrow=c(1, 2))
plot(x, y, xlab = "", ylab = "")
boxplot(y ~ bin, xlab = "", ylab = "")
residuals <- y - a - b*x
```
## Normality assumption not valid
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 4.5
#| fig-align: center
set.seed(456)
n <- 200
a <- 0
b <- 2
x <- seq(1, 10, length.out = n)[1:(n-1)]
e <- rnorm(n - 1)
scales <- seq(1, 2, length.out = n-1)^3
scales[e < 0] <- 1
e <- e*scales
y <- a + b*x + e
bin <- floor(x)
par(mfrow=c(1, 2))
plot(x, y, xlab = "", ylab = "")
boxplot(y ~ bin, xlab = "", ylab = "", side="right")
```
## QQ plot of residuals (Y|X normal)
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 4.5
#| fig-align: center
set.seed(456)
n <- 100
a <- 0
b <- 2
x <- seq(1, 10, length.out = n)[1:(n-1)]
y <- a + b*x + rnorm(n - 1)
bin <- floor(x)
par(mfrow=c(1, 2))
boxplot(y ~ bin, xlab = "", ylab = "", side="right")
residuals <- y - a - b*x
qqnorm(residuals, main = "Q-Q plot of Y|X")
qqline(residuals)
```
## QQ plot of residuals (Y|X not normal)
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 4.5
#| fig-align: center
set.seed(456)
n <- 200
a <- 0
b <- 2
x <- seq(1, 10, length.out = n)[1:(n-1)]
e <- rnorm(n - 1)
scales <- seq(1, 2, length.out = n-1)^3
scales[e < 0] <- 1
e <- e*scales
y <- a + b*x + e
bin <- floor(x)
par(mfrow=c(1, 2))
boxplot(y ~ bin, xlab = "", ylab = "", side="right")
residuals <- y - a - b*x
qqnorm(residuals, main = "Q-Q plot of Y|X")
qqline(residuals)
```
## Do not use QQ-plot of Y!
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 4.5
#| fig-align: center
set.seed(456)
n <- 100
a <- 0
b <- 2
x <- seq(1, 10, length.out = n)[1:(n-1)]
y <- a + b*x + rnorm(n - 1)
bin <- floor(x)
par(mfrow=c(1, 2))
residuals <- y - a - b*x
qqnorm(residuals, main = "Q-Q plot of Y|X")
qqline(residuals)
qqnorm(y, main = "Q-Q plot of Y")
qqline(y)
```
## Checking for normality with R diagnostic plots
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 5
#| fig-align: center
model <- lm(CWD.BASA ~ RIP.DENS, data = CWD)
par(mfrow=c(2,2))
plot(model)
```
## What if homoscedasticity or normality false?
- Transformation of **dependent variable** can help to obtain normality and homoscedasticity.
- Example transformations: $\sqrt{Y}, Y^2, 1/Y, \ln{Y}$.
- Transformation of **independent variable** does not change distribution of $Y$ for given $X$:
- does not help in obtaining normality or homoscedasticity.
- does help to obtain linearity if normality and homoscedasticity are ok.
## What if homoscedasticity or normality false?
- Often because outcome can only take on values in certain interval (e.g. $[0, 1]$, positive numbers, ...)
- **Solution**: transform outcome such that it can take on all real values
- Example: `CWD.BASA` is always positive: take $\ln$ to make outcome real-valued:
```{r eval=FALSE}
model2 <- lm(I(log(CWD.BASA)) ~ RIP.DENS, data = CWD)
summary(model2)
```
## Transforming the outcome
\footnotesize
```{r}
#| echo: true
model2 <- lm(log(CWD.BASA) ~ RIP.DENS, data = CWD)
summary(model2)
```
## Residual plots
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 5
#| fig-align: center
par(mfrow=c(2,2))
plot(model2)
```
# Higher-order regression
## What if linearity assumption is false?
- Transformation of dependent variable
- Transformation of independent variable
- If residuals reveal **quadratic association**, such that
$$
e_i\approx \delta_0+\delta_1 x_i+\delta_2 x_i^2
$$
then
$$
y_i=\hat{\alpha}+\hat{\beta}x_i+e_i\approx
(\hat{\alpha}+\delta_0)+(\hat{\beta}+\delta_1)x_i+\delta_2
x_i^2
$$
## Quadratic regression
- We assume
$$
E(Y|X)=\alpha+\beta X+\gamma X^2
$$
- Unknown parameters estimated by **least squares method**: minimize
$$
\sum_{i=1}^n (Y_i-\alpha-\beta X_i-\gamma X_i^2)^2
$$
## Quadratic regression
\footnotesize
```{r}
model3 <- lm(log(CWD.BASA) ~ RIP.DENS + I(RIP.DENS^2), data = CWD)
summary(model3)
```
## Residual plots
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 5
#| fig-align: center
par(mfrow=c(2,2))
plot(model3)
```
## Building model proceeds hierarchically
\small
```{r}
summary(model3)$coefficients
```
\normalsize
- Add terms to model and keep those as long as they are significantly associated with outcome
- Example: adding third order term is not significant contribution (p-value 0.26)
- Adding proceeds **hierarchically**: lower order terms are kept as long as higher order terms are in model
## Results
- We conclude
$$
E\{\ln(Y)|X\}=-9.69+0.017X-4.96 \ 10^{-6}X^2
$$
or equivalently that geometric mean CWD basal area for given tree density $X$ is equal to
$$
\exp(-9.69+0.017X-4.96 \ 10^{-6}X^2)
$$
- For $X=500$ we now find 0.086 m$^2$ (previously: -17 m$^2$)
- **How precise is this?**
# Interpreting the results of a regression model
## Inference for mean outcome
Given an input $X = x_h$, what do we expect the outcome $Y$ to be on average?
Use the regression equation:
- $\hat{y}_h = \hat{\alpha} + \hat{\beta} x_h$ is unbiased estimator of $E(Y|X = x_h) = \alpha+\beta x_h$.
What is the uncertainty of this estimator?
For this we need the standard error of $\hat{Y}_h$.
## Inference for mean outcome: uncertainty
- Standard error of $\hat{Y}_h$ is
$$
SE(\hat{Y}_h)=\sqrt{MSE\left\{\frac{1}{n}+\frac{(X_h-\bar X)^2}{\sum_i (X_i-\bar X)^2}\right\}}.
$$
- Tests and CI for $E(Y|X_h)$ based on
$$
\frac{\hat{Y}_h-E(Y|X_h)}{SE(\hat{Y}_h)}\sim t_{n-p}
$$
with $p$ number of unknown parameters in model.
## Inference for mean outcome: intuition
- Highest precision for predictions in $X_h=\bar X$: relative confidence in predictions for $X$ **close to mean**.
- Lower precision as predictions have $X$ further **away from the mean**.
## Prediction in R
\small
```{r}
#| echo: true
model3 <- lm(I(log(CWD.BASA)) ~ RIP.DENS + I(RIP.DENS^2),
data = CWD)
p <- predict(model3,
newdata = data.frame(RIP.DENS=c(1000, 1500)),
interval = "confidence")
p
```
## Prediction in R
Predictions and lower/upper bound of the CI are for `log(CWD.BASA)` and need to be transformed back:
```{r}
#| echo: true
exp(p)
```
## Expected outcome with 95\% CI
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 4.5
#| fig-align: center
p <- predict(model3, newdata = data.frame(RIP.DENS=800:2200),
interval = "confidence")
plot(CWD$RIP.DENS, CWD$CWD.BASA, xlab = "Riparian tree density", ylab = "CWD basal area", ylim = c(0, max(exp(p[,3]))))
abline(model, lty=3)
lines(800:2200, exp(p[,1]))
lines(800:2200, exp(p[,2]), lty=2)
lines(800:2200, exp(p[,3]), lty=2)
```
## Inference for slope $\beta$
- The regression coefficient $\hat{\beta}$ is an (unbiased) estimator of $\beta$, the population regression coefficient.
- It comes with a measure of uncertainty: standard error of $\hat{\beta}$:
$$
SE(\hat{\beta})=\sqrt{\frac{MSE}{\sum_i (X_i-\bar X)^2}}.
$$
with $MSE=\frac{1}{n-2}\sum_{i=1}^n(Y_i-\hat{Y}_i)^2$
- Large spread on $X$ improves precision.
## Spread and precision
```{r}
#| out-width: 4.5in
#| fig-width: 6
#| fig-height: 4.5
#| fig-align: center
set.seed(1234)
n <- 50
a <- 0
b <- 0.5
x1 <- seq(-3, 3, length.out = n)
y1 <- a + b*x1 + rnorm(n, sd = 0.5)
x2 <- x1 / 3
y2 <- a + b*x2 + rnorm(n, sd = 0.5)
m1 <- lm(y1 ~ x1)
m2 <- lm(y2 ~ x2)
se1 <- summary(m1)$coefficients["x1", "Std. Error"]
se2 <- summary(m2)$coefficients["x2", "Std. Error"]
par(mfrow=c(1, 2))
plot(x1, y1, xlim = c(-3, 3), xlab = "", ylab = "",
main = paste0("SE(beta) = ", round(se1, digits = 2)))
abline(m1)
plot(x2, y2, xlim = c(-3, 3), xlab = "", ylab = "",
main = paste0("SE(beta) = ", round(se2, digits = 2)))
abline(m2)
```
## Association tree density vs.\ CWD
Tests and confidence intervals for $\beta$ are based on
$$
\frac{\hat{\beta}-\beta}{SE(\hat{\beta})}\sim t_{n-2}
$$
\small
```{r}
#| echo: true
summary(model)$coefficients
```
## Association tree density vs.\ CWD
- 95\% CI for $\beta$ needs $t_{14,0.975} = 2.14$
- CI is given by
$$
[0.116 - 2.14\times 0.0234,0.116 + 2.14\times 0.0234]=[0.066,0.166]
$$
\vfill
\small
```{r}
#| echo: true
confint(model)
```