三维荧光平行因子学习记录–（三）使用DOMfluor工具箱进行平行因子分析–（二）

注 本文仅作为自己的学习记录以备以后复习查阅

本文的参考文献：
本文的参考博客链接：

https://zhuanlan.zhihu.com/p/407011475

大家对三维荧光数据分析有疑问的或者对平行因子分析感兴趣的可以去看看他的主页。

话不多说，上一篇我们用工具箱自带的数据集对平行因子分析的过程进行了初步的了解，这篇文章将记录如何对自己的数据集进行平行因子分析，数据集的制作参考我的另一篇文章，这里默认大家已经都做好了自己的数据集，这里在补充一下上一篇文章我没有对工具箱自带的数据集做太多的介绍，事实上自带的数据都是已经扣除了空白、瑞丽和拉曼散射，所以我们在做分析之前自己的数据也要扣除空白、瑞丽和拉曼散射，这样才能得到较好的分析效果，本文所用的数据集是采自本校理工湖土壤的DOM溶液，由于样本数不够严格来讲不具备分析的条件，但目前因为疫情也不能去取样，为了学习我将现有的样本复制了几份以达到应有的样本量，所以数据分析结果不具有参考价值。

三、分析自己的数据集

数据集基本情况：

第一步：绘制等高线图

现在我们将绘制EEM，以检查数据是否正确加载，我们输入：

PlotEEMby1(1:5,OriginalData,'R.U.')

将数据集中的前5个EEM一次绘制一个等高线图。按下键盘上的任何键都可以查看下一个图形。如果将1:5替换为1:65，则所有65个样本将一次绘制一个。通过在键盘上输入[Ctrl+C]，然后关闭图形窗口，可以随时停止打印过程 ：
注：‘R.U.’，表示你数据的荧光强度的单位，没有经过校正，即任意单位时，为’A.U.‘，经过拉曼校正后，为’R.U.’，经过硫酸奎宁校正后，为’QSE’。

②将一次绘制四个数据。这在处理大型数据集时更快。

PlotEEMby4(1,OriginalData,'R.U.')

③与上述操作相同，但使用从最小和最大测量数据自动导出的固定z轴（颜色条比例）绘制数据。

PlotEEMby4FixZ(1,OriginalData,'R.U.')

还可以尝试键入help PlotEEMby1、help PlotEEMby4，最后键入help PlotEEMby4FixZ。有关如何使用这些功能的一些说明会打印到命令窗口中。这适用于MATLAB中的所有函数（例如，尝试键入help load）。

也可以使用PlotSurfby1和PlotSurfby4功能以类似的方式绘制曲面图。

例如：

PlotSurfby1(1,OriginalData,'R.U') 

第二步：去除散射

这一步将创建一个新的数据副本，其中受散射峰影响的波长已被切割，并替换为缺失值或零（如不手动停止会自动进行到最后一个样本）。输入：

[CutData] = EEMCut(OriginalData,20,20,NaN,NaN,'No');

数据将被切割，然后绘制，以便比较切割前后的EEM。图表将从第一个样本自动绘制到最后一个样本。该函数删除无荧光区域（发射波长小于激发波长）和受一阶散射影响较大的区域（瑞利峰和拉曼峰主导信号）中的数据，并用缺失值（“NaN”（在MATLAB中不是数字）替换它们。此外，还插入了一个零区域，以辅助PARAFAC建模。
在继续下一步之前，我们输入：

[CutData] = EEMCut(OriginalData,20,20,NaN,NaN,'');

以便在本教程的其余部分对数据进行适当的处理。这将剪切数据，但不会绘制结果（因为最后的参数为’ ‘）。

第三步：初步探索性数据分析和异常值识别

这个步骤用于去除一些异常的样本以及确定模型最后的组分数。我们输入：

[Test1] = OutlierTest(CutData,2,1,7,'No','No');

进行第一次测试。按下任意键后，将运行从2个组件到7个组件的一系列模型。最初五个模型（2、3、4、5、6和7个组件模型）的结果可以通过多种方式进行评估，但在本教程中，我们将使用两种类型的图；负载和杠杆。

运行完之后我们输入：

PlotLL(Test1,2)

将创建一个图，显示模型的分数和载荷，绘制两个组件的发射和激发载荷。图a显示了两种成分的浓度在样品之间的变化。检查载荷是否合理，即它们是否平滑，光谱外观是否符合您对数据的理解。

这里我们发现有些组分的loadings存在负值，针对这种情况我们需要添加非负约束后重新运行，我们输入：

[Test2] = OutlierTest(CutData,2,1,7,'Yes','No');

我们再继续查看，输入:

PlotLL(Test2,4)

这时候发现样本7和19偏离的相对来说很厉害（大于0.3），所以我们对这两个样本进行去除，我们输入：

[removal] = RemoveOutliers(CutData,[7 19],[],[]);

再继续添加非负约束，我们输入：

