b2_regularization.R

#' ---
#' title: "Regularization methods"
#' author: "Jan-Philipp Kolb"
#' date: "`r format(Sys.time(), '%d %B, %Y')`"
#' fontsize: 10pt
#' output: 
#'   beamer_presentation: 
#'     theme: Dresden
#'     colortheme: dolphin
#'     fig_height: 3
#'     fig_width: 5
#'     fig_caption: no
#'     fonttheme: structuresmallcapsserif
#'     highlight: haddock
#' ---
#' 
## ----setup, include=FALSE------------------------------------------------
knitr::opts_chunk$set(echo = T,warning=F,message=F)

#' 
#' 
#' ## Insufficient Solution
#' 
#' - When the number of features exceed the number of observations ($p>n$), the OLS solution matrix is not invertible. 
#' - This causes significant issues because it means: 
#' 
#' (1) The least-squares estimates are not unique. There are an infinite set of solutions available and most of these solutions overfit the data. 
#' 
#' (2) In many instances the result will be computationally infeasible.
#' 
#' - To resolve this issue we can remove variables until $p<n$ and then fit an OLS regression model. 
#' - Although we can use pre-processing tools to apply this manual approach ([**Kuhn and Johnson**](http://appliedpredictivemodeling.com/), 2013, pp. 43-47), it can be cumbersome and prone to errors.
#' 
#' 
#' ## Regularized Regression
#' 
#' - When we experience these concerns, one alternative is to use regularized regression (also commonly referred to as penalized models or [**shrinkage methods**](https://gerardnico.com/lang/r/ridge_lasso)) to control the parameter estimates. 
#' 
#' - Regularized regression puts contraints on the magnitude of the coefficients and will progressively shrink them towards zero. This constraint helps to reduce the magnitude and fluctuations of the coefficients and will reduce the variance of our model.
#' 
#' 
#' 
#' ## [Regularization](https://elitedatascience.com/algorithm-selection)
#' 
#' ### elitedatascience.com definition
#' 
#' Regularization is a technique used to prevent overfitting by artificially penalizing model coefficients.
#' 
#' - It can discourage large coefficients (by dampening them).
#' -  It can also remove features entirely (by setting their coefficients to 0).
#' -  The "strength" of the penalty is tunable. 
#' 
#' ### Wikipedia definition of [Regularization](https://en.wikipedia.org/wiki/Regularization_(mathematics)) 
#' Regularization is the process of adding information in order to solve an ill-posed problem or to prevent overfitting. 
#' 
#' 
#' <!--
#' Rmarkdown tipps and tricks
#' no figure caption
#' https://stackoverflow.com/questions/38514954/removing-figure-text-in-rmarkdown
#' syntax highlighting
#' https://eranraviv.com/syntax-highlighting-style-in-rmarkdown/
#' https://latex-kurs.blogspot.com/2012/09/latex-plusminus.html
#' -->
#' 
#' 
#' ## Strenghts and weaknesses of [regularization](https://elitedatascience.com/machine-learning-algorithms)
#' 
#' <!--
#' Regularization is a technique for penalizing large coefficients in order to avoid overfitting, and the strength of the penalty should be tuned.
#' -->
#' 
#' ### Strengths: 
#' 
#' - Linear regression is straightforward to understand and explain, and can be regularized to avoid overfitting. 
#' - Linear models can be updated easily with new data
#' 
#' ### Weaknesses: 
#' 
#' - Linear regression in general performs poorly when there are non-linear relationships. - They are not naturally flexible enough to capture more complex patterns, and adding the right interaction terms or polynomials can be tricky and time-consuming.
#' 
#' ## [Three regularized regression algorithms](https://elitedatascience.com/algorithm-selection)
#' 
#' ![](figure/reg_3algos.PNG){ height=40% }
#' 
#' ### Lasso regression
#' 
#' - Absolute size of coefficients is penalized. 
#' - Coefficients can be exactly 0.
#' 
#' ### Ridge regression
#' 
#' - Squared size of coefficients is penalized.
#' - Smaller coefficients, but it doesn't force them to 0.
#' 
#' ### Elastic-net
#' 
#' - A mix of both absolute and squared size is penalzied.
#' <!--
#' - The ratio of the two penalty types should be tuned.
#' -->
#' 
#' 
#' ## The objective function of regularized regression methods...
#' 
#' - is very similar to OLS regression; 
#' - And a penalty parameter (P) is added.
#' 
#' $$
#' \text{minimize}\{SSE+P\}
#' $$
#'  
#' - There are two main penalty parameters, which have a similar effect. 
#' - They constrain the size of the coefficients such that the only way the coefficients can increase is if we experience a comparable decrease in the sum of squared errors (SSE). 
#' 
#' 
#' 
#' ## Preparations
#' 
#' - Most of the following slides are based on the [**UC Business Analytics R Programming Guide**](http://uc-r.github.io/regularized_regression)
#' 
## ----eval=F,echo=F-------------------------------------------------------
## install.packages("rsample")

#' 
#' ### Necessary packages
#' 
## ------------------------------------------------------------------------
library(rsample)  # data splitting 
library(glmnet)   # implementing regularized regression approaches
library(dplyr)    # basic data manipulation procedures
library(ggplot2)  # plotting
library(knitr)    # for tables

#' 
#' <!--
#' http://jse.amstat.org/v19n3/decock/DataDocumentation.txt
#' http://jse.amstat.org/v19n3/decock.pdf
#' -->
#' 
#' ## The example dataset
#' 
## ------------------------------------------------------------------------
library(AmesHousing) 
ames_data <- AmesHousing::make_ames()

