Statistical Methods with R

Simple Regression

Unit C · Chapter 07 · Lecture 07a

Developed by Jeffrey M. Girard

Roadmap: Regression

Introduction
Regression Inference
Centering
Regression Effect Sizes

Introduction

Regression as Fancy Correlation

A correlation quantifies how well a straight line represents the relationship between two continuous variables

Regression as Fancy Correlation

A correlation quantifies how well a straight line represents the relationship between two continuous variables

Regression as Fancy Correlation

But how did we know where to draw the line?
Many different lines could have been drawn…

Regression as Fancy Correlation

Regression is a way to select and quantify the “best” line

Defining a Line in General

y = \color{blue}{m}x + \color{blue}{b}

y and x are variables

\color{blue}m and \color{blue}b are coefficients

\color{blue}m is the slope
Change in y when x+1
\color{blue}b is the intercept
Value of y when x=0

Defining a Regression Line

If y is the outcome and x is the predictor…
We can define the regression line in the same way

\color{red}{\hat{y}_i} = \color{blue}{b_0} + \color{blue}{b_1}x_i

\color{red}{\hat{y}_i} is the model-predicted value of y for observation i

\color{blue}{b_0} is the intercept of the line (i.e., b on the previous slide)

\color{blue}{b_1} is the slope of the line (i.e., m on the previous slide)

Defining Regression Residuals

Our regression line is defined by the formula:

\color{red}{\hat{y}_i} = \color{blue}{b_0} + \color{blue}{b_1}x_i

But many points aren’t exactly on the line, so our predictions aren’t perfect

Defining Regression Residuals

Thus, predictions (\color{red}{\hat{y}_i}) differ from the observed values (y_i)
We refer to these differences as the model’s residuals

\begin{align} \color{red}{\hat{y}_i} &= \color{blue}{b_0} + \color{blue}{b_1}x_i\\ y_i &= \color{blue}{b_0} + \color{blue}{b_1}x_i + \color{magenta}{e_i} \\ \color{magenta}{e_i} &= y_i - \color{red}{\hat{y}_i} \end{align}

Visualizing Regression Residuals

Minimizing Regression Residuals

If we want the best possible estimates of y, we should…
…look for the line that minimizes the prediction errors
Residuals are the prediction errors
We thus seek to minimize the residuals

But residuals can be positive or negative…
So minimizing the sum of the residuals won’t work
We thus minimize the sum of the squared residuals

OLS Regression

Minimizing the sum of the squared residuals is called OLS
- This stands for “Ordinary Least Squares” Regression
So what are the OLS formulas for the slope and intercept?

b_1 = r_{xy}\frac{s_y}{s_x} \qquad\quad b_0 = \bar{y} - b_1\bar{x}

Regression and correlation are thus integrally connected

Example Dataset

The water dataset has information about mortality and water hardness in 61 UK towns

water <- read_csv("water.csv")
water

# A tibble: 61 × 4
   town        location mortality hardness
   <chr>       <chr>        <dbl>    <dbl>
 1 Bath        South         1247      105
 2 Birkenhead  North         1668       17
 3 Birmingham  South         1466        5
 4 Blackburn   North         1800       14
 5 Blackpool   North         1609       18
 6 Bolton      North         1558       10
 7 Bootle      North         1807       15
 8 Bournemouth South         1299       78
 9 Bradford    North         1637       10
10 Brighton    South         1359       84
# ℹ 51 more rows

Data Exploration

OLS Regression by Hand

x <- water$hardness
y <- water$mortality

b1 <- cor(x, y) * sd(y) / sd(x)
b1
## [1] -3.226092

b0 <- mean(y) - b1 * mean(x)
b0
## [1] 1676.356

\hat{y}_i = 1676.356 - 3.226 x_i

“We estimate that a town with a water hardness of zero would have a mortality of 1676.356. We estimate that towns with one more unit of water hardness have mortality that is lower by 3.226, on average.”

Visualizing our Result

