        table width = 900 border = 0><tr><td>
# <b>NMST539 | Lab Session 10</b>

## <b>Mixture models, correspondance analysis and variable selection techniques</b>
### <b><font color = "#880000">Summer Term 2019/2020 | May 11th - May 15th, 2020 (week 13)</font></b>
###### <a href = "individual_week13.Rmd">Rmd file (UTF8 coding)</a>

**************


The R-software is available for download from the official website: <a href = "https://www.r-project.org">https://www.r-project.org</a>

A user-friendly interface (one of many): <a href = "https://www.rstudio.com">RStudio</a><br><br>

<b>Some useful manuals and tutorials for <i>R</i>
(in Czech or English):</b><br>
    
<ul>
<li>Bína, V., Komárek, A. a Komárková, L.: <i>Jak na jazyk R.</i> (in Czech) | (<a href = "NMFM301/Rmanual2.pdf">PDF manual</a>)</li>
<li>Komárek, A.: <i>Základy práce s R.</i> (in Czech) | (<a href = "http://msekce.karlin.mff.cuni.cz/~komarek/vyuka/dataRko/Rmanual1.pdf">PDF manual</a>)</li>
<li>Kulich, M.: <i>Velmi stručný úvod do R.</i> (in Czech) | (<a href = "NMFM301/uvodr.pdf" >PDF manual</a>)</li>
<li>Venables, W.N., Smith, D.M., and the R Core Team: <i>An Introduction to R.</i> |  (<a href = "https://cran.r-project.org/doc/manuals/R-intro.pdf">PDF manual</a>)</li>
<li>De Vries, A. a Meys, J.: <i>R for Dummies.</i> (ISBN-13: 978-1119055808) |   (<a href = "http://www.ievbras.ru/ecostat/Kiril/R/Biblio/R_eng/R%20dummies.pdf">PDF book</a>)</li>    
</ul>
**************

  

### <font color = "#880000"><b>1. Mixture models</b></font>

We have already discussed some clustering techniques where the available observations are assigned to some groups with respec to their mutual similarities (depending on the 
similarity measure or the distance which is used to assess the differences among objetcs). 

However, a fundamental model for clustering -- a mixture model -- is more complex than that. In general, it is standardly assumed that the  observations are randomly drawn from K different populations but instead of providing just the corresponding labels for the observations (labeling observations accorging to their assignment to the groups) one estimates the whole population with the mixture model. Thus, a much more complex information is obtained as a results. On the other hand, more assumptions are needed as well. <br><br>

Let us consider the data on the consumption of automobiles in the United states in 80's (the dataset 'mtcars' available in R). The data consists of 32 observations (32 different cars on the US market available that time) and 11 covariates (some of them are continious, some of them are categorical). 
For an additional information on this dataset use the help session in R (type `?mtcars`). 

```{r, warning = F, message = F}
rm(list = ls())
attach(mtcars)
head(mtcars)
```
<br><br>

Using the given data, we can apriori assume three different populations -- for instance, the populations determined by the number of cylinders (four, six, and eight cylinder cars).
Each population has its own specific characteristics -- such as the mean, variance, distribution, or density. Considering, for instance, the car consumption (miles per gallon) 
we can estimate the proportion of each population in the given sample and, imposing the normality condition with the unit variance, we can estimate the overall mixture model by estimating the corresponding mean values and plugging everything all together.  The estimated mixture distribution is the following: 

```{r, fig.width=10}
### population mean estimates
mean1 <- mean(mpg[cyl == 4])
mean2 <- mean(mpg[cyl == 6])
mean3 <- mean(mpg[cyl == 8])

### estimated proportions for each population
prop <- table(cyl)/dim(mtcars)[1]

### the overall mixture distribution
mixture <- function(x, prop, means){
  mixValue <- 0
  for (i in 1:length(prop)){
    mixValue <- mixValue + prop[i] * dnorm(x, means[i])
  }
  return(mixValue)
}

xSeq <- seq(10, 34, length = 100)
density <- mixture(xSeq, prop, c(mean1, mean2, mean3))
plot(density ~ xSeq, col = "red", type = "l", xlab = "Miles per Gallon", ylab = "Density")
```

