37-Regression.Rmd

# Regression {#Regression}


<!-- Introductions; easier to separate by format -->
```{r, child = if (knitr::is_html_output()) {'./introductions/37-Regression-HTML.Rmd'} else {'./introductions/37-Regression-LaTeX.Rmd'}}
```


```{r}
data(RedDeer)
```


## Introduction {#Chap35-Intro}

The last chapter introduced *correlation*, which measures the *strength* of the *linear* relationship between two quantitative variables $x$ (an explanatory variable) and $y$ (a response variable).
We now study *regression*, which describes *what* the linear relationship is between $x$ and $y$.
The relationship is described using an *equation*, which allows us to:

1. *Predict* values of $y$ from given values of $x$ (Sect. \@ref(RegressionForPrediction)); and
2. *Understand* the relationship between $x$ and $y$ (Sect. \@ref(RegressionForUnderstanding)).

An example of a linear regression equation, describing the relationship between an explanatory variable $x$ and a response variable $y$, is    
\[
  \hat{y} = -4  + (2\times x), \qquad\text{usually written}\qquad \hat{y} = -4  + 2x.
\]
The notation $\hat{y}$ refers to the *mean predicted* value of the response variable given a value of $x$.


::: {.importantBox .important data-latex="{iconmonstr-warning-8-240.png}"}
$y$ refers to the observed data, whose values are known.
$\hat{y}$ refers to *mean predicted* (not observed) values of $y$ for some given value of $x$.
:::



:::: {.pronounceBox .pronounce data-latex="{iconmonstr-microphone-7-240.png}"}
$\hat{y}$ is pronounced as 'why-hat'; the 'caret' above the $y$ is called a 'hat', and designates a *predicted* value (of $y$).
:::


More generally, the equation of a straight line is  
\[
   \hat{y} = b_0 + (b_1 \times x), \qquad\text{usually written}\qquad \hat{y} = b_0 + b_1 x
\]
where $b_0$ and $b_1$ are (unknown) numbers estimated from sample data.
Notice that $b_1$ is the number multiplied by $x$.
In the example above, $b_0 = -4$ and $b_1 = 2$.


