Nothing
Multi-Block PLS Discriminant Analysis'
In rchemo: Dimension Reduction, Regression and Discrimination for Chemometrics
The Multi-Block PLS Discriminant Analysis method
Metabolomics is a powerful phenotyping tool, generating complex data that need
dedicated treatments to enrich knowledge of biological systems. In particular,
to investigate relations between experimental factors, phenotypes and
metabolism, discriminant statistical analyses are generally performed separately
on metabolomic datasets, complemented by associations with metadata. Another
relevant strategy is to simultaneously analyse thematic data blocks by a
multi-block partial least squares discriminant analysis (MBPLSDA) allowing
determining the importance of variables of the different data blocks in
discriminating groups of observations, taking into account data structure.
References:
-
Wold, S. (Ed.). (1984). Three PLS algorithms according to SW. In Report from
the symposium MULTDAST (multivariate data analysis in science and technology)
(pp. 26–30). Umeå, Sweden
-
Westerhuis, J. A., Kourti, T., & Macgregor, J. F. (1998). Analysis of multi
block and hierarchical PCA and PLS models. Journal of Chemometrics, 12, 301–321.
A selection of reduced non-targeted metabolomics datasets from an article of
Zhang et al. have been used to illustrate the method application.
(Zhang, Y.; Barupal, D.K.; Fan, S.; Gao, B.; Zhu, C.; Flenniken, A.M.;
McKerlie, C.; Nutter, L.M.J.; Lloyd, K.C.K.; Fiehn, O. Sexual Dimorphism of the
Mouse Plasma Metabolome Is Associated with Phenotypes of 30 Gene Knockout Lines.
Metabolites 2023, 13, 947. https://doi.org/10.3390/metabo13080947).
In this study, plasma were analyzed by several metabolomics platforms.
In our example, the 6 samples from 5 mutant groups (Dhfr, Gnpda1, Plk1, Sra1,
Ulk3) and the 40 controls were retained, with the exception of 2 animals that
had missing values in non-targeted metabolomics. Three HILIC-NEG and 2 HILIC-POS
variables were then removed due to missing or infinite values.
Lastly, we kept only 2 data blocks (GCTOF and HILIC POS) to reduce the
computation time.
Always in our example, the objective was to highlight metabolites discriminating
males and females.
R environment preparation
#install.packages("knitr")
library(knitr)
opts_chunk$set(echo = TRUE)
# To ensure reproducibility
set.seed(12)
pkgs <- c("rchemo")
sapply(pkgs, function(x) {
if (!requireNamespace(x, quietly = TRUE)) {
install.packages(x)
}
})
library(rchemo) # to load rchemo
Data preprocessing
A demonstration dataset used for this example is imported.
It contains 2 metabolomics datasets and 1 sample metadata block, with the
groups to be discriminated.
data(Zhang2023, package = "rchemo")
All data blocks should have exactly the same samples (rows).
The block dimension can be checked with the following command:
# Check dimension
BlockNames <- names(Zhang2023)
nbrBlocs <- length(BlockNames)
dims <- lapply(X=Zhang2023[BlockNames], FUN=dim)
names(dims) <- BlockNames
dims
# Remove unuseful object for the next steps
rm(nbrBlocs, dims)
The identical order of samples in the three blocks should be ensured.
# Check rows names in any order
row_names <- lapply(X=Zhang2023[BlockNames], FUN=rownames)
rns <- do.call(cbind, row_names)
rns.unique <- apply(rns, 1, function(x) length(unique(x)))
if (max(rns.unique) > 1) {
stop("Rows names are not identical between blocks.")
}
# Check order of samples
check_row_names <- all(sapply(X=row_names, FUN=identical, y = row_names[[1]]))
if (!check_row_names && max(rns.unique) == 1) {
print("Rows names are not in the same order for all blocks.")
}
# Remove unuseful object for the next steps
rm(row_names, rns, rns.unique, check_row_names)
GCTOF <- Zhang2023$GCTOF
HILICPOS <- Zhang2023$HILICPOS
sample_metadata <- Zhang2023$metadata
rm(Zhang2023)
Data visualization by PCA on separate data blocks
PCA on different datablocks is necessary to verify that there is no outlier,
and that the data scaling is suited.
In our example, unit variance scaling is suited: For both metabolomic dataset,
the scatterplot of loadings is homogeneous, and there are no outliers.
With other datasets, particularly those containing more noise, it would have
been necessary to apply a log transformation with Pareto scaling.
# GCTOF with unit variance scaling
pcaGCTOF <- pcanipalsna(X = scale(GCTOF[,1:dim(GCTOF)[2]]), nlv = nrow(GCTOF),
gs = TRUE,
tol = .Machine$double.eps^0.5, maxit = 200)
## diagram of explained variance
barplot(summary(pcaGCTOF,X = scale(GCTOF[,1:dim(GCTOF)[2]]))$explvar$pvar * 100, names.arg = 1:nrow(GCTOF), main = "diagram of explained variance - PCA GCTOF")
## score plot
plotxy(X= pcaGCTOF$T, group = sample_metadata$Gender,
asp = 0, col = 3:4, alpha.f = .8,
zeroes = TRUE, circle = FALSE, ellipse = FALSE,
labels = FALSE,
legend = TRUE, main = "components - PCA GCTOF", ncol = 1,
pch=16)
## loading plot
plotxy(X= pcaGCTOF$P, group = NULL,
asp = 0, col = NULL, alpha.f = .8,
zeroes = TRUE, circle = FALSE, ellipse = FALSE,
labels = TRUE,
legend = FALSE, main = "loadings - PCA GCTOF", ncol = 1,
cex=0.8)
#HILICPOS with unit variance scaling
pcaHILICPOS <- pcanipalsna(X = scale(HILICPOS[,1:dim(HILICPOS)[2]]), nlv = nrow(HILICPOS),
gs = TRUE,
tol = .Machine$double.eps^0.5, maxit = 200)
## diagram of explained variance
barplot(summary(pcaHILICPOS,X = scale(HILICPOS[,1:dim(HILICPOS)[2]]))$explvar$pvar * 100, names.arg = 1:nrow(HILICPOS), main = "diagram of explained variance - PCA HILIC POS")
## score plot
plotxy(X= pcaHILICPOS$T, group = sample_metadata$Gender,
asp = 0, col = 3:4, alpha.f = .8,
zeroes = TRUE, circle = FALSE, ellipse = FALSE,
labels = FALSE,
legend = TRUE, main = "components - PCA HILIC POS", ncol = 1,
pch=16)
## loading plot
plotxy(X= pcaHILICPOS$P, group = NULL,
asp = 0, col = NULL, alpha.f = .8,
zeroes = TRUE, circle = FALSE, ellipse = FALSE,
labels = TRUE,
legend = FALSE, main = "loadings - PCA HILIC POS", ncol = 1,
cex=0.8)
# Remove unuseful object for the next steps
rm(pcaGCTOF, pcaHILICPOS)
MB-PLS-DA model
In MB-PLS-DA, datablocks are separatly scaled, and weigthed by their Frobenius
norm. Then, they are concatenated before a PLS-DA to predict or explain a
categorical variable.
The main steps to apply this method are:
-
the determination of the optimal number of latent variables of the model,
proposed in the package by cross-validation
-
the model validation, proposed by permutation test
-
the model analysis, based on the scatter plot, the loading plot, the Variable
Importance in Projection
-
the calculation of the predicted categories for each observation.
In rchemo, all these steps can be performed with a single function. The arguments are:
-
Xlist: list of training X-data (n, p).
-
Xnames: names of the X-matrices
-
Xscaling: vector of Xlist length. X variable scaling among "none"
(mean-centering only), "pareto" (mean-centering and pareto scaling), "sd"
(mean-centering and unit variance scaling). If "pareto" or "sd", uncorrected
standard deviation is used.
-
Y: Training Y-data (n, 1)
-
Yscaling: Y variable scaling among "none" (mean-centering only), "pareto"
(mean-centering and pareto scaling), "sd" (mean-centering and unit variance
scaling). If "pareto" or "sd", uncorrected standard deviation is used.
-
weights Weights (n, 1) to apply to the training observations. Internally,
weights are "normalized" to sum to 1. Default to NULL (weights are set to 1/n).
-
newXlist: list of new X-data (m, p) to consider.
-
newXnames: names of the newX-matrices
-
method: method to apply among "mbplsrda","mbplslda","mbplsqda"
-
prior: for mbplslda or mbplsqda models : The prior probabilities of the
classes. Possible values are "unif" (default; probabilities are set equal for
all the classes) or "prop" (probabilities are set equal to the observed
proportions of the classes in y).
-
step: step of the analysis among "nlvtest" (cross-validation to help determine
the optimal number of latent variables), "permutation" (permutation test),
"model" (model calculation),"prediction" (prediction of newX-data or X-data if
any))
-
nlv: number of latent variables to test if step is "nlvtest"; number of latent
variables of the model if step is not "nlvtest".
-
modeloutput: if step is "model": outputs among "scores", "loadings", "coef"
(regression coefficients), "vip" (Variable Importance in Projection; the VIP
calculation being based on the proportion of Y-variance explained by the
components, as proposed by Mehmood et al (2012, 2020).)
-
cvmethod: if step is "nlvtest" or "permutation": "kfolds" for k-folds
cross-validation, or "loo" for leave-one-out.
-
nbrep: if step is "nlvtest" and cvmethod is "kfolds": An integer, setting the
number of CV repetitions. Default value is 30. Must me set to 1 if cvmethod is
"loo"
-
seed: if step is "nlvtest" and cvmethod is "kfolds", or if step is
"permutation: a numeric. Seed used for the repeated resampling
-
samplingk: A vector of length n. The elements are the values of a qualitative
variable used for stratified partition creation. If NULL, the first observation
is set in the first fold, the second observation in the second fold, etc...
-
nfolds: if cvmethod is "kfolds". An integer, setting the number of partitions
to create.Default value is 10.
-
npermut: if step is "permutation": An integer, setting the number of Y-Block
with permutated responses to create. Default value is 30.
-
criterion: if step is "nlvtest" or "permutation" and method is "mbplsrda",
"mbplslda" or "mbplsqda": optimisation criterion among "rmse" and "err" (for
classification error rate)))
-
selection: if step is "nlvtest": a character indicating the selection method
to use to choose the optimal combination of components, among "localmin",
"globalmin","1std". If "localmin": the optimal combination corresponds to the
first local minimum of the mean CV rmse or error rate. If "globalmin" : the
optimal combination corresponds to the minimum mean CV rmse or error rate. If
"1std" (one standard error rule) : it corresponds to the first combination after
which the mean cross-validated rmse or error rate does not decrease significantly.
determination by cross-validation of the optimal number of latent variables of the model
In this example,the mbplslda method with an uniform prior is performed: a linear
discriminant analysis is applied on the PLS scores to obtain the classifications.
In the mbplsrda method, the classification step is performed according to the
predicted value, and consequently, it doesn't take into account the variance of
the scatter plot. And the mbplsqda method is mode suited when the scatter plot
has an atypical shape.
The cvmethod is to choose according to the number of observations. If this
number is small, a leave-one-out cross-validation is more suited, but can lead
to optimistic results. In the other cases, it could be better to apply a k-fold
cross-validation, usually with 3, 5 or 10 folds.
In order to reduce the computational time, in this example, the number of cross-
validation repetition is set to 10. But it must be higher (at least 30).
Concerning the criterion, in a discriminant context, it is more logical to
choose the classification error rate. However, this can result in an error plot
depending on the number of latent variables that is not smooth, so some prefer
to use rmsecv instead.
In the example, the selection parameter is set to "localmin". Normally, the
cross-validated error decreases until it reaches the optimal number of latent
variables, before rising again. Sometimes it fluctuates before this rise, so
choosing “localmin” allows for a more parsimonious model. The “1std” option is
also very parsimonious because it only considers a larger number of latent
variables if doing so significantly improves the cross-validated error.
nlvtestmbplsda <- mbplsr_mbplsda_allsteps(Xlist = list(GCTOF = GCTOF[,1:dim(GCTOF)[2]],
HILICPOS = HILICPOS[,1:dim(HILICPOS)[2]]),
Xnames = c("GCTOF", "HILICPOS"),
Xscaling = c("none","pareto","sd")[3],
Y = sample_metadata[,"Gender", drop=FALSE],
Yscaling = c("none","pareto","sd")[1],
weights = NULL,
newXlist = NULL, newXnames = NULL,
method = c("mbplsrda","mbplslda","mbplsqda")[2],
prior = c("unif", "prop")[1],
step = c("nlvtest","permutation","model","prediction")[1],
nlv = 4,
#modeloutput = c("scores","loadings","coef","vip"),
cvmethod = c("kfolds","loo")[1],
nbrep = 10,
seed = 123,
samplingk = NULL,
nfolds = 3,
#npermut = 30,
criterion = c("err","rmse")[1],
selection = c("localmin","globalmin","1std")[1],
outputfilename = NULL)
nlvtestmbplsda
nlvoptmbplsda <- nlvtestmbplsda[nlvtestmbplsda$optimum==1,"nblv"] # to obtain the optimal number of LV.
# to plot the results of the cross-validation
plot(nlvtestmbplsda$nblv, nlvtestmbplsda$err_mean, xlab = "number of LV", ylab = "CV classification error rate", pch = 16, ylim = c(0,0.6))
segments(nlvtestmbplsda$nblv,nlvtestmbplsda$err_mean-nlvtestmbplsda$err_sd,nlvtestmbplsda$nblv,nlvtestmbplsda$err_mean+nlvtestmbplsda$err_sd)
segments(nlvtestmbplsda$nblv-0.1,nlvtestmbplsda$err_mean-nlvtestmbplsda$err_sd,nlvtestmbplsda$nblv+0.1,nlvtestmbplsda$err_mean-nlvtestmbplsda$err_sd)
segments(nlvtestmbplsda$nblv-0.1,nlvtestmbplsda$err_mean+nlvtestmbplsda$err_sd,nlvtestmbplsda$nblv+0.1,nlvtestmbplsda$err_mean+nlvtestmbplsda$err_sd)
# Remove unuseful object for the next steps
rm(nlvtestmbplsda)
model validation by permutation test
Usually, the cross-validation parameters are the same than for the
determination of the optimal number of latent variables.
In this example, as for the determination of the optimal number of latent
variables, in order to reduce the computational time, the number of cross-
validation repetition and the number of permuted responses are set to 10.
But they must be higher (at least 30).
To be valid, all the cross-validated errors obtained on permuted data must be
lower than the cross-validated error obtained on the original data.
permutmbplsda <- mbplsr_mbplsda_allsteps(Xlist = list(GCTOF = GCTOF[,1:dim(GCTOF)[2]],
HILICPOS = HILICPOS[,1:dim(HILICPOS)[2]]),
Xnames = c("GCTOF", "HILICPOS"),
Xscaling = c("none","pareto","sd")[3],
Y = sample_metadata[,"Gender",drop=FALSE],
Yscaling = c("none","pareto","sd")[1], weights = NULL,
newXlist = NULL, newXnames = NULL,
method = c("mbplsrda","mbplslda","mbplsqda")[2],
prior = c("unif", "prop")[1],
step = c("nlvtest","permutation","model","prediction")[2],
nlv = nlvoptmbplsda,
modeloutput = c("scores","loadings","coef","vip"),
cvmethod = c("kfolds","loo")[1],
nbrep = 10,
seed = 123,
samplingk = NULL,
nfolds = 3,
npermut = 10,
criterion = c("err","rmse")[1],
# selection = c("localmin","globalmin","1std")[1],
import = c("R","ChemFlow","W4M")[1],
outputfilename = NULL)
#plot of the results
plot(permutmbplsda, pch = 16, ylab = "CV classification error rate", xlab = "dyssimilarity Y-Ypermuted")
# Remove unuseful object for the next steps
rm(permutmbplsda)
model analysis, based on the scatter plot, the loading plot, the Variable Importance in Projection
The score plot allows to show the group discrimination.
Based on the VIP curve, the threshold to consider that a variable is important
can be set.
modelmbplsda <- mbplsr_mbplsda_allsteps(Xlist = list(GCTOF = GCTOF[,1:dim(GCTOF)[2]],
HILICPOS = HILICPOS[,1:dim(HILICPOS)[2]]),
Xnames = c("GCTOF", "HILICPOS"),
Xscaling = c("none","pareto","sd")[3],
Y = sample_metadata[,"Gender",drop=FALSE],
Yscaling = c("none","pareto","sd")[1],
weights = NULL,
newXlist = NULL, newXnames = NULL,
method = c("mbplsrda","mbplslda","mbplsqda")[2],
prior = c("unif", "prop")[1],
step = c("nlvtest","permutation","model","prediction")[3],
nlv = nlvoptmbplsda,
modeloutput = c("scores","loadings","coef","vip"),
cvmethod = c("kfolds","loo")[1],
# nbrep = 30,
# seed = 123,
# samplingk = NULL,
# nfolds = 5,
# npermut = 30,
# criterion = c("err","rmse")[1],
# selection = c("localmin","globalmin","1std")[1],
import = c("R","ChemFlow","W4M")[1],
outputfilename = NULL)
# score plot
plotxy(X= modelmbplsda$scores, group = sample_metadata$Gender,
asp = 0, col = 3:4, alpha.f = .8,
zeroes = TRUE, circle = FALSE, ellipse = FALSE,
labels = FALSE,
legend = TRUE, main = "scores - MB PLS DA", ncol = 1,
pch=16)
# loading plot
plotxy(X= modelmbplsda$loadings, group = substr(rownames(modelmbplsda$loadings),1,6),
asp = 0, col = NULL, alpha.f = .8,
zeroes = TRUE, circle = FALSE, ellipse = FALSE,
labels = FALSE,
legend = TRUE, main = "loadings - MB PLS DA", ncol = 1,
cex=0.8, pch = 16)
# VIP curve
plot(modelmbplsda$vip[order(modelmbplsda$vip[,nlvoptmbplsda], decreasing = TRUE),nlvoptmbplsda], pch = 16,cex = 0.8,
col = as.numeric(as.factor(substr(rownames(modelmbplsda$vip[order(modelmbplsda$vip[,nlvoptmbplsda], decreasing = TRUE),nlvoptmbplsda, drop=FALSE]),1,6))), ylab = "VIP value",
main = "VIP curve MB PLS DA")
legend("topright", legend = c("GCTOF", "HILICPOS"), pch = 16, col = 1:2)
# Remove unuseful object for the next steps
rm(modelmbplsda)
predicted values
The prediction step provides a table with:
-
subjects scores
-
predicted categories
-
probabilities to be predicted for each category.
predmbplsda <- mbplsr_mbplsda_allsteps(Xlist = list(GCTOF = GCTOF[,1:dim(GCTOF)[2]],
HILICPOS = HILICPOS[,1:dim(HILICPOS)[2]]),
Xnames = c("GCTOF", "HILICPOS"),
Xscaling = c("none","pareto","sd")[3],
Y = sample_metadata[,"Gender",drop=FALSE],
Yscaling = c("none","pareto","sd")[1],
weights = NULL,
newXlist = NULL, newXnames = NULL,
method = c("mbplsrda","mbplslda","mbplsqda")[2],
prior = c("unif", "prop")[1],
step = c("nlvtest","permutation","model","prediction")[4],
nlv = nlvoptmbplsda,
# modeloutput = c("scores","loadings","coef","vip"),
#
# cvmethod = c("kfolds","loo")[1],
# nbrep = 30,
# seed = 123,
# samplingk = NULL,
# nfolds = 5,
# npermut = 30,
#
# criterion = c("err","rmse")[1],
# selection = c("localmin","globalmin","1std")[1],
import = c("R","ChemFlow","W4M")[1],
outputfilename = NULL)
predmbplsda
# Remove unuseful object for the next steps
rm(predmbplsda, nlvoptmbplsda)
Reproducibility
This vignette was produced with the following R session configuration.
sessionInfo()
Try the rchemo package in your browser
Any scripts or data that you put into this service are public.
rchemo documentation built on June 30, 2026, 5:10 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
The Multi-Block PLS Discriminant Analysis method
Metabolomics is a powerful phenotyping tool, generating complex data that need dedicated treatments to enrich knowledge of biological systems. In particular, to investigate relations between experimental factors, phenotypes and metabolism, discriminant statistical analyses are generally performed separately on metabolomic datasets, complemented by associations with metadata. Another relevant strategy is to simultaneously analyse thematic data blocks by a multi-block partial least squares discriminant analysis (MBPLSDA) allowing determining the importance of variables of the different data blocks in discriminating groups of observations, taking into account data structure.
References:
-
Wold, S. (Ed.). (1984). Three PLS algorithms according to SW. In Report from the symposium MULTDAST (multivariate data analysis in science and technology) (pp. 26–30). Umeå, Sweden
-
Westerhuis, J. A., Kourti, T., & Macgregor, J. F. (1998). Analysis of multi block and hierarchical PCA and PLS models. Journal of Chemometrics, 12, 301–321.
A selection of reduced non-targeted metabolomics datasets from an article of Zhang et al. have been used to illustrate the method application. (Zhang, Y.; Barupal, D.K.; Fan, S.; Gao, B.; Zhu, C.; Flenniken, A.M.; McKerlie, C.; Nutter, L.M.J.; Lloyd, K.C.K.; Fiehn, O. Sexual Dimorphism of the Mouse Plasma Metabolome Is Associated with Phenotypes of 30 Gene Knockout Lines. Metabolites 2023, 13, 947. https://doi.org/10.3390/metabo13080947).
In this study, plasma were analyzed by several metabolomics platforms.
In our example, the 6 samples from 5 mutant groups (Dhfr, Gnpda1, Plk1, Sra1, Ulk3) and the 40 controls were retained, with the exception of 2 animals that had missing values in non-targeted metabolomics. Three HILIC-NEG and 2 HILIC-POS variables were then removed due to missing or infinite values. Lastly, we kept only 2 data blocks (GCTOF and HILIC POS) to reduce the computation time. Always in our example, the objective was to highlight metabolites discriminating males and females.
R environment preparation
#install.packages("knitr") library(knitr) opts_chunk$set(echo = TRUE) # To ensure reproducibility set.seed(12)
pkgs <- c("rchemo") sapply(pkgs, function(x) { if (!requireNamespace(x, quietly = TRUE)) { install.packages(x) } })
library(rchemo) # to load rchemo
Data preprocessing
A demonstration dataset used for this example is imported. It contains 2 metabolomics datasets and 1 sample metadata block, with the groups to be discriminated.
data(Zhang2023, package = "rchemo")
All data blocks should have exactly the same samples (rows). The block dimension can be checked with the following command:
# Check dimension BlockNames <- names(Zhang2023) nbrBlocs <- length(BlockNames) dims <- lapply(X=Zhang2023[BlockNames], FUN=dim) names(dims) <- BlockNames dims # Remove unuseful object for the next steps rm(nbrBlocs, dims)
The identical order of samples in the three blocks should be ensured.
# Check rows names in any order row_names <- lapply(X=Zhang2023[BlockNames], FUN=rownames) rns <- do.call(cbind, row_names) rns.unique <- apply(rns, 1, function(x) length(unique(x))) if (max(rns.unique) > 1) { stop("Rows names are not identical between blocks.") } # Check order of samples check_row_names <- all(sapply(X=row_names, FUN=identical, y = row_names[[1]])) if (!check_row_names && max(rns.unique) == 1) { print("Rows names are not in the same order for all blocks.") } # Remove unuseful object for the next steps rm(row_names, rns, rns.unique, check_row_names)
GCTOF <- Zhang2023$GCTOF HILICPOS <- Zhang2023$HILICPOS sample_metadata <- Zhang2023$metadata rm(Zhang2023)
Data visualization by PCA on separate data blocks
PCA on different datablocks is necessary to verify that there is no outlier, and that the data scaling is suited. In our example, unit variance scaling is suited: For both metabolomic dataset, the scatterplot of loadings is homogeneous, and there are no outliers. With other datasets, particularly those containing more noise, it would have been necessary to apply a log transformation with Pareto scaling.
# GCTOF with unit variance scaling pcaGCTOF <- pcanipalsna(X = scale(GCTOF[,1:dim(GCTOF)[2]]), nlv = nrow(GCTOF), gs = TRUE, tol = .Machine$double.eps^0.5, maxit = 200) ## diagram of explained variance barplot(summary(pcaGCTOF,X = scale(GCTOF[,1:dim(GCTOF)[2]]))$explvar$pvar * 100, names.arg = 1:nrow(GCTOF), main = "diagram of explained variance - PCA GCTOF") ## score plot plotxy(X= pcaGCTOF$T, group = sample_metadata$Gender, asp = 0, col = 3:4, alpha.f = .8, zeroes = TRUE, circle = FALSE, ellipse = FALSE, labels = FALSE, legend = TRUE, main = "components - PCA GCTOF", ncol = 1, pch=16) ## loading plot plotxy(X= pcaGCTOF$P, group = NULL, asp = 0, col = NULL, alpha.f = .8, zeroes = TRUE, circle = FALSE, ellipse = FALSE, labels = TRUE, legend = FALSE, main = "loadings - PCA GCTOF", ncol = 1, cex=0.8) #HILICPOS with unit variance scaling pcaHILICPOS <- pcanipalsna(X = scale(HILICPOS[,1:dim(HILICPOS)[2]]), nlv = nrow(HILICPOS), gs = TRUE, tol = .Machine$double.eps^0.5, maxit = 200) ## diagram of explained variance barplot(summary(pcaHILICPOS,X = scale(HILICPOS[,1:dim(HILICPOS)[2]]))$explvar$pvar * 100, names.arg = 1:nrow(HILICPOS), main = "diagram of explained variance - PCA HILIC POS") ## score plot plotxy(X= pcaHILICPOS$T, group = sample_metadata$Gender, asp = 0, col = 3:4, alpha.f = .8, zeroes = TRUE, circle = FALSE, ellipse = FALSE, labels = FALSE, legend = TRUE, main = "components - PCA HILIC POS", ncol = 1, pch=16) ## loading plot plotxy(X= pcaHILICPOS$P, group = NULL, asp = 0, col = NULL, alpha.f = .8, zeroes = TRUE, circle = FALSE, ellipse = FALSE, labels = TRUE, legend = FALSE, main = "loadings - PCA HILIC POS", ncol = 1, cex=0.8) # Remove unuseful object for the next steps rm(pcaGCTOF, pcaHILICPOS)
MB-PLS-DA model
In MB-PLS-DA, datablocks are separatly scaled, and weigthed by their Frobenius norm. Then, they are concatenated before a PLS-DA to predict or explain a categorical variable. The main steps to apply this method are:
-
the determination of the optimal number of latent variables of the model, proposed in the package by cross-validation
-
the model validation, proposed by permutation test
-
the model analysis, based on the scatter plot, the loading plot, the Variable Importance in Projection
-
the calculation of the predicted categories for each observation.
In rchemo, all these steps can be performed with a single function. The arguments are:
-
Xlist: list of training X-data (n, p).
-
Xnames: names of the X-matrices
-
Xscaling: vector of Xlist length. X variable scaling among "none" (mean-centering only), "pareto" (mean-centering and pareto scaling), "sd" (mean-centering and unit variance scaling). If "pareto" or "sd", uncorrected standard deviation is used.
-
Y: Training Y-data (n, 1)
-
Yscaling: Y variable scaling among "none" (mean-centering only), "pareto" (mean-centering and pareto scaling), "sd" (mean-centering and unit variance scaling). If "pareto" or "sd", uncorrected standard deviation is used.
-
weights Weights (n, 1) to apply to the training observations. Internally, weights are "normalized" to sum to 1. Default to NULL (weights are set to 1/n).
-
newXlist: list of new X-data (m, p) to consider.
-
newXnames: names of the newX-matrices
-
method: method to apply among "mbplsrda","mbplslda","mbplsqda"
-
prior: for mbplslda or mbplsqda models : The prior probabilities of the classes. Possible values are "unif" (default; probabilities are set equal for all the classes) or "prop" (probabilities are set equal to the observed proportions of the classes in y).
-
step: step of the analysis among "nlvtest" (cross-validation to help determine the optimal number of latent variables), "permutation" (permutation test), "model" (model calculation),"prediction" (prediction of newX-data or X-data if any))
-
nlv: number of latent variables to test if step is "nlvtest"; number of latent variables of the model if step is not "nlvtest".
-
modeloutput: if step is "model": outputs among "scores", "loadings", "coef" (regression coefficients), "vip" (Variable Importance in Projection; the VIP calculation being based on the proportion of Y-variance explained by the components, as proposed by Mehmood et al (2012, 2020).)
-
cvmethod: if step is "nlvtest" or "permutation": "kfolds" for k-folds cross-validation, or "loo" for leave-one-out.
-
nbrep: if step is "nlvtest" and cvmethod is "kfolds": An integer, setting the number of CV repetitions. Default value is 30. Must me set to 1 if cvmethod is "loo"
-
seed: if step is "nlvtest" and cvmethod is "kfolds", or if step is "permutation: a numeric. Seed used for the repeated resampling
-
samplingk: A vector of length n. The elements are the values of a qualitative variable used for stratified partition creation. If NULL, the first observation is set in the first fold, the second observation in the second fold, etc...
-
nfolds: if cvmethod is "kfolds". An integer, setting the number of partitions to create.Default value is 10.
-
npermut: if step is "permutation": An integer, setting the number of Y-Block with permutated responses to create. Default value is 30.
-
criterion: if step is "nlvtest" or "permutation" and method is "mbplsrda", "mbplslda" or "mbplsqda": optimisation criterion among "rmse" and "err" (for classification error rate)))
-
selection: if step is "nlvtest": a character indicating the selection method to use to choose the optimal combination of components, among "localmin", "globalmin","1std". If "localmin": the optimal combination corresponds to the first local minimum of the mean CV rmse or error rate. If "globalmin" : the optimal combination corresponds to the minimum mean CV rmse or error rate. If "1std" (one standard error rule) : it corresponds to the first combination after which the mean cross-validated rmse or error rate does not decrease significantly.
determination by cross-validation of the optimal number of latent variables of the model
In this example,the mbplslda method with an uniform prior is performed: a linear discriminant analysis is applied on the PLS scores to obtain the classifications.
In the mbplsrda method, the classification step is performed according to the predicted value, and consequently, it doesn't take into account the variance of the scatter plot. And the mbplsqda method is mode suited when the scatter plot has an atypical shape.
The cvmethod is to choose according to the number of observations. If this number is small, a leave-one-out cross-validation is more suited, but can lead to optimistic results. In the other cases, it could be better to apply a k-fold cross-validation, usually with 3, 5 or 10 folds. In order to reduce the computational time, in this example, the number of cross- validation repetition is set to 10. But it must be higher (at least 30).
Concerning the criterion, in a discriminant context, it is more logical to choose the classification error rate. However, this can result in an error plot depending on the number of latent variables that is not smooth, so some prefer to use rmsecv instead.
In the example, the selection parameter is set to "localmin". Normally, the cross-validated error decreases until it reaches the optimal number of latent variables, before rising again. Sometimes it fluctuates before this rise, so choosing “localmin” allows for a more parsimonious model. The “1std” option is also very parsimonious because it only considers a larger number of latent variables if doing so significantly improves the cross-validated error.
nlvtestmbplsda <- mbplsr_mbplsda_allsteps(Xlist = list(GCTOF = GCTOF[,1:dim(GCTOF)[2]], HILICPOS = HILICPOS[,1:dim(HILICPOS)[2]]), Xnames = c("GCTOF", "HILICPOS"), Xscaling = c("none","pareto","sd")[3], Y = sample_metadata[,"Gender", drop=FALSE], Yscaling = c("none","pareto","sd")[1], weights = NULL, newXlist = NULL, newXnames = NULL, method = c("mbplsrda","mbplslda","mbplsqda")[2], prior = c("unif", "prop")[1], step = c("nlvtest","permutation","model","prediction")[1], nlv = 4, #modeloutput = c("scores","loadings","coef","vip"), cvmethod = c("kfolds","loo")[1], nbrep = 10, seed = 123, samplingk = NULL, nfolds = 3, #npermut = 30, criterion = c("err","rmse")[1], selection = c("localmin","globalmin","1std")[1], outputfilename = NULL) nlvtestmbplsda nlvoptmbplsda <- nlvtestmbplsda[nlvtestmbplsda$optimum==1,"nblv"] # to obtain the optimal number of LV. # to plot the results of the cross-validation plot(nlvtestmbplsda$nblv, nlvtestmbplsda$err_mean, xlab = "number of LV", ylab = "CV classification error rate", pch = 16, ylim = c(0,0.6)) segments(nlvtestmbplsda$nblv,nlvtestmbplsda$err_mean-nlvtestmbplsda$err_sd,nlvtestmbplsda$nblv,nlvtestmbplsda$err_mean+nlvtestmbplsda$err_sd) segments(nlvtestmbplsda$nblv-0.1,nlvtestmbplsda$err_mean-nlvtestmbplsda$err_sd,nlvtestmbplsda$nblv+0.1,nlvtestmbplsda$err_mean-nlvtestmbplsda$err_sd) segments(nlvtestmbplsda$nblv-0.1,nlvtestmbplsda$err_mean+nlvtestmbplsda$err_sd,nlvtestmbplsda$nblv+0.1,nlvtestmbplsda$err_mean+nlvtestmbplsda$err_sd) # Remove unuseful object for the next steps rm(nlvtestmbplsda)
model validation by permutation test
Usually, the cross-validation parameters are the same than for the determination of the optimal number of latent variables. In this example, as for the determination of the optimal number of latent variables, in order to reduce the computational time, the number of cross- validation repetition and the number of permuted responses are set to 10. But they must be higher (at least 30). To be valid, all the cross-validated errors obtained on permuted data must be lower than the cross-validated error obtained on the original data.
permutmbplsda <- mbplsr_mbplsda_allsteps(Xlist = list(GCTOF = GCTOF[,1:dim(GCTOF)[2]], HILICPOS = HILICPOS[,1:dim(HILICPOS)[2]]), Xnames = c("GCTOF", "HILICPOS"), Xscaling = c("none","pareto","sd")[3], Y = sample_metadata[,"Gender",drop=FALSE], Yscaling = c("none","pareto","sd")[1], weights = NULL, newXlist = NULL, newXnames = NULL, method = c("mbplsrda","mbplslda","mbplsqda")[2], prior = c("unif", "prop")[1], step = c("nlvtest","permutation","model","prediction")[2], nlv = nlvoptmbplsda, modeloutput = c("scores","loadings","coef","vip"), cvmethod = c("kfolds","loo")[1], nbrep = 10, seed = 123, samplingk = NULL, nfolds = 3, npermut = 10, criterion = c("err","rmse")[1], # selection = c("localmin","globalmin","1std")[1], import = c("R","ChemFlow","W4M")[1], outputfilename = NULL) #plot of the results plot(permutmbplsda, pch = 16, ylab = "CV classification error rate", xlab = "dyssimilarity Y-Ypermuted") # Remove unuseful object for the next steps rm(permutmbplsda)
model analysis, based on the scatter plot, the loading plot, the Variable Importance in Projection
The score plot allows to show the group discrimination. Based on the VIP curve, the threshold to consider that a variable is important can be set.
modelmbplsda <- mbplsr_mbplsda_allsteps(Xlist = list(GCTOF = GCTOF[,1:dim(GCTOF)[2]], HILICPOS = HILICPOS[,1:dim(HILICPOS)[2]]), Xnames = c("GCTOF", "HILICPOS"), Xscaling = c("none","pareto","sd")[3], Y = sample_metadata[,"Gender",drop=FALSE], Yscaling = c("none","pareto","sd")[1], weights = NULL, newXlist = NULL, newXnames = NULL, method = c("mbplsrda","mbplslda","mbplsqda")[2], prior = c("unif", "prop")[1], step = c("nlvtest","permutation","model","prediction")[3], nlv = nlvoptmbplsda, modeloutput = c("scores","loadings","coef","vip"), cvmethod = c("kfolds","loo")[1], # nbrep = 30, # seed = 123, # samplingk = NULL, # nfolds = 5, # npermut = 30, # criterion = c("err","rmse")[1], # selection = c("localmin","globalmin","1std")[1], import = c("R","ChemFlow","W4M")[1], outputfilename = NULL) # score plot plotxy(X= modelmbplsda$scores, group = sample_metadata$Gender, asp = 0, col = 3:4, alpha.f = .8, zeroes = TRUE, circle = FALSE, ellipse = FALSE, labels = FALSE, legend = TRUE, main = "scores - MB PLS DA", ncol = 1, pch=16) # loading plot plotxy(X= modelmbplsda$loadings, group = substr(rownames(modelmbplsda$loadings),1,6), asp = 0, col = NULL, alpha.f = .8, zeroes = TRUE, circle = FALSE, ellipse = FALSE, labels = FALSE, legend = TRUE, main = "loadings - MB PLS DA", ncol = 1, cex=0.8, pch = 16) # VIP curve plot(modelmbplsda$vip[order(modelmbplsda$vip[,nlvoptmbplsda], decreasing = TRUE),nlvoptmbplsda], pch = 16,cex = 0.8, col = as.numeric(as.factor(substr(rownames(modelmbplsda$vip[order(modelmbplsda$vip[,nlvoptmbplsda], decreasing = TRUE),nlvoptmbplsda, drop=FALSE]),1,6))), ylab = "VIP value", main = "VIP curve MB PLS DA") legend("topright", legend = c("GCTOF", "HILICPOS"), pch = 16, col = 1:2) # Remove unuseful object for the next steps rm(modelmbplsda)
predicted values
The prediction step provides a table with:
-
subjects scores
-
predicted categories
-
probabilities to be predicted for each category.
predmbplsda <- mbplsr_mbplsda_allsteps(Xlist = list(GCTOF = GCTOF[,1:dim(GCTOF)[2]], HILICPOS = HILICPOS[,1:dim(HILICPOS)[2]]), Xnames = c("GCTOF", "HILICPOS"), Xscaling = c("none","pareto","sd")[3], Y = sample_metadata[,"Gender",drop=FALSE], Yscaling = c("none","pareto","sd")[1], weights = NULL, newXlist = NULL, newXnames = NULL, method = c("mbplsrda","mbplslda","mbplsqda")[2], prior = c("unif", "prop")[1], step = c("nlvtest","permutation","model","prediction")[4], nlv = nlvoptmbplsda, # modeloutput = c("scores","loadings","coef","vip"), # # cvmethod = c("kfolds","loo")[1], # nbrep = 30, # seed = 123, # samplingk = NULL, # nfolds = 5, # npermut = 30, # # criterion = c("err","rmse")[1], # selection = c("localmin","globalmin","1std")[1], import = c("R","ChemFlow","W4M")[1], outputfilename = NULL) predmbplsda # Remove unuseful object for the next steps rm(predmbplsda, nlvoptmbplsda)
Reproducibility
This vignette was produced with the following R session configuration.
sessionInfo()
Try the rchemo package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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.