The mixture is composed from three mixtures (normal distributions with the unit variance and some corresponding mean parameters). 

```{r, fig.width = 10}
plot(density ~ xSeq, col = "red", type = "l", lwd = 3, xlab = "Miles per Gallon", ylab = "Density")
d1 <- density(rnorm(10000, mean1, 1))
d2 <- density(rnorm(10000, mean2, 1))
d3 <- density(rnorm(10000, mean3, 1))
lines(d1$x, prop[1] * d1$y, col = "black", lwd = 2)
lines(d2$x, prop[2] * d2$y, col = "violet", lwd = 2)
lines(d3$x, prop[3] * d3$y, col = "gray", lwd = 2)
```

<br><br>

The joint density (mixture distribution) has multiple nodes -- one for each cluster mean $\mu_k$, for $k = 1, 2, 3$. However, for two centroids which are close enough, 
it is possible that only one node will be present in the final mixture. <br><br>

In practical problems the mixture estimation becomes way more complex. Usually, it is common that the number of populations included in the mixture is not known apriori 
or the variance can differ across different populations. <br><br>

Maximum likelihood methods are typically used to estimate the overall mixture models but some regularity conditions (restrictions) are required to obtain some reasonable solutions 
at the end. 

<br><br>

#### <font color = "#880000">Do by yourself</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
 <ul>
 <li>What will happen, if you assume the same approach as above (a normal mixture model with unit variance), however, with no restriction regarding the number of the underlying populations? The minimization problem can be formulated as 
$$
Maximize \prod_{k = 1}^{K} \prod_{i = 1}^{n_k} f_k(X_{k i}) \equiv Maximize \prod_{k = 1}^{K} \prod_{i = 1}^{n_k} \frac{1}{\sqrt{2 \pi}} \cdot exp \left\{\frac{-1}{2} (X_{k i} - \mu_k)^2\right\}
$$
 where the maximization takes place with respect to $K \in \{1, \dots, n\}$, and $\mu_1, \dots, \mu_K \in \mathbb{R}$ (with $n \in \mathbb{N}$ being the overall number of observations 
 and $n_k \in \mathbb{N}$ being the number of subjects from the $k$-th population -- observations $X_{k 1}, \dots, X_{k n_k}$). If no further restrictions are defined, the final mixture model 
 will look like this: 
 
```{r, fig.width = 10}
plot(density ~ xSeq, col = "gray", type = "l", lwd = 2, xlab = "Miles per Gallon", ylab = "Density", lty = 2, ylim = c(0,0.5))
for (i in 1:length(mpg)){
  lines(rep(mpg[i], 2), c(dnorm(0), 0), col = "red", lwd = 1, lty = 3)
  points(mpg[i], dnorm(0), pch = 21, cex = 0.5, bg = "red")
}
```
 
 
 </li>
 <li>Can you explain what you see in the figure above? What will happen if you run the maximization problem with respect to the unknown variance parameters $\sigma_k^2 > 0$, for $k = 1, \dots, K$? </li>
 
 <li>How would you utilize the overall mixture model for clustering? How would you define the decision algorithm if you need to assign some new observation into one of three clusters?</li>
 
 <li>Try to improve the mixture distribution, for instance, by considering heterogenous variances. Compare the results under homoscedasticity and heteroscedasticity. Define the overall likelihood and try to maximize it to obtain the final mixture model. </li>
 <li>The mixture can be further improved by considering some other than normal distribution. On the other hand, the normal distribution is quite convenient from some theoretical point of view. Some prametric distribution family is usually assumed as it allows to express the overall mixture structure in terms of the likelihood. If, in addition, the number of underlying populations is given as well, it improves the calculations significantly. 
 </li>
 
 <li>One needs to be carefull with the mixture estimation. Usually some regularization techniques are needed to avoid a so-called "Dirac catastrophe". The regularization approaches can be, however, implemented differently (for instance, instead of assuming the number of pupulations, one can assume some pre-defined shape for the final mixture -- such as unimodality, or log-concavity, etc.). </li>
 </ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

