---
title: "Confidence intervals (Notes)"
subtitle: "STAT 155"
author: "Your Name"
format:
  html:
    toc: true
    toc-depth: 2
    embed-resources: true
---


```{r setup}
#| include: false
knitr::opts_chunk$set(
  collapse = TRUE, 
  warning = FALSE,
  message = FALSE,
  error = TRUE,
  fig.height = 2.75, 
  fig.width = 4.25,
  fig.env = 'figure',
  fig.pos = 'h',
  fig.align = 'center')

# Use a color blind friendly color palette throughout doc
library(tidyverse)
cb_palette <- c("#0072B2", "#D55E00", "black", "#E69F00", "#56B4E9", "#009E73", "#F0E442",  "#CC79A7")
scale_colour_discrete <- function(...) scale_colour_manual(values = cb_palette, ...)
scale_fill_discrete   <- function(...) scale_fill_manual(values = cb_palette, ...)
theme_set(theme_bw())
```








## Exercises {-}

### Exercise 1: Confidence Interval Review

In the first set of exercises, we'll model the time it takes to complete a mountain hike.
To begin, let's explore the relationship of completion `time` (in hours) by hike `length` (in miles):

E[time | length] = $\beta_0$ + $\beta_1$ length

A sample *estimate* of this population model, obtained using our data on hiking trails in the Adirondack mountains, is below:

E[`time` | `length`] = $\hat{\beta}_0$ + $\hat{\beta}_1$ `length` = 2.048 + 0.684 `length`

```{r}
# Import the data & load important packages
library(tidyverse)
peaks <- read.csv("https://mac-stat.github.io/data/high_peaks.csv")

# Model the relationship
peaks_model_1 <- lm(time ~ length, data = peaks)
coef(summary(peaks_model_1))

# Visualize the relationship
peaks %>% 
  ggplot(aes(y = time, x = length)) + 
  geom_point() + 
  geom_smooth(method = "lm")
```

#### Part a

$\hat{\beta}_1$ simply provides a *point estimate*, or our single best guess, of $\beta_1$.
To also produce an *interval estimate*, use the 68-95-99.7 Rule to approximate a 95% CI for $\beta_1$.


#### Part b

We can calculate a more *accurate* CI by applying the `confint()` function to our *model*.
Your approximation from Part a should be close!

```{r}
confint(peaks_model_1, level = 0.95)
```    


### Exercise 2: Interpreting the CI

#### Part a

*Interpreting* the CI for $\beta_1$ in context requires that we can interpret $\beta_1$ itself!
So how can we interpret $\beta_1$ (in general, without assuming a specific value for the unknown $\beta_1$)? Choose 1.

- $\beta_1$ measures the expected completion time for hikes that are 0 miles long
- $\beta_1$ measures the difference in the expected completion time for hikes that long vs hikes that aren't long
- $\beta_1$ measures the change in the expected completion time for each additional 1 mile in length


#### Part b

Per the previous exercise: "We are 95% confident that $\beta_1$ is between 0.56 and 0.81".
Interpret this CI in *context*, drawing on your answer to Part a.


### Exercise 3: time versus elevation

Next, let's explore the relationship between the completion `time` and maximum `elevation` of a hike (in feet): 

E[`time` | `elevation`] = $\beta_0$ + $\beta_1$ `elevation`
        
We can gain insight into this relationship using our sample data:

```{r}
# Model the relationship
peaks_model_2 <- lm(time ~ elevation, data = peaks)

# Get confidence intervals
confint(peaks_model_2, level = 0.95)

# Plot the sample model
ggplot(peaks, aes(y = time, x = elevation)) + 
  geom_point() + 
  geom_smooth(method = "lm")
```

What can we conclude from *both* the CI for $\beta_1$ and the confidence bands in the plot above?    



### Exercise 4: time versus rating

Next, let's explore the relationship between the completion `time` and hike `rating` (difficult, easy, or moderate): 

E[`time` | `rating`] = $\beta_0$ + $\beta_1$ `ratingeasy` + $\beta_2$ `ratingmoderate`

**PAUSE:** Before analyzing this model, pause to reflect on what $\beta_0$, $\beta_1$, and $\beta_2$ *mean* in this context.
We can gain insight into these coefficients using our sample data:

```{r}
# Visualize the relationship
ggplot(peaks, aes(y = time, x = rating)) + 
  geom_boxplot()

# Model the relationship
peaks_model_3 <- lm(time ~ rating, data = peaks)

# Obtain CIs for the model coefficients
confint(peaks_model_3)
```


#### Part a

How can we interpret the confidence interval for $\beta_2$, the `ratingmoderate` coefficient:  (-5.92, -3.19)?
We're 95% confident that... (Choose 1)   

- the average completion time of a moderate hike is between 3.19 and 5.92 hours    
- the average completion time of a moderate hike is between 3.19 and 5.92 hours *less* than that of a difficult hike    
- the average completion time of moderate hikes in our sample was between 3.19 and 5.92 hours *less* than that of difficult hikes in our sample    

#### Part b

How can we interpret the confidence interval for the *intercept* $\beta_0$: (13.8, 16.2)?
We're 95% confident that... (Choose 1)      

- the average completion time of a hike with a 0 rating is between 13.8 and 16.2 hours.    
- the average completion time of a difficult hike is between 13.8 and 16.2 hours.
- the average completion time of a hike is between 13.8 and 16.2 hours.
- the average completion time of a difficult hike is between 13.8 and 16.2 hours *more* than that of an easy hike.

### Exercise 5: Multiple logistic regression REVIEW

