---
title: "Linear Regression"
author: 
- name: "Audrey Baguette"
  affiliation: Workshop Lead
- name: "Alex Diaz-Papkovich"
  affiliation: Original Material
- name: "Adrien Osakwe"
  affiliation: Original Material
date: "`r format(Sys.Date(), '%B %d, %Y')`"
---

# Loading libraries
```{r Loading libraries}
library(tidyverse)
```

# Verifying the dataset
Let's take the iris data set, and check if the petal length and width follow a linear correlation
```{r plot data}
ggplot(iris) + geom_point(aes(Petal.Length, Petal.Width, col = Species))
cor(iris$Petal.Length, iris$Petal.Width)
```
The data seems to be linearly correlated. The species have different ranges of sizes, but they seem to all follow the same trend.



# Fitting with MLE
## 1- Estimating beta1
$$\hat{\beta_1} = \frac{\sum^{N}_{i=1}(x_{i} - \bar{x})(y_{i} - \bar{y})}{\sum_{i=1}^{N}(x_{i} - \bar{x})^{2}} = \frac{SS_{xy}}{SS_{xx}}$$

```{r beta1 estimation}
mean_x = mean(iris$Petal.Length)
mean_y = mean(iris$Petal.Width)
numerator = sum((iris$Petal.Length - mean_x) * (iris$Petal.Width - mean_y))
denominator = sum((iris$Petal.Length - mean_x)**2)
beta1 = numerator / denominator
```
## 2- Estimating beta0
$$\hat{\beta_0} = \bar{y} - \hat{\beta_{1}}\bar{x}$$
```{r beta0 estimation}
beta0 = mean_y - beta1*mean_x
```
## 3- Visualize the regression
```{r visualization}
regression_df = data.frame(x = 1:7, y = beta0+c(1:7)*beta1)
ggplot(iris) + geom_point(aes(Petal.Length, Petal.Width, col = Species)) +
            geom_line(aes(x, y), data = regression_df)
```



# Build-in linear regression
lm is the build-in regression function
```{r lm}
linear_model <- lm(Petal.Width ~ Petal.Length, data = iris)
summary(linear_model)
```
Plot the results to compare with the MLE
```{r plot lm}
regression_df2 = data.frame(x = 1:7, y = linear_model$coefficients[1]+c(1:7)*linear_model$coefficients[2])
ggplot(iris) + geom_point(aes(Petal.Length, Petal.Width, col = Species)) +
  geom_line(aes(x, y), data = regression_df2)
```


# Predicting results
Say we remove the versicolor species during model construction,  can we predict the Petal Width in function of its Length?

1- We remove the versicolor data and compute the model
```{r pred 1}
iris_subset <- iris %>% filter(Species != 'versicolor')
iris_model <- lm(Petal.Width ~ Petal.Length, data = iris_subset)
```
2- We extract the versicolor Petal Length
```{r pred 2.1}
versicolor_p <- iris %>% filter(Species == 'versicolor') %>%
                          mutate(Predicted_pw = iris_model$coefficients[1] + Petal.Length*iris_model$coefficients[2])
```
Alternatively
```{r pred 2.2}
predict.lm(iris_model, versicolor_p, type = 'response')
```
3- Plot the results
```{r pred 3}
ggplot(iris) + geom_point(aes(Petal.Length, Petal.Width, col = Species)) +
  geom_point(aes(Petal.Length, Predicted_pw), data = versicolor_p, color = 'darkgreen', shape = '+', size = 5)
```
4- Compute the MSE
```{r pred 4}
MSE = mean((versicolor_p$Petal.Width - versicolor_p$Predicted_pw)**2)
print(MSE)
```


# Using multiple variables
What if we also use the Sepal values to predict the Petal Width?

1- We compute the model after removal of versicolor
```{r multi vars 1}
iris_model2 <- lm(Petal.Width ~ Petal.Length + Sepal.Length + Sepal.Width, data = iris_subset)
```
2- We extract the versicolor Petal Length
```{r multi vars 2}
versicolor_p2 <- iris %>% filter(Species == 'versicolor') %>%
  mutate(Predicted_pw = iris_model2$coefficients[1] +
                        Petal.Length*iris_model2$coefficients[2] +
                        Sepal.Length*iris_model2$coefficients[3] +
                        Sepal.Width*iris_model2$coefficients[4])
```
3- Plot the results
```{r multi vars 3}
ggplot(iris) + geom_point(aes(Petal.Length, Petal.Width, col = Species)) +
  geom_point(aes(Petal.Length, Predicted_pw), data = versicolor_p2, color = 'darkgreen', shape = '+', size = 5)
```

4- Compute the MSE
```{r multi vars 4}
MSE2 = mean((versicolor_p2$Petal.Width - versicolor_p2$Predicted_pw)**2)
print(MSE2)
```
In this case, the MSE is a lot larger than the previous model, so adding more variables did not help. WWhat are the correlations between the variables?
```{r multi vars 5}
ggplot(iris) + geom_point(aes(Sepal.Width, Petal.Width, col = Species))
cor(iris$Sepal.Width, iris$Petal.Width)
ggplot(iris) + geom_point(aes(Sepal.Length, Petal.Width, col = Species))
cor(iris$Sepal.Length, iris$Petal.Width)
```
Quick test with only the Petal Length and Sepal Length
```{r multi vars 6}
iris_model3 <- lm(Petal.Width ~ Petal.Length + Sepal.Length, data = iris_subset)
versicolor_p3 <- iris %>% filter(Species == 'versicolor') %>%
  mutate(Predicted_pw = iris_model3$coefficients[1] +
                        Petal.Length*iris_model3$coefficients[2] +
                        Sepal.Length*iris_model3$coefficients[3])
MSE3 = mean((versicolor_p3$Petal.Width - versicolor_p3$Predicted_pw)**2)
ggplot(iris) + geom_point(aes(Petal.Length, Petal.Width, col = Species)) +
  geom_point(aes(Petal.Length, Predicted_pw), data = versicolor_p3, color = 'darkgreen', shape = '+', size = 5)
print(MSE3)
```
This is the worst MSE of the 3, possibly because of overfitting.



# Hands-on
Create your own linear model

Using mtcars, find the best linear model to predict the miles per gallon (mpg)
Use the first 2/3 of the data as training data for the model and the last third as test data
```{r hands on}
mtcars_train <- mtcars[1:22,]
mtcars_test <- mtcars[23:32,]
```

Solution
```{r hands on solution}
# To help us select a set of variables to use, we can test the correlation of all variables with mpg
mtcars_train %>% apply(2, cor, y=mtcars_train$mpg)
# We have a good negative correlation with the number of cylinders, the displacement,
# horsepower and weight. We will construct a model with those

# 1- Compute the model
cars_model <- lm(mpg ~ cyl+disp+hp+wt, data = mtcars_train)

# 2- Compute the results
cars_p <- mtcars_test %>% 
  mutate(Predicted_mpg = cars_model$coefficients[1] +
                          cyl*cars_model$coefficients[2] +
                          disp*cars_model$coefficients[3] +
                          hp*cars_model$coefficients[4] +
                          wt*cars_model$coefficients[5])

# 3- Compute the MSE
MSE = mean((cars_p$mpg - cars_p$Predicted_mpg)**2)
```
Note: this may not be the overall best solution, it's only one example



#