
---

```{r}
library(dplyr)
library(ggplot2)
```


# 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).

Does the waiting time between eruptions have a linear relationship with the duration of the eruption?

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

$$
eruptions = \beta_0 + \beta_1*waiting + \epsilon
$$

## 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)
```

Looks generally linear (no curves). 

we can also see that it looks like there might be two groups of points.

### Constant Variance

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

if there was constant variance, we would expect the residuals to not have any patterns in the above plot. here, we can see that there is a slight fanning pattern (like < symbol, going from small to wide in terms of how the points are distributed about the line in the middle), so there is not constant variance.

### Normality of Residuals

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

It is not quite normal. it looks like there are almost 2 peaks.

### Independence

_What about the independence of the observations?_

No, it is not fulfilled because all of the observations are from the same geyser. 

### Interpretation of Results

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

$$
\hat{eruptions} = -1.87 + 0.0756*waiting
$$

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

$b_0$ = -1.87

expected y value when our x value is 0

We expect the duration of the eruption to be -1.87 minutes when the waiting time is 0 minutes. 

$b_1$ = 0.0756

expected increase in our y variable for every 1 increase in our x variable.

We expect that for every increase of 1 minute in waiting time, the average duration of the eruption will increase by 0.0756 minutes. 

Note: if b_1 < 0: we expect that for every increase in 1 waiting time, the average duration of the eruption will decrease by ____ minutes. 

_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?_

estimating for 120 minutes is extrapolation because our waiting variable only ranges from 40-100. 

```{r}
-1.87 + 0.0756*60
```


### 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}
# point estimate for beta1
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)
confint(lm_model, level = 0.9)
```


# 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 natrual scaled Ozone plots?  Which 
appears to be more linear?

For all three of the predictors of interest, the relationship with log-Ozone seems much more linear than with natural scale Ozone.

## 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)
```

## Checking Conditions for Estimates
1. Linearity
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 natrual 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 natrual 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)
```

The natural scale Ozone model certainly violates normality and constant variance conditions. The log scale Ozone model seems to more closely satisfy the conditions.

## 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.

In the natural scale, $b_3$ represents the expected change in Ozone for a one degree change in temperature. For the log scale, $b_3$ represents that for a one degree change in temperature, the expected Ozone is multiplied by $e^{b_3}$.

