--- title: 'Lab #3: Penalized Regression' output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) set.seed(24) library(MASS) ``` 1. We will be analyzing gene expression data and try to predict whether or not a patient has prostate cancer. Read in the prostate cancer gene expression data "Prostate_GSE6919_U95C.csv" into a data frame called `dat`. The "type" column indicates the patient cancer type. Convert the string to binary values. (The original dataset comes from https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE6919 and was pre-processed using the pipeline from https://sbcb.inf.ufrgs.br/cumida) ```{r} ``` ## Variable subset selection using forward selection Let's try to identify which gene expression values are important for determining whether or not a patient has prostate cancer. For simplicity, we will test out variable selection using the entire dataset. 2. Install the `leaps` package ```{r} ``` 3. Load the `leaps` package ```{r} ``` 4. How many predictors do you have available in the dataset? Are there categorical variables? ```{r} ``` 5. Perform forward selection using the `regsubsets` function. Use `type` as the outcome and all other variables as predictors. Look at a summary of the fitted models. Allow the method to fit models to use up to 5 predictors. To speed things up, try restricting to the first 50 variables. (It's very slow if you run on the entire data matrix. You will also run into issues of collinearity.) ```{r} ``` 6. You can grab the coefficients from variable selection procedure using the `coef` function. Print out the coefficients from the model with the best subset of size 5. ```{r} ``` ## Ridge regression Let's now fit a model that predicts cancer vs normal using ridge regression. 7. Install the package `glmnet` ```{r} ``` 8. Load the `glmnet` package ```{r} ``` 9. Convert the data into components `x` and `y` that are suitable for use in the `glmnet` package. ```{r} ``` 10. Generate a set of values to try for lambda from $\lambda = 10^{5}$ to $\lambda = 10^{-2} = 0.01$. Even space out these values on the log10 scale. ```{r} ``` 11. Fit ridge regression to the dataset for all of the lambda values you have created. ```{r} ``` 12. Look at the value of `lambda` associated with `ridgeMod[75]` ```{r} ``` 13. Plot a histogram with the fitted coefficients associated with this value of `lambda`. ```{r} ``` 14. And then the $l_2$ norm associated with this value of `lambda` ```{r} ``` 15. Repeat steps 17 to 19 for `lambda` associated with index 55. Is the difference in the $l_2$ norm in the direction that you would expect? Why/why not? ```{r} ``` ## Lasso We will now do the same as above but using the Lasso. 16. Fit lasso regression to the dataset for the same set of lambda values you considered in fitting ridge regression ```{r} ``` 17. How many of the estimated coefficients are nonzero for the 85th lambda value? Plot a histogram of the non-zero coefficients for the lasso fit. ```{r} ``` ## Cross-validation Up to now, we've been using the whole dataset to fit the models. Now let's be more formal when fitting and selecting our final model. To do this, we will perform cross-validation. 18. Install and load the `caret` package. Let's tune the lasso penalty parameter `lambda` using the `caret` package. ```{r} ``` 19. To run 5-fold cross-validation, use the `trainControl` method to specify that we'll be running 5-fold cv. ```{r} ``` 20. To perform 5-fold cross-validation using the `train` method in `caret`, we need to specify the grid of hyper-parameter values. Create a data.frame called `caret_grid` with two columns: `lambda` and `alpha`. We'll keep the value of the `alpha` column fixed at 1 because we want to use the lasso. Use the values in `grid` for the `lambda` column. ```{r} ``` 21. Perform 5-fold cross-validation using the `train` method in `caret` over the values in `caret_grid`. Which value of `lambda` did 5-fold CV pick? ```{r} ``` 22. The final model selected by `caret` is given by `cv_model$finalModel`. What are the coefficients in the selected model? How many genes are selected in the final model? ```{r} ```