Mixture models are also used in a regression modelling. As an illustration, consider a GLM model (with a logarithmic link) for Poisson counts (for instance, modelling the number of positive COVID tests in some subpopulation units). Due to some low prevalence, it can be assumed, that very likely no positive cases will be recorded while, occasionally, some occurrences will be observed. The overall model can be, for intance, expressed as 

$$
Y_i | \boldsymbol{x}_i \sim \left\{ 
\begin{array} 
 00 & \textrm{with probability } p \in (0,1);\\
Poiss(\lambda_i) & \textrm{with probability } 1 - p_i.
\end{array}   \right.
$$
Here, the subpopulation characteristics (independent covariates) are given as $\boldsymbol{x}_i$. 
In other words, the probability of observing no cases in a given subpopulation is $P[Y_i = 0 | \boldsymbol{x}_i] = p + (1 - p) \cdot e^{- \boldsymbol{x}_i^\top \boldsymbol{\beta}}$ and the probability of observing $k > 0$ cases can be expressed as $P[Y_i = k | \boldsymbol{x}_i] = (1 - p) \cdot \frac{(\boldsymbol{x}_i^\top \boldsymbol{\beta})^k}{k!}e^{- \boldsymbol{x}_i^\top \boldsymbol{\beta}}$ where we assume a standard Poisson (GLM) regression to model the mean parameter $\lambda_i = E[Y_i | \boldsymbol{x}_i] = \boldsymbol{x}_i^\top \boldsymbol{\beta}$. for some unknown vector of parameters $\boldsymbol{\beta}$. The final model is obtained by maximizing the corresponding likelihood. <br>
Such mixture models are very common in practical applications. The model above is ussually referred to as a zero-inflated Poisson regression model. <br><br>



### <font color = "#880000"><b>2. Canonical Correlation Analysis in R</b></font>

The canonical correlation analysis (CCA) is a multidimensional exploratory statistical method
which (roughly) operates on the same principle as the principal component analysis (PCA). However, unlike the PCA where one tries to explore the common variability within the given data set, the canonical correlation is used to explore the mutual 
variability (correlations) between two different sets of quantitative variables observed on the same experimental units. Thus, the PCA method deals with one dataset only and it tries to reduce the overall dimensionality of the dataset using some linear combinations of the initial variables with the maximum variace possible
while the CCA approach searhes for some linear combinations of the initial variables in both datasets in order to achieve maximum correlations amonng them instead. <br><br>

For instance, for two multivariate random vectors $\boldsymbol{X}$ and $\boldsymbol{Y}$ we assume that the common variance-covariance matrix of $(\boldsymbol{X}, \boldsymbol{Y})^\top$ is of the form 
$$
Var \left(\begin{array}{c}\boldsymbol{X}\\ \boldsymbol{Y}\end{array}\right) = \left(\begin{array}{cc} \Sigma_{XX} & \Sigma_{XY}\\ \Sigma_{YX} & \Sigma_{YY} \end{array}\right).
$$

If $\boldsymbol{X}$ is a $p$ dimensional random vector and $\boldsymbol{Y}$ is a $q$ dimensional random vecton, then the corresponding correlation of any linear combination of $\boldsymbol{X}$ and $\boldsymbol{Y}$ in terms
of $\boldsymbol{a}^\top \boldsymbol{X}$ and $\boldsymbol{b}^\top\boldsymbol{Y}$, for $\boldsymbol{a} \in \mathbb{R}^p$ and $\boldsymbol{b} \in \mathbb{R}^q$ can be expressed as 
$$
Cor(\boldsymbol{a}^\top \boldsymbol{X}, \boldsymbol{b}^\top\boldsymbol{Y}) = \boldsymbol{a}^\top \Sigma_{XY} \boldsymbol{b}.
$$

Therefore, the canonical correlations analysis searches for such vectors $\boldsymbol{a} \in \mathbb{R}^p$ and $\boldsymbol{b} \in \mathbb{R}^q$, that the previous correlation is as large as possible. However, 
the canonical correlation analysis does not search for just one pair of the canonical variables $\boldsymbol{a}^\top \boldsymbol{X}$ and $\boldsymbol{b}^\top\boldsymbol{Y}$ but, instead, it defines all pairs (up to $min(p, q)$), 
such that all canonical covariates within each dataset are uncorrelated among each other. <br><br>

In the R statistical software there is function `cancor()` available under the standard instalation. 
Some additional packages can be downloaded and installed as well (mainly the R packages 'CCA' and 'vegan'). 

In the following we will consider the same dataset we have already worked with before. The data consists 
of two datasets where each dataset represents measurements in some specific (64-65) river localities in the Czech republic. The first dataset contains measurements of biologial metrics for each locality (17 different metrics and taxons) and the secon dataset contains measurements on chemical concentrations and values at the same localities (7 covariates recored). 

The idea is to somehow relate both dataset in order to explain what bio metrics can correlate which which chemical concentrations. One way of answering this is to apply the canonical correlation approach as explained above. 


```{r, warning = F, message = F}
rm(list = ls())
bioData <- read.csv("http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/bioData.csv", header = T)
chemData <- read.csv("http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/chemData.csv", header = T)

head(bioData)
head(chemData)
```
<br><br>

Again, there is one locality which is missing the dataset of chemical measurements. We identify this locality and we do not consider it for further analysis. Both datasets are alligned with respect to he locality names - the first column in each dataset. 


```{r, warning = F, message = F}
ind <- match(chemData[,1], bioData[,1])
data <- data.frame(bioData[ind, ], chemData[, 2:8])
```


```{r}
X <- data[,2:9]
Y <- data[,19:25]
```

A nice graphical tool to view both datasets with respect to their overall (within and between) correlation structure is avialable within the R package 'CCA' (Canonical Correlation Analysis). The packages needs to be firstly installed (`install.packages("CCA")`) and later one can use functions `matcor()` and `img.matcor()` to graphically display the corration matrix 

$$ cor(\mathcal{X}, \mathcal{Y}) = \left(\begin{align}\Sigma_{XX} & \Sigma_{XY}\\\Sigma_{YX} & \Sigma_{YY}\end{align}\right)$$.

```{r, warning=F, message = F, fig.width=10, fig.height=8}
library("CCA")
correl <- matcor(X, Y )
img.matcor(correl, type = 2)
```


Next we can directly use `cancor()` function to obtain canonical correlations for the datasets $\mathcal{X}$ and $\mathcal{Y}$. Another option is to use function `cc()` (from the R package 'CCA'). Note, the the results are equivalent although not identical. 

```{r}
cc1 <- cancor(X, Y)  ### function from standard R instalation
cc2 <- cc(X, Y)      ### function for the R package 'CCA'
```

The canonical correlations are the same for both approaches (R commands used) 

```{r}
cc1$cor  ### function from standard R instalation
cc2$cor  ### function for the R package 'CCA'
```


The canonical covariates can be easily reconstructred as follows: 
```{r}
CanX <- as.matrix(X) %*% cc1$xcoef
CanY <- as.matrix(Y) %*% cc1$ycoef

cor(CanX, CanY)
```

Which corresponds with the output provided by the function `cancor()` and `cc()`above (see also the picture below
for a visual comparison).<br><br>



```{r, warning=F, message = F, fig.width=10, fig.height=8}
par(mfrow = c(1,2))
barplot(cc1$cor, main = "Canonical correlations for 'cancor()'", col = "gray")
barplot(cc2$cor, main = "Canonical correlations for 'cc()'", col = "gray")
```


However, the actual estimated linear combinations of the original covariates in both datasets $\mathcal{X}$ and $\mathcal{Y}$ are different. Compare the following: 

```{r}
cc1$xcoef  ### function from standard R instalation
cc2$xcoef  ### function for the R package 'CCA'

cc1$ycoef
cc2$ycoef
```

Using the 'CCA' package one can also take an advantage of nice graphical tools which are available for
the canonical correlation analysis. The idea is to display maximized correlations between transformed variables of the dataset $\mathcal{X}$ and the dataset $\mathcal{Y}$. 

```{r, warning=F, message = F, fig.width=10, fig.height=8}
plt.cc(cc2, var.label = TRUE, ind.names = data[,1])
```


#### <font color = "#880000">Questions</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
 <ul>
 <li>How would you interpret the left hand side of the figure?</li>
 
 <li>What is the relationship between the subfigure on the left hand side and the subfigure on the 
 right hand side? </li>
 
 <li>What would be your conclussion about the two available dataset and their mutual relationship? </li>
 </ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

A usefull extension to the classical canonical correlation approach is the regularized version - so called regularized canonical correlation analysis. It is meant for scenarios where the number of parameters $p \in \mathbb{N}$ is greater than the sample size $n \in \mathbb{n}$ (usually $p \gg n$). 
There is R function `rcc()` available in the R package 'CCA'. <br><br>

Another option how to perform canonical correlation in R is to install some additional R packages -- for instance, package 'vegan' (use `install.packages("vegan")` for instalation) and to apply function `cca()`.  Use the help session in R for more details about the function itself or the corresponding package respectively.<br><br>




#### <font color = "#880000">To Do</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
 <ul>
 <li>Use this markdown and consider (maybe) the whole dataset on biologial measurements. Compare how 
 much the results changes. </li>
 
 <li>Consider a lower dimensional profile and try to refit the whole canonical correlation analysis on
 a lower dimensional datasets (e.g. two covariates for bio data and two covariates for chemical data). </li>
 
 <li>Try to manually reconstruct the results from the above by using the correponding canonical correlation theory form the lecture. </li>
 
 
 </ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

#### <font color = "#880000">Note</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
 <ul>
 <li><b>Canonical Correspondance Analysis</b><br> 
 An analogous method that the canonical correlation analysis, however, mainly meant for discrete data.<br><br> </li>
 
 <li><b>Redundancy Analysis</b><br></li>
 A similar method than the canonical correlation analysis, however, not symetric. The idea is tu use linear combinations of one dataset 
 to maximize the corelation with the original covariates in the second dataset. Thus, it depends on which data set is used as the orginal one. 
 
 </ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>



### <font color = "#880000"><b>3. Variable selection via LASSO and other techniques</b></font>

Finally, we shortly return back to the principal component analysis and we briefly discuss some more popular alternative. The PCA approach 
is mostly meant as an exploratory tool. But, what can it be used for after the data exploration is performed? One possible direction is to 
use the recovered principal components to explain the covariate of interest -- a regression on the principal components. This is a classical 
regression techniques, however, instead of using the underlying covariates, one includes the principal components into the model instead. 

Sometimes, it is also convenient to use both -- the principal components for some dimensionality reduction and the original covariates for 
some better interpretation purposes. Indeed, regressing the dependent covariate with respect to the principal components introduces a huge challenge 
in terms of the following interpretation. Therefore, some alternative techniques for the variable selection and data dimensionality reduction 
are recently considered in statistics. 

A significant dimensional reduction and, at the same time a variable selection technique, is introdused by a reqularized 
regression based on aLASSO penalization. LASSO stands for the Least Absolute Selection and Shrinkage Operator
and it was introduced by Tibshirani (1996). The main idea of the LASSO reqularized (penalized) regression
is the following one: instead of minimizing 

$$
\|\boldsymbol{Y} - \mathbb{X}\boldsymbol{\beta}\|_{2}^{2},
$$

with respect to the vector of regression coefficients $\boldsymbol{\beta} \in \mathbb{R}^{p}$ for $p < n$ once considers an alternative minimization defined as 

$$
\|\boldsymbol{Y} - \mathbb{X}\boldsymbol{\beta}\|_{2}^{2} + \lambda \|\boldsymbol{\beta}\|_{1},
$$

where $\lambda > 0$ is some reqularization parameter and $\|\cdot\|_{1}$ stands for a classical $L_{1}$ norm. It is obvious that the larger the value of the regularization parameter $\lambda > 0$ the more dominant the second term in the minimization problem above. Thus, the elements of $\boldsymbol{\beta}$ are shrunk towards zero. Actually, for large values of $\lambda \to \infty$ the most elements of the vector of the regression parameters will be set to zero directly. Such solution is called a sparse solution (the vector of parameter estimates is sparse, it contains non-zero elements very rarely).

<img src = "contours.png" width = "900">

On the other hand, for $\lambda \to 0$ the first term in the minimization problem is dominant and thus, for $\lambda = 0$ the whole problem reduces to a classical least squares regression. 

The main advantage of the LASSO regularized regression approach is that it also allows for heavily overparametrized models (situations where $p \gg n$). Such problems commonly occur in genetic data, sociology surveys, econometrics and basically all other areas... 

Unlike the classical least squares problem where the vector of the parameter estimates is given as a solution of the normal equations and thus, is explicit, no explicit solution is generaly given for the LASSO regression. But... <br><br>


#### <font color = "#880000"><b>3.1. LARS-LASSO Algorithm</b></font>

A very elegant solution to the LASSO regression fitting problem was obtained once we realized that
the solution paths are piece-wise linear. Given this, one only needs to calculate the knot point positions and the rest of the path is a linear approximation only. 

Using the LARS-LASSO algorithm proposed by Efron et al. (2004) we can do that in a very effective way. Actually in $p$ linear steps only!!! 

<img src = "lasso_lars.png" width = "900">

The whole algorithm is based on a geometric interpretation of covariances between the current response estimate and the covariates still not included in the model yet. 

```{r}
source("http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/larsLASSO.r")
```

In the sourse code loaded above there are two functions for R: `solutionPath()` and `solutionPathPlot()`. 

<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li>`solutionPath(X, Y, noJumps)` - calculates the complete solution paths for the first 'noJumps' parameters entering the model. The value of 'noJumps' is an integer value between one and the second dimension of $X$ (the total number of predictors). The matrix $X$ stands for a classical regression matrix of the type $n \times p$ however, even situations for $p \gg n$ are allowed. Finally, the vector of $\boldsymbol{Y}$ is a response vector. As a result one obtains a list with two elements: 'betaVector' with augmented parameter estimates in each row and 'corrMatrix' which us a matrix with actual correlations of the actuall reponse estimate and the independent covariates.<br></li>
<li>`solutionPathPlot(solutionVector, corrMatrix, xSeq)` - a function which creates two plots: a complete solution path until 'noJump' number of nonzero predictors enter the model and the plot with actual correlation decrease at each step. </li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br>

```{r, warning = F, message = F}
Y <- data$Pcelk
X <- as.matrix(data[,5:18])
LARSsolution<- solutionPath(X, Y, 10)

betaVector <- LARSsolution[[1]]
corrMatrix <- LARSsolution[[2]]
```


And the graphical output can be obtained using this code: 

```{r fig.width=10, fig.height=10}
par(mfrow = c(2,1))
solutionPathPlot(betaVector, corrMatrix, Y, data[,1])
```


#### <font color = "#880000"><b>Cons and Prons</b></font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li>The main advantage of the LASSO regularized regression approach is that it allows for very high dimensions to be effectively taken care of. Indeed, with the LASSO reqularized regression there is no problem to solve regression problems with $p \gg n$. However, the number of nonzero parameters at the end (in the final model) is $n$ at most. Irrelevant regressors are set to zero.  </li>
<li>Instead of one set of estimates one can (with the same time efficiency) obtain the whole path of solutions. This makes the LASSO selection and estimation a very effective tool for high dimensional problems. <br></li>
<li>Probably the main disadvantage is that the final estimates are shrunked towards zero. There were different alternatives proposed to correct for this bias. The most popular approach is to use LASSO for the variable selection only and afterwards one should perform a classical linear regression on the selected covariates only. </li>
<li>Another disadvantage is a lack of inference for these models (see e.g. post-selection inference in the literature).</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>



#### <font color = "#880000"><b>3.2. LASSO Methods in R</b></font>
There are few usefull libraries in R which implement the LASSO shrinkage method and the LARS algorithm too. In addition, there were various modifications of the LASSO estimation approach proposed in the literature. Eventhough, the LASSO is a very effective variable selection tool it still suffers of some issues and disadvantages. Most of the proposed modifications aim to improve the LASSO performance especially with respect to the selection performance and consistency (for instance, to achieve the ORACLE properties or sign consistency instead). 

The most convenient package, which offers, beside the LASSO regression modeling, also the Ridge regression and the elastic net algorithms 
is implemented in the library `glmnet` (use the command `install.packages("glmnet")` for installation and `library("glmnent")` for initialization in the R environment). 


```{r, fig.width=10}
library("glmnet")

fit1=glmnet(X,Y)
plot(fit1)

coef(fit1,s=0.01)
```


Another library, which offers LARS fitting appraoch for various models (for instance, LASSO, Least Angle Regresssion, Forward Stagewise, or stepwise regression) is implemented in the package `lars` (run the command `install.packages("lars")` to install the library and `library("lars")` to initialize it in the R environment).

See the next examples for a brief illustration:

```{r, fig.width=10, message=F, warning=F}
library("lars")

object1 <- lars(X,Y, type = "lasso")
plot(object1)

object2 <- lars(X,Y, type = "lar")
plot(object2)

object3 <- lars(X,Y, type = "forward.stagewise")
plot(object3)

object4 <- lars(X,Y, type = "stepwise")
plot(object4)
```





<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<hr size="3" width="100%" noshade style="color:#440000" align="left" /><br><br> 



## <font color = "#880000"><b>Homework Assignment</b></font>
### (Deadline:  End of the summer term | 22.05.2020)

Use the command below to instal the `SMSdata` packages with various multivariate datasets: 

```{r, eval = F}
install.packages(pkgs="http://www.karlin.mff.cuni.cz/~hlavka/sms2/SMSdata_1.0.tar.gz", repos=NULL, type="source")
```

Chose one dataset of your choise from the list of all available datasets in the package:

```{r, eval = F}
data(package = "SMSdata") 
```
 
 There are 21 different datasets and you can load each of them by typing and running 
`data(dataset_name)`, for instance:

```{r, eval = F}
data(plasma) 
```

<ul>
<li>Define a reasonable mixture model and estimate the model by using the maximum likelihood approach; 
(you are free to use some apriori defined number of populations but try to use some reasonable mixture model -- for instance, 
use some prior knowledge or some exploratory analysis to determin the underlying number of populations)</li>

<li>Provide some grafical output for the estimated mixture model and discuss the results (interpret the model);</li>
</ul>

<font color = "red">The final PDF document with all partial tasks assigned during the term should be finalized and submitted 
by email to <b>maciak [AT] karlin.mff.cuni.cz</b> by Friday, May 22nd, 2020. </font> This only applies to students who intend 
to obtain the course credit. 



</td></tr></table>
<br><br><br>