## Lasso Regression

Let's now use lasso regression to select predictors for the Fish data.

```{r}
rm(list=ls())  # clean up workspace
fish.data <- read.csv("Fish.csv")
# before fitting your model, check the weight parameter values to see if all of them make sense.
# you might want to discard the row with weight = 0.
fish.data <- fish.data[fish.data[, "Weight"] > 0, ]
```

Now let's divide the data into a training set (~80% of original observations) and a validation set (~20% of original observations).

```{r}
library(glmnet)
library(caret)
set.seed(6040)

training_samples <- createDataPartition(fish.data$Weight, p = 0.8, list = FALSE)
```

Let's start from Model 8 from last lab session:
`log(Weight) ~ Species * (Length1 + Length2 + Length3 + Height + Width)`.
The `glmnet` package takes a matrix (of predictors) and a vector (of responses) as input.
We use `model.matrix` function to create them. 
`glmnet` will add intercept by default, so we drop intercept term when forming `x` matrix.

```{r}
# X and y from original data
x_all <- model.matrix(
  log(Weight) ~ Species * (Length1 + Length2 + Length3 + Height + Width), 
  data = fish.data)
y_all <- log(fish.data$Weight)
# training X and y
x_train <- x_all[training_samples[, 1], ]
y_train <- y_all[training_samples[, 1]]
# validation X and y
x_val <- x_all[-training_samples[, 1], ]
y_val <- y_all[-training_samples[, 1]]
```

### Now fit the lasso regression and plot solution path:

```{r}
# fill the code here
```


- Now we can evaluate the performance of the models (corresponding to different $\lambda$ values) on the validation set.

```{r}
# fill your code here
```

- Now the question is whether we should make decision on this single training-validation split? A common strategy is to use **k-fold cross validation**.
```{r}
set.seed(7260)
cv_lasso <- cv.glmnet(x_all, y_all, alpha = 1, type.measure = "deviance", nfolds = 5)  # check the documentation of the argument type.measure to find that "deviance" corresponds to squared-error for guassian models.
plot(cv_lasso)
```

The plot displays the cross-validation error (deviance by default) according to the log of lambda. The left dashed vertical line indicates that the log of the optimal value of log lambda is approximately -7.8, which is the one that minimizes the average deviance. This lambda value will give the most accurate model. The exact value of lambda and corresponding model can be viewed as follow:
```{r}
cv_lasso$lambda.min
coef(cv_lasso, cv_lasso$lambda.min)
```

### What is the difference between the above two procedures?



--------------------------------------------------------

## ANOVA 

Reference: [Two-way ANOVA test in R](http://www.sthda.com/english/wiki/two-way-anova-test-in-r)

Let's use R to perform a Two-way ANOVA test on the honeybee data to test what you have calculated in class.

Data preparation:
```{r}
library(tidyverse)
honey.bee <- matrix(c(20, 20, 3.1, 3.7, 4.7,
		20, 40, 5.5, 6.7, 7.3,
		20, 60, 7.9, 9.2, 9.3,
		30, 20, 6, 6.9, 7.5,
		30, 40, 11.5, 12.9, 13.4,
		30, 60, 17.5, 15.8, 14.7,
		40, 20, 7.7, 8.3, 9.5,
		40, 40, 15.7, 14.3, 15.9,
		40, 60, 19.1, 18.0, 19.9), 9, 5, byrow = TRUE)
colnames(honey.bee) <- c("Temp", "Suc", "rep.1", "rep.2", "rep.3")

honey.bee.df <- as_tibble(honey.bee) %>%  # %>% is the pipe operator, read more about it here: https://r4ds.had.co.nz/pipes.html
  pivot_longer(cols = c("rep.1", "rep.2", "rep.3"), names_to = NULL, values_to = "Consumption") %>%
  mutate(Temp = as.factor(Temp), Suc = as.factor(Suc))  # change variable types to factor
```

### Use the function `aov` to perform ANOVA test.  You may refer to the reference link above.

```{r}
res.aov <- aov(Consumption ~ Temp * Suc, data = honey.bee.df)
summary(res.aov)
```

### Use the function `interaction.plot` to generate an interaction plot with `Temp` as the x factor and `Suc` as the trace factor.


### Use the function `ggline` from `ggpubr` packge to generate the same interaction plot.

```{r}
library("ggpubr")
```


--------------------------------------------------------------------------------
## [Optional practice] Implementing the Metropolis-Hastings algorithm
Requested by Sang-Eun (last year).

You may practice implementing the Metropolis-Hastings Algorithm and apply it to sample from some distributions.

Dr. Vladimir Minin has a great lab session on doing so.

- You should first read about the lab session [here](https://github.com/vnminin/SISMID_MCMC_I/blob/master/2020/labs/mh-lab.pdf).

- You can then fill out the code in this R script [here for sampling the standard normal distribution](https://github.com/vnminin/SISMID_MCMC_I/blob/master/2020/code/norm_mh_reduced.R) and [here for sampling from the distribution of the time of infection](https://github.com/vnminin/SISMID_MCMC_I/blob/master/2020/code/infect_time_reduced.R)

- Check your implementation and results with the filled version [here for sampling the standard normal distribution](https://github.com/vnminin/SISMID_MCMC_I/blob/master/2020/code/norm_mh.R) and [here for sampling from the distribution of the time of infection](https://github.com/vnminin/SISMID_MCMC_I/blob/master/2020/code/infect_time.R)

- For more materials, you can visit their course site [here](http://vnminin.github.io/SISMID_MCMC_I/)



