参考信息：
Sparse PLS-DA的简单实现：http://www.bioconductor.org/packages/release/bioc/vignettes/mixOmics/inst/doc/vignette.html#principle-of-sparse-pls-da
Sparse PLS-DA分析的参数选择：http://www.bioconductor.org/packages/release/bioc/vignettes/mixOmics/inst/doc/vignette.html#tuning:sPLSDA
Sparse PLS-DA以上信息的中文翻译信息：https://www.jianshu.com/p/52eae84670c1/
本文的重点参考信息：
http://mixomics.org/case-studies/splsda-srbct/
正文：
数据：小圆形蓝细胞肿瘤 (SRBCT) 数据集包括在 63 个样本上测量的 2,308 个基因的表达水平。样品分为以下四类：8 伯基特淋巴瘤 (BL)、23 尤文肉瘤 (EWS)、12 神经母细胞瘤 (NB) 和 20 横纹肌肉瘤 (RMS)。
srbct 数据集包含以下内容：
$gene：一个有 63 行 2308 列的数据框。63 名受试者 2,308 个基因的表达水平。
$class：包含每个个体的肿瘤类的类向量（共 4 个类）。
$gene.name：一个包含 2,308 行和 2 列的数据框，其中包含有关基因的更多信息。
library(knitr)
library(mixOmics)
knitr::opts_chunk$set(dpi = 100, echo= TRUE, warning=FALSE, message=FALSE, fig.align = 'center', fig.show=TRUE, fig.keep = 'all', out.width = '50%')
data(srbct)
X = srbct$gene #the gene expression datadim(X)
Y = srbct$class
一、主成分分析
可以先进行非监督的主成分分析，查看在非监督的情况下，不同肿瘤样本之间的是否可以基于当前的基因表达信息区分开。
pca.srbct = pca(X, ncomp = 10, center = TRUE, scale = TRUE)
plotIndiv(pca.srbct, group = srbct$class, ind.names = FALSE, legend = TRUE, title = 'PCA on SRBCT')
二、PLS-DA分析
进行PLS-DA分析，初步判定利用多少个主成分comp（ncomp：the number of components to include in the model）可以较好的进行分类。PLS-DA配备了10个comp，用于评估获得能够得到最佳判定效果的comp个数。
对于判别分析，需要设置表示每个样本的分组因子 Y。在 PLS-DA 过程中，Y 因子将被转换为虚拟矩阵。
srbct.plsda <- plsda(X, Y, ncomp = 10) # set ncomp to 10 for performance assessment
plotIndiv(srbct.plsda , comp = 1:2,group = srbct$class, ind.names = FALSE, ellipse = TRUE, legend = TRUE, title = 'PLSDA on SRBCT')
与无监督的 PCA 样本图相比，观察到四种肿瘤类别的清晰分离。绘制每个类别的置信椭圆以突出区分的强度（置信水平默认设置为 95%）。
对PLS-DA 模型的分类性能进行评估
PLS-DA 模型的分类性能通过perf函数使用重复 10 次的 5 倍交叉验证进行评估。重复次数对于确保分类错误率的良好估计是必要的（因为 CV -folds是以随机方式确定的）。根据性能结果，我们可以决定最终 PLS 模型的ncomp值。
set.seed(2543) # for reproducibility, only when the `cpus' argument is not used
perf.plsda.srbct <- perf(srbct.plsda, validation = "Mfold", folds = 5, progressBar = FALSE, auc = TRUE, nrepeat = 10)
plot(perf.plsda.srbct, col = color.mixo(5:7), sd = TRUE, legend.position = "horizontal")
从性能图中，我们观察到 overall error rate和the Balanced Error Rate (BER) 相似，从 1 个comp到 3 个comp过程中急剧下降。在 6 个comp后稳定。
perf.plsda.srbct$choice.ncomp #查看得到最佳结果的comp数
显示的分类性能表明 ncomp = 6 时 BER 和最大距离足以实现良好的判别性能（0.06 错误率，也就意味着并不是必须选择“最佳结果的comp数”，在分类性能良好的基础上，可以选择使用最少的comp）。
auc.plsda = auroc(srbct.plsda, roc.comp = 6)
也可以使用函数auroc获得 AUC 图，其中 AUC 是根据训练集的交叉验证结果计算并取平均值。
三、sPLS-DA 分析
PLS-DA 模型建立在 X 中的所有基因上，其中许多基因在进行分类上可能并未发挥作用。因此可以选择使用sPLS-DA 分析，识别最能区分分类类别的一小部分基因。
sPLS-DA分析最佳参数测定
首先使用函数tune.splsda 估计与模型中所选变量数量相关的分类性能（错误率）。根据 PLS-DA 性能评估的建议设置最大值ncomp = 6。选择重复 10 次的5 折交叉验证（folds = 5），并指定预测距离为最大距离 。
list.keepX <- c(1:10, seq(20, 300, 10)) #每个组件进行评估时使用的基因个数
tune.splsda.srbct <- tune.splsda(X, Y, ncomp = 6, validation = 'Mfold', folds = 5,progressBar = TRUE, dist = 'max.dist', measure = "BER",test.keepX = list.keepX, nrepeat = 10, cpus = 2)
error <- tune.splsda.srbct$error.rate # 指定的keepX下每个comp的错误率
ncomp <- tune.splsda.srbct$choice.ncomp$ncomp # 基于t检验找到最佳的分析comp个数
要在每个comp上选择的特征值数量在$choice.keepX 中输出：
select.keepX <- tune.splsda.srbct$choice.keepX[1:ncomp] # optimal number of variables to select
select.keepX
对于函数tune.splsda中指定的所有comp，以最后一个comp为条件的每个comp的分类错误率如下所示。
plot(tune.splsda.srbct, col = color.jet(6))
当 sPLS-DA 中包含更多comp时，分类错误率会降低（预测精度越低越好）。3 个comp足以让最终的 sPLS-DA 模型达到最佳性能。
最终模型包括 3 个comp和前 3 个comp上的 9、280、30 个选定特征值。
进行 sPLS-DA分析并绘图
splsda.srbct <- splsda(X, Y, ncomp = ncomp, keepX = select.keepX)
plotIndiv(splsda.srbct, comp = c(1,2),group = srbct$class, ind.names = FALSE,ellipse = TRUE, legend = TRUE,title = 'sPLS-DA on SRBCT, comp 1 & 2')
添加第三个comp进一步区分了 NB 和 RMS：
plotIndiv(splsda.srbct, comp = c(1,3),group = srbct$class, ind.names = FALSE,ellipse = TRUE, legend = TRUE,title = 'sPLS-DA on SRBCT, comp 1 & 3')
对于 PLS-DA 分析，也可以使用函数auroc获得 AUC 图。
仅包含 2 个comp的AUROC ：
auc.splsda = auroc(splsda.srbct, roc.comp = 2)
包含最佳个数的comp的AUROC
auc.splsda = auroc(splsda.srbct, roc.comp = ncomp)
PLS-DA分类性能评估
通过指定预测距离算法，使用perf函数评估最终sPLS-DA 模型的分类性能。
set.seed(40) # for reproducibility, only when the `cpus' argument is not used# takes about 1 min to run
perf.srbct <- perf(splsda.srbct, validation = "Mfold", folds = 5,dist = 'max.dist', nrepeat = 10,progressBar = FALSE)
可以观察到，随着模型中添加更多comp，每类的整体和平衡错误率 (BER) 会降低。
perf.srbct$error.rate #可查看添加comp后的错误率变化情况，可对该结果进行呈图
plot(perf.srbct, col = color.mixo(5))
还可以查看每个comp下的变量及其在CV-folds验证过程中的稳定性。
par(mfrow=c(1,3))
plot(perf.srbct$features$stable[[1]], type = 'h', ylab = 'Stability',xlab = 'Features', main = 'Comp 1', las =2)
plot(perf.srbct$features$stable[[2]], type = 'h', ylab = 'Stability',xlab = 'Features', main = 'Comp 2', las =2)
plot(perf.srbct$features$stable[[3]], type = 'h', ylab = 'Stability',xlab = 'Features', main = 'Comp 3', las =2)
par(mfrow=c(1,1))
利用selectVar函数可以输出所选择的变量及其 loading weight value，从最重要到最不重要的顺序。
下面的代码还提取了所选变量的稳定性信息，以辅助读者更好的对变量的性能进行判定
ind.match = match(selectVar(splsda.srbct, comp =1)$name, names(perf.srbct$features$stable[[1]])) #extract the frequency of selection of those selected variables
Freq = as.numeric(perf.srbct$features$stable[[1][ind.match])
data.frame(selectVar(splsda.srbct, comp = 1)$value, Freq)
还有一些其他可输出的图形，具体可详见代码信息：
链接：https://pan.baidu.com/s/1Hs48f_eQ8reYV0d8L4fWuQ
提取码：e7sc
