MATLAB: REPEATED MEASURES ANOVA MATLAB

anovacovariateMATLABranovarepeated measuresstatistics

Hi everyone!
I am having troubles understanding how to perform the correct anova test on my data. I have 16 subjects who are being treated with 2 different treatments. For each subject and each treatment measures are taken at 6 different time points (attached is the excel file). Is repeated measures anova the best way to go? If so, what is the most apporpriate command to use?
Here is the code I have so far. The results are kind of strange so I am not sure if it is wrong or I am just interpreting it the wrong way.
Thanks,
Giulia
 Q2 = xlsread('Q2.xls');
 Treatment = categorical(Q2(:,2)); 
 ID = Q2(:,1); 
 t = table(Treatment, Q2(:,3),Q2(:,4),Q2(:,5),Q2(:,6),Q2(:,7), Q2(:,8),'VariableNames', {'Treatment', 't1', 't2', 't3', 't4', 't5', 't6'});
 Time = [1 2 3 4 5 6]';
 rm = fitrm( t, 't1-t6~Treatment', 'WithinDesign', Time); 
 ranovatbl = ranova(rm); 
 disp(ranovatbl);
 multcompare(rm, 'Time', 'By', 'Treatment')

Best Answer

Is repeated measures anova the best way to go?
Your design is certainly applicable to a repeated measures design but the answer to that question depends on what you're testing.
It looks like you want to perform a repeated measures ANOVA with a covariate. In addition to patient number, the treatment type is the covariate, and the 6 time points is the within-subject repeated measure.
There are 4 null hypotheses to look at:
  • Patient number has no effect on the population mean
  • Treatment type has no effect on the population mean
  • There is no interaction between patient number and treatment type (usually the most important question)
  • There is no 3-way interaction between patient number, treatment type, and treatment time.
Setting up a RM Anova with a convariate in Matlab
Your data are attached in a file named Q2unlocked.xlsx. The following lines set up a RM Anova with a covariate (written in r2019b).
% Read in the table
Q2 = readtable('Q2unlocked.xlsx');
% Replace the 6 measurment headers 
Q2.Properties.VariableNames(3:8) = {'t1', 't2', 't3', 't4', 't5', 't6'}; 
% The within-subjects design
withinDsgn = table((1:6)','VariableNames',{'Time'});  
% Run the RM Anova (note the interaction between PATIENT and TREAT!)
rm = fitrm(Q2, 't1-t6~PATIENT*TREAT', 'WithinDesign', withinDsgn);
Also see a similar example provided by Matlab.
Testing for assumptions
Before we look at the results, you must make sure your data are appropriate for a RM Anova by confirming the following assumptions.
  1. Independent observations: You collected the data randomly and without bias
  2. Normality: The measurments have an approximately normal distribution
  3. Sphericity: we'll use the Mauchly test.
% Are the data approximately normal? This should ideally be done for each factor
histogram(reshape(Q2{:,3:end},[],1))
A small rightward tail but not too bad. 
% Does the data pass the sphericity test?
rm.mauchly
ans = 1×4 table
       W       ChiStat    DF    pValue 
    _______    _______    __    _______

    0.48842    11.537     14    0.64344
Yes: a pValue < 0.05 would fail
https://www.mathworks.com/help/stats/mauchlys-test-of-sphericity.html
RM Anova results
ranova(rm) or rm.ranova produce the results table showing the p-value and then 3 adjusted p-values depending on whether the the response variables have different variances (symmetry assumption). 
% Produce ranova table
rm.ranova
ans = 5×8 table
                            SumSq       DF      MeanSq        F        pValue     pValueGG    pValueHF    pValueLB
                          __________    __    __________    ______    ________    ________    ________    ________

    (Intercept):Time      3.9654e+06     5    7.9308e+05    2.2798    0.053221    0.071668    0.053221    0.14843 
    PATIENT:Time          4.0894e+06     5    8.1788e+05    2.3511    0.047012    0.064753    0.047012    0.14259 
    TREAT:Time            3.2657e+06     5    6.5314e+05    1.8775     0.10606     0.12652     0.10606    0.18747 
    PATIENT:TREAT:Time    3.9524e+06     5    7.9048e+05    2.2723    0.053916    0.072433    0.053916    0.14905 
    Error(Time)           3.1309e+07    90    3.4788e+05
Row 1 representing all differences across the within-subjects factors (treatment times). It has a p value of 0.053 which is just beyond the socially accepted threshold of 0.05. A boxplot will show us if this value seems reasonable.
figure()
boxplot(Q2{:,3:end})
xlabel('Treatment times')
ylabel('measurement value')
title('Data pooled between all Patients and treatment factors')
It doesn't look like there's much of an effect of treatment time when factors are combined
Row 2 shows the interaction between Patient and treatment times and is p=0.047
bpdata = [];  
for i = 1:max(Q2.PATIENT) %assuming patient numbers are 1:max
    bpdata = [bpdata, Q2{Q2.PATIENT==i,3:8},nan(size(unique(Q2.TREAT)))]; 
end
figure()
boxplot(bpdata)
arrayfun(@xline,7:7:size(bpdata,2))
xlabel('6 treatment times across 11 patients')
ylabel('measurement value')
title('Data pooled between treatment factors')
set(gca,'XTick', [])
Clearly the difference between the 6 treatment times differs across subjects (Treatment type combined)
Row 3 shows the interaction between treatment type and treatment times and is p=0.106.
bpdata = [];  
for i = 1:max(Q2.TREAT) %assuming TREAT numbers are 1:max
    bpdata = [bpdata, Q2{Q2.TREAT==i,3:8},nan(size(unique(Q2.PATIENT)))]; 
end
figure()
boxplot(bpdata)
xline(7)
xlabel('6 treatment times across 2 TREAT factors')
ylabel('measurement value')
title('Data pooled between patients')
set(gca,'XTick', [])
Other than the 2nd treatment factor there is no difference between the 2 Treatment types.
Row 4 shows the 3-way interaction between patient, treatment type, and treatment times.
Related Question