```{r}
#| echo: false
#| message: false
#| warning: false

library(dplyr)
library(ggplot2)
library(broom)
```

# Simple Linear Regression

## Dataset

The dataset **faithful** describes the waiting time between eruptions _(waiting)_ and the duration of the eruption _(eruptions)_ for the Old Faithful geyser in Yellowstone National Park, Wyoming, USA. Both time are in minutes.

```{r}
glimpse(faithful)
```

## Linear Regression Setup

We would like to explore the linear relationship between the duration of the eruption (response variable) and the waiting time between eruptions (explanatory variable).

_Write out the linear equation between response variable and explanatory variable._

## Check the Conditions of Least Squares Regression


```{r}
lm_model <- lm(eruptions ~ waiting, data = faithful)

# The following lines stores residuals and predictions (y_hat) in the faithful data frame
faithful <- faithful |>
  modelr::add_residuals(lm_model) |> 
  modelr::add_predictions(lm_model)

glimpse(faithful)
```

### Linearity

```{r}
ggplot(faithful, aes(x = waiting, y = eruptions)) +
  geom_point() +
  geom_smooth(method = 'lm', se = FALSE)
```

### Constant Variance

```{r}
faithful |> 
  ggplot(aes(x = waiting, y = resid)) +
  geom_point() +
  geom_hline(yintercept = 0)
```



### Normality of Residuals

```{r}
faithful |> 
  ggplot(aes(x = resid)) +
  geom_density()
```

### Independence

_What about the independence of the observations?_

### Interpretation of Results

```{r}
lm_model_summary <- broom::tidy(lm_model)
lm_model_summary
```

_What are the point estimates for slope and intercept? How to interpret them?_

$b_0$ = ?  
$b_1$ = ?

_What is the estimated eruption duration when the waiting time is 60 minutes? What is the estimated eruption duration when the waiting time is 120 minutes?_

### Inference

_Perform the hypothesis testing on whether the linear relationship between the duration of the eruption and waiting time between eruptions exists._

_Do we have significant evidence that the waiting time between eruptions affects the duration of the eruption?_

_Find the $95\%$ confidence interval for the slope. How to interpret this interval?_



```{r}
b1 <- as.numeric(lm_model_summary[2,2]) # This gets the second row and second column of the lm_model_summary table

# b0 <- as.numeric(lm_model_summary[1, 2]) # you can also use this if needed

critical_val <- qt(0.975, df = 271) 
se_b1 <- as.numeric(lm_model_summary[2,3])

me <- critical_val * se_b1
c(b1 - me, b1 + me)
```

```{r}
confint(lm_model)
```


# Multiple Linear Regression
```{r}
#| echo: false
#| message: false
#| warning: false
library(tidyverse)
library(patchwork)
library(broom)
library(modelr)
```
For this example, we explore the power of the log transformation, along with a downside.  
We start by looking at air quality data from New York City, May to September 1973.  
The response of interest here is the Ozone, so we also add the log transformation in the dataset.
Our predictors will be solar.R, wind, and temp.  So the models we will look at are
$$
y_i = \beta_0 + \beta_1x_{1,i} + \beta_2x_{2,i} + \beta_3x_{2,i} + \epsilon,
$$
and 
$$
\log y_i = \beta_0 + \beta_1x_{1,i} + \beta_2x_{2,i} + \beta_3x_{2,i} + \epsilon.
$$
```{r}
air_quility_df = airquality |>
  as_tibble() |>
  mutate(log_ozone = log(Ozone)) |>
  na.omit()
glimpse(air_quility_df)
```

