Clear Filters
Clear Filters

Set up Repeated Measures Anova function MATLAB

17 views (last 30 days)
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.
  3 Comments
Jacob Jacobo
Jacob Jacobo on 25 Jul 2022
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
Scott MacKenzie on 25 Jul 2022
@Jacob Jacobo, thanks for the clarification. I just posted an answer. Good luck.

Sign in to comment.

Accepted Answer

Scott MacKenzie
Scott MacKenzie on 25 Jul 2022
Edited: Scott MacKenzie on 25 Jul 2022
@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

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!