#' 
## ----eval=F,echo=F-------------------------------------------------------
## kable(head(ames_data[,1:4]))

#' 
## ----eval=F,echo=F-------------------------------------------------------
## library(DT)
## DT::datatable(ames_data)

#' 
#' ![](figure/dt_amesdata.PNG)
#' 
#' ## Create training (70%) and test (30%) sets 
#' 
#' - `set.seed` is used for reproducibility
#' - `initial_split` is used to split data in training and test data
#' 
## ------------------------------------------------------------------------
set.seed(123)
ames_split <- rsample::initial_split(ames_data, prop = .7, 
                            strata = "Sale_Price")
ames_train <- rsample::training(ames_split)
ames_test  <- rsample::testing(ames_split)

#' 
#' <!--
#' ## Multicollinearity
#' 
#' - As p increases we are more likely to capture multiple features that have some multicollinearity. 
#' - When multicollinearity exists, we often see high variability in our coefficient terms. 
#' - E.g. we have a correlation of 0.801 between `Gr_Liv_Area` and `TotRms_AbvGrd` 
#' - Both variables are strongly correlated to the response variable (`Sale_Price`).
#' 
## ------------------------------------------------------------------------
cor(ames_data[,c("Sale_Price","Gr_Liv_Area","TotRms_AbvGrd")])

#' 
#' 
#' ## Multicollinearity
#' 
## ------------------------------------------------------------------------
lm(Sale_Price ~ Gr_Liv_Area + TotRms_AbvGrd, data = ames_train)

#' 
#' - When we fit a model with both these variables we get a positive coefficient for `Gr_Liv_Area` but a negative coefficient for `TotRms_AbvGrd`, suggesting one has a positive impact to Sale_Price and the other a negative impact.
#' 
#' ## Seperated models
#' 
#' - If we refit the model with each variable independently, they both show a positive impact. 
#' - The `Gr_Liv_Area` effect is now smaller and the `TotRms_AbvGrd` is positive with a much larger magnitude.
#' 
## ------------------------------------------------------------------------
lm(Sale_Price ~ Gr_Liv_Area, data = ames_train)$coefficients

#' 
## ------------------------------------------------------------------------
lm(Sale_Price ~ TotRms_AbvGrd, data = ames_train)$coefficients

#' 
#' - This is a common result when collinearity exists. 
#' - Coefficients for correlated features become over-inflated and can fluctuate significantly. 
#' 
#' 
#' ## Consequences
#' 
#' - One consequence of these large fluctuations in the coefficient terms is [**overfitting**](https://en.wikipedia.org/wiki/Overfitting), which means we have high variance in the bias-variance tradeoff space. 
#' - We can use tools such as variance inflaction factors (Myers, 1994) to identify and remove those strongly correlated variables, but it is not always clear which variable(s) to remove. 
#' - Nor do we always wish to remove variables as this may be removing signal in our data.
#' -->
#' 
#' <!--
#' ## Interpretability
#' 
#' - With a large number of features, we often would like to identify a smaller subset that has the strongest effects. 
#' - We sometimes prefer techniques that provide feature selection. One approach to this is called hard threshholding feature selection, which can be performed with linear model selection approaches. 
#' - Model selection approaches can be computationally inefficient, do not scale well, and they simply assume a feature as in or out. 
#' - We may wish to use a soft threshholding approach that slowly pushes a feature’s effect towards zero. As will be demonstrated, this can provide additional understanding regarding predictive signals.
#' -->
#' 
#' ## Ridge Regression
#' 
#' - Ridge regression (Hoerl, 1970) controls the coefficients by adding
#' $\lambda\sum_{j=1}^p\beta_j^2$ to the objective function. 
#' 
#' - This penalty parameter is referred to as "$L_2$" as it signifies a second-order penalty being used on the coefficients.
#' 
#' $$
#' \text{minimize}\{\text{SSE}+\lambda\sum_{j=1}^p\beta_j^2 \}
#' $$
#' 
#' - This penalty parameter can take on a wide range of values, which is controlled by the tuning parameter $\lambda$. 
#' 
#' - When $\lambda=0$, there is no effect and our objective function equals the normal OLS regression objective function of simply minimizing SSE. 
#' 
#' - As $\lambda \rightarrow \infty$, the penalty becomes large and forces our coefficients to zero. 
#' 
#' 
#' ## Exemplar coefficients
#' 
#' Exemplar coefficients have been regularized with $\lambda$ ranging from 0 to over 8,000 (log(8103)=9).
#' 
#' ![](figure/ridge_coef.png)
#' 
#' <!--
#' Für was stehen die unterschiedlichen Linien
#' 
#' Jede Linie müsste für einen Koeffizienten im Regressionsmodell stehen.
#' 
#' Siehe hier: https://beta.vu.nl/nl/Images/werkstuk-fonti_tcm235-836234.pdf
#' -->
#' 
#' 
#' ## How to choose the right $\lambda$ 
#' 
#' - Although these coefficients were scaled and centered prior to the analysis, you will notice that some are extremely large when $\lambda\rightarrow 0$. 
#' 
#' - We have a large negative parameter that fluctuates until $log(\lambda)\approx 2$ where it then continuously shrinks to zero. 
#' 
#' - This is indicates multicollinearity and likely illustrates that constraining our coefficients with $log(\lambda)>2$ may reduce the variance, and therefore the error, in our model.
#' 
#' - But how do we find the amount of shrinkage (or $\lambda$) that minimizes our error? 
#' 
#' ## Implementation in `glmnet`
#' 
## ----fig.height=2.5------------------------------------------------------
plot(density(ames_data$Sale_Price),main="")

