Set up Repeated Measures Anova function MATLAB

Hello Everyone,

I want to run an ANOVA comparing 3 treatments (labeled as 1, 2, and 3) and then run a post-hoc comparison between treatments if the ANOVA shows a difference. I have 16 individual patients who all receive 3 different treatments and each treatment has 3 independent measurements taken at the same time. Here is the code I have so far:

data = readtable('Data.xlsx') ;

% Removing the 'Patients' column 
t = data(:,2:end) ;

% Within design
WithinDesign = table((1:3)','VariableNames',{'Measurements'}) ;

% Repeated measures model
rm = fitrm(t,'AA-FE~Treatment','WithinDesign',WithinDesign) ;

% Sphericity test
rm.mauchly

ans = 1×4 table
       W       ChiStat    DF      pValue  
    _______    _______    __    __________

    0.70763    15.563     2     0.00041744

% Anova
ranova(rm)

ans = 3×8 table
                                  SumSq       DF    MeanSq      F         pValue       pValueGG      pValueHF      pValueLB 
                                __________    __    ______    ______    __________    __________    __________    __________

    (Intercept):Measurements    1.9978e+05     2     99890     245.2    1.3625e-37    1.1519e-29    1.9743e-30    4.6883e-20
    Treatment:Measurements           27529     2     13764    33.787    9.9539e-12    1.3755e-09    8.5385e-10    5.5235e-07
    Error(Measurements)              37480    92    407.39                                                                  

I understand that the data does not pass the sphericity test caclulated by rm.mauchly, but I would still like to know whether or not my ANOVA set-up represents what I wanted to acquire from ANOVA since I would like to do this in the future.

[Scott MacKenzie]

@Jacob Jacobo I can help with the ANOVA and post hoc comparisons, but first I'd like to verify that I understand your design. Are AA, IE, and FE measurements taken for each treatment? You refer to them as "independent measurements", but they sound like dependent variables to me. That is, AA is a measurement that depends on the treatment (1, 2, or 3). The same can be said for IE and FE, if I understand correctly. Just a terminology thing, perhaps. Seems to me your independent variable is treatment, with levels 1, 2, and 3. So, I would characterize your design as a single-factor design with one within-subjects independent variable (treatment with levels 1, 2, and 3) and three dependent variables (AA, IE, and FE). Is this correct?

[Jacob Jacobo]

Hi @Scott MacKenzie,

Seems like I was wrong in my original post. To answer your questions; yes, AA, IE, and FE should be dependent variables. These 3 measurements are taken for each treatment and you are correct that the measurement value depends on the treatment. I apologize for the confusion. The goal would be to compare how a single measurement type compares among the three treatments, i.e. how the AA measurements in treatment 1 compares to the AA measurements in treatments 2 and 3 and so on (and also being aware that each patient is given each treatment). On that note, your suggested design seems correct, although I am not sure what needs to change with the lines of code that I currently have. Would I need to include each variable in the within design or maybe change the format of my data table? I have looked in the documentation and I cannot seem to find an answer. Thank you very much for helping.

[Scott MacKenzie]

@Jacob Jacobo, thanks for the clarification. I just posted an answer. Good luck.

Answer 1

[Scott MacKenzie]

@Jacob Jacobo, your setup for fitrm is slightly wrong, since you only have a single within-subjects factor. Below is what I put together for the AA dependent variable. The effect of treatment on AA was statistically significant, F(2,30) = 31.3, p < .0001. All six of the pairwise comparisons are also significant. You'll get similar results for the IE and FE dependent variables.

M = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1073380/Data.xlsx');

% extract and reorganize data for the AA dependent variable (1 row per subject, 1 column per treatment)
AA = reshape(M(:,3), [], 3)

AA = 16×3

100.6000 79.3000 66.8000 99.9000 73.0000 68.8000 108.7000 79.3000 68.7000 118.4000 75.5000 64.9000 84.4000 105.9000 84.0000 92.0000 86.3000 92.1000 87.7000 97.6000 76.9000 92.6000 81.7000 73.9000 117.2000 67.7000 62.6000 113.0000 80.2000 69.3000

% put the AA data into a table
T = array2table(AA, 'VariableNames', {'T1', 'T2', 'T3'});

withinDesign = table([1 2 3]', 'VariableNames', {'Treatment'});
withinDesign.Treatment = categorical(withinDesign.Treatment);
rm = fitrm(T, 'T1-T3 ~ 1', 'WithinDesign', withinDesign);
AT = ranova(rm, 'WithinModel', 'Treatment');

% output a conventional ANOVA table
disp(anovaTable(AT, 'AA'));

ANOVA table for AA
===============================================================================
Effect                   df          SS             MS         F         p
-------------------------------------------------------------------------------
Participant              15      907.06583       60.47106                      
Treatment                 2     7155.41792     3577.70896     31.272     0.0000
Participant(Treatment)   30     3432.15542      114.40518                      
===============================================================================

% do the pairwise comparisons (3 treatments, therefore 6 comparisons)
multcompare(rm, 'Treatment')

ans = 6×7 table
    Treatment_1    Treatment_2    Difference    StdErr      pValue       Lower      Upper 
    ___________    ___________    __________    ______    __________    _______    _______

         1              2           19.038      4.6778     0.0027196      6.887     31.188
         1              3           29.494      4.3018    1.5408e-05      18.32     40.667
         2              1          -19.038      4.6778     0.0027196    -31.188     -6.887
         2              3           10.456      1.5858    2.4072e-05     6.3371     14.575
         3              1          -29.494      4.3018    1.5408e-05    -40.667     -18.32
         3              2          -10.456      1.5858    2.4072e-05    -14.575    -6.3371

% -------------------------------------------------------------------------
% Function to create a conventional ANOVA table from the overly-complicated 
% and confusing anova table created by the ranova function.
function [s] = anovaTable(AT, dvName)

c = table2cell(AT);

% remove erroneous entries in F and p columns
for i=1:size(c,1)
    if c{i,4} == 1
        c(i,4) = {''};
    end
    if c{i,5} == .5
        c(i,5) = {''};
    end
end

% use conventional labels in Effect column
effect = AT.Properties.RowNames;
for i=1:length(effect)
    tmp = effect{i};
    tmp = erase(tmp, '(Intercept):');
    tmp = strrep(tmp, 'Error', 'Participant');
    effect(i) = {tmp};
end

% determine the required width of the table
fieldWidth1 = max(cellfun('length', effect)); % width of Effect column
fieldWidth2 = 57; % width for df, SS, MS, F, and p columns
barDouble = repmat('=', 1, fieldWidth1 + fieldWidth2);
barSingle = repmat('-', 1, fieldWidth1 + fieldWidth2);

% re-organize the data
c = c(2:end,[2 1 3 4 5]);
c = [num2cell(repmat(fieldWidth1, size(c,1), 1)), effect(2:end), c]';

% create the ANOVA table
s = sprintf('ANOVA table for %s\n', dvName);
s = [s sprintf('%s\n', barDouble)];
s = [s sprintf('%-*s %4s %11s %14s %9s %9s\n', fieldWidth1, 'Effect', 'df', 'SS', 'MS', 'F', 'p')];
s = [s sprintf('%s\n', barSingle)];
s = [s sprintf('%-*s %4d %14.5f %14.5f %10.3f %10.4f\n', c{:})];
s = [s sprintf('%s\n', barDouble)];
end

[Jacob Jacobo]

Thank you very much for your help!

[Scott MacKenzie]

You're welcome. Glad to help. Good luck.