在现代微生物组学分析中，高通量的测试方法使得研究者可以一次性获取大量的数据信息，这时候所获得的数据里可能存在大量“冗余”；此外，在实际操作中，研究人员为避免遗漏重要的系统特征，往往倾向于较周到的选取测试指标，这些变量之间也很可能存在多重共线性。因此，在大数据量的多个数据集之间进行分析时，常常难以有效的进行数据挖掘。
基本原理
为了消除冗余数据，选取系统的主要特征，可以使用排序方法进行降维并挑选主要变化因子（应用到生态中就是挑选主要物种或环境因子）。然而约束排序只能使用一个解释变量数据集对一个响应变量数据集进行分析，当有多个数据集时，就需要多种方法结合进行分析。例如大豆根瘤微生物群落、土壤理化性质、大豆种子质量（含油量、粒重、数目等）这三个数据集，我们倾向于用土壤理化因子去解释另外两者，而对于小鼠肠道微生物、食物参数（蛋白质、纤维、油脂含量）、身体状况（体重、血糖等）这三个数据集，我们倾向于用食物参数和肠道微生物去解释身体状况。
当研究认定多个数据集有依次的解释关系时，可以使用连续的解释变量-响应变量模式进行分析，这里介绍一种约束排序-PLS回归模型分析方法。在OLS回归分析中，多重共线性会造成模型的回归系数置信区间过分夸大，造成模型的解释能力大大降低，这时可以采用偏最小二乘（Partialleast squares，PLS）回归的方法。偏最小方差回归是在对因变量y进行拟合的同时使用y对自变量进行X进行约束降维，使用X的主成分来对y进行拟合，具体方法如下所示。
假设有以下自变量X和因变量y：
第一步，计算y与X的协方差向量w1（若因变量是多变量矩阵Y则为协方差矩阵W，这里为简化计算以一元因变量y代替）：
根据协方差向量提取X的主成分t1：
第二步，分别构建y、x1、x2…xp的关于t1的回归：
第三步，使用残差向量y(1)、xj(1)代替前面的y、X重复第一步、第二步，依此类推不断迭代，直到预测平方和（PRESS）与残差平方和（SS）之比较大（或不再增长）为止，也即要选择尽可能小的主成分个数A，来使得交差：
尽可能的小。接下来以实例进行分析。
冗余分析
现有三组处理的小鼠分别为正常食物饮食（NCD）、高脂肪酸饮食（HFD）、牛磺熊去氧胆酸（TUDCA），我们以这个因子变量对其肠道微生物群落进行约束排序也即RDA分析，筛选受不同处理影响较大的物种，来预测小鼠生理特征也即表型（phenotype）数据。首先读取并处理数据，如下所示：
代码语言：javascript复制#读取处理、表型、微生物群落信息
treatment=read.table(file="1treatment.txt", header=TRUE, check.names=FALSE)
rownames(treatment)=treatment[,1]
tret=treatment
phenotype=read.table(file="2phenotype.txt",header=TRUE,check.names=FALSE)
rownames(phenotype)=phenotype[,1]
pheno=as.matrix(phenotype[,-1])
pheno=pheno[rownames(tret),]
species=read.table(file="3otu_table.txt",header=TRUE,check.names=FALSE)
rownames(species)=species[,1]
speci=as.matrix(species[,-1])
#计算物种相对丰度
library(vegan)
speci=decostand(speci, MARGIN=2, "total")*100
speci=t(speci)
speci=speci[rownames(tret),]因为PLS模型为线性模型，因此我们首先选用同样基于线性模型的RDA分析筛选受不同处理影响大的关键物种（keyotu），具体方法如下：
代码语言：javascript复制#RDA分析
library("maptools")
rda=rda(speci~Diet_, data=tret, scale=FALSE)
#这里的data为数据框格式
plot(rda, display=c("sp", "si", "cn"), scaling=3)
#scaling=1关注物种之间的关系，scaling=2关注样本之间的关系，scaling=3 关注样本与物种之间的关系
#“cn”为小组样品形心初步作图如下所示：
接下来我们进一步进行显著性检验，使用蒙特卡罗置换检验：
代码语言：javascript复制#蒙特卡罗置换检验
permutest(rda, permu=999)检验结果如下：
检验结果p=0.001，说明不同处理的小鼠肠道微生物变化很显著。
接下来提取RDA分析结果，并筛选主坐标RDA1和RDA2解释量最大的100个otu，也即根据主坐标得分以及其解释量筛选物种：
代码语言：javascript复制#提取RDA分析结果并筛选主坐标得分高的OTU
rda_sum=summary(rda, scaling=2)
#scaling=2表示保留两个主约束坐标
sp=as.matrix(rda_sum$species[,1:2])
si=as.matrix(rda_sum$sites[,1:2])
cn=as.matrix(rda_sum$centroids[,1:2])
rdap=as.vector(rda_sum$cont$importance[2,1:2], mode="any")
rdapp=paste(rdap*100, "%",sep="")
#计算物种加权得分并排序
tscore=numeric(length(rownames(sp)))
for (i in 1:length(tscore)) {
tscore[i]=abs(rdap[1]*sp[i,1]) abs(rdap[2]*sp[i,2])
}
spt=cbind(sp, tscore)
spt=spt[rev(order(spt[,3])),]
sph=spt[1:100,]根据物种在两个约束排序轴的得分（即系数绝对值）以及两个约束排序轴的方差解释量计算加权得分，从而筛选出受不同饮食处理影响较大的物种，以便最后进行PLS回归分析，最终的图形绘制如下所示：
代码语言：javascript复制#最终绘图
group=as.vector(tret$Diet_)
sicolor=character(length(group))
sicolor[which(group=="HFD")]="red"
sicolor[which(group=="NCD")]="blue"
sicolor[which(group=="TUDCA")]="green4"
cncolor=c("red", "blue", "green4")
rdap1=paste("RDA1(", rdapp[1], ")", sep="")
rdap2=paste("RDA2(", rdapp[2], ")", sep="")
limx=max(abs(si[,1]/10))
limy=max(abs(si[,2]/10))
plot(si/10, pch=19, col=sicolor, ylim=c(-limy 0.2, limy), xlim=c(-limx, limx),
xlab=rdap1, ylab=rdap2, main="RDA analysis", cex=1, las=1)
pointLabel(x=cn[,1]/10, y=cn[,2]/10, labels=rownames(cn), cex=1.2, col=cncolor)
abline(h=0, lty=2)
abline(v=0, lty=2)
arrows(0, 0, sph[,1]/4, sph[,2]/4, col="black", angle=20, length=0.1)
text(sph[1:20,1]/4 0.01, sph[1:20,2]/4 0.04, rownames(sph)[1:20],col="black", cex=0.6,adj=0.2,pos=1)绘图结果如下所示：
PLS回归
接下来我们关注的是肠道微生物与宿主也即小鼠生理状况的关系，也即通过其肠道微生物主要物种的丰度变化能否预测小鼠的生理状况，能较好的预测那些生理参数，这时物种为自变量，生理参数为响应变量。首先我们需要提取通过RDA分析筛选的100个物种及其丰度，然后与小鼠生理数据构建PLS回归预测模型，具体如下：
代码语言：javascript复制#PLS回归分析
#提取100个物种丰度
sph=t(sph)
specih=as.matrix(speci[,colnames(sph)])
#将生理状况数据标准化
pheno=decostand(pheno, method='standardize')
library(pls)
rdapls=plsr(pheno~specih, validation="LOO", jackknife=TRUE)
summary(rdapls, what="all")其中jackknife=TRUE表示对系数使用刀切法交叉验证，初步模型拟合结果如下：
结果中给出了每个响应变量交叉验证后的均方根误差（Root-mean-squareerror of cross-validation，RMSECV）以及方差解释量。接下来要选择合适的主成分个数，对CV验证的可以进行作图展示（篇幅限制这里只展示9个）：
代码语言：javascript复制#对交叉验证结果进行作图
validationplot(rdapls)
我们需要在选择的成分个数尽可能小的前提下，选择使均方根误差最小或几乎不变以及training方差解释量尽可能大的主成分个数。由上面结果我们选择ncomp=3，继续进行优化建模：
代码语言：javascript复制#根据初步拟合结果选取最佳主成分数目
rdapls2=plsr(pheno~specih, ncomp=3, validation="LOO", jackknife=TRUE)
summary(rdapls2,what="all")
#对最终预测模型进行作图
predplot(rdapls2, line=TRUE, validation="CV")由结果可以看出Body_weight、ALT拟合效果最好，预测准确度最高。最后我们提取分析结果：
代码语言：javascript复制#提取最终的RMSECV
rmse=RMSEP(rdapls2)
rmsecv=rmse$val
#提取相关系数R
r=R2(rdapls2)
#提取回归系数
coef=coef(rdapls2)
#提取检验结果
test=jack.test(rdapls2)
#解释量
summary(rdapls2, what="training")三个主成分的最终解释量为：
END