When doing regression is it is always good to start by creating plots of our 
predictors with the response.  This results in the following for the natural scale
Ozone (first row) and log scale Ozone (second row).
```{r}
#| echo: false
#| message: false
#| warning: false
solar_ozone_plot = ggplot(data = air_quility_df, aes(x = Solar.R, y = Ozone)) +
  geom_point() +
  theme_bw()
wind_ozone_plot = ggplot(data = air_quility_df, aes(x = Wind, y = Ozone)) +
  geom_point() +
  theme_bw()
temp_ozone_plot = ggplot(data = air_quility_df, aes(x = Temp, y = Ozone)) +
  geom_point() +
  theme_bw()
solar_log_ozone_plot = ggplot(data = air_quility_df, aes(x = Solar.R, y = log_ozone)) +
  geom_point() +
  theme_bw()
wind_log_ozone_plot = ggplot(data = air_quility_df, aes(x = Wind, y = log_ozone)) +
  geom_point() +
  theme_bw()
temp_log_ozone_plot = ggplot(data = air_quility_df, aes(x = Temp, y = log_ozone)) +
  geom_point() +
  theme_bw()

(solar_ozone_plot + wind_ozone_plot + temp_ozone_plot) /
  (solar_log_ozone_plot + wind_log_ozone_plot + temp_log_ozone_plot)
```

### Question :
What can you notice when be compare the log and natural scaled Ozone plots?  Which appears to be more linear?

## Fitting MLR to Both :
```{r}
model_no_log = lm(Ozone ~ Solar.R + Wind + Temp, data = air_quility_df)
model_log = lm(log_ozone ~ Solar.R + Wind + Temp, data = air_quility_df)
```

Lets look at what our estimates are for regression coefficients with corresponding
p-values.

```{r}
tidy(model_no_log)
tidy(model_log)
```

In our log model, what is the notation that would represent 0.0492?


## Checking Conditions for Estimates
1. Linearity
   a. addressed above
2. Normality of residuals
3. Constant Variance
4. Independence

### Question :
Which conditions does the natural scale Ozone model violate?  Which conditions 
does the log scale Ozone model violate?  *Use the plots below.*
```{r}
#| echo: false
#| message: false
#| warning: false
residuals_no_log_df = tibble(
  fitted_values = fitted.values(model_no_log),
  residuals = residuals(model_no_log)
)
residuals_log_df = tibble(
  fitted_values = fitted.values(model_log),
  residuals = residuals(model_log)
)
res_dens_no_log_plot = ggplot(residuals_no_log_df, aes(x = residuals)) +
  geom_density(color = "dodgerblue", fill = "dodgerblue", alpha = .1) +
  labs(title = "Distribution of residuals",
       subtitle = "for Model with Ozone on natural scale.") +
  theme_bw()
res_dens_log_plot = ggplot(residuals_log_df, aes(x = residuals)) +
  geom_density(color = "dodgerblue", fill = "dodgerblue", alpha = .1) +
  labs(title = "Distribution of residuals",
       subtitle = "for Model with Ozone on log scale.") +
  theme_bw()
res_no_log_plot = ggplot(residuals_no_log_df, aes(x = fitted_values, y = residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "dodgerblue") +
  labs(title = "Residuals vs Fitted Values", 
       subtitle = "for Model with Ozone on natural scale.") +
  theme_bw()
res_log_plot = ggplot(residuals_log_df, aes(x = fitted_values, y = residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "dodgerblue") +
  labs(title = "Residuals vs Fitted Values", 
       subtitle = "for Model with Ozone on log scale.") +
  theme_bw()

(res_dens_no_log_plot + res_dens_log_plot) /
  (res_no_log_plot + res_log_plot)
```

## Interpretation
When interpreting coefficients for a model, we must remember to convert transformed
variables back into their original scale.  So, in the log scale, we need to use
$$
\hat{y_i} = \exp\left[ b_0 + b_1x_{1,i} + b_2x_{2,i} + b_3x_{2,i}\right]
$$
as a basis for interpreting coefficients.

### Question 
Give an interpretation for $b_3$ in both the natural scale and log scale models.


Using notation, write out the hypotheses to test for the slope coefficient for wind. Then, using the output from the log scale model, determine the results of this hypothesis test.

```{r}
confint(model_no_log)
```

Interpret the 95% CI for the Wind slope coefficient.

