---
title: "Model objects"
author: 
  name: Dr Gavin Rhys Lloyd
  affiliation: Phenome Centre Birmingham, University of Birmingham, UK
output:
  BiocStyle::html_document:
    toc_float: true
package: structToolbox
abstract: Introduction to model objects
vignette: >
  %\VignetteIndexEntry{Model objects}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
    dpi=72
)
library(structToolbox)
library(gridExtra)
```

</br></br>

# Using `struct` model objects
PCA (Principal Component Analysis) is a commonly applied method for exploring multivariate datasets. We will use the iris DatasetExperiment as an example, which is included in the package and already prepared as a DatasetExperiment object.

```{r}
D = iris_DatasetExperiment()
head(D$data)
```
 
</br></br>

## PCA model
Before we apply PCA we first need to create a PCA object. This object contains all the inputs, outputs and methods needed to apply PCA. We can set parameters such as the number of components when the PCA model is created, but we can also use dollar notation to change/view it later. 

```{r}
P = PCA(number_components=15)
P$number_components=5
P$number_components
```
  
The inputs for a model can be listed using `param_ids(object)`:

```{r}
param_ids(P)
```
</br></br>

## Model sequences
Unless you have very good reason not to, it is usally sensible to mean centre the columns of the data before PCA. Using the `STRUCT` framework we can create a model sequence that will mean centre and then apply PCA to the mean centred data.

```{r}
M = mean_centre() + PCA(number_components = 4)
```
  
In `STRUCT` mean centring and PCA are both model objects, and therefore joining them creates a model_sequence object. The objects in the sequence can be accessed by indexing, and we can combine this with dollar notation. For example, the PCA object is the second object in our sequence and we can access the number of components like this:

```{r}
M[2]$number_components
```
</br></br>

## Training/testing models
Model and model_sequence objects need to be trained using a training DatasetExperiment.

```{r}
M = model_train(M,D)
```
  
Model objects can be used to generate predictions for test datasets. For this example we will just use the training data (sometimes called autoprediction).

```{r}
M = model_predict(M,D)
```

The available outputs for an object can be listed and accessed using dollar notation:
  
```{r}
output_ids(M[2])
M[2]$scores
```
</br></br>

## Model charts
The struct framework includes charts. Charts associated with a model object can be listed.

```{r}
chart_names(M[2])
```
  
Like model objects, chart objects need to be created before they can be used. Here we will plot the PCA scores plot for our mean centred PCA model_

```{r}
C = pca_scores_plot(factor_name='Species') # colour by Species
chart_plot(C,M[2])
```
  
If we makes changes to our chart object, we must call `chart_plot` again.

```{r}
# add petal width to emta data of pca scores
M[2]$scores$sample_meta$Petal.Width=D$data$Petal.Width
# update plot
C$factor_name='Petal.Width'
chart_plot(C,M[2])
```
  
The `chart_plot` method can return e.g. a ggplot object so that you can easily combine it with other plots using the gridExtra package for example.

```{r,fig.width=10}
C1 = pca_scores_plot(factor_name='Species') # colour by Species
g1 = chart_plot(C1,M[2])
C2 = pca_scree_plot()
g2 = chart_plot(C2,M[2])
grid.arrange(grobs=list(g1,g2),nrow=1)
```
</br></br>

## STATO Integration
Some model objects are also STATO objects. STATO is a general purpose statistics ontology (http://stato-ontology.org/). In the `STRUCT` framework we use it to provide standarded definitions for objects. The PCA model object is also a STATO object.

```{r}
is(PCA(),'stato')
```

We can access the STATO ontology using some methods specific to stato objects.

```{r}
# this is the stato id for PCA
stato_id(P)

# this is the stato name
stato_name(P)

# this is the stato definition
stato_definition(P)
```

This information is more succinctly displayed using `stato_summary`. This method also scans over all inputs and outputs for those with STATO definitions and displays those as well. For PCA the number of components is present, but none of the outputs are STATO objects and therefore no definition is provided.

```{r}
stato_summary(P)
```

# Validating models
Validation is an important aspect of chemometric modelling. The `STRUCT` framework enables this kind of iterative model testing through `iterator` objects. In order to demonstrate this we will first load the iris data set, which as been pre-prepared as a DatasetExperiment object as part of the `STRUCT` package.

```{r}
D = iris_DatasetExperiment()
D
```

</br></br>

## Cross-validation
Cross validation is a common technique for assessing the performance of classification models. For this example we will use a PLSDA model_ Data should be mean centred prior to PLS, so we will build a model sequence first.

```{r}
M = mean_centre() + PLSDA(number_components=2,factor_name='Species')
M
```

Iterators objects like the k-fold cross-validation object can be created just like any other struct object. Parameters can be set at creation =, and accessed/changed later using dollar notation.

```{r}
XCV = kfold_xval(folds=5,factor_name='Species')
# change the number of folds
XCV$folds=10
XCV$folds
```

The model to be cross-validated can be set/accessed used the `models` method.

```{r}
models(XCV)=M
models(XCV)
```

Alternatively,  iterators can be combined with models using the multiplication symbol:
```{r}
XCV = kfold_xval(folds=5,method='venetian',factor_name='Species') * 
      (mean_centre()+PLSDA(number_components = 2,factor_name='Species'))
```

The `run` method can be used with any iterator object. The iterator will then run the model sequence multiple times. In our case we will run cross-validation 5 times splitting the data into different training and test sets each time. The `run` method also needs a `metric` to be specified. This metric may be calculated once after all iterations, or after each iteration, depending on the iterator type (resampling, permutation etc). For cross-validation we will calculate balanced accuracy after all iterations.

```{r}
XCV = run(XCV,D,balanced_accuracy())
XCV$metric
```

</br></br>

Like other `STRUCT` objects, iterators can have chart objects associated with them. The `chart_names` function will list them for an object.

```{r}
chart_names(XCV)
```

Charts for iterator objects can be plotted in the same way as charts for any other object.

```{r,warning=FALSE}
C = kfoldxcv_grid(factor_name='Species')
chart_plot(C,XCV)[[2]] # produces multiple figures. only plot second one.
```

It is possible to combine multiple iterators by multiplying them together. This is equivalent to nesting one iterator inside the other. For example, we can repeat our cross-validation multiple times by permuting the sample order.


```{r}
P = permute_sample_order(number_of_permutations = 10) * 
    kfold_xval(folds=5,factor_name='Species')*
    (mean_centre() + PLSDA(factor_name='Species',number_components=2))
P = run(P,D,balanced_accuracy())
P$metric
```


# Session Info
```{r}
sessionInfo()
```
