mixOmics
Omics Data Integration Project

Package index
Search the mixOmics package
Vignettes
Functions
445
Source code
131
Man pages
86
Home
/
Bioconductor
/
mixOmics
/
mixOmics

mixOmics
In mixOmics: Omics Data Integration Project 

BiocStyle::markdown()

library(knitr)
# global options
knitr::opts_chunk$set(dpi = 100, echo=TRUE, warning=FALSE, message=FALSE, eval = TRUE,
                      fig.show=TRUE, fig.width= 7,fig.height= 6,fig.align='center', out.width = '80%'#, 
                      #fig.path= 'Figures/'
                      )

# knitr::opts_chunk$set(
#   collapse = TRUE,
#   comment = "#>"
# )

Bookdown version

This vignette is also available as a bookdown format on our gitHub page https://mixomicsteam.github.io/Bookdown/.

Introduction {#intro}

mixOmics is an R toolkit dedicated to the exploration and integration of biological data sets with a specific focus on variable selection. The package currently includes nineteen multivariate methodologies, mostly developed by the mixOmics team (see some of our references in \@ref(intro:pubs)). Originally, all methods were designed for omics data, however, their application is not limited to biological data only. Other applications where integration is required can be considered, but mostly for the case where the predictor variables are continuous (see also \@ref(intro:datatypes)). In mixOmics, a strong focus is given to graphical representation to better translate and understand the relationships between the different data types and visualize the correlation structure at both sample and variable levels.

Input data {#intro:datatypes}

Note the data pre-processing requirements before analysing data with mixOmics:

Methods

Some background knowledge {#intro:background}

We list here the main methodological or theoretical concepts you need to know to be able to efficiently apply mixOmics:

Overview {#intro:overview}

Here is an overview of the most widely used methods in mixOmics that will be further detailed in this vignette, with the exception of rCCA and MINT. We depict them along with the type of data set they can handle. 

library(mixOmics)
coul <- color.mixo(1:3)

plot(0, type="n", xlim=c(0,100), ylim=c(0,100), axes=FALSE,
     xlab="",ylab="", main="mixOmics overview")
box()

# PCA
rect(xleft = 20, ybottom = 75, xright = 40, ytop = 95, col=coul[1])
text(5, 85, "PCA")

# PLS-DA
rect(xleft = 20, ybottom = 50, xright = 40, ytop = 70, col=coul[1])
rect(xleft = 43, ybottom = 50, xright = 45, ytop = 70, col=coul[2])
text(5, 60, "PLS-DA")

# PLS
rect(xleft = 20, ybottom = 25, xright = 40, ytop = 45, col=coul[1])
rect(xleft = 43, ybottom = 25, xright = 60, ytop = 45, col=coul[1])
text(5, 35, "PLS")

# DIABLO
rect(xleft = 20, ybottom = 0, xright = 40, ytop = 20, col=coul[1])
rect(xleft = 43, ybottom = 0, xright = 60, ytop = 20, col=coul[1])
points(x=61, y=10, pch=16, col=coul[3], cex=0.5)
points(x=62.5, y=10, pch=16, col=coul[3], cex=0.5)
points(x=64, y=10, pch=16, col=coul[3], cex=0.5)
rect(xleft = 65, ybottom = 0, xright = 80, ytop = 20, col=coul[1])
rect(xleft = 85, ybottom = 0, xright = 88, ytop = 20, col=coul[2])
text(5, 10, "DIABLO")

# legend
rect(xleft = 75, ybottom = 95, xright = 77, ytop = 97, col=coul[1])
text(90, 96, "Quantitative", cex=0.75)
rect(xleft = 75, ybottom = 90, xright = 77, ytop = 92, col=coul[2])
text(90, 91, "Qualitative", cex=0.75)

knitr::include_graphics("XtraFigs/Methods.png")

knitr::include_graphics("XtraFigs/cheatsheet.png")

Key publications {#intro:pubs}

Methods papers

The methods implemented in mixOmics are described in detail in the following publications. A more extensive list can be found at this link.

Biological application papers

For those interested in how these methods have been applied to omics biological problems:

Outline of this Vignette

Each of the methods chapter has the following outline:

  1. Type of biological question to be answered
  2. Brief description of an illustrative data set
  3. Principle of the method
  4. Quick start of the method with the main functions and arguments
  5. To go further: customized plots, additional graphical outputs and tuning parameters
  6. FAQ

Other methods not covered in this vignette

Other methods not covered in this document are described on our website and the following references:

Let's get started {#start}

Installation

First, download the latest \texttt{mixOmics} version from Bioconductor.:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
 BiocManager::install("mixOmics")

Alternatively, you can install the latest development version of the package from Github.

install.packages("devtools")
# then load
library(devtools)
install_github("mixOmicsTeam/mixOmics")

The mixOmics package will directly import the following packages: r Rpackage("igraph"), r Rpackage("rgl"), r Rpackage("ellipse"), r Rpackage("corpcor"), r Rpackage("RColorBrewer"), r Rpackage("parallel"), r Rpackage("dplyr"), r Rpackage("tidyr"), r Rpackage("reshape2"), r Rpackage("methods"), r Rpackage("matrixStats"), r Rpackage("rARPACK"), r Rpackage("gridExtra").

For apple mac users, if you are unable to install the imported library rgl, you will need to install the XQuartz software first.

Load the package {#start:upload}

library(mixOmics)

Check that there is no error when loading the package, especially for the rgl library (see above).

Upload data

The examples we give in this vignette use data that are already part of the package, which means that the objects are stored as list. However, you do not need to store your data as list, just data frames will do. To upload your own data, check first that your working directory is set, then read your data from a .txt or .csv format, either by using File > Import Dataset in RStudio (although be mindful this option may pose problem with row and col name headers) or via one of these command lines:

# from csv file
data <- read.csv("your_data.csv", row.names = 1, header = TRUE)

# from txt file
data <- read.table("your_data.txt", header = TRUE)

# then in the argument in the mixOmics functions, just fill with:
# X = data

For more details about the arguments used to modify those functions, type ?read.csv or ?read.table in the R console. In this vignette we assume you are sufficiently proficient in R to input your data to carry on with the mixOmics analyses!

Quick start in mixOmics {#start:PCA}

Each analysis should follow this workflow:

  1. Run the method
  2. Graphical representation of the samples
  3. Graphical representation of the variables

Then use your critical thinking and additional functions and visual tools to make sense of your data! (some of which are listed in \@ref(intro:overview)) and will be described in the next Chapters.

For instance, for Principal Components Analysis, we first load the data:

data(nutrimouse)
X <- nutrimouse$gene

Then use the following steps:

MyResult.pca <- pca(X)  # 1 Run the method
plotIndiv(MyResult.pca) # 2 Plot the samples
plotVar(MyResult.pca)   # 3 Plot the variables

This is only a first quick-start, there will be many avenues you can take to deepen your exploratory and integrative analyses. The package proposes several methods to perform variable, or feature selection to identify the relevant information from rather large omics data sets. The sparse methods are listed in the Table in \@ref(intro:overview).

Following our example here, sparse PCA can be applied to select the top 5 variables contributing to each of the two components in PCA. The user specifies the number of variables to selected on each component, for example here 5 variables are selected on each of the first two components (keepX=c(5,5)): 

MyResult.spca <- spca(X, keepX=c(5,5)) # 1 Run the method
plotIndiv(MyResult.spca)               # 2 Plot the samples
plotVar(MyResult.spca)                 # 3 Plot the variables

You can see know that we have considerably reduced the number of genes in the plotVar correlation circle plot.

Do not stop here! We are not done yet. You can enhance your analyses with the following:

Principal Component Analysis (PCA) {#pca}

library(mixOmics)
coul <- color.mixo(1:3)

plot(0, type="n", xlim=c(0,100), ylim=c(83,97), axes=FALSE,
     xlab="",ylab="", main="PCA overview")
box()

# PCA
rect(xleft = 20, ybottom = 85, xright = 60, ytop = 95, col=coul[1])
text(5, 90, "PCA")

# legend
rect(xleft = 85, ybottom = 90, xright = 87, ytop = 92, col=coul[1])
text(90, 94, "Quantitative", cex=0.75)
#rect(xleft = 75, ybottom = 90, xright = 77, ytop = 92, col=coul[2])
#text(90, 91, "Qualitative", cex=0.75)

Biological question

I would like to identify the major sources of variation in my data and itdentify whether such sources of variation correspond to biological conditions, or experimental bias. I would like to visualise trends or patterns between samples,whether they 'naturally' cluster according to known biological conditions.

The liver.toxicity study

The liver.toxicity is a list in the package that contains:

More details are available at ?liver.toxicity.

To illustrate PCA, we focus on the expression levels of the genes in the data frame liver.toxicity$gene. Some of the terms mentioned below are listed in \@ref(intro:background).

Principle of PCA

The aim of PCA [@Jol05] is to reduce the dimensionality of the data whilst retaining as much information as possible. 'Information' is referred here as variance. The idea is to create uncorrelated artificial variables called principal components (PCs) that combine in a linear manner the original (possibly correlated) variables (e.g. genes, metabolites, etc.).

Dimension reduction is achieved by projecting the data into the space spanned by the principal components (PC). In practice, it means that each sample is assigned a score on each new PC dimension - this score is calculated as a linear combination of the original variables to which a weight is applied. The weights of each of the original variables are stored in the so-called loading vectors associated to each PC. The dimension of the data is reduced by projecting the data into the smaller subspace spanned by the PCs, while capturing the largest sources of variation between samples.

The principal components are obtained so that their variance is maximised. To that end, we calculate the eigenvectors/eigenvalues of the variance-covariance matrix, often via singular value decomposition when the number of variables is very large. The data are usually centred (center = TRUE), and sometimes scaled (scale = TRUE) in the method. The latter is especially advised in the case where the variance is not homogeneous across variables.

The first PC is defined as the linear combination of the original variables that explains the greatest amount of variation. The second PC is then defined as the linear combination of the original variables that accounts for the greatest amount of the remaining variation subject of being orthogonal (uncorrelated) to the first component. Subsequent components are defined likewise for the other PCA dimensions. The user must therefore report how much information is explained by the first PCs as these are used to graphically represent the PCA outputs.

Load the data

We first load the data from the package. See \@ref(start:upload) to upload your own data. 

library(mixOmics)
data(liver.toxicity)
X <- liver.toxicity$gene

Quick start

MyResult.pca <- pca(X)     # 1 Run the method
plotIndiv(MyResult.pca)    # 2 Plot the samples
plotVar(MyResult.pca)      # 3 Plot the variables

If you were to run pca with this minimal code, you would be using the following default values:

Other arguments can also be chosen, see ?pca.

This example was shown in Chapter \@ref(start:PCA). The two plots are not extremely meaningful as specific sample patterns should be further investigated and the variable correlation circle plot contains too many variables to be easily interpreted. Let's improve those graphics as shown below to improve interpretation.

To go further

Customize plots

Plots can be customized using numerous options in plotIndiv and plotVar. For instance, even if PCA does not take into account any information regarding the known group membership of each sample, we can include such information on the sample plot to visualize any `natural' cluster that may corresponds to biological conditions.

Here is an example where we include the sample groups information with the argument group:

plotIndiv(MyResult.pca, group = liver.toxicity$treatment$Dose.Group, 
          legend = TRUE)

Additionally, two factors can be displayed using both colours (argument group) and symbols (argument pch). For example here we display both Dose and Time of exposure and improve the title and legend:

plotIndiv(MyResult.pca, ind.names = FALSE,
          group = liver.toxicity$treatment$Dose.Group,
          pch = as.factor(liver.toxicity$treatment$Time.Group),
          legend = TRUE, title = 'Liver toxicity: genes, PCA comp 1 - 2',
          legend.title = 'Dose', legend.title.pch = 'Exposure')

By including information related to the dose of acetaminophen and time of exposure enables us to see a cluster of low dose samples (blue and orange, top left at 50 and 100mg respectively), whereas samples with high doses (1500 and 2000mg in grey and green respectively) are more scattered, but highlight an exposure effect.

To display the results on other components, we can change the comp argument provided we have requested enough components to be calculated. Here is our second PCA with 3 components:

MyResult.pca2 <- pca(X, ncomp = 3)
plotIndiv(MyResult.pca2, comp = c(1,3), legend = TRUE,
          group = liver.toxicity$treatment$Time.Group,
          title = 'Multidrug transporter, PCA comp 1 - 3')

Here, the 3rd component on the y-axis clearly highlights a time of exposure effect.

Amount of variance explained and choice of number of components

The amount of variance explained can be extracted with the following: a screeplot or the actual numerical proportions of explained variance, and cumulative proportion.

plot(MyResult.pca2)
MyResult.pca2

There are no clear guidelines on how many components should be included in PCA: it is data dependent and level of noise dependent. We often look at the `elbow' on the screeplot above as an indicator that the addition of PCs does not drastically contribute to explain the remainder variance.

Other useful plots

We can also have a look at the variable coefficients in each component with the loading vectors. The loading weights are represented in decreasing order from bottom to top in plotLoadings. Their absolute value indicates the importance of each variable to define each PC, as represented by the length of each bar. See ?plotLoadings to change the arguments.

# a minimal example
plotLoadings(MyResult.pca)

# a customized example to only show the top 100 genes 
# and their gene name
plotLoadings(MyResult.pca, ndisplay = 100, 
             name.var = liver.toxicity$gene.ID[, "geneBank"],
             size.name = rel(0.3))

Such representation will be more informative once we select a few variables in the next section \@ref(sPCA).

Plots can also be displayed in 3 dimensions using the option style="3d", and interactively (we use the rgl package for this).

plotIndiv(MyResult.pca2,
          group = liver.toxicity$treatment$Dose.Group, style="3d",
          legend = TRUE, title = 'Liver toxicity: genes, PCA comp 1 - 2 - 3')

Variable selection with sparse PCA {#sPCA}

Biological question

I would like to apply PCA but also be able to identify the key variables that contribute to the explanation of most variance in the data set.

Variable selection can be performed using the sparse version of PCA implemented in spca [@She08]. The user needs to provide the number of variables to select on each PC. Here for example we ask to select the top 15 genes contributing to the definition of PC1, the top 10 genes contributing to PC2 and the top 5 genes for PC3 (keepX=c(15,10,5)).

MyResult.spca <- spca(X, ncomp = 3, keepX = c(15,10,5))                 # 1 Run the method
plotIndiv(MyResult.spca, group = liver.toxicity$treatment$Dose.Group,   # 2 Plot the samples
          pch = as.factor(liver.toxicity$treatment$Time.Group),
          legend = TRUE, title = 'Liver toxicity: genes, sPCA comp 1 - 2',
          legend.title = 'Dose', legend.title.pch = 'Exposure')
plotVar(MyResult.spca, cex = 1)                                        # 3 Plot the variables
# cex is used to reduce the size of the labels on the plot

Selected variables can be identified on each component with the selectVar function. Here the coefficient values are extracted, but there are other outputs, see ?selectVar:

selectVar(MyResult.spca, comp = 1)$value

Those values correspond to the loading weights that are used to define each component. A large absolute value indicates the importance of the variable in this PC. Selected variables are ranked from the most important (top) to the least important.

We can complement this output with plotLoadings. We can see here that all coefficients are negative.

plotLoadings(MyResult.spca)

If we look at component two, we can see a mix of positive and negative weights (also see in the plotVar), those correspond to variables that oppose the low and high doses (see from the `plotIndiv):

selectVar(MyResult.spca, comp=2)$value
plotLoadings(MyResult.spca, comp = 2)

Tuning parameters

For this set of methods, two parameters need to be chosen:

The function tune.pca calculates the percentage of variance explained for each component, up to the minimum between the number of rows, or column in the data set. The `optimal' number of components can be identified if an elbow appears on the screeplot. In the example below the cut-off is not very clear, we could choose 2 components.

tune.pca(X)

Regarding the number of variables to select in sparse PCA, there is not clear criterion at this stage. As PCA is an exploration method, we prefer to set arbitrary thresholds that will pinpoint the key variables to focus on during the interpretation stage.

Additional resources

Additional examples are provided in example(pca) and in our case studies on our website in the Methods and Case studies sections.

Additional reading in [@She08].

FAQ

PLS - Discriminant Analysis (PLS-DA) {#plsda}

library(mixOmics)
coul <- color.mixo(1:3)

plot(0, type="n", xlim=c(0,100), ylim=c(83,97), axes=FALSE,
     xlab="",ylab="", main="PLSDA overview")
box()

# PLS-DA
rect(xleft = 20, ybottom = 85, xright = 60, ytop = 95, col=coul[1])
rect(xleft = 63, ybottom = 85, xright = 65, ytop = 95, col=coul[2])
text(5, 90, "PLS-DA")


# legend
rect(xleft = 85, ybottom = 92, xright = 87, ytop = 94, col=coul[1])
text(90, 95, "Quantitative", cex=0.75)
rect(xleft = 85, ybottom = 87, xright = 87, ytop = 89, col=coul[2])
text(90, 90, "Qualitative", cex=0.75)

Biological question

I am analysing a single data set (e.g. transcriptomics data) and I would like to classify my samples into known groups and predict the class of new samples. In addition, I am interested in identifying the key variables that drive such discrimination.

The srbct study

The data are directly available in a processed and normalised format from the package. The Small Round Blue Cell Tumours (SRBCT) dataset from [@Kha01] includes the expression levels of 2,308 genes measured on 63 samples. The samples are classified into four classes as follows: 8 Burkitt Lymphoma (BL), 23 Ewing Sarcoma (EWS), 12 neuroblastoma (NB), and 20 rhabdomyosarcoma (RMS).

The srbct dataset contains the following:

$gene: a data frame with 63 rows and 2308 columns. The expression levels of 2,308 genes in 63 subjects.

$class: a class vector containing the class tumour of each individual (4 classes in total).

$gene.name: a data frame with 2,308 rows and 2 columns containing further information on the genes.

More details can be found in ?srbct.

To illustrate PLS-DA, we will analyse the gene expression levels of srbct$gene to discriminate the 4 groups of tumours.

Principle of sparse PLS-DA

Although Partial Least Squares was not originally designed for classification and discrimination problems, it has often been used for that purpose [@Ngu02a;@Tan04]. The response matrix Y is qualitative and is internally recoded as a dummy block matrix that records the membership of each observation, i.e. each of the response categories are coded via an indicator variable (see [@mixomics] Suppl. Information S1 for an illustration). The PLS regression (now PLS-DA) is then run as if Y was a continuous matrix. This PLS classification trick works well in practice, as demonstrated in many references [@Bar03;@Ngu02a;@Bou07;@Chung10].

Sparse PLS-DA [@Lec11] performs variable selection and classification in a one step procedure. sPLS-DA is a special case of sparse PLS described later in \@ref(pls), where $\ell_1$ penalization is applied on the loading vectors associated to the X data set.

Inputs and outputs

We use the following data input matrices: X is a $n \times p$ data matrix, Y is a factor vector of length $n$ that indicates the class of each sample, and $Y^*$ is the associated dummy matrix ($n \times K$) with $n$ the number of samples (individuals), $p$ the number of variables and $K$ the number of classes. PLS-DA main outputs are:

Set up the data

We first load the data from the package. See \@ref(start:upload) to upload your own data.

We will mainly focus on sparse PLS-DA that is more suited for large biological data sets where the aim is to identify molecular signatures, as well as classifying samples. We first set up the data as X expression matrix and Y as a factor indicating the class membership of each sample. We also check that the dimensions are correct and match:

library(mixOmics)
data(srbct)
X <- srbct$gene
Y <- srbct$class 
summary(Y)
dim(X); length(Y)

Quick start

For a quick start we arbitrarily set the number of variables to select to 50 on each of the 3 components of PLS-DA (see section \@ref(tuning:sPLSDA) for tuning these values).

MyResult.splsda <- splsda(X, Y, keepX = c(50,50)) # 1 Run the method
plotIndiv(MyResult.splsda)                          # 2 Plot the samples
plotVar(MyResult.splsda)                            # 3 Plot the variables
selectVar(MyResult.splsda, comp=1)$name             # Selected variables on component 1

As PLS-DA is a supervised method, the sample plot automatically displays the group membership of each sample. We can observe a clear discrimination between the BL samples and the others on the first component (x-axis), and EWS vs the others on the second component (y-axis). Remember that this discrimination spanned by the first two PLS-DA components is obtained based on a subset of 100 variables (50 selected on each component).

From the plotIndiv the axis labels indicate the amount of variation explained per component. Note that the interpretation of this amount is not the same as in PCA. In PLS-DA, the aim is to maximise the covariance between X and Y, not only the variance of X as it is the case in PCA!

If you were to run splsda with this minimal code, you would be using the following default values:

It is important to note that all Discriminat Analyses (e.g. plsda, mint.splsda etc.) use mode='regression' for deflation.

PLS-DA without variable selection can be performed as:

MyResult.plsda <- plsda(X,Y) # 1 Run the method
plotIndiv(MyResult.plsda)    # 2 Plot the samples
plotVar(MyResult.plsda)      # 3 Plot the variables

To go further {#plsda-tgf}

Customize sample plots {#splsda:plotIndiv}

The sample plots can be improved in various ways. First, if the names of the samples are not meaningful at this stage, they can be replaced by symbols (ind.names=TRUE). Confidence ellipses can be plotted for each sample (ellipse = TRUE, confidence level set to 95\% by default, see the argument ellipse.level), Additionally, a star plot displays arrows from each group centroid towards each individual sample (star = TRUE). A 3D plot is also available, see plotIndiv for more details.

plotIndiv(MyResult.splsda, ind.names = FALSE, legend=TRUE,
          ellipse = TRUE, star = TRUE, title = 'sPLS-DA on SRBCT',
          X.label = 'PLS-DA 1', Y.label = 'PLS-DA 2')

Customize variable plots

The name of the variables can be set to FALSE (var.names=FALSE):

plotVar(MyResult.splsda, var.names=FALSE)

In addition, if we had used the non sparse version of PLS-DA, a cut-off can be set to display only the variables that mostly contribute to the definition of each component. Those variables should be located towards the circle of radius 1, far from the centre. 

plotVar(MyResult.plsda, cutoff=0.7)

In this particular case, no variable selection was performed. Only the display was altered to show a subset of variables.

Other useful plots

Background prediction

A 'prediction' background can be added to the sample plot by calculating a background surface first, before overlaying the sample plot. See ?background.predict for more details. More details about prediction, prediction distances can be found in [@mixomics] in the Suppl. Information.

background <- background.predict(MyResult.splsda, comp.predicted=2,
                                dist = "max.dist") 
plotIndiv(MyResult.splsda, comp = 1:2, group = srbct$class,
          ind.names = FALSE, title = "Maximum distance",
          legend = TRUE,  background = background)

ROC

As PLS-DA acts as a classifier, we can plot a ROC Curve to complement the sPLS-DA classification performance results detailed in \@ref(tuning:sPLSDA). The AUC is calculated from training cross-validation sets and averaged. Note however that ROC and AUC criteria may not be particularly insightful, or may not be in full agreement with the PLSDA performance, as the prediction threshold in PLS-DA is based on specified distance as described in [@mixomics].

auc.plsda <- auroc(MyResult.splsda)

Variable selection outputs

First, note that the number of variables to select on each component does not need to be identical on each component, for example: 

MyResult.splsda2 <- splsda(X,Y, ncomp=3, keepX=c(15,10,5))

Selected variables are listed in the selectVar function: 

selectVar(MyResult.splsda2, comp=1)$value

and can be visualised in plotLoadings with the arguments contrib = 'max' that is going to assign to each variable bar the sample group colour for which the mean (method = 'mean') is maximum. See example(plotLoadings) for other options (e.g. min, median)

plotLoadings(MyResult.splsda2, contrib = 'max', method = 'mean')

Interestingly from this plot, we can see that all selected variables on component 1 are highly expressed in the BL (orange) class. Setting contrib = 'min' would highlight that those variables are lowly expressed in the NB grey class, which makes sense when we look at the sample plot.

Since 4 classes are being discriminated here, samples plots in 3d may help interpretation: 

plotIndiv(MyResult.splsda2, style="3d")

Tuning parameters and numerical outputs {#tuning:sPLSDA}

For this set of methods, three parameters need to be chosen:

1 - The number of components to retain ncomp. The rule of thumb is usually $K - 1$ where $K$ is the number of classes, but it is worth testing a few extra components.

2 - The number of variables keepX to select on each component for sparse PLS-DA,

3 - The prediction distance to evaluate the classification and prediction performance of PLS-DA.

For item 1, the perf evaluates the performance of PLS-DA for a large number of components, using repeated k-fold cross-validation. For example here we use 3-fold CV repeated 10 times (note that we advise to use at least 50 repeats, and choose the number of folds that are appropriate for the sample size of the data set):

MyResult.plsda2 <- plsda(X,Y, ncomp=10)
set.seed(30) # for reproducibility in this vignette, otherwise increase nrepeat
MyPerf.plsda <- perf(MyResult.plsda2, validation = "Mfold", folds = 3, 
                  progressBar = FALSE, nrepeat = 10) # we suggest nrepeat = 50

# type attributes(MyPerf.plsda) to see the different outputs
# slight bug in our output function currently see the quick fix below
#plot(MyPerf.plsda, col = color.mixo(5:7), sd = TRUE, legend.position = "horizontal")

# quick fix
matplot(MyPerf.plsda$error.rate$BER, type = 'l', lty = 1, 
        col = color.mixo(1:3), 
        main = 'Balanced Error rate')
legend('topright', 
       c('max.dist', 'centroids.dist', 'mahalanobis.dist'), 
       lty = 1,
       col = color.mixo(5:7))

The plot outputs the classification error rate, or Balanced classification error rate when the number of samples per group is unbalanced, the standard deviation according to three prediction distances. Here we can see that for the BER and the maximum distance, the best performance (i.e. low error rate) seems to be achieved for ncomp = 3.

In addition (item 3 for PLS-DA), the numerical outputs listed here can be reported as performance measures: 

MyPerf.plsda

Regarding item 2, we now use tune.splsda to assess the optimal number of variables to select on each component. We first set up a grid of keepX values that will be assessed on each component, one component at a time. Similar to above we run 3-fold CV repeated 10 times with a maximum distance prediction defined as above. 

list.keepX <- c(5:10,  seq(20, 100, 10))
list.keepX # to output the grid of values tested
set.seed(30) # for reproducbility in this vignette, otherwise increase nrepeat
tune.splsda.srbct <- tune.splsda(X, Y, ncomp = 3, # we suggest to push ncomp a bit more, e.g. 4
                                 validation = 'Mfold',
                                 folds = 3, dist = 'max.dist', progressBar = FALSE,
                                 measure = "BER", test.keepX = list.keepX,
                                 nrepeat = 10)   # we suggest nrepeat = 50

We can then extract the classification error rate averaged across all folds and repeats for each tested keepX value, the optimal number of components (see ?tune.splsda for more details), the optimal number of variables to select per component which is summarised in a plot where the diamond indicated the optimal keepX value:

error <- tune.splsda.srbct$error.rate
ncomp <- tune.splsda.srbct$choice.ncomp$ncomp # optimal number of components based on t-tests on the error rate
ncomp
select.keepX <- tune.splsda.srbct$choice.keepX[1:ncomp]  # optimal number of variables to select
select.keepX
plot(tune.splsda.srbct, col = color.jet(ncomp))

Based on those tuning results, we can run our final and tuned sPLS-DA model:

MyResult.splsda.final <- splsda(X, Y, ncomp = ncomp, keepX = select.keepX)
plotIndiv(MyResult.splsda.final, ind.names = FALSE, legend=TRUE,
          ellipse = TRUE, title="SPLS-DA, Final result")

Additionally we can run perf for the final performance of the sPLS-DA model. Also note that perf will output features that lists the frequency of selection of the variables across the different folds and different repeats. This is a useful output to assess the confidence of your final variable selection, see a more detailed example here.

Additional resources

Additional examples are provided in example(splsda) and in our case studies on our website in the Methods and Case studies sections, and in particular here. Also have a look at [@Lec11]

FAQ

Projection to Latent Structure (PLS) {#pls}

library(mixOmics)
coul <- color.mixo(1:3)

plot(0, type="n", xlim=c(0,100), ylim=c(83,97), axes=FALSE,
     xlab="",ylab="", main="PLS overview")
box()

# PLS
rect(xleft = 20, ybottom = 85, xright = 50, ytop = 95, col=coul[1])
rect(xleft = 52, ybottom = 85, xright = 73, ytop = 95, col=coul[1])
text(5, 90, "PLS")


# legend
rect(xleft = 85, ybottom = 92, xright = 87, ytop = 94, col=coul[1])
text(90, 95, "Quantitative", cex=0.75)

Biological question

I would like to integrate two data sets measured on the same samples by extracting correlated information, or by highlighing commonalities between data sets.

The nutrimouse study

The nutrimouse study contains the expression levels of genes potentially involved in nutritional problems and the concentrations of hepatic fatty acids for forty mice. The data sets come from a nutrigenomic study in the mouse from our collaborator [@Mar07], in which the effects of five regimens with contrasted fatty acid compositions on liver lipids and hepatic gene expression in mice were considered. Two sets of variables were measured on 40 mice:

More details can be found in ?nutrimouse.

To illustrate sparse PLS, we will integrate the gene expression levels (gene) with the concentrations of hepatic fatty acids (lipid).

Principle of PLS

Partial Least Squares (PLS) regression [@Wol66; @Wol01] is a multivariate methodology which relates (\textit{integrates}) two data matrices X (e.g. transcriptomics) and Y (e.g. lipids). PLS goes beyond traditional multiple regression by modelling the structure of both matrices. Unlike traditional multiple regression models, it is not limited to uncorrelated variables. One of the many advantages of PLS is that it can handle many noisy, collinear (correlated) and missing variables and can also simultaneously model several response variables in Y.

PLS is a multivariate projection-based method that can address different types of integration problems. Its flexibility is the reason why it is the backbone of most methods in mixOmics. PLS is computationally very efficient when the number of variables $p + q >> n$ the number of samples. It performs successive local regressions that avoid computational issues due to the inversion of large singular covariance matrices. Unlike PCA which maximizes the variance of components from a single data set, PLS maximizes the covariance between components from two data sets. The mathematical concepts of covariance and correlation are similar, but the covariance is an unbounded measure and covariance has a unit measure (see \@ref(intro:background)). In PLS, linear combination of variables are called latent variables or latent components. The weight vectors used to calculate the linear combinations are called the loading vectors. Latent variables and loading vectors are thus associated, and come in pairs from each of the two data sets being integrated.

Principle of sparse PLS

Even though PLS is highly efficient in a high dimensional context, the interpretability of PLS needed to be improved. sPLS has been recently developed by our team to perform simultaneous variable selection in both data sets X and Y data sets, by including LASSO $\ell_1$ penalizations in PLS on each pair of loading vectors [@Lec08].

Inputs and outputs

We consider the data input matrices: X is a $n \times p$ data matrix and Y a $n \times q$ data matrix, where $n$ the number of samples (individuals), $p$ and $q$ are the number of variables in each data set. PLS main outputs are:

Set up the data

We first set up the data as X expression matrix and Y as the lipid abundance matrix. We also check that the dimensions are correct and match:

library(mixOmics)
data(nutrimouse)
X <- nutrimouse$gene  
Y <- nutrimouse$lipid
dim(X); dim(Y)

Quick start

We will mainly focus on sparse PLS for large biological data sets where variable selection can help the interpretation of the results. See ?pls for a model with no variable selection. Here we arbitrarily set the number of variables to select to 50 on each of the 2 components of PLS (see section \@ref(tuning:PLS) for tuning these values).

MyResult.spls <- spls(X,Y, keepX = c(25, 25), keepY = c(5,5))  
plotIndiv(MyResult.spls)                                      
plotVar(MyResult.spls)

If you were to run spls with this minimal code, you would be using the following default values:

Because PLS generates a pair of components, each associated to each data set, the function plotIndiv produces 2 plots that represent the same samples projected in either the space spanned by the X-components, or the Y-components. A single plot can also be displayed, see section \@ref(pls:plotIndiv).

To go further {#pls-tgf}

Customize sample plots {#pls:plotIndiv}

Some of the sample plot additional arguments were described in \@ref(splsda:plotIndiv). In addition, you can choose the representation space to be either the components from the X-data set, the Y- data set, or an average between both components rep.space = 'XY-variate'. See more examples in examples(plotIndiv) and on our website. Here are two examples with colours indicating genotype or diet:

plotIndiv(MyResult.spls, group = nutrimouse$genotype,
          rep.space = "XY-variate", legend = TRUE,
          legend.title = 'Genotype',
          ind.names = nutrimouse$diet,
          title = 'Nutrimouse: sPLS')

plotIndiv(MyResult.spls, group=nutrimouse$diet,
          pch = nutrimouse$genotype,
          rep.space = "XY-variate",  legend = TRUE,
          legend.title = 'Diet', legend.title.pch = 'Genotype',
          ind.names = FALSE, 
          title = 'Nutrimouse: sPLS')

Customize variable plots {#pls:plotVar}

See (example(plotVar)) for more examples. Here we change the size of the labels. By default the colours are assigned to each type of variable. The coordinates of the variables can also be saved as follows:

plotVar(MyResult.spls, cex=c(3,2), legend = TRUE)
coordinates <- plotVar(MyResult.spls, plot = FALSE)

Other useful plots for data integration

We extended other types of plots, based on clustered image maps and relevance networks to ease the interpretation of the relationships between two types of variables. A similarity matrix is calculated from the outputs of PLS and represented with those graphics, see [@Gon12] for more details, and our website

Clustered Image Maps

A clustered image map can be produced using the cim function. You may experience figures margin issues in RStudio. Best is to either use X11() or save the plot as an external file. For example to show the correlation structure between the X and Y variables selected on component 1:

X11()
cim(MyResult.spls, comp = 1)
cim(MyResult.spls, comp = 1, save = 'jpeg', name.save = 'PLScim')

Relevance networks {#pls:network}

Using the same similarity matrix input in CIM, we can also represent relevance bipartite networks. Those networks only represent edges between on type of variable from X and the other type of variable, from Y. Whilst we use sPLS to narrow down to a few key correlated variables, our keepX and keepY values might still be very high for this kind of output. A cut-off can be set based on the correlation coefficient between the different types of variables.

Other arguments such as interactive = TRUE enables a scrollbar to change the cut-off value interactively, see other options in ?network. Additionally, the graph object can be saved to be input into Cytoscape for an improved visualisation.

X11()
network(MyResult.spls, comp = 1)
network(MyResult.spls, comp = 1, cutoff = 0.6, save = 'jpeg', name.save = 'PLSnetwork')
# save as graph object for cytoscape
myNetwork <- network(MyResult.spls, comp = 1)$gR

Arrow plots

Instead of projecting the samples into the combined XY representation space, as shown in \@ref(pls:plotIndiv), we can overlap the X- and Y- representation plots. One arrow joins the same sample from the X- space to the Y- space. Short arrows indicate a good agreement found by the PLS between both data sets.

plotArrow(MyResult.spls,group=nutrimouse$diet, legend = TRUE,
          X.label = 'PLS comp 1', Y.label = 'PLS comp 2', legend.title = 'Diet')

Variable selection outputs

The selected variables can be extracted using the selectVar function for further analysis.

MySelectedVariables <- selectVar(MyResult.spls, comp = 1)
MySelectedVariables$X$name # Selected genes on component 1
MySelectedVariables$Y$name # Selected lipids on component 1

The loading plots help visualise the coefficients assigned to each selected variable on each component:

plotLoadings(MyResult.spls, comp = 1, size.name = rel(0.5))

Tuning parameters and numerical outputs {#tuning:PLS}

For PLS and sPLS, two types of parameters need to be chosen:

1 - The number of components to retain ncomp,

2 - The number of variables to select on each component and on each data set keepX and keepY for sparse PLS.

For item 1 we use the perf function and repeated k-fold cross-validation to calculate the Q$^2$ criterion used in the SIMCA-P software [@Ume96]. The rule of thumbs is that a PLS component should be included in the model if its value is $\leq 0.0975$. Here we use 3-fold CV repeated 10 times (note that we advise to use at least 50 repeats, and choose the number of folds that are appropriate for the sample size of the data set).

We run a PLS model with a sufficient number of components first, then run perf on the object.

MyResult.pls <- pls(X,Y, ncomp = 4)  
set.seed(30) # for reproducibility in this vignette, otherwise increase nrepeat
perf.pls <- perf(MyResult.pls, validation = "Mfold", folds = 5,
                  progressBar = FALSE, nrepeat = 10)
plot(perf.pls, criterion = 'Q2.total')

We can see that the values decrease as we add new components which generally indicates overfitting so a model with high number of components would not be suitable for prediction.

Item 2 can be quite difficult to tune. Here is a minimal example where we only tune keepX based on the correlation between the cross-validated components and those from the full model with all the data. Other measures proposed are Residual Sum of Squares (RSS) (see ?tune.spls):

list.keepX <- c(2:10, 15, 20)
# tuning based on correlations
set.seed(30) # for reproducibility in this vignette, otherwise increase nrepeat
tune.spls.cor <- tune.spls(X, Y, ncomp = 3,
                           test.keepX = list.keepX,
                           validation = "Mfold", folds = 5,
                           nrepeat = 10, progressBar = FALSE,
                           measure = 'cor')
plot(tune.spls.cor, measure = 'cor')

Based on the correlation of the cross-validated components, the optimal number of variables to select in the X data set, including all variables in the Y data set would be:

tune.spls.cor$choice.keepX

Tuning keepX and keepY conjointly is also possible:

# tuning both X and Y
set.seed(30) # for reproducibility in this vignette, otherwise increase nrepeat
tune.spls.cor.XY <- tune.spls(X, Y, ncomp = 3,
                           test.keepX = c(8, 20, 50),
                           test.keepY = c(4, 8, 16),
                           validation = "Mfold", folds = 5,
                           nrepeat = 10, progressBar = FALSE,
                           measure = 'cor')
## visualise correlations
plot(tune.spls.cor.XY, measure = 'cor')
## visualise RSS
plot(tune.spls.cor.XY, measure = 'RSS')

PLS modes {#PLS:details}

You may have noticed the mode argument in PLS. We can calculate the residual matrices at each PLS iteration differently. Note: this is for advanced users.

Regression mode The PLS regression mode models uni-directional (or 'causal') relationship between two data sets. The Y matrix is deflated with respect to the information extracted/modelled from the local regression on X. Here the goal is to predict Y from X (Y and X play an assymmetric role). Consequently the latent variables computed to predict Y from X are different from those computed to predict X from Y. More details about the model can be found in[ Appendix [@Lec08].

PLS regression mode, also called PLS2, is commonly applied for the analysis of biological data [@Bou05; @Byl07] due to the biological assumptions or the biological dogma. In general, the number of variables in Y to predict are fewer in number than the predictors in X.

Canonical mode Similar to a Canonical Correlation Analysis (CCA) framework, this mode is used to model a bi-directional (or symmetrical) relationship between the two data sets. The Y matrix is deflated with respect to the information extracted or modelled from the local regression on Y. Here X and Y play a symmetric role and the goal is similar to CCA. More details about the model can be found in [@Lec09a].

PLS canonical mode is not well known (yet), but is applicable when there is no a priori relationship between the two data sets, or in place of CCA but when variable selection is required in large data sets. In [@Lec09a], we compared the measures of the same biological samples on different types of microarrays, cDNA and Affymetrix arrays, to highlight complementary information at the transcripts levels. Note however that for this mode we do not provide any tuning function.

Other modes The 'invariant' mode performs a redundancy analysis, where the Y matrix is not deflated. The 'classic' mode is similar to a regression mode. It gives identical results for the variates and loadings associated to the X data set, but differences for the loadings vectors associated to the Y data set (different normalisations are used). Classic mode is the PLS2 model as defined by [@Ten98], Chap 9.

Difference between PLS modes For the first PLS dimension, all PLS modes will output the same results in terms of latent variables and loading vectors. After the first dimension, these vectors will differ, as the matrices are deflated differently.

Additional resources

Additional examples are provided in example(spls) and in our case studies on our website in the Methods and Case studies sections, see also [@Lec08; @Lec09a].

FAQ

Multi-block Discriminant Analysis with DIABLO {#diablo}

library(mixOmics)
coul <- color.mixo(1:3)

plot(0, type="n", xlim=c(0,100), ylim=c(83,97), axes=FALSE,
     xlab="",ylab="", main="DIABLO overview")
box()


rect(xleft = 12, ybottom = 85, xright = 30, ytop = 95, col=coul[1])
rect(xleft = 32, ybottom = 85, xright = 45, ytop = 95, col=coul[1])
points(x=47, y=90, pch=16, col='black', cex=0.5)
points(x=48.5, y=90, pch=16, col='black', cex=0.5)
points(x=50, y=90, pch=16, col='black', cex=0.5)
rect(xleft = 52, ybottom = 85, xright = 61, ytop = 95, col=coul[1])
rect(xleft = 63, ybottom = 85, xright = 65, ytop = 95, col=coul[2])
text(5, 90, "DIABLO")


# legend
rect(xleft = 85, ybottom = 92, xright = 87, ytop = 94, col=coul[1])
text(90, 95, "Quantitative", cex=0.75)
rect(xleft = 85, ybottom = 87, xright = 87, ytop = 89, col=coul[2])
text(90, 90, "Qualitative", cex=0.75)

DIABLO is our new framework in mixOmics that extends PLS for multiple data sets integration and PLS-discriminant analysis. The acronyms stands for Data Integration Analysis for Biomarker discovery using a Latent cOmponents [@Sin16].

Biological question

I would like to identify a highly correlated multi-omics signature discriminating known groups of samples.

The breast.TCGA study

Human breast cancer is a heterogeneous disease in terms of molecular alterations, cellular composition, and clinical outcome. Breast tumours can be classified into several subtypes, according to levels of mRNA expression [@Sor01. Here we consider a subset of data generated by The Cancer Genome Atlas Network [@TCGA12]. For the package, data were normalised and drastically prefiltered for illustrative purposes but DIABLO can handle larger data sets, see [@mixomics] Table 2. The data were divided into a training set with a subset of 150 samples from the mRNA, miRNA and proteomics data, and a test set including 70 samples, but only with mRNA and miRNA data (proteomics missing). The aim of this integrative analysis is to identify a highly correlated multi-omics signature discriminating the breast cancer subtypes Basal, Her2 and LumA.

The breast.TCGA is a list containing training and test sets of omics data data.train and data.test which both include:

More details can be found in ?breast.TCGA.

To illustrate DIABLO, we will integrate the expression levels of miRNA, mRNA and the abundance of proteins while discriminating the subtypes of breast cancer, then predict the subtypes of the test samples in the test set.

Principle of DIABLO

The core DIABLO method extends Generalised Canonical Correlation Analysis [@Ten11], which contrary to what its name suggests, generalises PLS for multiple matching datasets, and the sparse sGCCA method [@Ten14]. Starting from the R package RGCCA, we extended these methods for different types of analyses, including unsupervised N-integration (block.pls, block.spls) and supervised analyses (block.plsda, block.splsda).

The aim of N-integration with our sparse methods is to identify correlated (or co-expressed) variables measured on heterogeneous data sets which also explain the categorical outcome of interest (supervised analysis). The multiple data integration task is not trivial, as the analysis can be strongly affected by the variation between manufacturers or omics technological platforms despite being measured on the same biological samples. Before you embark on data integration, we strongly suggest individual or paired analyses with sPLS-DA and PLS to first understand the major sources of variation in each data set and to guide the integration process.

More methodological details can be found in [@Sin16].

Inputs and outputs

We consider as input a list X of data frames with $n$ rows (the number of samples) and different number of variations in each data frame. Y is a factor vector of length $n$ that indicates the class of each sample. Internally and similar to PLS-DA in Chapter \@ref(plsda) it will be coded as a dummy matrix.

DIABLO main outputs are:

Set up the data

We first set up the input data as a list of data frames X expression matrix and Y as a factor indicating the class membership of each sample. Each data frame in X should be named consistently to match with the keepX parameter.

We check that the dimensions are correct and match. We then set up arbitrarily the number of variables keepX that we wish to select in each data set and each component.

library(mixOmics)
data(breast.TCGA)
# extract training data and name each data frame
X <- list(mRNA = breast.TCGA$data.train$mrna, 
          miRNA = breast.TCGA$data.train$mirna, 
          protein = breast.TCGA$data.train$protein)
Y <- breast.TCGA$data.train$subtype
summary(Y)

list.keepX <- list(mRNA = c(16, 17), miRNA = c(18,5), protein = c(5, 5))

Quick start

MyResult.diablo <- block.splsda(X, Y, keepX=list.keepX)
plotIndiv(MyResult.diablo)
plotVar(MyResult.diablo)

Similar to PLS (Chapter \@ref(pls)), DIABLO generates a pair of components, each associated to each data set. This is why we can visualise here 3 sample plots. As DIABLO is a supervised method, samples are represented with different colours depending on their known class.

The variable plot suggest some correlation structure between proteins, mRNA and miRNA. We will further customize these plots in Sections \@ref(diablo:plotIndiv) and \@ref(diablo:plotVar).

If you were to run block.splsda with this minimal code, you would be using the following default values:

We focused here on the sparse version as would like to identify a minimal multi-omics signature, however, the non sparse version could also be run with block.plsda:

MyResult.diablo2 <- block.plsda(X, Y)

To go further

Customize sample plots {#diablo:plotIndiv}

Here is an example of an improved plot, see also Section \@ref(splsda:plotIndiv) for additional sources of inspiration.

plotIndiv(MyResult.diablo, 
          ind.names = FALSE, 
          legend=TRUE, cex=c(1,2,3),
          title = 'BRCA with DIABLO')

Customize variable plots {#diablo:plotVar}

Labels can be omitted in some data sets to improve readability. For example her we only show the name of the proteins:

plotVar(MyResult.diablo, var.names = c(FALSE, FALSE, TRUE),
        legend=TRUE, pch=c(16,16,1))

Other useful plots for data integration

Several plots were added for the DIABLO framework.

plotDiablo

A global overview of the correlation structure at the component level can be represented with the plotDiablo function. It plots the components across the different data sets for a given dimension. Colours indicate the class of each sample.

plotDiablo(MyResult.diablo, ncomp = 1)

Here, we can see that a strong correlation is extracted by DIABLO between the mRNA and protein data sets. Other dimensions can be plotted with the argument comp.

circosPlot

The circos plot represents the correlations between variables of different types, represented on the side quadrants. Several display options are possible, to show within and between connexions between blocks, expression levels of each variable according to each class (argument line = TRUE). The circos plot is built based on a similarity matrix, which was extended to the case of multiple data sets from [@Gon12]. A cutoff argument can be included to visualise correlation coefficients above this threshold in the multi-omics signature.

circosPlot(MyResult.diablo, cutoff=0.7)

cimDiablo

The cimDiablo function is a clustered image map specifically implemented to represent the multi-omics molecular signature expression for each sample. It is very similar to a classic hierarchical clustering:

# minimal example with margins improved:
# cimDiablo(MyResult.diablo, margin=c(8,20))

# extended example:
cimDiablo(MyResult.diablo, color.blocks = c('darkorchid', 'brown1', 'lightgreen'), comp = 1, margin=c(8,20), legend.position = "right")

plotLoadings

The plotLoadings function visualises the loading weights of each selected variables on each component (default is comp = 1) and each data set. The color indicates the class in which the variable has the maximum level of expression (contrib = "max") or minimum (contrib ="min"), on average (method="mean") or using the median (method ="median"). We only show the last plot here: 

#plotLoadings(MyResult.diablo, contrib = "max")
plotLoadings(MyResult.diablo, comp = 2, contrib = "max")

Relevance networks

Another visualisation of the correlation between the different types of variables is the relevance network, which is also built on the similarity matrix [@Gon12]. Each color represents a type of variable. A threshold can also be set using the argument cutoff.

See also Section \@ref(pls:network) to save the graph and the different options, or ?network.

network(MyResult.diablo, blocks = c(1,2,3),
        color.node = c('darkorchid', 'brown1', 'lightgreen'), 
        cutoff = 0.6, save = 'jpeg', name.save = 'DIABLOnetwork')

Numerical outputs

Classification performance

Similar to what is described in Section \@ref(tuning:sPLSDA) we use repeated cross-validation with perf to assess the prediction of the model. For this complex classification problems, often a centroid distance is suitable, see details in [@mixomics] Suppl. Material S1. 

set.seed(123) # for reproducibility in this vignette
MyPerf.diablo <- perf(MyResult.diablo, validation = 'Mfold', folds = 5, 
                   nrepeat = 10, 
                   dist = 'centroids.dist')

#MyPerf.diablo  # lists the different outputs

# Performance with Majority vote
#MyPerf.diablo$MajorityVote.error.rate

AUC

An AUC plot per block is plotted using the function auroc see [@mixomics] for the interpretation of such output as the ROC and AUC criteria are not particularly insightful in relation to the performance evaluation of our methods, but can complement the statistical analysis.

Here we evaluate the AUC for the model that includes 2 components in the miRNA data set.

Myauc.diablo <- auroc(MyResult.diablo, roc.block = "miRNA", roc.comp = 2)

Prediction on an external test set

The predict function predicts the class of samples from a test set. In our specific case, one data set is missing in the test set but the method can still be applied. Make sure the name of the blocks correspond exactly.

# prepare test set data: here one block (proteins) is missing
X.test <- list(mRNA = breast.TCGA$data.test$mrna, 
                      miRNA = breast.TCGA$data.test$mirna)

Mypredict.diablo <- predict(MyResult.diablo, newdata = X.test)
# the warning message will inform us that one block is missing
#Mypredict.diablo # list the different outputs

The confusion table compares the real subtypes with the predicted subtypes for a 2 component model, for the distance of interest:

confusion.mat <- get.confusion_matrix(
                truth = breast.TCGA$data.test$subtype, 
                predicted = Mypredict.diablo$MajorityVote$centroids.dist[,2])
kable(confusion.mat)
get.BER(confusion.mat)

Tuning parameters

For DIABLO, the parameters to tune are:

1 - The design matrix design indicates which data sets, or blocks should be connected to maximise the covariance between components, and to which extend. A compromise needs to be achieved between maximising the correlation between data sets (design value between 0.5 and 1) and maximising the discrimination with the outcome Y (design value between 0 and 0.5), see [@Sin16] for more details.

2 - The number of components to retain ncomp. The rule of thumb is usually $K-1$ where $K$ is the number of classes, but it is worth testing a few extra components.

3 - The number of variables to select on each component and on each data set in the list keepX.

For item 1, by default all data sets are linked as follows: 

MyResult.diablo$design

The design can be changed as follows. By default each data set will be linked to the Y outcome.

MyDesign <- matrix(c(0, 0.1, 0.3,
                     0.1, 0, 0.9,
                     0.3, 0.9, 0),
                   byrow=TRUE,
                   ncol = length(X), nrow = length(X),
                 dimnames = list(names(X), names(X)))
MyDesign
MyResult.diablo.design <- block.splsda(X, Y, keepX=list.keepX, design=MyDesign)

Items 2 and 3 can be tuned using repeated cross-validation, as we described in Chapter \@ref(plsda). A detailed tutorial is provided on our website in the different DIABLO tabs.

Additional resources

Additional examples are provided in example(block.splsda) and in our DIABLO tab in http://www.mixomics.org. Also have a look at [@Sin16]

FAQ

Session Information

sessionInfo()

References



Try the mixOmics package in your browser

Any scripts or data that you put into this service are public.

mixOmics documentation built on April 15, 2021, 6:01 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close