[Test3] = OutlierTest(removal,1,1,7,'Yes','No');

可以看到现在几乎所有的样本都在理想的范围内（小于0.3），再次尝试绘制所有模型的载荷和杠杆（例如PlotLL（Test3,3））。现在，载荷似乎更符合逻辑（即，它们不显示负荧光），接着我们进行下一步。

第四步 确定组分数

查看残差图（不应该出现大面积的负值）：

EvalModel(Test3,4)

绘制曲面的残差图：

EvalModelSurf(Test3,4)

接下来，我们可以使用Compare2Models（或Compare2ModelsUrf）和CompareSpescse函数比较两个不同模型的结果，我们输入：

Compare2Models(Test3,3,4)

如果要绘制曲面图可以输入：

Compare2ModelsSurf(Test3,3,4)

接下来，可以使用比较函数检查平方误差的谱和，我们输入：

CompareSpecSSE(Test3,3,4,5)

第五步：拆半检验

现在数据将被分成两半。不只是进行一种类型的拆分，而是进行两种不同的拆分，每个拆分以不同的方式对数据进行拆分。这使我们能够运行两个对半分析，并使用性能最好的一个。有时，如果样本类型在两半之间分布不均，则拆分可能是次优的。理想情况下，分割应尽可能“随机”。下图显示了如何使用工具箱中的SplitData函数拆分数据。

我们输入：

[AnalysisData]=SplitData(Test3);

接下来，使用SplitHalfAnalysis功能执行SplitHalfAnalysis，我们输入：

[AnalysisData]=SplitHalfAnalysis(AnalysisData,(3:7),'MyData.mat');

现在，可以使用SplitHalfValidation函数检查分析结果（输入help SplitHalfValidation以获取解释）。该函数使用（Lorenzo Seva和Berge，2006）描述的Tucker同余系数，在数学上比较在单独的数据分割上运行的模型的激励和发射载荷，并在命令窗口中说明模型是否得到验证。此外，还绘制了每一半的发射和激发载荷，以便直观地进行比较，我们输入：

SplitHalfValidation(AnalysisData,'1-2',4)

SplitHalfValidation(AnalysisData,'3-4',4)

检验结果：

第六步:随机初始化

接下来，将使用模型的随机初始化对整个数据拟合一系列模型。我们需要确保我们导出的模型实际上是最小二乘结果，而不是局部极小值。PARAFAC模型使用交替最小二乘算法进行拟合，与其他迭代拟合程序（例如简单非线性回归）一样，它可能会受到初始估计值的影响。在这一步中，使用随机初始化过程来拟合具有相同数量组件的多个模型，以便我们找到真正的最小二乘结果（最佳拟合）。

使用随机初始化运行10个四组件模型。这需要15-20分钟，我们输入：

[Random]=RandInitAnal(AnalysisData,4,10);

一旦建模完成，将绘制一个显示每个模型的误差平方和的图。具有最小二乘结果的模型用绿色圆圈高亮显示。（无需使用RandInitResult函数重新运行分析，即可重新绘制曲线图）。

最后，检查发射和激发载荷是否与之前的对半验证中发现的相同。应始终如此，但建议进行检查。这可以通过重新绘制对半分析结果来实现 ，我们输入：

SplitHalfValidation(Random,'1-2',4)

如果没有问题我们进行下一步。

第七步：绘制组分图

要创建每个组件的曲面或等高线图，可以使用ComponentEEM或ComponentSurf函数。每个组件都将在单独的窗口中标绘，我们输入：

ComponentEEM(Random,4)

绘制曲面的组分图：

ComponentSurf(Random,4)

注：在使用ComponentEEM或者ComponentSurf这两个函数时，经常会出现报错的情况。

我在网上查了很多，现有的解决资料比较少，找到了一种解决的方法，我们打开这个函数的位置，将这一行注释掉，如图：

原因可能是这一行代码会将激发大于发射的部分替换成NaN并且全部删除掉，所以导致矩阵的行数和列数与之前矩阵的维度不一样导致报错，但是这里直接删除掉出来的组分图就不会在出现本应该扣除的部分（我也不知道解释的清不清楚，反正大概是这个意思，这两个函数都有这一行，都要注释掉）。

第八步：数据导出

如果需要，可以使用ModelOut函数将结果导出到MATLAB中。这将创建一个excel文件，其中包含每个样品中每个成分的荧光强度以及每个成分的发射和激发载荷，我们输入：

[FMax,B,C]=ModelOut(Random,4,'D:\MyParafacResults.xls')

导出路径可以根据自己的需要进行修改。

如果想要在MATLAB中继续工作，这些数据也可以在MATLAB工作区中分别作为FMax、B和C获得。

四、总结

学习平行因子分析法的道路确实很坎坷，但学习不能停止，如果之后又学习到了新的知识也会再进行补充，有问题也欢迎大家私信我一起探讨。

以上！！