---
title: "Quick Start"
author: "Martin Vincent"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Quick Start}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, echo = FALSE}
set.seed(150)
```

## Quick Start (for msgl version `r packageVersion("msgl")`)

### 1. Load the msgl library in R
```{r results="hide"}
library(msgl)
```

### 2. Load your data
Load data containing N samples and p features (covariates):

```{r eval = FALSE}
x <- # load design matrix (of size N x p)
classes <- # load class labels (a vector of size N)
```

For the purpose of this tutorial we will load a data set consisting of microRNA normalized expression measurements of primary cancer samples.
```{r}
data(PrimaryCancers)
x[1:5,1:5]
dim(x)
table(classes)
```
Hence, p = 384, N = 165 and the number of classes K = 9, this implies that the multinomial classification model has 9*(384+1) = 3465 parameters.

Let us take out a small test set:
```{r}
idx <- 1:10
x.test <- x[idx,]
x <- x[-idx,]
classes.test <- classes[idx]
classes <- classes[-idx]
```

### 3. Estimate error using cross validation
Choose `lambda` (fraction of lambda.max) and `alpha`, with `alpha = 1` for lasso, `alpha = 0` for group lasso and `alpha` in the range (0,1) for spares group lasso.

Use `msgl::cv` to estimate the error for each lambda in a sequence decreasing from the data derived *lambda max* to `lambda` * *lambda max*.
Lambda max is the lambda at which the first penalized parameter becomes non-zero.
A smaller `lambda` will take longer to fit and include more features.
The following command will run a 10 fold cross validation for each lambda value in the lambda sequence using 2 parallel units (using the [foreach](https://CRAN.R-project.org/package=foreach) and [doParallel](https://CRAN.R-project.org/package=doParallel) packages.

```{r}
cl <- makeCluster(2)
registerDoParallel(cl)

fit.cv <- msgl::cv(x, classes, fold = 10, alpha = 0.5, lambda = 0.1, use_parallel = TRUE)

stopCluster(cl)
```

(for the current version *no progress bar will be shown*)

**Get a summery of the validated models.**
We have now cross validated the models corresponding to the lambda values, one model for each lambda value. We may get a summery of this validation by doing:
```{r}
fit.cv
```

Hence, the best model is obtained using lambda index `r best_model(fit.cv)` and it has a cross validation error of `r round(Err(fit.cv)[best_model(fit.cv)],2)`. The expected number of selected features is `r colMeans(features_stat(fit.cv))[best_model(fit.cv)]` and the expected number of parameters is `r colMeans(parameters_stat(fit.cv))[best_model(fit.cv)]`.

### 4. Fit the final model

Use msgl to fit a final model.
```{r}
fit <- msgl::fit(x, classes, alpha = 0.5, lambda = 0.1)
```
**Get a summery of the estimated models**
```{r}
fit
```

**Take a look at the estimated models.**
As we saw in the previous step the model with index `r best_model(fit.cv)` had the best cross validation error, we may take a look at the included features using the command:
```{r}
features(fit)[[best_model(fit.cv)]] # Non-zero features in best model
```

Hence `r length(features(fit)[[best_model(fit.cv)]])` features are included in the model, this is close to the expected number based on the cross validation estimate.

The sparsity structure of the parameters belonging to these `r length(features(fit)[[best_model(fit.cv)]])` features may be viewed using  
```{r}
parameters(fit)[[best_model(fit.cv)]]
```
We may also take a look at the estimate parameters (or coefficients)

```{r}
coef(fit, best_model(fit.cv))[,1:5] # First 5 non-zero parameters of best model
```

If we count the total number of non-zero parameters in the model we get, in this case `r sum(parameters(fit)[[best_model(fit.cv)]])` which is close to the expected based on the cross validation estimate.

### 6. Use your model for predictions

**Load test data** containing M samples and p features.
```R
x.test <- # load matrix with test data (of size M x p)
```

Use the final model to predict the classes of the M samples in x.test.
```{r}
res <- predict(fit, x.test)

res$classes[,best_model(fit.cv)] # Classes predicted by best model
classes.test # True classes
```

We may also get the estimated probabilities for each of the classes
```{r}
res$response[[best_model(fit.cv)]]
```
