三维荧光平行因子学习记录–（三）使用DOMfluor工具箱进行平行因子分析–（一）
注 本文仅作为自己的学习记录以备以后复习查阅
本文的参考文献：
大家对三维荧光数据分析有疑问的或者对平行因子分析感兴趣的可以去看看他的主页。
一、前言
本文讲分为两个大模块来展示DOMfluor工具箱的使用，第一个部分是用官方提供的数据（PARAFACexample.mat）进行分析，第二个部分使用我们自己的数据集进行分析，其中也会穿插一些报错的解决办法，这里默认大家都已经安装了工具箱，没安装的可以去看看我之前发的博客，下面话不多说进入正题。
二、官方自带数据
第一步：加载数据
load PARAFACexample.mat加载数据并绘制EEM，在命令窗口中，把教程数据加载到MATLAB工作空间中，它由一个名为OriginalData的数据结构组成。
第二步：绘制等高线图
现在我们将绘制EEM，以检查数据是否正确加载。工具箱中有几个函数可用于绘制EEM。
①将数据集中的前5个EEM一次绘制一个等高线图。按下键盘上的任何键都可以查看下一个图形。如果将1:5替换为1:65，则所有65个样本将一次绘制一个。通过在键盘上输入[Ctrl+C]，然后关闭图形窗口，可以随时停止打印过程 ：
注：‘R.U.’，表示你数据的荧光强度的单位，没有经过校正，即任意单位时，为’A.U.‘，经过拉曼校正后，为’R.U.’，经过硫酸奎宁校正后，为’QSE’
PlotEEMby1(1:5,OriginalData,'R.U.')
②将一次绘制四个数据。这在处理大型数据集时更快。
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')
PlotSurfby4(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个组件模型）的结果可以通过多种方式进行评估，但在本教程中，我们将使用两种类型的图；负载和杠杆。这些测试大约需要3分钟，您的计算机才能运行。
运行完之后我们输入：
PlotLoadings(Test1,2)
将创建一个图，显示模型的分数和载荷，绘制两个组件的发射和激发载荷。图a显示了两种成分的浓度在样品之间的变化。检查载荷是否合理，即它们是否平滑，光谱外观是否符合您对数据的理解。
同时也可以输入以下代码，查看函数相关信息：
PlotLeverage(Test1,2)
将创建一个显示杠杆的图形，在图中，我们去寻找具有极端杠杆的样本，这些杠杆指示极端和潜在的边缘样本。
使用函数PlotLL可以创建杠杆图和加载图的组合，我们输入：
PlotLL(Test1,2)
现在尝试为其他模型（3、4、5、6和7个组件）创建相同的图。在检查加载曲线图之后，我们显然需要开始约束模型。有些负荷是负的。例如，在4组件模型中，其中一个组件的a值为负值，似乎抵消了其他组件中的一个。这些“错误”荷载可能是因为使用了错误数量的组件，但在假设四个组件有效的情况下进行分析，通过强制参数为非负来避免问题是合理的。在一个更为最终的模型中，关注为什么模型需要这样一个附加约束来提供具有化学意义的结果是有用的，但这在现阶段是不必要的。
重新运行测试，但输入非负性约束 ：
[Test2] = OutlierTest(CutData,2,1,7,'Yes','No');现在绘制负载和杠杆。（可以通过绘制两个测试的结果来比较应用非负性约束的效果，例如，对于四组件模型类型PlotLL（Test1,4），然后是PlotLL（Test2,4））。然后绘制其他模型的载荷和杠杆。，下面我们来检查一下Test2，输入：
PlotLL(Test2,4)
杠杆图似乎表明，样本5、21、30、40和53可能存在问题。因此，我们现在应该尝试绘制这些样本的EEM，并将它们与其他样本进行比较，看看它们是否包含一些测量误差。（例如，使用PlotEEMby4（1，CutData，'R.U.），查看这些样本是否与其他样本不同）。只有样本5和30是带有测量误差的明显异常值，所以我们输入：
[Test3]=RemoveOutliers(CutData,[5 30],[],[])再输入：
[Test3]=OutlierTest(Test3,1,1,7,'No','No')然后绘制每个模型的杠杆和载荷数据，例如PlotLL（Test2,3）。杠杆似乎略有改善。然而，正如前面提到的，我们仍然存在负浓度问题。这可以通过对模型应用非负性约束来解决。我们可以重新运行模型，并使用非负性约束覆盖Test3，我们输入：
[Test3] = OutlierTest(removal,1,1,7,'Yes','No');此时我们再检查以下Test3：
PlotLL(Test3,4)
再次尝试绘制所有模型的载荷和杠杆（例如PlotLL（Test3,3））。现在，载荷似乎更符合逻辑（即，它们不显示负荧光）。
第五步：确定组分数
我们将使用的下一个函数是EvalModel，它创建了一系列图表，我们可以通过查看残差（测量值减去建模数据）来评估模型拟合度，我们输入：
EvalModel(Test3,4)
输入下面的代码来绘制曲面图（曲面图会一直绘制到最后一个样本）：
EvalModelSurf(Test3,4)
并检查为四分量模型绘制的图。特别是检查剩余EEM。理想情况下，这不应包含任何系统信号，且主要由仪器噪声组成。可以看出，四组分模型通常适用于大多数样本。然而，有少数样品在325 nm激发和400 nm发射区域具有无法解释的荧光。接下来，使用相同的功能评估其他模型（2、3、5和6个组件）。还可以尝试在去除异常值之前绘制拟合模型的结果（例如Test2），并检查已识别的异常值样本的残差。
接下来，我们可以使用Compare2Models（或Compare2ModelsUrf）和CompareSpescse函数比较两个不同模型的结果，我们输入：
Compare2Models(Test3,3,4)如果想要绘制曲面图可以输入：
Compare2ModelsSurf(Test3,3,4)比较3个和4个组件模型。同样地，比较4和5组件模型以及5和6、6和7组件模型。从这些图来看，似乎至少需要四个组件。
接下来，通过键入，比较4、5、6和7组件模型的载荷 ：
PlotLoadings(Test3,4)
PlotLoadings(Test3,5)
PlotLoadings(Test3,6)
PlotLoadings(Test3,7)
在这里，我们可以看到随着模型复杂性的增加，组件是如何演变的。每个模型的组成部分通常非常相似，每一步，一个组成部分被分成两个，以便更好地建模数据，并消除残差中留下的系统可变性。这是一个积极的结果，表明模型之间的一致性，而不是随机的局部极小值结果。
接下来，可以使用比较函数检查平方误差的谱和，我们输入：
CompareSpecSSE(Test3,3,4,5)
绘制了三种不同模型在激发和发射方向上的残差平方和。还可以尝试比较4、5和6组件模型，最后比较5、6和7组件模型。从这些结果来看，似乎很明显，从6到7组分模型的步骤对拟合的改善很小，这表明6个或更少的组分足以获得这些数据。
第六步：拆半检验
在这一阶段的分析之后，我们已经识别并移除了可疑的异常样本，并发现为该数据建模所需的正确分量数量介于4到6个分量之间。现在是时候尝试验证这些模型了。
现在数据将被分成两半。不只是进行一种类型的拆分，而是进行两种不同的拆分，每个拆分以不同的方式对数据进行拆分。这使我们能够运行两个对半分析，并使用性能最好的一个。有时，如果样本类型在两半之间分布不均，则拆分可能是次优的。理想情况下，分割应尽可能“随机”。下图显示了如何使用工具箱中的SplitData函数拆分数据。
我们输入：
[AnalysisData]=SplitData(Test3);
可以从图中看到，所有的样本随机分配到四个组，第一半：（1-2）；第二半：（3-4）。进入命令窗口。绘制了四张图，显示了激发和发射方向上每一半数据的平方值之和。每对两半的曲线应该相似（1-2,3-4）。
接下来，使用SplitHalfAnalysis功能执行SplitHalfAnalysis，我们输入：
[AnalysisData]=SplitHalfAnalysis(AnalysisData,(3:7),'MyData.mat');进行分析。这将使带有3到7个组件的模型与数据相匹配。这一步需要10-20分钟，具体取决于计算机的速度。
现在，可以使用SplitHalfValidation函数检查分析结果（输入help SplitHalfValidation以获取解释）。该函数使用（Lorenzo Seva和Berge，2006）描述的Tucker同余系数，在数学上比较在单独的数据分割上运行的模型的激励和发射载荷，并在命令窗口中说明模型是否得到验证。此外，还绘制了每一半的发射和激发载荷，以便直观地进行比较，我们输入：
SplitHalfValidation(AnalysisData,'1-2',4)
SplitHalfValidation(AnalysisData,'3-4',4)
查看第一次拆分（1-2）时的三组件模型结果。你会发现它是经过验证的。现在尝试其他组件模型（3-7）和其他拆分（3-4）。你会发现四个组件可以被分成两部分。
第七步：随机初始化
接下来，将使用模型的随机初始化对整个数据拟合一系列模型。我们需要确保我们导出的模型实际上是最小二乘结果，而不是局部极小值。PARAFAC模型使用交替最小二乘算法进行拟合，与其他迭代拟合程序（例如简单非线性回归）一样，它可能会受到初始估计值的影响。在这一步中，使用随机初始化过程来拟合具有相同数量组件的多个模型，以便我们找到真正的最小二乘结果（最佳拟合）。
使用随机初始化运行10个四组件模型。这需要15-20分钟，我们输入：
[Random]=RandInitAnal(AnalysisData,4,10);一旦建模完成，将绘制一个显示每个模型的误差平方和的图。具有最小二乘结果的模型用绿色圆圈高亮显示。（无需使用RandInitResult函数重新运行分析，即可重新绘制曲线图）。
使用之前使用的PlotLL和EvalModel函数（例如PlotLoadings（AnalysisData，4）、EvalModel（AnalysisData，4））检查模型的加载和拟合。
最后，检查发射和激发载荷是否与之前的对半验证中发现的相同。应始终如此，但建议进行检查。这可以通过重新绘制对半分析结果来实现 ，我们输入：
SplitHalfValidation(Random,'1-2',4)如果没有问题我们进行下一步。
第八步：绘制组分图
要创建每个组件的曲面或等高线图，可以使用ComponentEEM或ComponentSurf函数。每个组件都将在单独的窗口中标绘，我们输入：
ComponentEEM(Random,4)
绘制曲面的组分图：
ComponentSurf(Random,4)
第九步：数据导出
现在，一个六分量模型已经拟合并验证了数据，如果需要，可以使用ModelOut函数将结果导出到MATLAB中。这将创建一个excel文件，其中包含每个样品中每个成分的荧光强度以及每个成分的发射和激发载荷，我们输入：
[FMax,B,C]=ModelOut(Random,4,'D:\MyParafacResults.xls')如果想要在MATLAB中继续工作，这些数据也可以在MATLAB工作区中分别作为FMax、B和C获得。
下一篇将重点说怎么用DOMfluor处理自己的数据集。
以上！！