#' 
#' - `glmnet` does not use the formula method (y ~ x) so prior to modeling we need to create our feature and target set.
#' - The `model.matrix` function is used on our feature set, which will automatically dummy encode qualitative variables 
#' - We also log transform our response variable due to its skeweness.
#' 
#' 
#' <!--
#' infinity symbol in LateX
#' https://praxistipps.chip.de/latex-unendlich-zeichen-eingeben-so-gehts_92332
#' -->
#' 
#' ## Training and testing feature model matrices and response vectors.
#' 
#' - We use `model.matrix(...)[, -1]` to discard the intercept
#' 
## ------------------------------------------------------------------------
ames_train_x <- model.matrix(Sale_Price ~ ., ames_train)[, -1]
ames_train_y <- log(ames_train$Sale_Price)

ames_test_x <- model.matrix(Sale_Price ~ ., ames_test)[, -1]
ames_test_y <- log(ames_test$Sale_Price)

# What is the dimension of of your feature matrix?
dim(ames_train_x)

#' 
#' 
#' ## Behind the scenes 
#' 
#' - The alpha parameter tells `glmnet` to perform a Ridge ($\alpha = 0$), Lasso ($\alpha = 1$), or Elastic Net $(0\leq \alpha \leq 1)$ model. 
#' - Behind the scenes, `glmnet` is doing two things that you should be aware of:
#' 
#' (1.) It is essential that predictor variables are standardized when performing regularized regression. `glmnet` performs this for you. If you standardize your predictors prior to `glmnet` you can turn this argument off with `standardize=FALSE`.
#' 
#' (2.) `glmnet` will perform Ridge models across a wide range of $\lambda$ parameters, which are illustrated in the figure on the next slide.
#' 
## ------------------------------------------------------------------------
ames_ridge <- glmnet(x = ames_train_x,y = ames_train_y,
  alpha = 0)

#' 
#' 
#' ## A wide range of $\lambda$ parameters
#' 
## ------------------------------------------------------------------------
plot(ames_ridge, xvar = "lambda")

#' 
#' 
#' ## $\lambda$ values in `glmnet`
#' 
#' - We can see the exact $\lambda$ values applied with `ames_ridge$lambda`. 
#' - You can specify your own $\lambda$ values, 
#' - By default `glmnet` applies 100 $\lambda$ values that are data derived. 
#' - Normally you will have little need to adjust the default $\lambda$ values.
#' 
## ------------------------------------------------------------------------
head(ames_ridge$lambda)

#' 
#' ## Access the coefficients with `coef`.
#' 
#' - The coefficients are stored for each model in order of largest to smallest $\lambda$. 
#' - The coefficients for the `Gr_Liv_Area` and `TotRms_AbvGrd` features for the largest $\lambda$ (279.1035) and smallest $\lambda$ (0.02791035) are visible. 
#' - The largest $\lambda$ value has pushed these coefficients to nearly 0.
#' 
#' <!--
#' # coefficients for the largest and smallest lambda parameters
#' ## Coefficients of the `ames_ridge` model 
#' -->
#' 
## ------------------------------------------------------------------------
coef(ames_ridge)[c("Gr_Liv_Area", "TotRms_AbvGrd"),100]

#' 
## ------------------------------------------------------------------------
coef(ames_ridge)[c("Gr_Liv_Area", "TotRms_AbvGrd"), 1] 

#' 
#' - But how much improvement we are experiencing in our model.
#' 
#' ## Tuning
#' 
#' - Recall that $\lambda$ is a tuning parameter that helps to control our model from over-fitting to the training data. 
#' - To identify the optimal $\lambda$ value we need to perform cross-validation (CV). 
#' - `cv.glmnet` provides a built-in option to perform k-fold CV, and by default, performs 10-fold CV.
#'  
## ------------------------------------------------------------------------
ames_ridge <- cv.glmnet(x = ames_train_x,y = ames_train_y,
  alpha = 0)

#' 
#' ## Results of cv Ridge regression
#' 
## ------------------------------------------------------------------------
plot(ames_ridge)

#' 
#' - The plot illustrates the 10-fold CV mean squared error (MSE) across the $\lambda$ values. 
#' - We see no substantial improvement; 
#' 
#' ## The plot explained (I)
#' 
#' - As we constrain our coefficients with $log(\lambda)\leq 0$ penalty, the MSE rises considerably. 
#' - The numbers at the top of the plot (301) just refer to the number of variables in the model. 
#' - Ridge regression does not force any variables to exactly zero so all features will remain in the model.
#' 
#' ## The plot explained (II)
#' 
#' <!--
#' - The first and second vertical dashed lines represent the $\lambda$  value with the minimum MSE and the largest $\lambda$ value within one standard error of the minimum MSE.
#' -->
#' 
#' 
## ------------------------------------------------------------------------
min(ames_ridge$cvm)       # minimum MSE
ames_ridge$lambda.min     # lambda for this min MSE
# 1 st.error of min MSE
ames_ridge$cvm[ames_ridge$lambda == ames_ridge$lambda.1se]  
ames_ridge$lambda.1se  # lambda for this MSE

#' 
#' 
#' ## The plot explained (III)
#' 
#' - The advantage of identifying the $\lambda$ with an MSE within one standard error becomes more obvious with the Lasso and Elastic Net models. 
#' - For now we can assess this visually. 
#' - We plot the coefficients across the $\lambda$ values and the dashed red line represents the largest $\lambda$ that falls within one standard error of the minimum MSE. 
#' - This shows you how much we can constrain the coefficients while still maximizing predictive accuracy.
#' 
#' 
#' 
## ------------------------------------------------------------------------
ames_ridge_min <- glmnet(x = ames_train_x,y = ames_train_y,
  alpha = 0)