Let's turn our attention to weather in 3 Australia locations, Hobart, Uluru, and Wollongong: When controlling for location, is today's 9am humidity level (0-100%) a "significant" predictor of whether or not it will rain tomorrow?
The following population model represents this relationship:

log(odds of `rain_tmrw`) = $\beta_0$ + $\beta_1$ `locationUluru` + $\beta_2$ `locationWollongong` + $\beta_3$ `humidity9am`

Below we build a sample model using a sample of 500 data points:


```{r}
# NOTE: We'll take a sub-sample of 500 days for TEACHING PURPOSES ONLY
# In practice, we'd use all data in our sample :)
set.seed(155)
weather <- read_csv("https://mac-stat.github.io/data/weather_3_locations.csv") %>% 
  mutate(rain_tmrw = (raintomorrow == "Yes")) %>% 
  sample_n(size = 500, replace = FALSE)

# Build the sample model
weather_model_1 <- glm(rain_tmrw ~ location + humidity9am, data = weather, family = "binomial")
coef(summary(weather_model_1))

# Exponentiate the coefficients
exp(coef(weather_model_1))
```

#### Part a

Why are we using logistic regression?


#### Part b

Interpret $\hat{\beta}_0$, the estimated intercept, on the odds scale.
Circle any option that's correct!
**On a day with 0 humidity at 9am in Hobart...**

- The odds of rain are 0.015.
- The chance of rain is 0.015.
- The chance of rain is 1.5% as high as the chance of no rain.
- The odds of rain increase by 0.15%.

#### Part c

Interpret $\hat{\beta}_1$, the estimated `locationUluru` coefficient, on the odds scale.
**When controlling for 9am humidity levels...**

- the odds of rain are 0.60 in Uluru.
- the odds of rain are 60% higher in Uluru than in Hobart.
- the odds of rain in Uluru are 60% as high as (hence 40% lower than) the odds in Hobart.
- the chance of rain in Uluru is 60% as high as (hence 40% lower than) the chance that it doesn't rain.


#### Part d

Interpret $\hat{\beta}_3$, the estimated `humidity9am` coefficient, on the odds scale.
**When controlling for location...**

- the odds of rain increase by 4% for every additional 1-percentage point increase in 9am humidity.
- the odds of rain increase by 1.04 for every additional 1-percentage point increase in 9am humidity.
- the odds of rain are 1.04.


\
\
\
\



### Exercise 6: Inference for logistic regression

Let's focus on $\beta_1$, the "true" `locationUluru` coefficient.
In the previous exercise, we *estimated* $\beta_1$ to be -0.504 (or 0.604 on the odds scale).

#### Part a

*On the log(odds) scale*, the 95% CI for $\beta_1$ is roughly

$$
-0.504 \pm 2*0.414 = (-1.332, 0.324)
$$

More accurately:

```{r}
confint(weather_model_1)
```

What can we conclude?

- We do *not* have evidence that the chance of rain significantly differs in Uluru and Hobart.
- We do *not* have evidence that the chance of rain significantly differs in Uluru and Hobart when controlling for 9am humidity (i.e. when the 9am humidity is the same in both locations).
- We *do* have evidence that the chance of rain significantly differs in Uluru and Hobart.
- We *do* have evidence that the chance of rain significantly differs in Uluru and Hobart when controlling for 9am humidity.


#### Part b

We can also convert the 95% CI to the *odds scale* by exponentiating the limits:

$$
(e^{-1.332}, e^{0.324}) = (0.264, 1.38)
$$

More accurately:

```{r}
exp(confint(weather_model_1))
```

What can we conclude?

- We do *not* have evidence that the chance of rain significantly differs in Uluru and Hobart when controlling for 9am humidity.
- We *do* have evidence that the chance of rain significantly differs in Uluru and Hobart.


#### Part c

Your answers to Parts a and b should agree!
Changing between the log(odds) and odds scales *does* change interpretations, but *doesn't* change our conclusions.
To answer Part a, it was important to check whether **0** was in the CI for $\beta_1$ (on the log(odds) scale).
Equivalently, in Part b, what value did you check for in the CI for $e^{\beta_1}$ (on the odds scale)?




\
\
\
\


### Exercise 7: Comparing models

#### Part a

As shown in Part b of the previous exercise, a 95% CI for the *exponentiated* `humidity9am` coefficient is (1.025, 1.058).
So, is `humidity9am` a "significant" predictor of rain in this model?
REMEMBER: Check for "1" not "0"!


#### Part b

Let's add another predictor to our model: `humidity3pm`, the humidity level at 3pm.

```{r}
weather_model_2 <- glm(rain_tmrw ~ location + humidity9am + humidity3pm, data = weather, family = "binomial")
exp(confint(weather_model_2))
```

Using the 95% CI for the `humidity9am` coefficient (exponentiated), is `humidity9am` a "significant" predictor of rain in this model?


#### Part c

In Parts a and b, you should have concluded that `humidity9am` was a "significant" predictor of rain in `weather_model_1` but *not* in `weather_model_2`.
How could this be?!?




## Done!

- Finalize your notes: (1) Render your notes to an HTML file; (2) Inspect this HTML in your Viewer -- check that your work translated correctly; and (3) Outside RStudio, navigate to your `inclass_activities` subfolder within your `stat155` folder and locate the HTML file -- you can open it again in your browser.
- Clean up your RStudio session: End the rendering process by clicking the 'Stop' button in the 'Background Jobs' pane.
- Check the solutions in the course website, at the bottom of the corresponding chapter.
- Work on homework and/or any extra practice exercises!