qplot(x = hardness, y = mortality, data = water, geom = "point") + 
  geom_abline(slope = b1, intercept = b0, color = "salmon")

Regression using lm()

fit <- lm(formula = mortality ~ hardness, data = water)
fit
## 
## Call:
## lm(formula = mortality ~ hardness, data = water)
## 
## Coefficients:
## (Intercept)     hardness  
##    1676.356       -3.226

Just like before, the formula is outcome ~ predictor
R refers to b_0 as (Intercept) and b_1 as hardness
This will be helpful when we add multiple predictors

Regression Inference

Uncertainty in the Population

OLS finds the best coefficients for the sample
But we usually care more about the population

Thus, we estimate the population intercept and slope
There will be some inferential uncertainty in this process

As before, we can use NHST and/or estimation approaches
Let’s calculate 95% CIs, test statistics, and p-values

Slope Confidence Interval

b_1 \pm \text{SE}(b_1)\times t_{crit}

\text{SE}(b_1) = \frac{s_y}{s_x}\sqrt{\frac{1-r_{xy}^2}{n-2}}

t_{crit}\text{ comes from the } t(n-2, 0, 1)\text{ distribution}

Intercept Confidence Interval

b_0 \pm \text{SE}(b_0)\times t_{crit}

\text{SE}(b_0) = \sqrt{\frac{\sum_i(y_i-\hat{y}_i)^2}{n-2}} \sqrt{\frac{1}{n} + \frac{\bar{x}^2}{(n-1)s_x^2}}

t_{crit}\text{ comes from the } t(n-2, 0, 1)\text{ distribution}

Confidence Intervals in R

fit
## 
## Call:
## lm(formula = mortality ~ hardness, data = water)
## 
## Coefficients:
## (Intercept)     hardness  
##    1676.356       -3.226

confint(fit)
##                  2.5 %      97.5 %
## (Intercept) 1617.73030 1734.980903
## hardness      -4.19602   -2.256163

I’ll show you the {easystats} way after we cover NHST

Plotting Regression Results

library(modelbased)
est <- estimate_relation(fit, by = "hardness")
plot(est, show_data = TRUE)

Hypothesis Testing

Null: The true coefficient (in the population) is zero

H_0: b = 0

Alternative: The true coefficient (in the population) is not zero

H_1: b \neq 0

Note that this applies to both intercepts and slopes

Testing Individual Coefficients

We can re-use the logic/procedure of the one-sample t-test
We assume the null sampling distribution for \hat{b} is normal

t_{obs} = \frac{b}{\text{SE}(b)}

t_{obs} \sim t(\mathit{df}=n-k-1, \mu=0, \sigma=1)

Statistical Inference in R

library(easystats)
model_parameters(fit)
## Parameter   | Coefficient |    SE |             95% CI | t(59) |      p
## -----------------------------------------------------------------------
## (Intercept) |     1676.36 | 29.30 | [1617.73, 1734.98] | 57.22 | < .001
## hardness    |       -3.23 |  0.48 | [  -4.20,   -2.26] | -6.66 | < .001

“A town with zero water hardness is expected to have a mortality rate of 1676.36, which is significantly different from zero, 95% CI: [1617.73, 1734.98], p<.001. For each one-unit increase in water hardness, the model predicts a significant decrease in mortality of 3.23, 95% CI: [–4.20, –2.26], p<.001.”

Centering

Centering a Predictor

The intercept (b_0) is the value of \hat{y} when x=0
This isn’t so useful if x never equals 0 in the sample

In such cases, it can be useful to “center” x before modeling

x_i^c = x_i - \bar{x}

The slope will stay the same but the intercept will change
The intercept will become the value of \hat{y} when x=\bar{x}

Centering a Predictor

b_0 = 1676.356

“A town with zero water hardness
would have 1676.356 mortality.”

b_0 = 1524.148

“A town with average water hardness
would have 1524.148 mortality.”

Centering in R

# Two-step process
water$hardness_c <- water$hardness - mean(water$hardness)