#' 
#' ## Coefficients across the $\lambda$ values
#' 
## ------------------------------------------------------------------------
plot(ames_ridge_min, xvar = "lambda")
abline(v = log(ames_ridge$lambda.1se), col = "red", 
       lty = "dashed")

#' 
#' 
#' 
#' ## Advantages and Disadvantages
#' 
#' - The Ridge regression model has pushed many of the correlated features towards each other rather than allowing for one to be wildly positive and the other wildly negative. 
#' - Many of the non-important features have been pushed closer to zero. 
#' - We have reduced the noise in our data $\Rightarrow$ more clarity in identifying the true signals.
#' 
#' 
## ----eval=F--------------------------------------------------------------
## coef(ames_ridge, s = "lambda.1se") %>%
##   filter(row != "(Intercept)") %>%
##   top_n(25, wt = abs(value)) %>%
##   ggplot(aes(value, reorder(row, value))) +
##   geom_point() +
##   ggtitle("Top 25 influential variables") +
##   xlab("Coefficient") +
##   ylab(NULL)

#' 
#' 
#' ## Top 25 influential variables
#' 
#' ![](figure/ridgeTop25influentalVars.PNG)
#' 
#' ## Exercise: ridge regression (I)
#' 
#' 1) Load the `lars` package and the `diabetes` dataset 
#' <!--
#' (Efron, Hastie, Johnstone and Tibshirani (2003) “Least Angle Regression” (with discussion) Annals of Statistics). 
#' This is the same dataset from the LASSO exercise set and has patient level data on the progression of diabetes. 
#' -->
#' 2) Load the `glmnet` package to implement ridge regression.
#' 
#' The dataset has three matrices x, x2 and y. x has a smaller set of independent variables while x2 contains the full set with quadratic and interaction terms. y is the dependent variable which is a quantitative measure of the progression of diabetes.
#' 
#' 3) Generate separate scatterplots with the line of best fit for all the predictors in x with y on the vertical axis.
#' 
#' 4) Regress y on the predictors in x using OLS. We will use this result as benchmark for comparison.
#' 
#' ## Exercise: ridge regression (II)
#' <!-- 
#' Exercise 2
#' -->
#' 
#' 5) Fit the ridge regression model using the `glmnet` function and plot the trace of the estimated coefficients against lambdas. Note that coefficients are shrunk closer to zero for higher values of lambda.
#' 
#' <!--
#' Exercise 3
#' -->
#' 
#' 6) Use the cv.glmnet function to get the cross validation curve and the value of lambda that minimizes the mean cross validation error.
#' 
#' <!--
#' Exercise 4
#' -->
#' 
#' 7) Using the minimum value of lambda from the previous exercise, get the estimated beta matrix. Note that coefficients are lower than least squares estimates.
#' 
#' <!--
#' Exercise 5
#' -->
#' 
#' 8) To get a more parsimonious model we can use a higher value of lambda that is within one standard error of the minimum. Use this value of lambda to get the beta coefficients. Note the shrinkage effect on the estimates.
#' 
#' ## Exercise: ridge regression (III)
#' 
#' <!--
#' Exercise 6
#' -->
#' 
#' 9) Split the data randomly between a training set (80%) and test set (20%). We will use these to get the prediction standard error for least squares and ridge regression models.
#' 
#' <!--
#' Exercise 7
#' -->
#' 
#' 10) Fit the ridge regression model on the training and get the estimated beta coefficients for both the minimum lambda and the higher lambda within 1-standard error of the minimum.
#' 
#' <!--
#' Exercise 8
#' -->
#' 
#' 11) Get predictions from the ridge regression model for the test set and calculate the prediction standard error. Do this for both the minimum lambda and the higher lambda within 1-standard error of the minimum.
#' 
#' <!--
#' Exercise 9
#' -->
#' 
#' 12) Fit the least squares model on the training set.
#' 
#' 
#' <!-- 
#' Exercise 10
#' -->
#' 
#' 13) Get predictions from the least squares model for the test set and calculate the prediction standard error.
#' 
#' 
#' 
#' 
#' ## Ridge and Lasso
#' 
#' <!--
#' A ridge model will retain all variables.
#' 
#' -->
#' 
#' ### A Ridge model...
#' 
#' - ... is good if we need to retain all features, yet reduce the noise that less influential variables may create and minimize multicollinearity. 
#' - ... does not perform feature selection. If greater interpretation is necessary where you need to reduce the signal in your data to a smaller subset then a Lasso model may be preferable.
#' 
#' <!--
#' - LASSO is a feature selection method.
#' 
#' https://eight2late.wordpress.com/2017/07/11/a-gentle-introduction-to-logistic-regression-and-lasso-regularisation-using-r/
#' 
#' - LASSO regression has inbuilt penalization functions to reduce overfitting.
#' 
#' https://www.analyticsvidhya.com/blog/2016/12/introduction-to-feature-selection-methods-with-an-example-or-how-to-select-the-right-variables/
#' -->
#' 
#' 
#' <!--
#' ## [Lasso regression](https://elitedatascience.com/algorithm-selection)
#' 
#' ### LASSO, stands for least absolute shrinkage and selection operator 
#' 
#' -  Lasso regression penalizes the absolute size of coefficients.
#' -   Practically, this leads to coefficients that can be exactly 0.
#' -   Thus, Lasso offers automatic feature selection because it can completely remove some features.
#' -   The "strength" of the penalty should be tuned.
#' -  A stronger penalty leads to more coefficients pushed to zero.
#' 
#' 
#' ## [Lasso](https://en.wikipedia.org/wiki/Lasso_(statistics)) regression overview
#' 
#' - Lasso is a regression analysis method that performs variable selection and regularization (reduce overfitting)
#' - We want to enhance prediction accuracy and interpretability of the statistical model.
#' 
#' <!--
#' https://eight2late.wordpress.com/2017/07/11/a-gentle-introduction-to-logistic-regression-and-lasso-regularisation-using-r/
#' -->
#' 
#' - We could remove less important variables.
#' - We can do that manually by examining p-values of coefficients and discarding those variables whose coefficients are not significant.
#' - But this can become tedious for classification problems with many independent variables
#' 
#' 
#' 
#' <!--
#' ## History of Lasso
#' 
#' - Originally introduced in geophysics literature in 1986
#' - Independently rediscovered and popularized in 1996 by Robert Tibshirani, who coined the term and provided further insights into the observed performance.
#' 
#' Lasso was originally formulated for least squares models and this simple case reveals a substantial amount about the behavior of the estimator, including its relationship to ridge regression and best subset selection and the connections between lasso coefficient estimates and so-called soft thresholding. It also reveals that (like standard linear regression) the coefficient estimates need not be unique if covariates are collinear.
#' -->
#' 
#' <!--
#' ## What is [lasso regression](http://www.statisticshowto.com/lasso-regression/)
#' 
#' - Lasso regression uses shrinkage
#' - Data values are shrunk towards a central point
#' 
#' - [Ridge and lasso regularization work by adding a penalty term to the log likelihood function.](https://eight2late.wordpress.com/2017/07/11/a-gentle-introduction-to-logistic-regression-and-lasso-regularisation-using-r/)
#' 
#' - A tuning parameter, $\lambda$ controls the strength of the L1 penalty.
#' 
#' $$
#' \sum\limits_{i=1}^n \big( y_i -\beta_0 - \sum\limits_{j=1}^p \beta_jx_{ij} \big)^2 + \lambda \sum\limits_{j=1}^p |\beta_j| = RSS + \lambda\sum\limits_{j=1}^p |\beta_j|.
#' $$
#' 
#' -->
#' <!--
#' wir haben einen penalty term, der hoch ist, wenn die Parameterschätzwerte hoch sind.
#' 
#' Youtube Video zu Lasso
#' https://www.youtube.com/watch?v=A5I1G1MfUmA
#' -->
#' 
#' ## Lasso Regression
#' 
#' - Originally introduced in geophysics literature in 1986
#' - The least absolute shrinkage and selection operator (Lasso) model was rediscovered and popularized in [**1996 by Robert Tibshirani**](https://www.jstor.org/stable/2346178?seq=1#page_scan_tab_contents) 
#' - It is an alternative to Ridge regression that has a small modification to the penalty in the objective function.
#' - Rather than the $L_2$ penalty we use the following $L_1$ penalty $\lambda\sum_{j=1}^p |\beta_j|$ in the objective function.
#' 
#' $$
#' \text{minimize} \{\text{SSE}+\lambda\sum_{j=1}^p |\beta_j|\}
#' $$
#' 
#' 
#' ## Lasso penalty pushes coefficients to zero
#' 
#' ![Lasso regression coefficients as $\lambda$ grows from $0 \rightarrow \infty$.](figure/fig3_loglambda.PNG){ heigth=60% }
#' 
#' Lasso improves the model with regularization and conducts automated feature selection.
#' 
#' ## The reduction of coefficients
#' 
#' - 15 variables for $\text{log}(\lambda)=-5$
#' - 12 variables for $\text{log}(\lambda)=-1$
#' - 3 variables for $\text{log}(\lambda)=1$
#' 
#' When a data set has many features, Lasso can be used to identify and extract those features with the largest (and most consistent) signal.
#' 
#' 
#' ### Implementation Lasso regression to ames data
#' 
#' - Implementing Lasso follows the same logic as implementing the Ridge model, we just need to switch $\alpha = 1$ within `glmnet`.
#' 
## ------------------------------------------------------------------------
ames_lasso<-glmnet(x=ames_train_x,y=ames_train_y,alpha=1)

