如何用MATLAB对多组被试重复测量数据进行rANOVA分析

好的，我来一步步指导你用MATLAB完成这个重复测量方差分析（rANOVA）。首先咱们先明确你的数据结构：你已经把3组受试者（青年/中年/老年）的6项测量值整理成了MATLAB Table，其中meas1/meas2是静息心率的两次重复测试，meas3/meas4是运动后心率的两次重复测试，meas5是年龄，meas6是血压。接下来分两种常见的分析场景来讲解：

场景1：单独分析某一类心率的重复测量（比如仅静息心率） #

如果你只想针对静息心率的两次测试块差异，结合不同年龄组的组间差异做分析，可以按以下步骤操作：

1. 预处理数据，明确变量类型
   首先确保分组变量是分类类型（MATLAB的重复测量函数对分类变量支持更好），然后提取需要的列：

   % 把分组变量转成分类类型
   Meas.groups = categorical(Meas.groups);
   % 提取静息心率相关数据，重命名变量方便理解
   rest_data = table(Meas.groups, Meas.meas1, Meas.meas2, ...
       'VariableNames', {'AgeGroup', 'Rest_Block1', 'Rest_Block2'});
   

2. 创建重复测量模型并执行rANOVA
   使用MATLAB的fitrm函数创建重复测量模型，再用ranova执行方差分析：

   % 创建重复测量模型：组内因子是两次测试块，组间因子是年龄组
   rm_rest = fitrm(rest_data, 'Rest_Block1-Rest_Block2 ~ AgeGroup');
   % 执行重复测量方差分析
   ranova_results = ranova(rm_rest);
   % 查看结果
   disp(ranova_results);
   

3. 结果解读
   ranova_results表格会输出以下关键信息：

   • (Intercept)：整体心率的均值效应（一般不用重点关注）
   • AgeGroup：不同年龄组之间的静息心率组间差异（p值<0.05说明组间有显著差异）
   • Time：两次测试块之间的静息心率组内差异（也就是重复测量的主效应）
   • AgeGroup:Time：年龄组和测试块的交互效应（p值<0.05说明不同年龄组在两次测试块中的心率变化趋势不同）

场景2：同时分析静息+运动后心率的双重复测量 #

如果想同时分析**心率状态（静息/运动后）和测试块（1/2）**两个组内因子，结合年龄组的组间差异，需要先把宽格式数据转成长格式：

1. 重塑数据为长格式
   把每个受试者的4次心率测量拆分成单独的行，同时标记对应的心率状态和测试块：

   % 把宽格式的4项心率测量转成长格式
   long_data = stack(Meas, {'meas1','meas2','meas3','meas4'}, ...
       'NewDataVariableName', 'HeartRate', ...
       'IndicatorVariableName', 'Condition');
   % 给每个测量标记心率状态和测试块
   long_data.HeartState = repmat({'Rest','Rest','PostExercise','PostExercise'}, height(Meas), 1);
   long_data.Block = repmat({'Block1','Block2','Block1','Block2'}, height(Meas), 1);
   % 转成分类变量
   long_data.HeartState = categorical(long_data.HeartState);
   long_data.Block = categorical(long_data.Block);
   long_data.groups = categorical(long_data.groups);
   

2. 创建双组内因子的重复测量模型
   这里的模型包含组间因子（年龄组）、两个组内因子（心率状态、测试块），以及它们的交互效应：

   % 创建重复测量模型，指定组内因子为HeartState和Block
   rm_combined = fitrm(long_data, ...
       'HeartRate ~ groups + HeartState*Block + groups*HeartState*Block', ...
       'WithinDesign', {'HeartState','Block'});
   % 执行方差分析
   ranova_combined = ranova(rm_combined);
   disp(ranova_combined);
   

3. 结果解读
   除了场景1中的组间/组内效应，你还会看到：

   • HeartState：静息和运动后心率的差异主效应
   • Block：测试块的主效应
   • HeartState:Block：心率状态和测试块的交互效应（比如运动后心率在两次测试块中的变化是否和静息不同）
   • groups:HeartState：年龄组和心率状态的交互效应（比如不同年龄组的运动后心率提升幅度是否有差异）

关键注意事项 #

• 球形性检验：ranova的结果中会包含Mauchly球形性检验的p值，如果p<0.05说明球形性假设不成立，此时需要看表格中GG Correction或HF Correction列的校正p值，避免第一类错误。
• 事后检验：如果发现显著的主效应或交互效应，可以用multcompare函数做两两比较，比如：

  % 比较不同年龄组的心率差异
  mc_groups = multcompare(rm_rest, 'AgeGroup');
  % 比较两次测试块的心率差异
  mc_block = multcompare(rm_rest, 'Time');
  

• 协变量分析：如果想把年龄（meas5）或血压（meas6）作为协变量，可以在fitrm的模型公式中加入协变量，比如：'Rest_Block1-Rest_Block2 ~ AgeGroup + meas6'

内容的提问来源于stack exchange，提问作者Juliette