::: {.example #RegressionEquations name="Regression equations"}
If $\hat{y} = 15  - 102x$, then $b_0 = 15$ and $b_1 = -102$.
:::


::: {.thinkBox .think data-latex="{iconmonstr-light-bulb-2-240.png}"}
Consider the regression equation $\hat{y} = -0.0047x + 2.1$.\label{thinkBox:RegressionCoefs}
What are the values of $b_0$ and $b_1$?
(Be careful!)

`r if (knitr::is_latex_output()) '<!--'`
`r webexercises::hide()`
$b_0 = 2.1$ and $b_1 = -0.0047$.
Remember: $b_1$ is the number multiplied by the $x$, and that $b_0$ is the number by itself.

Some software and some journal papers write regression equations with the value of $b_0$ first, and some with the value of $b_1$ first. 
The *order* is not important.
`r webexercises::unhide()`
`r if (knitr::is_latex_output()) '-->'`
:::


## Linear equations: review {#RegressionEquationsReview}

To introduce, or revise, the idea of a linear equation, consider the (artificial) data in 
`r if (knitr::is_latex_output()) {
   'Fig. \\@ref(fig:ExampleScatterLATEX) (left panel),'
} else {
   'Fig. \\@ref(fig:ExampleScatterHTML),'
}`
with an explanatory variable $x$ and a response variable $y$.
In the regression equation $\hat{y} = b_0 + b_1 x$, the numbers $b_0$ and $b_1$ are called *regression coefficients*, where

* $b_0$ is called the *intercept*.
  Its value corresponds to the *predicted* mean value of $y$ when $x = 0$.
* $b_1$ is the *slope*.
  It measures how much the value of $y$ changes, on average, when the value of $x$ *increases* by 1.

We will use software to find the values of $b_0$ and $b_1$ (the formulas are tedious to use).

A rough approximation of the value of the *intercept* can be made by drawing a sensible-looking straight line through the scatterplot of the data, and determining what that line predicts for the value of $y$ when $x = 0$.
Also, a rough approximation of the *slope* can be made by computing  
\[
   \text{approximate slope} 
   = \frac{\text{Change in $y$}}{\text{Corresponding change in $x$}}
   = \frac{\text{rise}}{\text{run}}.
\]
This approximation of the slope is the *change* in the value of $y$ (the 'rise') divided by the corresponding *change* in the value of $x$ (the 'run').
In
`r if (knitr::is_latex_output()) {
   'Fig. \\@ref(fig:ExampleScatterLATEX) (left panel),'
} else {
   'Fig. \\@ref(fig:ExampleScatterHTML),'
}`
I have drawn a sensible line on the graph to capture the relationship (your line may look a little different).
When $x = 0$, the regression line predicts the value of $y$ to be about $2$, so $b_0$ is approximately $2$.

To approximate the slope, use the 'rise over run'
`r if (knitr::is_latex_output()) {
   'idea (Fig. \\@ref(fig:ExampleScatterLATEX), right panel; the online version has an animation).'
} else {
   'idea; the animation below may help explain.'
}`
When the value of $x$ increases from $1$ to $5$ (a change of $5 - 1 = 4$), the corresponding value of $y$ changes from 5 to 17 (a change of $17 - 5 = 12$).
So:  
\[
   \frac{\text{rise}}{\text{run}} 
   = \frac{17 - 5}{5 - 1}
   = \frac{12}{4} 
   = 3.
\]
The value of $b_1$ is about $3$.
The regression line is approximately $\hat{y} = 2 + (3\times x)$, usually written  
\[
  \hat{y} = 2 + 3x.
\]

::: {.importantBox .important data-latex="{iconmonstr-warning-8-240.png}"}
The regression equation has $\hat{y}$ (not $y$) on the left-hand side.
That is, the equation *predicts* the mean values of $y$, which may not be equal to the observed values (which are denoted by $y$).

A 'good' regression equation would produce predicted values $\hat{y}$ close to the observed values $y$; that is, the lines passes close to each point on the scatterplot.
:::


The *intercept* has the same measurement units as the response variable.
The measurement unit for the *slope* is the 'measurement units of the response variable', per 'measurement units of the explanatory variable'.



:::{.example name="Measurement units for intercept and slope"}
For the red-deer data, explanatory variable is age (in years), and the response variable is molar weight (in g).
Hence, the intercept $b_0$ is measured in grams, and the slope $b_1$ is measured in grams per year.
:::

  
<!-- FOR LaTeX: Put this example data plot and the static rise/run in the same plot -->
<!-- FOT HTML, separate plots -->

```{r ExampleScatterHTML, fig.cap="An example scatterplot", fig.align="center", fig.width=5, fig.height=4}
if (knitr::is_html_output()){
  set.seed(5000)
x <- seq(0.5, 5, 
        length = 10)
mu <- 2 + 3 * x
y <- mu + rnorm(length(x), 0, 0.5) 
m1 <- lm(y ~ x)

plot(y ~ x, 
     xlim = c(0, 5), 
     ylim = c(0, 20), 
     las = 1, 
     xlab = expression(paste("Explanatory, ", italic(x)) ), 
     ylab = expression(paste("Response, ", italic(y)) ), 
     pch = 19)
grid()

abline( coef(m1), 
        col = "grey",
        lwd = 2)
}
```

```{r AnimateRiseRun, animation.hook="gifski", interval=0.85, dev=if (is_latex_output()){"pdf"}else{"png"}}
if (knitr::is_html_output()){
  NumPlots <- 9

  set.seed(5000)

  x <- seq(0.5, 5,  
           length = 10)
  b0 <- 2
  b1 <- 3
    
  mu <- b0 + b1 * x
  y <- mu + rnorm(length(x), 0, 0.5)

  x1 <- 1
  x2 <- 5
  y1 <- b0 + (b1 * x1) 
  y2 <- b0 + (b1 * x2)
    

  for (i in (1:NumPlots)){
    par( mar = c(5, 4, 4, 2) + 0.1)
    # DATA
    
    plot(y ~ x, 
         xlim = c(0, 5), 
         ylim = c(0, 20), 
         type = "n",
         las = 1, 
         xlab = expression(paste("Explanatory, ", italic(x)) ), 
         ylab = expression(paste("Response, ", italic(y)) ), 
         pch = 19)
    rug(1:20, 
        side = 2, 
        ticksize = -0.02)
    grid()
    points(y ~ x, 
           pch = 19) 
    
    if ( i >=2 ){ # DATA + LINE
      abline(b0, b1,
             col = "red", 
             lwd = 2)
    }
    if (i >= 3) { # DATA + LINE + x1 
      lines( c(x1, x1), 
             c(0, y1),
             lwd = 2, 
             col = "darkgrey",
             lty = 2)
    }
    
    if ( i >= 4){  # DATA + LINE + x1 + y1
      lines( c(0, x1), 
             c(y1, y1),
             lwd = 2, 
             col = "darkgrey",
             lty = 2)
    }
    
    if ( i >= 5 ){  # DATA + LINE + x1 + y1 + x2
      lines( c(x2, x2), 
             c(0, y2),
             lwd = 2, 
             col = "darkgrey",
             lty = 2)
    }
    
    if ( i >= 6) {  # DATA + LINE + x1 + y1 + x2 + y2
      lines( c(0, x2), 
             c(y2, y2),
             lwd = 2, 
             col = "darkgrey",
             lty = 2)
    }
    
    if (i >= 7 ) {  # DATA + LINE + x1 + y1 + x2 + y2 + RUN
      run <- x2 - x1
      arrows(x1, 2, 
             x2, 2,
             lwd = 2,
             length = 0.15,
             col = "darkgrey",
             code = 2)
      text( mean( c(x1, x2)), 2.5, 
            labels = bquote("Run: "~ .(x2)-.(x1) == .(run)),
            pos = 2 )  
    }
    
    if (i >= 8 ){ # DATA + LINE + x1 + y1 + x2 + y2 + RUN + RISE
      rise <- y2 - y1
      arrows(0.5, y1, 
             0.5, y2,
             lwd = 2,
             length = 0.15,
             col = "darkgrey",
             code = 2)
      text( 0.5, 
            mean( c(y1, y2)), 
            labels = bquote("Rise: "~ .(y2)-.(y1) == .(rise)),
            pos = 4 )  
    }
    
    
    if (i >= 9 ) {  # DATA + LINE + x1 + y1 + x2 + y2 + RUN + RISE + SLOPE
      slope <- rise / run
      polygon( c(2.15, 2.15, 4.85, 4.85), 
               c(6, 8, 8, 6), 
               border = NA, # No border
               col = "antiquewhite1")
      text(3.5, 
           7, 
            labels = bquote("Slope: "~ .(rise) %/%  .(run) == .(slope)) )
    }
  }
}
```


<!-- FOR LaTeX: Put this example data plot and the static rise/run in the same plot -->
```{r ExampleScatterLATEX, fig.cap="An example scatterplot. Left: Approximating $b_0$; right: approximating $b_1$ using rise-over-run", fig.align="center", fig.width=9, fig.height=4, out.width='100%'}
set.seed(5000)
x <- seq(0.5, 5, 
         length = 10)
mu <- 2 + 3 * x
y <- mu + rnorm(length(x), 0, 0.5)
m1 <- lm(y ~ x)

if (knitr::is_latex_output()){
  
  par(mfrow = c(1, 2))
  
  plot(y ~ x, 
       xlim = c(0, 5), 
       ylim = c(0, 20), 
       las = 1, 
       main = expression(Approximating~italic(b)[0]),
       xlab = expression(Explanatory~italic(x) ), 
       ylab = expression(Response~italic(y) ), 
       pch = 19,
       col = "grey")
  grid()
   rug(1:20, 
      side = 2, 
      ticksize = -0.02)

  
  
  abline( coef(m1), 
          col = "black",
          lwd = 2)
  
  # Beta 0
  abline(v = 0,
         lwd = 2,
         col = "grey")
  lines( x = c(0, 0),
         y = c(0, coef(m1)[1]),
         lty = 1,
         lwd = 2,
         col = "grey")
  lines( x = c(-1, 0.45),
         y = c(coef(m1)[1], coef(m1)[1]),
         lty = 1,
         lwd = 2,
         col = "grey")
  
  arrows(x0 = 2,
         x1 = 0.5,
         y0 = 2,
         y1 = coef(m1)[1],
         angle = 15,
         length = 0.25,
         lwd = 2)
  text(x = 2,
       y = 2,
       pos = 4,
       labels = expression(Approximately~italic(b)[0]))
  
  ######
  
  b0 <- 2
  b1 <- 3
  
  x1 <- 1
  x2 <- 5
  y1 <- b0 + (b1 * x1) 
  y2 <- b0 + (b1 * x2)
  
  NumPlots <- i <- 9
  #  for (i in (NumPlots:NumPlots)){
  
  # DATA
  plot(y ~ x, 
       xlim = c(0, 5), 
       ylim = c(0, 20), 
       type = "n",
       las = 1, 
       main = expression(Approximating~italic(b)[1]),
       xlab = expression(Explanatory~italic(x) ), 
       ylab = expression(Response~italic(y) ), 
       pch = 19)
  rug(1:20, 
      side = 2, 
      ticksize = -0.02)
  grid()
  
  
  points(y ~ x, 
         pch = 19) 
  
  if ( i >=2 ){ # DATA + LINE
    abline(b0, b1,
           col = "red", 
           lwd = 2)
  }
  if (i >= 3) { # DATA + LINE + x1 
    lines( c(x1, x1), 
           c(0, y1),
           lwd = 2, 
           col = "darkgrey",
           lty = 2)
  }
  
  if ( i >= 4){  # DATA + LINE + x1 + y1
    lines( c(0, x1), 
           c(y1, y1),
           lwd = 2, 
           col = "darkgrey",
           lty = 2)
  }
  
  if ( i >= 5 ){  # DATA + LINE + x1 + y1 + x2
    lines( c(x2, x2), 
           c(0, y2),
           lwd = 2, 
           col = "darkgrey",
           lty = 2)
  }
  
  if ( i >= 6) {  # DATA + LINE + x1 + y1 + x2 + y2
    lines( c(0, x2), 
           c(y2, y2),
           lwd = 2, 
           col = "darkgrey",
           lty = 2)
  }
  
  if (i >= 7 ) {  # DATA + LINE + x1 + y1 + x2 + y2 + RUN
    run <- x2 - x1
    arrows(x1, 2, 
           x2, 2,
           lwd = 2,
           length = 0.15,
           col = "darkgrey",
           code = 2)
    text( mean( c(x1, x2)), 
          2, 
          labels = bquote("Run: "~ .(x2)-.(x1) == .(run)),
          pos = 3 )  
  }
  
  if (i >= 8 ){ # DATA + LINE + x1 + y1 + x2 + y2 + RUN + RISE
    arrows(0.5, y1, 
           0.5, y2,
           lwd = 2,
           length = 0.15,
           col = "darkgrey",
           code = 2)
    rise <- y2 - y1
    text( 0.5, 
          mean( c(y1, y2)) + 0.5, 
          labels = bquote("Rise: "~ .(y2)-.(y1) == .(rise)),
          pos = 4 )  
  }
  #  }
}
```


`r if (knitr::is_html_output()){
  'You may like the following interactive activity, which explores slopes and intercepts.'}`

<iframe src="https://phet.colorado.edu/sims/html/graphing-slope-intercept/latest/graphing-slope-intercept_en.html" width="600" height="450" scrolling="no" allowfullscreen="allowfullscreen"></iframe>


::: {.example #CycloneRegressionGuesses name="Estimating regression parameters"}
The relationship [@mypapers:dunnsmyth:glms] between the number of cyclones $y$ in the Australian region each year from 1969 to 2005, and a climatological index called the *Ocean `r readr::parse_character( c("Niño"), locale=locale(encoding="UTF-8"))` Index* (ONI, $x$), is shown in Fig. \@ref(fig:ONIcyclones).

When the value of $x$ is zero, the predicted value of $y$ is about $12$; $b_0$ is about $12$.
(You may get something slightly different.)
Notice that the intercept is the *predicted* value of $y$ when $x = 0$, which is *not* at the left of the graph.

To approximate the value of $b_1$, use the 'rise over run' idea.
When $x$ is about $-2$, the predicted value of $y$ is about $17$; when $x$ is about $2$, the predicted value of $y$ is about $8$.
The value of $x$ increases by $2 - (-2) = 4$,  while the value of $y$ changes by $8 - 17 = -9$ (a *decrease* of about $9$).
Hence, $b_1$ is approximately $-9/4 = -2.25$.
(You may get something slightly different.)
The relationship has a *negative* direction, so the slope must be *negative*.
The regression line is approximately $\hat{y} = 12 - 2.25x$.
:::

```{r ONIcyclones, fig.width=4, fig.align="center", fig.cap="The number of cyclones in the Australian region each year from 1969 to 2005, and the ONI for October, November, December", fig.width=5, fig.height=4}
data(Cyclones)

plot(Total ~ OND, 
     data = Cyclones,
     xlim = c(-2, 3),
     ylim = c(0, 20),
     pch = 19,
     ylab = "Total number of cyclones",
     xlab = "ONI averaged over Oct., Nov., Dec.",
     las = 1)

clm <- lm( Total ~ OND, 
           data = Cyclones)
abline( coef(clm), 
        lwd = 2, 
        col = "grey")
grid(lty = 2)
```


The above method gives a crude approximation to the values of the intercept $b_0$ and the slope $b_1$.
In practice, *many* reasonable lines could be drawn through a scatterplot of data.
However, one of those lines is the 'best' line in some sense^[For those interested: The 'line of best fit' is the line for which the sum of the *squared* vertical distances between the observations $y$ and the predicted values $\hat{y}$ (i.e, the regression line) is as small as possible.].
Software calculates this 'line of best fit' for us.


<iframe src="https://learningapps.org/watch?v=p84zyy2ot22" style="border:0px;width:100%;height:800px" allowfullscreen="true" webkitallowfullscreen="true" mozallowfullscreen="true"></iframe>



## Regression using software {#Regression-Software}

The regression line, explained above, is computed from the observed sample data.
This assumes that a regression equation exists in the population that is being estimated using the sample information.
In the population, the intercept is denoted by $\beta_0$ and the slope by $\beta_1$.
These population values are unknown, and are estimated by the statistics $b_0$ and $b_1$ respectively.


:::: {.pronounceBox .pronounce data-latex="{iconmonstr-microphone-7-240.png}"}

::: {style="display: flex;"}
The symbol $\beta$ is the Greek letter 'beta', pronounced 'beater' (as in 'egg beater').  
So $\beta_0$ is said as 'beater-zero', and $\beta_1$ as 'beater-one'.
:::

::: {}
```{r}
htmltools::tags$video(src = "./Movies/beta.mp4", 
                      width = "121", 
                      loop = "FALSE", 
                      controls = "controls", 
                      loop = "loop", 
                      style = "padding:5px; border: 2px solid gray;")
```
:::

::::


Every sample is likely to produce a slightly different value for both $b_0$ and $b_1$ (*sampling variation*), so both $b_0$ and $b_1$ have a standard error.
The formulas for computing the values of $b_0$ and $b_1$ are intimidating, so we will allow software to do the calculations. 

```{r}
RD.lm <- lm(Weight ~ Age, 
            data = RedDeer)
RD.b0 <- round( coef(RD.lm)[1], 3)
RD.b1 <- round( coef(RD.lm)[2], 3)
```


<div style="float:right; width: 222x; border: 1px; padding:10px">
<img src="Illustrations/diana-parkhouse-Qf4eZXC3t2Y-unsplash.jpg" width="200px"/>
</div>



For the red-deer data again
`r if( knitr::is_latex_output() ) {
  '(Fig. \\@ref(fig:RedDeerScatter)),'
} else {
  '(Fig. \\@ref(fig:RedDeerScatterHTML)),'
}
`
part of the relevant output is shown in Fig.&nbsp;\@ref(fig:RedDeerRegressionjamoviSPSSraw).
From the output, the sample *slope* is $b_1 = `r RD.b1`$, and the sample $y$-intercept is $b_0 = `r RD.b0`$: the values of $b_0$ and $b_1$ are in the column labelled  `Estimate` in jamovi, and the column labelled `B` in SPSS.
These are the values of the two *regression coefficients*; so  
\[
	\hat{y} = `r RD.b0` + (`r RD.b1`\times x),
	\qquad\text{usually written as}\qquad  
	\hat{y} = `r RD.b0` - `r abs(RD.b1)` x.
\]         


::: {.tipBox .tip data-latex="{iconmonstr-info-6-240.png}"}
The *sign* of the slope $b_1$ and the sign of correlation coefficient $r$ are always the same.
For example, if the slope is negative, the correlation coefficient will also be negative.
:::


```{r RedDeerRegressionjamoviSPSSraw, fig.cap="Output for the red-deer data,Left: jamovi; right: SPSS", fig.align="center", out.width=c("45%", "54%"), fig.show="hold"}
knitr::include_graphics( "jamovi/RedDeer/RedDeer-Regression.png")
knitr::include_graphics( "SPSS/RedDeer/RedDeer-Regression.png")
```



```{r}
cyclone.lm <- lm(Total ~ OND, 
                 data = Cyclones)
```

::: {.example #RegressionCoefficients name="Regression coefficients"}
The regression equation for the cyclone data (Fig. \@ref(fig:ONIcyclones)) is found from the jamovi output (Fig. \@ref(fig:CyclonesRegressionjamoviSPSS)):  
\[
  \hat{y} = 12.14 - 2.23x,
\]
where $x$ is the ONI (averaged over October, November, December) and $y$ is the number of cyclones.
These values are close the approximations made in Example \@ref(exm:CycloneRegressionGuesses) ($12$ and $-2.25$ respectively).
:::


```{r CyclonesRegressionjamoviSPSS, fig.cap="Output for the cyclone data. Left: jamovi; right: SPSS", fig.align="center", out.width=c("43%", "55%"), fig.show="hold"}
knitr::include_graphics( "jamovi/Cyclones/Cyclones-Regression.png")
knitr::include_graphics( "SPSS/Cyclones/SPSS-Cycones-Regression.png")
```


`r if (knitr::is_html_output()){
  'You may like the following interactive activity, which explores regression equations.'}`

<div style="text-align:center;">
<iframe src="https://phet.colorado.edu/sims/html/least-squares-regression/latest/least-squares-regression_en.html" width="600" height="450" scrolling="no" allowfullscreen="true"></iframe>
</div>



## Regression for predictions {#RegressionForPrediction}


<div style="float:right; width: 222x; border: 1px; padding:10px">
<img src="Illustrations/diana-parkhouse-Qf4eZXC3t2Y-unsplash.jpg" width="200px"/>
</div>


The regression equation for the red-deer data, estimated from one of the many possible samples, is $\hat{y} = `r RD.b0` - `r abs(RD.b1)` x$, and can be used to make *predictions* of the mean value of $y$ for a given value of $x$; for example, predicting the *average* molar weight for deer $10$ years old.
Since $x$ represents the age, use $x = 10$ in the regression equation:  
\begin{eqnarray*}
         \hat{y} 
         &=& `r RD.b0` - (`r abs(RD.b1)`\times 10)\\
         &=& `r RD.b0` - `r abs(RD.b1 * 10)` = `r RD.b0 + RD.b1*10`.
\end{eqnarray*}
Male red deer aged $10$ years are predicted to have a *mean* molar weight of $`r  RD.b0 + RD.b1*10`$ grams.
Some individual deer aged $10$ will have molars weighing *more* than this, and some weighing *less* than this; the model predicts that the *mean* molar weight for male red deer aged 10 will be about $`r RD.b0 + RD.b1*10`$ grams.

(Note: The value of $\hat{y}$ is computed using the estimates $b_0$ and $b_1$, which are in turn computed from sample data. 
Hence, the value of $\hat{y}$ is also based on the sample informatiom, so also has a standard error.)


::: {.thinkBox .think data-latex="{iconmonstr-light-bulb-2-240.png}"}
What is the predicted mean molar weight for for male red deer $12$ years of age?\label{thinkBox:Predictions}

`r if (knitr::is_latex_output()) '<!--'`
`r webexercises::hide()`
Prediction: $4.398 - (0.181 \times 12) = 2.226$, or about 2.23 grams.
`r webexercises::unhide()`
`r if (knitr::is_latex_output()) '-->'`
:::

Suppose we were interested in male red deer $20$ years of age; the mean predicted weight of the molars would be $4.398 - (0.181 \times 20) = 0.778$, or about $0.78$ grams.
However, while this prediction *may* be a useful prediction... it may be rubbish.

In the data, the oldest deer is aged 14.4 years, so the regression line may not even apply for deer aged over 14.4 years of age.
For example, the relationship may be non-linear after 14 years of age, or red deer may not even live to 20 years of age.
The prediction *may* be sensible... but it *may not* be either.

We don't know whether the prediction is sensible or not, because we have no data for deer aged over 14.4 years to inform us.
Making predictions outside the range of the available data is called *extrapolation*, and *extrapolation* beyond the data may lead to nonsense predictions.


::: {.definition #Extrapolation name="Extrapolation"}
*Extrapolation* refers to making predictions outside the range of the available data.
Extrapolation beyond the data may lead to nonsense predictions.
:::


## Regression for understanding {#RegressionForUnderstanding}

The regression equation can be used to *understand* the relationship between the two variables.
Consider again the red-deer regression equation:  
\begin{equation}
   \hat{y}
   = 
   `r RD.b0` 
   - 
    `r abs(RD.b1)` x.
   (\#eq:RedDeerEquation)
\end{equation}
What does this equation reveal about the relationship between $x$ and $y$?


### The meaning of $b_0$

$b_0$ is the *predicted* value of $y$ when $x = 0$.
Using $x = 0$ in Eq. \@ref(eq:RedDeerEquation) predicts a mean molar weight of  
\[
  \hat{y} = 4.398 - (0.181\times 0) = 4.398
\]
for deer zero years of age (i.e., newborn male red deer).
This prediction is the predicted *mean* molar weight; some individual deer will have molars weights greater than this, and some less than this.
But in any case, this predicted mean may be nonsense: it is *extrapolating* beyond the data (the youngest deer in the sample is aged $4.4$ years).


::: {.tipBox .tip data-latex="{iconmonstr-info-6-240.png}"}
The value of the intercept $b_0$ is sometimes (but not always) meaningless.
The value of the slope $b_1$ is usually of greater interest, as it explains the *relationship* between the two variables.
:::


### The meaning of $b_1$ {#SlopeMeaning}

$b_1$ tells us how much the value of $y$ changes (on average) when the value of $x$ *increases* by one.
For the red-deer data, $b_1$ tells us how much the molar weight changes (on average) when age increases by one year.

Each extra year older is associated with an average change of $`r RD.b1`$ grams in molar weight; that is, a *decrease* in molar weight by a mean of $`r abs(RD.b1)`$g.
The molars of some individual deer will lose more weight than this in some years, and some will lose less... the value is a *mean* weight loss per year.

To demonstrate: When $x = 10$, $y$ is predicted to be $\hat{y}= `r RD.b0 + 10 * RD.b1`$g.
For deer one year older than this (i.e.\ $x = 11$) we predict $y$ to increase by $b_1 = `r RD.b1`$g higher (or, equivalently, $0.181$ *lower*).
That is, we would predict $\hat{y}= `r RD.b0 + 10*RD.b1`  - `r abs(RD.b1)` = `r RD.b0 + 11*RD.b1`$g.
This is the same prediction made by using $x = 11$ in Eq. \@ref(eq:RedDeerEquation).


::: {.tipBox .tip data-latex="{iconmonstr-info-6-240.png}"}
If the value of $b_1$ is *positive*, then the predicted values of $y$ *increase* as the values of $x$ *increase*.
If the value of $b_1$ is *negative*, then the predicted values of $y$ *decrease* as the values of $x$ *increase*.
:::


This interpretation of $b_1$ explains the relationship: Each extra year of age reduces the weight of the molars by 0.181 grams, on average, in male red deer.
Recall that the units of measurement of the slope here are 'grams per year'.

Observe what happens if the slope is *zero*.
Since $b_1$ is the change in $y$ (on average) when $x$ increase by one, $b_1 = 0$ means that the predicted mean value of $y$ changes by *zero* if the value of $x$ changes by one.
In other words, if the value of $x$ changes, the predicted mean value of $y$ doesn't change.
This is equivalent to saying that there is *no relationship* between the variables.
(We would also find $r = 0$.)


::: {.importantBox .important data-latex="{iconmonstr-warning-8-240.png}"}
If the value of the slope is zero, there is *no linear relationship* between $x$ and $y$.
In this case, the correlation coefficient is also zero.
:::


## Confidence intervals for the regression parameters {#RegressionCI}

### Introduction {#RegressionTests-Intro}

The regression line is computed from the *sample* data, assuming a linear relationship actually exists in the *population*.
The values of $b_0$ and $b_1$ are estimates of the population parameters $\beta_0$ and $\beta_1$.

The *sample* estimates can be used to ask questions about the unknown *population* regression coefficients.
As usual, the sample estimates can vary across all possible samples (and so have a *sampling distribution*).

Usually the slope of more interest than the intercept, because the slope explains the  *relationship* between the two variables (Sect. \@ref(RegressionForUnderstanding)); however the same ideas apply for the intercept.
The CIs can be obtained from software or computed manually.


::: {.thinkBox .think data-latex="{iconmonstr-light-bulb-2-240.png}"}
Using the output (Fig. \@ref(fig:RedDeerRegressionjamoviSPSSraw)) for the red-deer data, can you determine the *approximate* 95% CI for $\beta_1$?\label{thinkBoxCI}
:::

### Describing the sampling distribution {#SamplingDistributionSlopeCI}

The observed value of the slope will vary across all possible samples; *sampling variation* is present.
However, under certain conditions (Sect. \@ref(RegressionStatValidity)), the sampling distribution can be described as:

* an approximate normal distribution,
* with a mean of $\beta_1 = 0$, and 
* a standard deviation, called the *standard error of the slope* and denoted $\text{s.e.}(b_1)$.

(A formula exists for finding $\text{s.e.}(b_1)$, but it is tedious to use and we will not give it.)
This distribution describes the possible values of the sample slope, through *sampling variation*.
For the red-deer data then, the values of the sample slope across all possible samples is described (Fig. \@ref(fig:RedDeerSlopeSampDistCI)) as:

* an approximate normal distribution,
* with a sampling mean whose values is $\beta_1 = 0$, and 
* with a standard deviation, called the *standard error of the slope* $\text{s.e.}(b_1) = 0.029$ (from software; Fig. \@ref(fig:RedDeerRegressionjamoviSPSSraw)).


```{r RedDeerSlopeSampDistCI, fig.cap="The distribution of sample slopes for the red-deer data, around the true slope $\\beta_1$", fig.align="center", fig.height=3.5, fig.width=8.25, out.width='85%'}
mn <- 0
se <- 0.0289

out <- plotNormal(mu = mn, 
                  sd = se,
                  round.dec = 4,
                  cex.axis = 0.95,
                  ylim = c(0, 27),
                  showXlabels = c( expression(beta[1] - 0.087),
                                   expression(beta[1] - 0.058),
                                   expression(beta[1] - 0.029),
                                   expression(beta[1]),
                                   expression(beta[1] + 0.029),
                                   expression(beta[1] + 0.058),
                                   expression(beta[1] + 0.087)),
                  main = "", # Main title
                  xlab = "Sample regression slopes", # horizontal axis labels
                  showZ = TRUE) # Whether to show z = -3:3 or not

arrows( x0 = 0,
        x1 = 0,
        y0 = 22,
        y1 = 14,
        angle = 15,
        length = 0.15)
text(x = 0,
     y = 22,
     pos = 3,
     labels = "Sampling mean")

arrows(x0 = 0,
       x1 = 0.029,
       y0 = 3,
       y1 = 3,
       angle = 15,
       length = 0.15)
text(x = 0.029/2,
     y = 3,
     pos = 3,
     labels = "Std. error")
```


Most CIs have the form  
\[
  \text{statistic} \pm ( \text{multiplier} \times \text{standard error}),
\]
where the multiplier is two for an approximate 95% CI (from the 68--95--99.7 rule),.
Using the standard error reported by the software, an approximate 95% CI is $-0.181 \pm (2\times 0.029)$, or $-0.181 \pm 0.058$, or from $-0.239$ to $-0.123$.

Software can be used to produce *exact* CIs too (Fig. \@ref(fig:RedDeerRegressionjamoviSPSSCI)).
The *approximate* and *exact* 95% CIs are very similar.
We could write:

> The approximate 95% confidence, for the slope in the regression line explaining molar weight of male red deer from age, is between $-0.239$ to $-0.123$ grams per year.


```{r RedDeerRegressionjamoviSPSSCI, fig.cap="Output for the red-deer data, including the CIs for the regression parameters. Top: jamovi; bottom: SPSS", fig.align="center", out.width=c("75%", "90%"), fig.show="hold"}
knitr::include_graphics( "jamovi/RedDeer/RedDeer-Regression-CI.png")
knitr::include_graphics( "SPSS/RedDeer/RedDeer-Regression-CI.png")
```


::: {.example #EDpatientsCI name="Emergency department patients"}
A study examined the relationship between the number of emergency department (ED) patients and the number of days following the distribution of monthly welfare monies [@brunette1991correlation] from 1986 to 1988 in Minneapolis, USA.
The data (extracted from Fig. 2 of @brunette1991correlation) are displayed in Fig. \@ref(fig:EDScatterjamovi) (left panel), and the jamovi output in Fig. \@ref(fig:EDScatterjamovi) (right panel)).
From the jamovi output (Fig. \@ref(fig:EDScatterjamovi)), the regression line is estimated as 
\[
  \hat{y} = 150.19 - 0.348x, 
\]
where $y$ represents the mean number of ED patients, and $x$ the number of days since welfare distribution.

A 95% CI for the slope is not shown in the output (though it could have been requested).
However, since most CIs have the form
\[
  \text{statistic} \pm ( \text{multiplier} \times \text{standard error}),
\]
an approximate 95% CI is easily computed:
\[
  -0.34790 \pm (2\times 0.046672), 
\]
or $-0.34790 \pm 0.093344$, equivalent to $-0.441$ to $-0.255$ patients per day.
:::


```{r EDScatterjamovi, fig.align = "center", fig.cap = "The number of emergency department patients, and the number of days since distribution of welfare. Left: scatterplot; right: jamovi output.", fig.width=5, fig.height=3.5, out.width=c("45%", "54%"), fig.show="hold"}
data(EDpatients)

plot( ED ~ Days,
      data = EDpatients,
      xlab = "Days after distribution",
      ylab = "Mean number of ED patients",
      main = "Mean number of ED patients vs \ndays after distribution of welfare",
      pch = 19,
      las = 1,
      xlim = c(0, 30),
      ylim = c(135, 155) )
knitr::include_graphics( "jamovi/ED/EDRegression-crop.png")
```



## Hypothesis testing for the regression parameters {#RegressionHT}

### Statistical hypotheses {#RegressionHypotheses}

The null hypothesis is the usual 'no relationship' hypothesis.
In this context, 'no relationship' means that the slope is zero (Sect. \@ref(SlopeMeaning)), so the null hypotheses (about the *population*) is $H_0$: $\beta_1 = 0$.

This hypothesis proposes that $\beta_1$ is zero, but $b_1$ is not zero due to sampling variation.
As part of the [decision-making process](#DecisionMaking), the null hypothesis is initially *assumed* to be true.
For the red-deer data (Sect. \@ref(RedDeerData)), testing if a relationship exists between the age of the deer and the weight of their molars implies these hypotheses:

\[
   \text{$H_0$: } \beta_1 = 0;\quad\text{and}\quad\text{$H_1$: } \beta_1 \ne 0.
\]
The parameter is $\beta_1$, the population slope for the regression equation predicting molar weight from age.
The alternative hypothesis is two-tailed, based on the RQ.


```{r}
RD.fitinfo <- summary(RD.lm)$coefficients
RD.t <- RD.fitinfo[, "t value"]
RD.P <- RD.fitinfo[, "Pr(>|t|)"]
```



### Describing the sampling distribution {#SamplingDistributionSlopeHT}

*Assuming* the null hypothesis is true (that $\beta_1 = 0$), the possible values of the sample slope $b_1$ can be described, due to *sampling variation*.
The variation in the sample slope across all possible samples is described (Fig. \@ref(fig:RedDeerSlopeSampDist)) using:

* an approximate normal distribution,
* with a sampling mean whose value is $\beta_1 = 0$ (from $H_0$), and 
* a standard deviation, called the *standard error of the slope* and denoted $\text{s.e.}(b_1)$, with a value of $0.0289$ (from software; Fig. \@ref(fig:RedDeerRegressionjamoviSPSSCI)).


```{r RedDeerSlopeSampDist, fig.cap="The distribution of sample slopes for the red-deer data, if the population slope is 0", fig.align="center", fig.height=3.5, fig.width=8.25, out.width='85%'}
par( mar = c(5, 2, 1, 2))
mn <- 0
se <- 0.0289

out <- plotNormal(mu = mn, 
                  sd = se,
                  ylim = c(0, 27),
                  round.dec = 4,
                  main = "", # Main title
                  xlab = "Sample regression slopes", # horizontal axis labels
                  showZ = TRUE) # Whether to show z = -3:3 or not

arrows( x0 = 0,
        x1 = 0,
        y0 = 22,
        y1 = 14,
        angle = 15,
        length = 0.15)
text(x = 0,
     y = 22,
     pos = 3,
     labels = "Sampling mean")

arrows(x0 = 0,
       x1 = 0.029,
       y0 = 3,
       y1 = 3,
       angle = 15,
       length = 0.15)
text(x = 0.029/2,
     y = 3,
     pos = 3,
     labels = "Std. error")
```


### Computing the test statistic

The *observed* sample slope is $b_1 = -0.181$.
Locating this on Fig. \@ref(fig:RedDeerSlopeSampDist) (way to the left) shows that it is very unlikely that one of the many possible samples would produce such a slope, just through sampling variation (if the population slope really was $\beta_1 = 0$).
The *test statistic* is found using the usual approach when the sampling distribution has an approximate normal distribution:  
\begin{align*}
   t 
   &= \frac{\text{observed sample slope} - \text{slope sampling mean}}{\text{standard error of the slope}}\\
   &= \frac{ b_1 - \beta_1}{\text{s.e.}(b_1)} \\
   &= \frac{-0.181 - 0}{0.0289} = -6.27,
\end{align*}
where the values of $b_1$ and $\text{s.e.}(b_1)$ are taken from the software output (Fig. \@ref(fig:RedDeerRegressionjamoviSPSSCI)).
This $t$-score is also reported by the software.


### Determining the $P$-value

To determine if the statistic is *consistent* with the null hypothesis, the $P$-value can be approximated using the [68--95--99.7 rule](#def:EmpiricalRule), approximated using tables, or taken from software output (Fig. \@ref(fig:RedDeerRegressionjamoviSPSSCI)).
Using software, the two-tailed $P$-value is $P < 0.001$.


### Writing conclusions

We write:

> The sample presents very strong evidence ($t = -6.27$; two-tailed $P < 0.001$) that the slope in the population between age of the deer and molar weight is not zero (slope: $-0.181$; 95% CI from $-0.239$ to $-0.124$).

Notice the three features of writing conclusions again:
An *answer to the RQ*; evidence to support the conclusion ('$t = -6.27$; two-tailed $P < 0.001$); and some *sample summary information* ('slope: $-0.181$; 95% CI from $-0.239$ to $-0.124$').


<!-- ```{r RedDeerRegressionjamoviSPSS, fig.cap="Output for the red-deer data. Top: jamovi; bottom: SPSS", fig.align="center", out.width=c("60%", "80%"), fig.show="hold"} -->
<!-- knitr::include_graphics( "jamovi/RedDeer/RedDeer-Regression.png") -->
<!-- knitr::include_graphics( "SPSS/RedDeer/RedDeer-Regression.png") -->
<!-- ``` -->


::: {.example #EDpatientsRegression name="Emergency department patients"}
For the emergency department data used in Example \@ref(exm:EDpatientsCI), the scatterplot of the data and the jamovi output are shown in Fig. \@ref(fig:EDScatterjamovi).
The regression line is estimated as  
\[
  \hat{y} = 150.19 - 0.348x, 
\]
where $y$ represents the mean number of ED patients, and $x$ the number of days since welfare distribution.
This regression equation suggests that each extra day after welfare distribution is associated with a *decrease* in the mean number of ED patients of about 0.35.
It may be easier to understand this way: 

> Each 10 extra days after welfare distribution is associated with a *decrease* in the number of ED patients of about $10\times 0.35 = 3.5$.

The scatterplot and the regression equation suggests a negative relationship between the number of ED patients and the days after distribution.
However, every sample is likely to be different, so the relationship may not actually be present in the population (the slope may be non-zero due to sampling variation).
We can test this hypothesis:  
\[
   \text{$H_0$: } \beta_1 = 0\quad\text{and}\quad\text{$H_1$: } \beta_1 \ne 0,
\]
where $\beta_1$ is the *population* slope.
The test is two-tailed, based on the authors' aim.
From the output, the test statistic is $t = -7.45$, which is very large; unsurprisingly, the two-tailed $P$-value is very small: $P < 0.001$.
We write:

> Very strong evidence exists ($t = -7.45$; two-tailed $P < 0.001$) of a relationship between the mean number of ED patients and the number of days after welfare distribution (slope: $-0.348$).
:::


<iframe src="https://learningapps.org/watch?v=phxny7fec22" style="border:0px;width:100%;height:500px" allowfullscreen="true" webkitallowfullscreen="true" mozallowfullscreen="true"></iframe>


## Statistical validity conditions {#RegressionStatValidity}

The results for the CI and the hypothesis test hold under [certain conditions to be met](#exm:StatisticalValidityAnalogy).
The conditions for which the test is statistically valid are the same as for correlation (Sect. \@ref(ValidityCorrelation)):

1. The relationship is approximately linear.
1. The variation in the response variable is approximately constant for all values of the explanatory variable.
1. The sample size is at least 25.

The sample size of 25 is a rough figure here, and some books give other values.


::: {.example #StatisticalValidityDeerReg name="Statistical validity"}
For the red-deer data, the relationship is approximately linear, and the variation in molar weight appears to be somewhat constant for various ages
`r if( knitr::is_latex_output() ) {
  '(Fig. \\@ref(fig:RedDeerScatter)),'
} else {
  '(Fig. \\@ref(fig:RedDeerScatterHTML)),'
}
`
so regression is appropriate.
The sample size is $n = 78$.
The results from the hypothesis test should be statistically valid.
:::


::: {.thinkBox .think data-latex="{iconmonstr-light-bulb-2-240.png}"}
Are the conditions for statistically validity met for the cyclones data (Fig. \@ref(fig:ONIcyclones))?\label{thinkBoxRegressionConditions}

`r if (knitr::is_latex_output()) '<!--'`
`r webexercises::hide()`
Looking at all three conditions:

* Linearity: The relationship looks approximately linear.
* Variation: The variations seems somewhat constant for all values of $x$ (the ONI), but it is a little hard to be sure.
* Sample size: The sample size is larger than 25.

The conditions for statistical validity appear to have been met.
`r webexercises::unhide()`
`r if (knitr::is_latex_output()) '-->'`
:::


::: {.example #EDpatientsValidity name="Emergency department patients"}
In Example \@ref(exm:EDpatientsRegression), the scatterplot (Fig. \@ref(fig:EDScatterjamovi)) shows the relationship is approximately linear, that the variation in ED patients seems reasonably constant for different numbers of days after distribution, and the sample size is larger than 25.
The CI (Example \@ref(exm:EDpatientsCI)) and the test (Example \@ref(exm:EDpatientsRegression)) should be statistically valid.
:::
                    



## Example: obstructive sleep apnoea {#Apnoea}

```{r}
OSAoriginal <- data(OSA)
OSA <- OSA[ -c(1, 7) ] # Remove ID (col1 ) and SAOS (col 6)
OSA$Gender <- factor(OSA$Gender, 
                     levels = 1:2, 
                     labels = c("Male", 
                                "Female"))
useOSA <- array( dim = c(11, 5))

useOSA[1:5, ] <- as.matrix(head(OSA, 5))
useOSA[7:11, ] <- as.matrix(tail(OSA, 5))
useOSA[6, ] <- rep("$\\vdots$", 5)
```

In a study of obstructive sleep apnoea (OSA) in adults with Down Syndrome [@carvalho2020stop], $n = 60$ adults underwent a sleep study and had various information recorded.
The main response variable of interest was OSA severity (measured using the Respiratory Event Index, REI): the average number of episodes of sleep disruption (according to specific criteria) per hour of sleep.
One research question is

> Among Down Syndrome adults, is there a relationship between the REI and neck size?

Here, $x$ is the neck size (in cm), and $y$ is the REI value.
`r if( knitr::is_html_output() ) {
  'The data are shown in Fig. \\@ref(fig:OSADT).'
} else
{
  'Part of the data are shown in Table \\@ref(tab:OSAKable).'
}
`


\begin{figure}
\centering
\begin{minipage}{0.45\textwidth}%
\captionof{table}{Part of the obstructive sleep apnoea dataset ($n = 60$)\label{tab:OSAKable}}
\fontsize{10}{12}\selectfont
```{r OSAKable}

knitr::kable( useOSA,
              format = "latex",
              longtable = FALSE,
              booktabs = TRUE,
              escape = FALSE,
              table.env = "@empty",
              linesep = "",
              align = "c",
              col.names = c("Age",
                            "Gender",
                            "BMI",
                            "Neck (cm)",
                            "REI")) %>%
  row_spec(0, bold = TRUE) %>%
  column_spec(1, width = "10mm") %>%
  column_spec(2, width = "13mm") %>%
  column_spec(3, width = "10mm") %>%
  column_spec(4, width = "10mm") %>%
  column_spec(4, width = "10mm")
#kable_styling(font_size = 10) # CANNOT USE THIS IN THE MINIPAGE
```
\end{minipage}
\hspace{0.05\textwidth}
\begin{minipage}{0.45\textwidth}%
\centering
```{r OSAscatter, fig.width=4, fig.height = 3.25, out.width='99%'}
plot( REI ~ Neck,
      data = OSA,
      xlab = "Neck circumference (in cm)",
      ylab = "REI",
      main = "REI vs neck circumference\nfor Down Syndrome adults",
      pch = 19,
      las = 1,
      ylim = c(0, 100),
      xlim = c(30, 50) )
```
\caption{Scatterplot of the neck circumference vs REI for Down Syndrome adults\label{fig:OSAscatter}}
\end{minipage}
\end{figure}



```{r OSADT, fig.cap="The obstructive sleep apnoea dataset", fig.align="center"}
if( knitr::is_html_output() ) {
  DT::datatable(OSA,
              list(scrollX = TRUE, 
                   scrollY = TRUE, 
                   ordering = FALSE),
              caption = "The Obstructive sleep apnoea dataset")
}
```

<!-- The figure for LaTeX is in the minipage (combined with data table), so only need show it for the HTML -->
`r if (knitr::is_latex_output()) '<!--'`
```{r OSAscatter2, fig.align = "center", fig.cap = "Scatterplot of the neck circumference vs REI for Down Syndrome adults", fig.width=5, fig.height=4}
plot( REI ~ Neck,
      data = OSA,
      xlab = "Neck circumference (in cm)",
      ylab = "REI",
      main = "REI vs neck circumference\nfor Down Syndrome adults",
      pch = 19,
      las = 1,
      ylim = c(0, 100),
      xlim = c(30, 50) )
```
`r if (knitr::is_latex_output()) '-->'`


Using the jamovi output (Fig. \@ref(fig:OSARegressionjamovi)), the correlation coefficientis $r = 0.264$ (and $R$-squared is about 7%), so the relationship is not very strong.
The value of the slope and $y$-intercept for the regression line in the *sample* are $b_0 = -24.47$ and $b_1 = 1.36$, so the regression equation is   
\[ 
  \hat{y} = -24.47 + 1.36x.
\]
The slope means that for each one centimetre increase in neck circumference, the number of sleep disruptions per hour increase (on average) by about $1.36$.


```{r OSARegressionjamovi, fig.cap="jamovi output for the Obstructive sleep apnoea data", fig.align="center", out.width="60%"}
knitr::include_graphics( "jamovi/OSA/OSARegression.png") 
```


Each sample will produce slightly different *sample* slopes, so we can test to see if the slope in the *population* is non-zero due to sampling variation, using a hypothesis test:  
\[
   \text{$H_0$: } \beta_1 = 0\quad\text{and}\quad\text{$H_1$: } \beta_1 \ne 0\quad \text{(that is, two-tailed)}.
\]
The parameter is $\beta_1$.
From the software output, $t = 2.09$ and the two-tailed $P$-value is $P = 0.041$, so moderate evidence exists that the neck circumference is associated with greater REI.

An *approximate* 95% CI for the population slope $\beta_1$ is  
\[
   1.3663 \pm (2\times 0.65441),
\]
or from $0.057$ to $2.68$.
From the scatterplot 
`r if (knitr::is_latex_output()) {
   '(Fig. \\@ref(fig:OSAscatter))'
} else {
   '(Fig. \\@ref(fig:OSAscatter2))'
}`
the results appear statistically valid.
We write:

> The sample presents moderate evidence ($t = 2.09$; two-tailed $P = 0.041$) of a relationship between neck circumference and REI (slope: $1.36$; $n = 60$; 95% CI from $0.057$ to $2.68$) in the population.


## Example: removal efficiency {#RegressionRemovalEfficiency}

In wastewater treatment facilities, air from biofiltration is passed through a membrane and dissolved in water, and is transformed into harmless byproducts.
The removal efficiency $y$ (in %) may depend on the inlet temperature (in $^\circ$C; $x$).
A RQ is

> In treating biofiltation wastewater, how does the removal efficiency depend on the inlet temperature?

The scatterplot of the data is shown in Fig. \@ref(fig:CorrelationRemovalEfficiency) and software output in Fig. \@ref(fig:RemovalEfficiencyjamoviSPSSRegression).
From the output, $b_0 = 97.5$ and $b_1 = 0.0757$, so the regression equation is  
\[
  \hat{y} = 97.5 + 0.0757x
\]
for $x$ and $y$ defined above.
The slope quantifies the relationship, so we can test  
\[
   \text{$H_0$: } \beta_1 = 0 \qquad\text{and}\qquad \text{$H_1$: } \beta_1 \ne 0.
\]
From the output, $t = 10.7$ whicis huge; the $P$#-value is small as expected: $P < 0.001$. The output does not include the CI, but since $\text{s.e.}(b_1) = 0.00705$, the approximate 95% CI is  
\[
  0.0757 \pm (2 \times 0.00705), \quad\text{ or }\quad 0.0757 \pm 0.0141.
\]
We can write:

> Very strong evidence exists ($t = 10.7$; $P < 0.001$) of a relationship between the inlet temperature and removal efficiency (slope: $0.0757$; approximate 95% CI: $0.0616$ to $0.0898$).


```{r RegressionRemovalEfficiency, fig.cap="The relationship between removal efficiency and inlet temperature", fig.align="center", fig.width=5, fig.height=3.5}
data(Removal)

plot(Removal ~ Temp, 
     data = Removal,
     xlab = "Temperature (degrees C)",
     ylab = "Removal efficiency (%)",
     pch = 19,
     las = 1
     )
```


```{r RemovalEfficiencyjamoviSPSSRegression, fig.cap="Regression output for the removal-efficiency data. Top: jamovi; bottom SPSS", fig.align="center", out.width=c("55%", "75%"), fig.show="hold"}
knitr::include_graphics("jamovi/Removal/RemovalRegression.png")
knitr::include_graphics("SPSS/Removal/RemovalRegression-crop.png")
```


<iframe src="https://learningapps.org/watch?v=pbefaxc9t22" style="border:0px;width:100%;height:500px" allowfullscreen="true" webkitallowfullscreen="true" mozallowfullscreen="true"></iframe>



## Summary {#Chap35-Summary}

In this chapter, we have learnt about **regression**, which mathematically describes the relationship between two *quantitative* variables.
The response variable is denoted by $y$, and the explanatory variable by $x$.
The linear relationship between them (the **regression equation**), in the sample, is

\[
   \hat{y} = b_0 + b_1 x,
\]
where $b_0$ is a number (the **intercept**), $b_1$ is a number (the **slope**), and the 'hat' above the $y$ indicates that the equation gives an *predicted* mean value of $y$ for the given $x$ value.

The *intercept* is the predicted mean value of $y$ when the value of $x$ is zero.
The *slope* is how much the predicted mean value of $y$ changes, on average, when the value of $x$ *increases* by 1.

The regression equation can be used to make *predictions* or to *understand* the relationship between the two variables.
Predictions made with values of $x$ outside the values of $x$ used to create the regression equation (called *extrapolation*) may not be reliable.

In the population, the regression equation is  
\[
   \hat{y} = \beta_0 + \beta_1 x.
\]
To test a hypothesis about a population slope $\beta_1$, based on the value of the sample slope $b_1$, **assume** the value of $\beta_1$ in the null hypothesis (usually zero) to be true.
Then, the sample slope varies from sample to sample and, under certain statistical validity conditions, varies with an approximate normal distribution centred around the hypothesised value of $\beta_1$, with a standard deviation of $\text{s.e.}(b_1)$.
This distribution describes what values of the sample slope could be **expected** in the sample if the value of $\beta_1$ in the null hypothesis was true.

The *test statistic* is  
\[
   t = \frac{ b_1 - \beta_1}{\text{s.e.}(b_1)},
\]
where $\beta_1$ is the hypothesised value given in the null hypothesis (usually zero).
The $t$-value is like a $z$-score, and so an approximate **$P$-value** can be estimated using the  [68--95--99.7 rule](#def:EmpiricalRule).


`r if (knitr::is_html_output()){
  'The following short video may help explain some of these concepts:'
}`

<div style="text-align:center;">
<iframe width="560" height="315" src="https://www.youtube.com/embed/4sl8xhLG6Ls" frameborder="0" allow="accelerometer; encrypted-media; gyroscope; picture-in-picture"></iframe>
</div>


## Quick review questions {#Chap35-QuickReview}


```{r, results='hide'}
data(AISsub)
AIS.out <- lm(Wt ~ Ht, 
              data = subset(AISsub, Sport == "Rowing") )
summary(AIS.out)
```


::: {.webex-check .webex-box}
A study of athletes [@data:Telford1991:sexsportsize] examined the relationship between the height and weight of $n = 37$ rowers at the Australian Institute of Sport (AIS), as shown in Fig. \@ref(fig:ScatterAISRowers).


```{r ScatterAISRowers, results ='hide', fig.cap="Scatterplot of Weight against Height rowers at the AIS", fig.width=5, fig.height=3.5, fig.align="center"}
plot( Wt ~ Ht,
      data = subset(AISsub, Sport == "Rowing"),
      las = 1,
      xlim = c(150, 200),
      ylim = c(50, 100),
      main = "Weight vs Height: AIS rowers",
      xlab = "Height (in cm)",
      ylab = "Weight (in kg)",
      pch = 19)
```


1. The $y$-variable is \tightlist
`r if( knitr::is_html_output() ) {mcq( c(
  answer="Weight",
  "Height") )} else {"________________."}`
1. Using the 'rise over run' idea, the *slope* is approximately 
`r if( knitr::is_html_output() ) {mcq( c(
  -140,
  answer = 1,
  -1,
  50 ))} else {"________________."}`
1. The regression equation is $\hat{y} = -138 + 1.2 x$.
   What does $x$ represent?
`r if( knitr::is_html_output() ) {longmcq( c(
  "The weight of the rowers",
  answer = "The height of the rowers") )} else {"________________."}`
1. What does the 'hat' above the $y$ mean?
`r if( knitr::is_html_output() ) {longmcq( c(
  "That the weights are not measured accurately",
  "That the weights are population values",
  "That the regression model gives poor estimates",
  "That the regression model gives good estimates", 
  answer = "That the regression model gives an estimate of the weight") )}`
1. To two decimal places, what weight would be predicted for a rower who is 180cm tall?
`r if( knitr::is_html_output() ) {fitb( answer = 78, num=TRUE, tol=0.1 )}`
1. The standard error of the slope is $0.112$. What is the value of the *test statistic* (to one decimal place) to test if the population slope is zero?
`r if( knitr::is_html_output() ) {fitb( answer = 10.7, num=TRUE, tol=0.1 )}`
1. True or false? The $P$-value for this test will be *very small*.
`r if( knitr::is_html_output() ) {torf( answer = TRUE)}`
1. True or false? The *units* of the slope are kg/cm.
`r if( knitr::is_html_output() ) {torf( answer = TRUE)}`
1. True or false? Making a prediction for the weight of a rower weighing 220 kg would be an example of *extrapolation*.
`r if( knitr::is_html_output() ) {torf( answer = TRUE)}`
:::




## Exercises {#RegressionExercises}

Selected answers are available in Sect. \@ref(RegressionAnswer).



::: {.exercise #RegressionExerciseManifold}
A study [@amin2019robust] of gas engines measured the throttle angle ($x$) and the manifold air pressure  ($y$) as a fraction of the maximum value.

1. The value of $r$ is given in the paper as $0.972986604$. 
   Comment on this value, and what it means.
1. Comment on the use of a regression model, based on the scatterplot (Fig. \@ref(fig:ThrottleAngle)). 
1. The authors fitted the following regression model: $y = 0.009 + 0.458x$.
   Identify errors that the researchers have made when giving this regression equation (there are more than one).
1. Critique the researchers' approach.
:::



```{r ThrottleAngle, fig.cap="Manifold air pressure plotted against throttle angle for an internal-combustion gas engine", fig.align="center", fig.width=5, fig.height=3.5}

par(mar = c(5.1, 5.1, 4.1, 2.1))

data(Throttle)

plot(Throttle, 
     pch = 19, 
     xlab = "Throttle angle (degrees)", 
     ylab = "Manifold air pressure\n(as fraction of the max.)", 
     ylim = c(0, 1),
     xlim = c(0, 90),
     las = 1)
par(mar = c(5.1, 4.1, 4.1, 2.1))

TAlm <- lm(MAPvalue ~ ThrottleAngle, 
           data = Throttle)
abline(coef(TAlm), 
       lwd = 2, 
       col = "grey")
```


::: {.exercise #RegressionExercisePunting}
(These data were also seen in Exercise \@ref(exr:CorrelationExercisePunting).)
A study (@MyersRegression, p. 75) of American footballers measured the right-leg strengths $x$ of 13 players (using a weight lifting test), and the distance $y$ they punt a football (with their right leg).

1. Use the plot (Fig. \@ref(fig:PuntPlot), left panel) to approximate of the values of the intercept and slope.
1. Using the jamovi output (Fig. \@ref(fig:PuntPlot), right panel), write down the value of the slope ($b_1$) and $y$-intercept ($b_0$).
1. Hence write down the regression equation.
1. Interpret the slope ($b_1$).
1. Write the hypotheses for testing for a relationship in the population
1. Write down the $t$-score and $P$-value from the output.
1. Determine an *approximate* 95% CI for the population slope $\beta_1$.
1. Write a conclusion.
:::


```{r PuntPlot, fig.align="center", fig.cap="Punting distance and right leg strength", fig.width=4, fig.height=3.5, out.width=c("45%", "50%"), fig.show="hold"}
data(Punting)

plot(Punt ~ Right, 
     data = Punting,
     las = 1,
     xlab = "Right leg strength (in pounds)",
     ylab = "Punting distance (in feet)",
     pch = 19,
     xlim = c(110, 180),
     ylim = c(100, 200)
)
abline( coef( Punt.lm <- lm( Punt ~ Right, 
                             data = Punting)), 
        lwd = 2, 
        col = "grey")

knitr::include_graphics("jamovi/Punting/Punting-Regression.png")
```


::: {.exercise #RegressionExercisesGorillas}
(These data were also seen in Exercise \@ref(exr:TwoQuantExercisesGorillas).)
A study [@wright2021chest] examined 25 gorillas are recorded information about their chest beating and their size (measured by the breadth of the gorillas' backs).
The relationship is shown in a scatterplot in Fig. \@ref(fig:GorillaWindmillPlot) (left panel).

Use the output (Fig. \@ref(fig:GorillajamoviHT)) to find the regression line, and perform a hypothesis tests for the slope.
Write a conclusion.

:::

```{r GorillajamoviHT, fig.cap="jamovi regression output for the gorilla data", fig.align="center", out.width="55%"}
knitr::include_graphics("jamovi/Gorillas/GorillasRegression-jamovi.png")
```


```{r}
data(Bitumen)
```


::: {.exercise #RegressionExerciseBitumen}
(These data were also seen in Exercise \@ref(exr:CorrelationExerciseBitumen).)
A study of hot mix asphalt [@data:Panda2018:Bitumen] created $n=42$ samples of asphalt and measured the volume of air voids and the bitumen content by weight (Fig. \@ref(fig:BitumenPlotReg)).
Use the software output (Fig. \@ref(fig:Bitumenjamovi)) to answer these questions.

1. Write down the regression equation.
1. Interpret what the regression equation means.
1. Perform a test to determine is there is a relationship between the variables.
1. Predict the mean percentage of air voids by volume when the percentage bitumen is 5.0%. 
   Do you expect this to be a good prediction? Why or why not?
1. Predict the mean percentage of air voids by volume when the percentage bitumen is 6.0%. 
   Do you expect this to be a good prediction? Why or why not?
:::


```{r BitumenPlotReg, fig.cap="Air voids in bitumen samples", fig.align="center", fig.width=5, fig.height=3.5}
plot(AirVoids ~ Bitumen, 
     data = Bitumen,
     pch = 19,
     xlab = "Actual bitumen content\nby weight (%)",
     ylab = "Air voids by volume (%)",
     las = 1,
     xlim = c(4.7, 5.4),
     ylim = c(3.5, 5.5)
     )
```


```{r BitumenOutputR, results='hide'}
Bitumen.lm <- lm( AirVoids ~ Bitumen, 
                  data = Bitumen)

if( knitr::is_latex_output() ) {
  kable(summary(Bitumen.lm)$coefficient,
        format = "latex",
        booktabs = TRUE,
        longtable = FALSE,
        digits = 4,
        caption = "Regression output for the bitumen data")
}

if( knitr::is_html_output() ) {
  kable(summary(Bitumen.lm)$coefficient,
        format = "html",
        booktabs = TRUE,
        longtable = FALSE,
        digits = 4,
        caption = "Regression output for the bitumen data")
}
```



```{r Bitumenjamovi, fig.cap="jamovi regression output for the bitumen data", fig.align="center", out.width="55%"}
knitr::include_graphics("jamovi/Bitumen/Bitumen-Regression.png")
```




::: {.exercise #RegressionExercisePossums}
(These data were also seen in Exercise \@ref(exr:CorrelationExercisePossums).)
A study of Leadbeater's possums in the Victorian Central Highlands [@data:Williams2022:Possums] recorded, among other information, the body weight of the possums and their location, including the elevation of the location.
A scatterplot of the data is shown in Fig. \@ref(fig:PossumsScatterjamovi).
The jamovi output for fitting the regression line is shown in Fig. \@ref(fig:PossumRegressionjamovi).

1. Write down the regression equation.
1. Interpret what the regression equation means.
1. Perform a test to determine is there is a relationship between the possum weight and teh elevation.
1. Interpret the meaning of the slope.
1. Predict the mean weight of male possums at an elevation of 1000m.
   Do you expect this to be a good prediction? Why or why not?
1. Predict the mean weight of male possums at an elevation of 200m.
   Do you expect this to be a good prediction? Why or why not?
:::

```{r PossumRegressionjamovi, fig.cap="The relationship between weight of possums and the elevation of their location", fig.align="center", fig.width=5, fig.height=3.5, out.width=c("45%", "45%"), fig.show="hold"}
knitr::include_graphics("jamovi/Possums/Possums-Regression-jamovi.png")
```







::: {.exercise #RegressionExerciseSoil}
(These data wee also seen in Exercise \@ref(exr:CorrelationExerciseSoil).)
The *California Bearing Ratio* (CBR) value is used to describe soil-sub grade for flexible pavements (such as in the design of air field runways).
One study [@talukdar2014study] examined the relationship between CBR and other properties of soil, including the plasticity index (PI, a measure of the plasticity of the soil).
The scatterplot from 16 different soil samples from Assam, India, is shown in  Fig. \@ref(fig:SoilPlotCor).

Estimate the regression equation, using the 'rise-over-run' idea.
:::


::: {.exercise #RegressionExerciseSunscreen}
A study of $n = 31$ people on the use of sunscreen [@data:Heerfordt2018:sunscreen] explored the relationship between the time (in minutes) spent on sunscreen application $x$ and the amount (in grams) of sunscreen applied ($y$).
The fitted regression equation was $\hat{y} = 0.27 + 2.21x$.

1. Interpret the meaning of $b_0$ and $b_1$. Do they seems sensible?
1. According to the article, a hypothesis for testing $\beta_0$ produced a $P$-value much larger than $0.05$.
   What does this mean?
1. If someone spent 8 minutes applying sunscreen, how much sunscreen would you predict that they used?
1. The article reports that $R^2=0.64$. 
   Interpret this value.
1. What is the value of the correlation coefficient?
:::


::: {.exercise #RegressionPredictBirthWeight}
One study [@bhargava1985mid] stated:

> In developing countries [...] logistic problems prevent the weighing of every newborn child.
> A study was performed to see whether other simpler measurements could be substituted for weight to identify neonates of low birth weight and those at risk.
>
> --- @bhargava1985mid, p. 1617

The relationship between infant chest circumference (in cm) $x$ and birth weight (in grams) $y$ was given as:  
\[
   \hat{y} = -3440.2403 + 199.2987x.
\]
The correlation coefficient was $r = 0.8696$ with $P < 0.001$.

1. Based on the regression equation only, could chest circumference be used as a useful predictor of birth weight?
   Explain.
1. Based on the correlation information only, could chest circumference be used as a useful predictor of birth weight?
   Explain.
1. Interpret the intercept and the slope of the regression equation.
1. What units of measurement are the intercept and slope measured in?
1. Predict the birth weight of an infant with a chest circumference of 30cm.
1. Critique the way in which the regression equation and correlation coefficient are reported. 
:::



<!-- QUICK REVIEW ANSWERS -->
`r if (knitr::is_html_output()) '<!--'`
::: {.EOCanswerBox .EOCanswer data-latex="{iconmonstr-check-mark-14-240.png}"}
**Answers to in-chapter questions:**

* Sect. \ref{thinkBox:Predictions}: $2.23$g
* Sect. \ref{thinkBox:RegressionCoefs}: $b_0 = 2.1$ and $b_1 = -0.0047$.
  Remember: $b_1$ is the number multiplied by the $x$, and that $b_0$ is the number by itself.
* Sect. \ref{thinkBoxCI}: Recall that all CIs have the same general form.
* Sect. \ref{thinkBoxRegressionConditions}: The relationship looks approximately linear; the variations seems somewhat constant for all values of $x$ (the ONI), but it is a little hard to be sure. 
  The sample size is larger than 25. 
  The conditions for statistical validity appear to have been met.
* \textbf{Answers to \textit{Quick Revision} questions:}
**1.** 1.
**2.** Weight.
**3.** The height of the rowers.
**4.** That the regression model gives an estimate of the weight.
**5.** 78.
**6.** 10.7.
**7.** True.
**8.** True.
**9.** True.
:::
`r if (knitr::is_html_output()) '-->'`