#' 
#' ## A quick drop in number of features
#' 
#' - Very large coefficients for ols (highly correlated)
#' - As model is constrained - these noisy features are pushed to 0.
#' - CV is necessary to determine right value for $\lambda$
#' 
## ----fig.height=3--------------------------------------------------------
plot(ames_lasso, xvar = "lambda")

#' 
#' ## Tuning with `cv.glmnet`
#' 
#' - `cv.glmnet` with `alpha=1` is used to perform cv.
#' 
## ------------------------------------------------------------------------
ames_lasso<-cv.glmnet(x=ames_train_x,y=ames_train_y,alpha=1)
names(ames_lasso)

#' 
#' ## MSE for cross validation
#' 
#' - MSE can be minimized with $-6\leq log (\lambda) \leq -4$
#' - Also the number of features can be reduced ($156 \leq p \leq 58$)
#' 
## ------------------------------------------------------------------------
plot(ames_lasso)

#' 
#' ## Minimum and one standard error MSE and $\lambda$ values.
#' 
## ------------------------------------------------------------------------
min(ames_lasso$cvm) # minimum MSE
ames_lasso$lambda.min # lambda for this min MSE

# 1 st.error of min MSE
ames_lasso$cvm[ames_lasso$lambda == ames_lasso$lambda.1se]  
ames_lasso$lambda.1se  # lambda for this MSE

