方差分析，即analysis of variance，简写为ANOVA，指的是对样本的总方差进行分解分析，当自变量的因子大于等于三个类别情况下，检验其各类别间平均数是否相等的统计模式。还记得当年上心理统计和SPSS操作与应用时，老师对方差分析的平方和分解极为细致的讲解，也是到了研究生阶段才领悟到方差分析中对因变量变异来源分解的思想。
其实用Matlab进行方差分析大可不必，用SPSS或者R language进行分析应该都会比较Matlab更友好。只是最近在做时序数据的分析，需要看随着时间的变化某个变量的效应如何改变。例如EEG数据就是一个时序数据，当我们想看不同实验条件在某个脑电指标上的效应是如何随着时间的改变而改变（例如，什么时候开始达到显著性水平）的时候，就需要在每个时间点进行一次比较。由于每次计算和比较的数据结构其实就是矩阵或者向量，前期的数据分段、根据时间点进行累积到后面的t检验，一般就用Matlab一次完成了。
我在数据分析过程中遇到的问题在于，当我将数据成功分段并累积后，发现无法通过方差分析检验某个自变量的效应。因为Matlab中自带的方差分析函数只有anova、anova2以及anovan，但是这几个函数只能进行普通的方差分析，对于重复测量设计的实验并不适用，特别是我的实验设计中每个组内自变量都有3个因素。想着实在不行，需要自己写重复测量方差分析的函数了，还是很沮丧的，因为这对我肯定是又一个大工程。不过天无绝人之路，我在MathWorks的论坛上看到有用户（Aaron Schurger）上传了他自己写的两因素重复测量方差分析的函数代码，我就下载下来开始研究。
这个代码不长，最终的函数结构长这样：
rm_anova2(Y,S,F1,F2,FACTNAMES)一共有5个参数，分别是：①Y代表因变量；②S代表样本编号；③F1代表第一个自变量的水平；④F2代表第二个自变量的水平；⑤ FACTNAMES 代表两个自变量的名称。
除了 FACTNAMES 时字符串构成的元胞之外，其余四个参数都必须是一个一维的列向量，并且需要一一对应。举个例子，我们常见的两因素重复测量方差分析（以2×2组内设计为例）数据结构为：
| participants_id | a1b1 | a1b2 | a2b1 | a2b2 |
| 1 | 2 | 4 | 3 | 6 |
| 2 | 2 | 3 | 4 | 6 |
| 3 | 1 | 2 | 3 | 5 |
在使用这个函数之前，就需要将常见的数据结构转换为该函数所要求的数据结构，在rm_anova2函数中，每一个参数的行数都是总数据点个数，即被试数×自变量1水平数×自变量2水平数。将上例转换后的数据结构为：
| Y | S | F1 | F2 |
| 2 | 1 | 1 | 1 |
| 2 | 2 | 1 | 1 |
| 1 | 3 | 1 | 1 |
| 4 | 1 | 1 | 2 |
| 3 | 2 | 1 | 2 |
| 2 | 3 | 1 | 2 |
| 3 | 1 | 2 | 1 |
| 4 | 2 | 2 | 1 |
| 3 | 3 | 2 | 1 |
| 6 | 1 | 2 | 2 |
| 6 | 2 | 2 | 2 |
| 5 | 3 | 2 | 2 |
函数最终输出的结果是一个7×6的元胞：
| source | SS | df | MS | F | p |
| F1 | |||||
| F2 | |||||
| F1×F2 | |||||
| F1×Sub | |||||
| F2×Sub | |||||
| F1×F2×Sub |
该函数最终可以输出不同变异来源的F值和对应的p值（显著性水平），但是不能进一步进行多重比较和简单效应分析。但是如果只是想看不同时间段的p值如何随着时间的变化而变化，却是完全够用的，该函数的一个主要缺点在于适用其的数据结构和常见的两因素重复测量方差分析的数据结构有很大区别。基于此，我对该代码进行了简单改良，在用该进行方差分析之前，可以先用几行代码将常见的数据结构转化为适合该函数的数据结构，并能够用for循环来将多个时间点下p值一一输出。
%设置参数
sub = 30; % 被试数量
F1_con = 3; %F1自变量水平数
F2_con = 3; %F2自变量水平数
total_con = F1_con * F2_con ; %实验条件数目，这里是3×3的实验设计。
%构建参数矩阵
%Y为因变量，需要整合为一维列向量
Y = [此处根据原始结构的数据一列列的排列]
%S为被试，和读取文件数目有关
S = repmat([1:sub], 1, total_con).';
% F1和F2 为自变量的水平
F1 = sort(repmat([1:F1_con], 1, sub*F2_con)).';
F2 = repmat(sort(repmat([1:F2_con], 1, sub)), 1, F1_con).';
%FACTNAMES为变量名称
FACTNAMES = {'wenti', 'bianshu'};
%重复测量方差分析
tab = rm_anova2(Y, S, F1, F2, FACTNAMES);
%输出p值
P = tab{2, 6};
以上是一个可以考虑的解决方案，但需要注意的是，在进行多次方差分析比较时，一定要进行多重比较校正（FDR），校正的次数和比较的次数相关，该模块需要和for循环的模块配合使用。
Citation
Aaron Schurger (2021). Two-way repeated measures ANOVA (https://www.mathworks.com/matlabcentral/fileexchange/6874-two-way-repeated-measures-anova), MATLAB Central File Exchange. Retrieved July 31, 2021.
如果有同仁还有其他解决方案，或者有在此基础上进行多重比较和简单分析的改良代码，欢迎交流与分享。可从联系页面联系到我。