---
title: "Linear Regression, Single Predictor"
subtitle: |
  | IMS2 Ch. 7 
  | Math 115
author: "Yurk"
format: 
  revealjs:
    theme: beige
    embed-resources: true
editor: source
metadata:
    license: CC BY-SA
---

## Body Measurements

```{r}
#| include: false
#| echo: false

library(tidyverse)
library(openintro)
library(kableExtra)
```

- `bdims` body measurement dataset, available [here](https://www.openintro.org/data/csv/bdims.csv)
- 507 physically active individuals (247 men, 260 women)
- `age`, weight (`wgt`), height (`hgt`), `sex`, 21 body girth variables (e.g., hip girth)

## Weight vs. Height

```{r}
#| include: true
#| echo: false
#| fig-cap: "Scatter plot of weight vs. height."

bdims |> ggplot(aes(x = hgt, y = wgt)) +
  geom_point() +
  theme_minimal() +
  theme(text = element_text(size = 18))
```

It appears that the data fall roughly along a line.

## Linear Model

```{r}
#| include: true
#| echo: false
#| fig-cap: "Scatter plot of weight vs. height with line of best fit."

bdims |> ggplot(aes(x = hgt, y = wgt)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal() +
  theme(text = element_text(size = 18))
```

We can add a *line of best fit* to the scatter plot.

---

- Equation for line: 
$$y = b_0 + b_1 x$$
- $b_0$ and $b_1$ are **coefficients**
   - $b_0$ = **intercept**
   - $b_1$ = **slope**
- $b_0$ and $b_1$ are *statistics* (fit using sample)
- $\beta_0$ and $\beta_1$ are the corresponding *parameters*
- The fitted values are $b_0=-105.0$, $b_1=1.018$

## Variable Roles

- `wgt` = **outcome/response** (dependent variable, $y$)
- `hgt` = **predictor** (independent variable, $x$) 
- We use a *hat* to indicate an estimate or prediction
$$\widehat{wgt} = -105.0 + 1.018 \times hgt$$

## Using a Model to Make Predictions

- Use the model to predict the weight of a person with a given height
- The predicted weight of a 170 cm tall individual is
$$\begin{array}{rcl}\widehat{wgt} &=& -105.0 + 1.018 \times hgt\\ &=& -105.0 + 1.018 \times 170 \\ &=& 68.06\, kg\end{array}$$

## Correlation

- The **correlation coefficient** describes *strength* and *direction* of a linear relationship
- Denoted $r$ for a sample, $\rho$ for a population
- $-1\leq r\leq1$

---

- Direction of linear relationship
   - $r>0$ indicates a *positive association*
  - $r<0$ indicates a *negative association*.
- Strength of linear relationship
   - Values close to 0 indicate a weak linear association
   - Values close to -1 or 1 indicate a strong linear association
   
---

![Some scatter plots and their correlations. IMS 2 Figure 7.10.](https://openintrostat.github.io/ims/model-slr_files/figure-html/fig-posNegCorPlots-1.png)

---

- Let $(x_i,y_i)$ be the $i$th observation of the numeric variables $x$ and $y$
- Then $r$ is
$$r=\frac{1}{n-1}\sum_{i=1}^n\frac{x_i-\bar{x}}{s_x}\cdot\frac{y_i-\bar{y}}{s_y}$$
- Here $\bar{x}$ and $\bar{y}$ are the sample means, and $s_x$ and $s_y$ are the sample standard deviations of the $x$ and $y$
- $r$ is independent of the units of measurement of $x$ and $y$

---

```{r}
#| include: true
#| echo: false
#| fig-cap: "Scatter plot of weight vs. height with line of best fit."

bdims |> ggplot(aes(x = hgt, y = wgt)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal() +
  theme(text = element_text(size = 18))
```

Correlation between height and weight: $r=0.717$

```{r}
#| include: false
#| echo: false

bdims |>
  summarize(n = n(),
            r = cor(hgt, wgt)) |>
    kbl(booktabs = TRUE, escape = FALSE) |>
    kable_styling(font_size = 22)
```

## Interpretation of coefficients

$$\widehat{wgt} = -105.0 + 1.018 \times hgt$$

- Slope: for each additional centimeter of height, we expect weight to increase by 1.018 kg
- Intercept: we would predict a 0 cm tall individual to weigh -105.0 kg
- In many cases, this intercept interpretation is not useful
- Better to think of intercept as positioning line vertically so it passes through the data cloud

## Extrapolation

- Predicting weight for individual with height outside of the range of the observed data is an example of **extrapolation**
- We should not expect the model to apply outside of this range
- Extrapolation can lead to nonsensical predictions (0 cm tall individuals with negative weight) or inaccurate ones

## Least Squares Regression

- How is the best fit line determined? 
- Slope and intercept chosen to minimize the error between the observed and predicted response

___

![Plot highlighting three residuals. IMS2 Figure 7.8. ](https://openintrostat.github.io/ims/model-slr_files/figure-html/fig-scattHeadLTotalLLine-highlighted-1.png)

The **residual** (error) for the $i$th observation $(x_i,y_i)$ is 
$$e_i = y_i - \hat{y}_i$$

## Least Squares Line

- The **least squares regression line** minimizes the sum of the squared residuals,
$$e_1^2+e_2^2+\cdots+e_n^2$$
- Properties of least squares line
  
  - The line passes through the point $(\bar{x},\bar{y})$
  - The slope is $b_1=\frac{s_y}{s_x}r$
- Note: We can use these properties to compute the slope and intercept if we know the means, SDs, and correlation

## Using Software 

- Typically we will use Jamovi or other Statistical software to compute the coefficients of the least squares line
- We will learn to do this in *J Lab 3*
- Results are usually presented in a **regression table**

```{r}
#| include: true
#| echo: false

library(jmv)

jmv::linReg(
    data = bdims,
    dep = wgt,
    covs = hgt,
    blocks = list(
        list(
            "hgt")),
    refLevels = list()) |>
    pluck("models")
```

## Coefficient of determination ($R^2$)

- The **coefficient of determination**, also known as **R-squared** ($R^2$) is used to measure how well a model describes the data
- $R^2$ is the proportion of variation in the outcome/response variable that is explained by the model
- For simple linear regression, $R^2 = r^2$

---

- $R^2$ will always have values between 0 and 1
- Value close to 1: linear model fits the data well (describes nearly 100% of the variability in outcomes)
- Value close to 0 indicates that it does not fit well

<!-- ## Total Sum of Squares -->

<!-- - **total sum of squares**, denoted *SST*, describes the total variation in the outcome -->
<!-- $$SST = (y_1-\bar{y})^2 + (y_2-\bar{y})^2 + \cdots + (y_n-\bar{y})^2$$ -->
<!-- - Note that *SST* does not involve the model at all -->
<!-- - However, can think of a **null model** that uses the sample mean as the prediction  -->
<!-- - *SST* is the sum of the squared residuals for the null model -->

<!-- ## Sum of Squared Errors -->

<!-- - **sum of squared errors**, denoted *SSE*, quantifies the variation in outcomes that the model *fails* to describe -->
<!-- $$\begin{array}{rcl}SSE &=& (y_1-\hat{y}_1)^2 + (y_2-\hat{y}_2)^2 + \cdots + (y_n-\hat{y}_n)^2 \\ &=& e_1^2 + e_2^2 + \cdots + e_n^2\end{array}$$ -->
<!-- - Given by the sum of the squared residuals, which we have encountered before -->

<!-- ## Regression Sum of Squares -->

<!-- - **regression sum of squares**, denoted **SSR**, measures the variation that *is* accounted for by the model  -->
<!-- $$SSR = SST - SSE$$ -->
<!-- - Hence, the proportion of variation in the outcome that is described by the model is -->
<!-- $$R^2 = \frac{SST - SSE}{SST} = 1 - \frac{SSE}{SST}$$ -->

---

- We can compute $R^2$ using Jamovi or other Statistical Software
- Height explains about 51.5% of the variability in weights

```{r}
#| include: true
#| echo: false
#| 
jmv::linReg(
    data = bdims,
    dep = wgt,
    covs = hgt,
    blocks = list(
        list(
            "hgt")),
    refLevels = list()) |>
    pluck("modelFit")
```


## Residual plots

- **residual plot** is a plot of residuals vs. predicted values (scatter plot with points $(\hat{y}_i,e_i)$
- Useful for diagnosing problems with the linear models
- If there is a pattern in the residual plot, then a more complicated model (e.g., a nonlinear model or a model that includes more predictors) may be more appropriate

---

There are no obvious patterns in the height vs. weight residual plot.

```{r}
#| include: true
#| echo: false
#| fig-cap: "Residual plot for weight vs. height with horizontal line at $e=0$ for reference."

library(broom)

lm1 <- lm(wgt ~ hgt, data = bdims)

bdims_aug <- augment(lm1, bdims)

bdims_aug |>
  ggplot(aes(x = .fitted, y = .resid)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  theme_minimal() +
  theme(text = element_text(size = 18))
```
---

More residual plots

![Some scatter plots (top) and corresponding residual plots (bottom). From IMS2 Example. ](https://openintrostat.github.io/ims/model-slr_files/figure-html/sampleLinesAndResPlots-1.png)
---

More residual plots

![More residual plots. From IMS2 Exercise. 7.2](https://openintrostat.github.io/ims/model-slr_files/figure-html/unnamed-chunk-39-1.png)

## Outliers

- *outliers* are observations that fall far from the point cloud
- **high leverage points** fall horizontally far from the center of the point cloud
- high leverage points have more pull on the regression line
- **influential points** have a strong influence on the slope of the regression line
- influential points can be identified by fitting a line with the point removed. If the slope is very different than when the point is included, then the point is influential.

---

Each of the following plots has an outlier. Which are high leverage? Influential?

![Scatter plots with outliers. From IMS2 Fig. 7.16.](https://openintrostat.github.io/ims/model-slr_files/figure-html/fig-outlier-plots-1.png)