#' 
#' 
#' ## MSE within one standard error
#' 
#' - The advantage of identifying the $\lambda$ with an MSE within one standard error becomes more obvious. 
#' - If we use the $\lambda$ that drives the minimum MSE we can reduce our feature set from 307 down to less than 160. 
#' - There is some variability with this MSE and we can assume that we can achieve a similar MSE with a slightly more constrained model (only 63 features). 
#' - If describing and interpreting the predictors is an important outcome of your analysis, this will help.
#' 
#' <!--
#' may significantly aid your endeavor.
#' -->
#' 
#' ## Model with minimum MSE
#' 
## ------------------------------------------------------------------------
plot(ames_lasso, xvar = "lambda")
abline(v=log(ames_lasso$lambda.min),col="red",lty="dashed")
abline(v=log(ames_lasso$lambda.1se),col="red",lty="dashed")

#' 
#' ## Lasso for other models than least squares
#' 
#' - Though originally defined for least squares, Lasso regularization is easily extended to a wide variety of statistical models including generalized linear models, generalized estimating equations, proportional hazards models, and M-estimators, in a straightforward fashion.
#' 
#' - Lasso’s ability to perform subset selection relies on the form of the constraint and has a variety of interpretations including in terms of geometry, Bayesian statistics, and convex analysis.
#' 
#' - The Lasso is closely related to basis pursuit denoising.
#' 
#' 
#' ## Advantages and Disadvantages
#' 
#' - Similar to Ridge, the Lasso pushes many of the collinear features towards each other rather than allowing for one to be wildly positive and the other wildly negative. 
#' - Unlike Ridge, the Lasso will actually push coefficients to zero and perform feature selection. 
#' - This simplifies and automates the process of identifying those feature most influential to predictive accuracy.
#' 
#' ## Rcode for plotting influential variables
#' 
## ----eval=F,echo=F-------------------------------------------------------
## datcoef <- coef(ames_lasso, s = "lambda.1se")
## dfcoef <- data.frame(coefficient=names(datcoef[,1]),
##                      value=as.vector(datcoef[,1]))
## dfcoef <- dfcoef[order(dfcoef[,2]),]
## ggplot(aes(value, color = value > 0,data=dfcoef)) +
##   geom_point(show.legend = FALSE) +
##   ggtitle("Influential variables") +
##   xlab("Coefficient") +
##   ylab(NULL)

#' 
#' 
## ----eval=F--------------------------------------------------------------
## coef(ames_lasso, s = "lambda.1se") %>%
##   tidy() %>%
##   filter(row != "(Intercept)") %>%
##   ggplot(aes(value, reorder(row, value), color = value > 0)) +
##   geom_point(show.legend = FALSE) +
##   ggtitle("Influential variables") +
##   xlab("Coefficient") +
##   ylab(NULL)

#' 
#' ## Plot Influential variables
#' 
#' ![](figure/influentalValues_lasso.PNG)
#' 
#' <!--
#' ## Plot Influential variables explained
#' 
#' - Often when we remove features we give up some accuracy. 
#' 
#' - To gain the refined clarity and simplicity that lasso provides, we sometimes reduce the level of accuracy. 
#' - We do not see large differences in the minimum errors.
#' - This may not be significant but if you are purely competing on minimizing error (i.e. Kaggle competitions) this may make all the difference!
#' -->
#' 
#' ## MSE for Ridge and Lasso
#' 
## ------------------------------------------------------------------------
# minimum Ridge MSE
min(ames_ridge$cvm)
# minimum Lasso MSE
min(ames_lasso$cvm)

#' 
#' ## [Elastic net](https://elitedatascience.com/algorithm-selection)
#' 
#' - Elastic-Net is a compromise between Lasso and Ridge.
#' 
#' - Elastic-Net penalizes a mix of both absolute and squared size.
#'   - The ratio of the two penalty types should be tuned.
#'   - The overall strength should also be tuned.
#' 
#' 
#' ## Elastic Nets
#' 
#' A generalization of the Ridge and Lasso models is the Elastic Net ([**Zou and Hastie, 2005**](https://rss.onlinelibrary.wiley.com/doi/full/10.1111/j.1467-9868.2005.00503.x)), which combines the two penalties.
#' 
#' $$
#' minimize \{SSE+\lambda \sum^p_{j=1} \beta^2_j+\lambda_2\sum_{j=1}^p |\beta_j|\}
#' $$
#' 
#' ## Summary overview 
#' 
#' - A result of Lasso is that typically when two strongly correlated features are pushed towards zero, one may be pushed fully to zero while the other remains in the model. 
#' - The process of one being in and one being out is not very systematic. 
#' - In contrast, the Ridge regression penalty is a little more effective in systematically reducing correlated features together. 
#' - The advantage of the Elastic Net model is that it enables effective regularization via the Ridge penalty with the feature selection characteristics of the Lasso penalty.
#' 
#' ## Implementation
#' 
#' - `alpha=.5` performs an equal combination of penalties
#' 
## ------------------------------------------------------------------------
lasso    <- glmnet(ames_train_x, ames_train_y, alpha = 1.0) 
elastic1 <- glmnet(ames_train_x, ames_train_y, alpha = 0.25) 
elastic2 <- glmnet(ames_train_x, ames_train_y, alpha = 0.75) 
ridge    <- glmnet(ames_train_x, ames_train_y, alpha = 0.0)