fit2a <- lm(formula = mortality ~ hardness_c, data = water)
model_parameters(fit2a)
## Parameter   | Coefficient |    SE |             95% CI | t(59) |      p
## -----------------------------------------------------------------------
## (Intercept) |     1524.15 | 18.31 | [1487.50, 1560.79] | 83.23 | < .001
## hardness c  |       -3.23 |  0.48 | [  -4.20,   -2.26] | -6.66 | < .001

# Or equivalent one-step process
fit2b <- lm(formula = mortality ~ center(hardness), data = water)
model_parameters(fit2b)
## Parameter        | Coefficient |    SE |             95% CI | t(59) |      p
## ----------------------------------------------------------------------------
## (Intercept)      |     1524.15 | 18.31 | [1487.50, 1560.79] | 83.23 | < .001
## center(hardness) |       -3.23 |  0.48 | [  -4.20,   -2.26] | -6.66 | < .001

Regression Effect Sizes

Standardized Slopes

It is common to standardize regression slopes (called \betas)
We can do this by standardizing both y and x and refitting

x_i^z = \frac{x_i - \bar{x}}{s_x} \qquad \quad y_i^z = \frac{y_i - \bar{y}}{s_y}

\hat{y}_i^z = \beta_0 + \beta_1x_i^z

The standardized slope (\beta_1) is interpreted in SD units
The standardized intercept (\beta_0) will equal/approach zero

Standardizing in R

# Method 1: Fit normally then refit
fit <- lm(mortality ~ hardness, data = water)
model_parameters(fit, standardize = "refit")
## Parameter   | Coefficient |   SE |         95% CI |     t(59) |      p
## ----------------------------------------------------------------------
## (Intercept) |   -3.99e-16 | 0.10 | [-0.20,  0.20] | -4.09e-15 | > .999
## hardness    |       -0.65 | 0.10 | [-0.85, -0.46] |     -6.66 | < .001

# Method 2: Standardize data then fit
fit2 <- lm(mortality ~ hardness, data = standardize(water))
model_parameters(fit2)
## Parameter   | Coefficient |   SE |         95% CI |     t(59) |      p
## ----------------------------------------------------------------------
## (Intercept) |   -3.99e-16 | 0.10 | [-0.20,  0.20] | -4.09e-15 | > .999
## hardness    |       -0.65 | 0.10 | [-0.85, -0.46] |     -6.66 | < .001

“For each one SD increase in water hardness, we would expect a 0.65 SD
decrease in mortality, \beta=-0.65, 95% CI: [-0.85, -0.46], p<.001.”

R-Squared

R^2 is the variance in y explained by the regression model
- Ranges from 0.0 (0% explained) to 1.0 (100% explained)
- It is equivalent to squaring the correlation of y and \hat{y}
- It is also called the “coefficient of determination”

R^2 = 1 - \frac{\text{SS}_{res}}{\text{SS}_{tot}} = 1 - \frac{\sum_i(y_i - \hat{y}_i)^2}{\sum_i(y_i - \bar{y})^2}

Adjusted R-Squared

R^2 will never decrease when adding more predictors
So we often adjust it to penalize adding more predictors

R_{adj}^2 = 1 - \left(\frac{\text{SS}_{res}}{\text{SS}_{tot}} \times \frac{n - 1}{n - k - 1}\right)

Negative values of R_{adj}^2 are possible
Its interpretation is less straightforward

Model Performance in R

# By hand (with the residuals function shortcut)
R2a <- 1 - sum(residuals(fit)^2) / sum((y - mean(y))^2)
R2a
## [1] 0.4288267

# By function (recommended)
R2b <- r2(fit, ci = 0.95)
R2b
##        R2: 0.429 [0.230, 0.601]
##   adj. R2: 0.419 [0.220, 0.593]

“The regression model in which mortality was regressed on water hardness explained 42.9% of the variance in mortality, R² = 0.429, 95% CI: [0.230, 0.601], Adjusted R² = 0.419.”