#' 
#' ## The four model results plottet
#' 
## ----eval=F,echo=F-------------------------------------------------------
## plot(lasso, xvar = "lambda",main="Lasso (Alpha = 1)\n\n\n")
## plot(elastic1,xvar="lambda",main="Elastic Net (Alpha = .25)")
## plot(elastic2, xvar = "lambda", main = "Elastic Net (Alpha = .75)")
## plot(ridge, xvar = "lambda", main = "Ridge (Alpha = 0)")

#' 
#' ![](figure/four_regmods.PNG)
#' 
#' ## Tuning the Elastic Net model
#' 
#' - $\lambda$  is the primary tuning parameter in Ridge and Lasso models. 
#' - With elastic nets, we want to tune the $\lambda$ and the $\alpha$ parameters. 
#' - To set up our tuning, we create a common `fold_id`, which just allows us to apply the same CV folds to each model. 
#' 
## ------------------------------------------------------------------------
# maintain the same folds across all models
fold_id <- sample(1:10, size = length(ames_train_y), 
                  replace=TRUE)

#' 
#' ## Creation of a tuning grid
#' 
#' - We then create a tuning grid that searches across a range of alphas from 0-1, and empty columns where we’ll dump our model results into.
#' 
## ------------------------------------------------------------------------
# search across a range of alphas
tuning_grid <- tibble::tibble(
  alpha      = seq(0, 1, by = .1),
  mse_min    = NA,
  mse_1se    = NA,
  lambda_min = NA,
  lambda_1se = NA
)

#' 
#' ## Iteration over $\alpha$ values - Elastic Net
#' 
#' Now we can iterate over each $\alpha$ value, apply a CV Elastic Net, and extract the minimum and one standard error MSE values and their respective $\lambda$ values.
#' 
## ------------------------------------------------------------------------
for(i in seq_along(tuning_grid$alpha)) {
  # fit CV model for each alpha value
  fit <- cv.glmnet(ames_train_x, ames_train_y, 
                   alpha = tuning_grid$alpha[i], 
                   foldid = fold_id)
  # extract MSE and lambda values
tuning_grid$mse_min[i]<-fit$cvm[fit$lambda==fit$lambda.min]
tuning_grid$mse_1se[i]<-fit$cvm[fit$lambda==fit$lambda.1se]
tuning_grid$lambda_min[i]<-fit$lambda.min
tuning_grid$lambda_1se[i]<-fit$lambda.1se
}

#' 
#' ## The resulting tuning grid
#' 
## ------------------------------------------------------------------------
tuning_grid

#' 
#' 
#' 
#' ## Plot the MSE
#' 
#' - If we plot the MSE $\pm$  one standard error for the optimal $\lambda$ value for each alpha setting, we see that they all fall within the same level of accuracy. 
#' - We could select a full Lasso model with $\lambda=0.02062776$, gain the benefits of its feature selection capability and reasonably assume no loss in accuracy.
#' 
## ----eval=F--------------------------------------------------------------
## tuning_grid %>%
##   mutate(se = mse_1se - mse_min) %>%
##   ggplot(aes(alpha, mse_min)) +
##   geom_line(size = 2) +
##   geom_ribbon(aes(ymax = mse_min + se, ymin = mse_min - se),
##               alpha = .25) +
##   ggtitle("MSE +/- one standard error")

#' 
#' 
#' ## MSE +/- one standard error
#' 
## ----echo=F--------------------------------------------------------------
tuning_grid %>%
  mutate(se = mse_1se - mse_min) %>%
  ggplot(aes(alpha, mse_min)) +
  geom_line(size = 2) +
  geom_ribbon(aes(ymax = mse_min + se, ymin = mse_min - se), alpha = .25)

#' 
#' 
#' ## Predicting
#' 
#' - With the preferred model, you can `predict` the same model on a new data set. 
#' - The only caveat is you need to supply predict an s parameter with the preferred models $\lambda$  value. 
#' - E.g., here we create a Lasso model, with a minimum MSE of 0.022. 
## ------------------------------------------------------------------------
# some best model
cv_lasso   <- cv.glmnet(ames_train_x, ames_train_y, alpha = 1.0)
min(cv_lasso$cvm)

#' 
#' - I use the minimum $\lambda$ value to predict on the unseen test set and obtain a slightly lower MSE of 0.015.
#' 
## ------------------------------------------------------------------------
# predict
pred <- predict(cv_lasso, s = cv_lasso$lambda.min, ames_test_x)
mean((ames_test_y - pred)^2)

#' 
#' 
#' ## The package `caret` - Classification and Regression Training
#' 
## ----eval=F,echo=F-------------------------------------------------------
## install.packages("caret")

#' 
#' ### [**Vignette for the `caret` package **](https://cran.r-project.org/web/packages/caret/vignettes/caret.html)
#' 
#' 
## ------------------------------------------------------------------------
library(caret)
train_control <- trainControl(method = "cv", number = 10)
caret_mod <- train(x = ames_train_x,y = ames_train_y,
                   method = "glmnet",
                   preProc = c("center", "scale", "zv", "nzv"),
                   trControl = train_control,
                   tuneLength = 10)

#' 
#' 
#' 
#' ## Output for `caret` model 
#' 
## ------------------------------------------------------------------------
caret_mod

#' 
#' <!--
#' https://www.amazon.com/Applied-Predictive-Modeling-Max-Kuhn/dp/1461468485/ref=sr_1_1?ie=UTF8&qid=1522246635&sr=8-1&keywords=applied+predictive+modelling
#' 
#' 
#' 
#' -->
#' 
#' 
#' 
#' ## Which regularization method should we choose?
#' 
#' - There’s no "best" type of penalty. It depends on the dataset and the problem. 
#' - We recommend trying different algorithms that use a range of penalty strengths as part of the tuning process
#' 
#' 
#' 
#' ## Advantages and Disadvantages
#' 
#' - The advantage of the Elastic Net model is that it enables effective regularization via the Ridge penalty with the feature selection characteristics of the Lasso penalty. 
#' - Elastic Nets allow us to control multicollinearity concerns, perform regression when $p>n$, and reduce excessive noise in our data so that we can isolate the most influential variables while balancing prediction accuracy.
#' 
#' - Elastic Nets, and regularization models in general, still assume linear relationships between the features and the target variable. 
#' - We can incorporate non-additive models with interactions, but it is tedious and difficult for a large number of features. 
#' - When non-linear relationships exist, its beneficial to start exploring non-linear regression approaches.
#' 
#' 
#' <!--
#' ## [The L1 norm explained](https://stats.stackexchange.com/questions/347257/geometrical-interpretation-of-l1-regression)
#' 
#' ![](figure/BBRXC.png)
#' -->
#' 
#' <!--
#' ## [Ridge regression](https://elitedatascience.com/algorithm-selection)
#' 
#' 
#' -    Ridge regression penalizes the squared size of coefficients.
#' -    Practically, this leads to smaller coefficients, but it doesn't force them to 0.
#' -    In other words, Ridge offers feature shrinkage.
#' -    Again, the "strength" of the penalty should be tuned.
#' -    A stronger penalty leads to coefficients pushed closer to zero.
#' -->
#' 
#' <!--
#' ## Lasso regression with package `glmnet`
#' 
## ----eval=F--------------------------------------------------------------
## install.packages("glmnet")

#' 
## ------------------------------------------------------------------------
library(glmnet)

#' 
## ------------------------------------------------------------------------
x=matrix(rnorm(100*20),100,20)
g2=sample(1:2,100,replace=TRUE)
fit2=glmnet(x,g2,family="binomial")

#' 
## ----eval=T--------------------------------------------------------------
caret::varImp(fit2,lambda=0.0007567)

#' 
#' 
#' 
#' ## 
#' 
#' - The logarithmic function is used for the link between probability and logits
#' 
#' - The Logit function is used to [linearize sigmoid curves](https://de.wikipedia.org/wiki/Logit).
#' -->
#' 
#' <!--
#' Die Logit-Funktion wird zur Linearisierung von sigmoiden Kurven verwendet.
#' -->
#' 
#' ## Further packages 
#' 
## ----eval=F--------------------------------------------------------------
## # https://cran.rstudio.com/web/packages/biglasso/biglasso.pdf
## install.packages("biglasso")

#' 
#' ![](figure/biglasso.PNG)
#' 
#' 
#' ## Resources and Links  
#' 
#' - Myers (1994) [**Classical and Modern Regression with Applications**](https://www.amazon.com/Classical-Regression-Applications-Duxbury-Classic/dp/0534380166)
#' 
#' ### Links
#' 
#' 
#' [A comprehensive beginners guide for Linear, Ridge and Lasso Regression](https://www.analyticsvidhya.com/blog/2017/06/a-comprehensive-guide-for-linear-ridge-and-lasso-regression/)
#' 
#' <!--
#' Sehr ausführlich aber mit Python
#' -->
#' 
#' - Course for statistical learning - [Youtube - Videos](https://www.r-bloggers.com/in-depth-introduction-to-machine-learning-in-15-hours-of-expert-videos/)
#' 
#' - [pcLasso: a new method for sparse regression](https://www.r-bloggers.com/pclasso-a-new-method-for-sparse-regression/)
#' 
#' - [Youtube - Lasso regression - clearly explained](https://www.youtube.com/watch?v=NGf0voTMlcs) 
#' 
#' - [`glmnet` Vignette](https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html)
#' 
#' - [Regularization Methods in R](https://www.geo.fu-berlin.de/en/v/soga/Geodata-analysis/multiple-regression/Regularization-Methods/Regularization-Methods-in-R/index.html)
#' 
#' - [A gentle introduction to logistic regression and Lasso regularisation using R](https://eight2late.wordpress.com/2017/07/11/a-gentle-introduction-to-logistic-regression-and-lasso-regularisation-using-r/)
#' 
#' - [Penalized Regression in R](https://machinelearningmastery.com/penalized-regression-in-r/)
#' 
#' <!--
#' - [Penalized Logistic Regression Essentials in R](http://www.sthda.com/english/articles/36-classification-methods-essentials/149-penalized-logistic-regression-essentials-in-r-ridge-lasso-and-elastic-net/)
#' 
#' - [All you need to know about Regularization](https://towardsdatascience.com/all-you-need-to-know-about-regularization-b04fc4300369)
#' --
#' <!--
#' Evtl könnte ich hier noch eine Aufgabe draus machen:
#' https://www.r-bloggers.com/ridge-regression-and-the-lasso/
#' -->
#' 
#' <!--
#' http://www.sthda.com/english/articles/36-classification-methods-essentials/149-penalized-logistic-regression-essentials-in-r-ridge-lasso-and-elastic-net/
#' 
#' https://eight2late.wordpress.com/2017/07/11/a-gentle-introduction-to-logistic-regression-and-lasso-regularisation-using-r/
#' 
#' https://www.r-bloggers.com/ridge-regression-and-the-lasso/
#' 
#' http://www.science.smith.edu/~jcrouser/SDS293/labs/lab10-r.html
#' 
#' https://gigadom.in/2017/10/20/practical-machine-learning-with-r-and-python-part-3/
#' 
#' https://beta.vu.nl/nl/Images/werkstuk-fonti_tcm235-836234.pdf
#' -->