Pairwise post-hoc testing using coefTest
20 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hello,
I found this post from back in 2015 but it doesn't seem to answer the question (and if it does, I am sorry that I am missing it!). I am trying to understand how I can test for differences within a level of an interaction. For this purpose I have created a random dataset with known properties. In my mock dataset, I have a group of subjects who participate in something-or-other. They can belong to one of three Groups (ABC) at the time of test. Then, we have three factors (ABC) under which we have three levels for each (A-DEF, B-GHI, C-JKL). I arranged the data such that of the groupings, only Group C has any effect (i.e., 2) and then there is a significant interaction factor C, in that if you have level L, you get a big boost.
Thus, any analysis should hopefully detect a signficiant main effect of Group and a significant Group x Factor C interaction effect, but no other obvious effects. Within the effect of Group, it should be level C that stands out. Within the interaction, it should be level L that stands out.
My mock dataset is attached.
To test my understanding, I generated the below script. In it, I am (very sure) I am correctly using coefTest to detect the significant main and interaction effects. Within the main effect, I am (somewhat sure) I am correctly using coefTest to perform pairwise comparisons among Groups ABC; as predicted, I find that C is different from both A and B, while A and B do not differ.
However, within the interaction effect, I would like to test for differences among levels JKL. Can this be done? If so, can someone please help me to do so?
Thank you!!!
load table.mat
mdl = fitlme(d, ...
'Value ~ Group*FactorA + Group*FactorB + Group*FactorC + (1|Subject)', ...
'DummyVarCoding', 'effects', 'FitMethod', 'REML');
disp(mdl)
%% Generate predictions for ALL data
g = unique(d.Group);
a = unique(d.FactorA);
b = unique(d.FactorB);
c = unique(d.FactorC);
cv = sortrows(combvec(1:numel(g), 1:numel(a), 1:numel(b), 1:numel(c))');
predTable = table();
predTable.Subject(1:size(cv, 1),1) = categorical({'Generic'});
predTable.Group = g(cv(:,1),1);
predTable.FactorA = a(cv(:,2),1);
predTable.FactorB = b(cv(:,3),1);
predTable.FactorC = c(cv(:,4),1);
predTable.Value = predict(mdl, predTable, 'Conditional', false);
% Trim to only Groups and then only groupxfactor C interaction
[~, ia, ~] = unique(predTable.Group);
groupTable = removevars(predTable(ia,:), ...
{'Subject', 'FactorA', 'FactorB', 'FactorC'});
groupTable.Value = splitapply(@mean, predTable.Value, ...
findgroups(predTable.Group));
[~, ia, ~] = unique(predTable.FactorC);
AxFactorCTable = removevars(predTable(ia,:), ...
{'Subject', 'Group', 'FactorA', 'FactorB'});
AxFactorCTable.Value = splitapply(@mean, ...
predTable.Value(predTable.Group == 'A'), ...
findgroups(predTable.FactorC(predTable.Group == 'A')));
BxFactorCTable = removevars(predTable(ia,:), ...
{'Subject', 'Group', 'FactorA', 'FactorB'});
BxFactorCTable.Value = splitapply(@mean, ...
predTable.Value(predTable.Group == 'B'), ...
findgroups(predTable.FactorC(predTable.Group == 'B')));
CxFactorCTable = removevars(predTable(ia,:), ...
{'Subject', 'Group', 'FactorA', 'FactorB'});
CxFactorCTable.Value = splitapply(@mean, ...
predTable.Value(predTable.Group == 'C'), ...
findgroups(predTable.FactorC(predTable.Group == 'C')));
figure('units', 'normalized', 'outerposition', [0; 0; 1; 1])
subplot(2, 3, 2)
bar(groupTable.Group, groupTable. Value)
title('Group Means')
ylim([-3, 6])
axis square
subplot(2, 3, 4)
bar(AxFactorCTable.FactorC, AxFactorCTable. Value)
title('Factor C Means for Group A')
ylim([-3, 6])
axis square
subplot(2, 3, 5)
bar(BxFactorCTable.FactorC, BxFactorCTable. Value)
title('Factor C Means for Group B')
ylim([-3, 6])
axis square
subplot(2, 3, 6)
bar(CxFactorCTable.FactorC, CxFactorCTable. Value)
title('Factor C Means for Group C')
ylim([-3, 6])
axis square
% Perform pairwise testing among the main effects for each group.
% Because Group C is not assigned to an effect, it can be tricky to do
% pairwise comparisons with it. However, following the comparison
% approach found at
% https://www.mathworks.com/help/stats/generalizedlinearmixedmodel.coeftest.html,
% we can determine what the post hoc arrangements should be:
% syms A B C
% f1 = A + B + C == 0;
% C = solve(f1, C)
% fAB = A - B
% fAC = A - C
% fBC = B - C
% The results are:
% C = -A - B
% fAB = A - B
% fAC = 2*A + B
% fBC = A + 2*B
% Thus:
% [Int, G_A, G_B, FA_D, FA_E, FB_G, FB_H, FC_J, FC_K, G_A:FA_D, G_B:FA_D, G_A:FA_E, G_B:FA_E, G_A:FB_G, G_B:FB_G, G_A:FB_H, G_B:FB_H, G_A:FC_J, G_B:FC_J, G_A:FC_K, G_B:FC_K]
HAB = [ 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
HAC = [ 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
HBC = [ 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
pG = coefTest(mdl, [HAB; HAC; HBC]);
pAB = coefTest(mdl, HAB);
pAC = coefTest(mdl, HAC);
pBC = coefTest(mdl, HBC);
disp(['Significance of Group effect: ', num2str(pG, 3)])
disp([' * Significance of Group A vs. Group B difference: ', ...
num2str(pAB, 3)])
disp([' * Significance of Group A vs. Group C difference: ', ...
num2str(pAC, 3)])
disp([' * Significance of Group B vs. Group C difference: ', ...
num2str(pBC, 3)])
% Test for significant Group x Factor interactions:
% [Int, G_A, G_B, FA_D, FA_E, FB_G, FB_H, FC_J, FC_K, G_A:FA_D, G_B:FA_D, G_A:FA_E, G_B:FA_E, G_A:FB_G, G_B:FB_G, G_A:FB_H, G_B:FB_H, G_A:FC_J, G_B:FC_J, G_A:FC_K, G_B:FC_K]
HG_FA = [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]; ...
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]; ...
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]; ...
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0]];
HG_FB = [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0]; ...
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0]; ...
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]; ...
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0]];
HG_FC = [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0]; ...
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0]; ...
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0]; ...
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]];
pG_FA = coefTest(mdl, HG_FA);
pG_FB = coefTest(mdl, HG_FB);
pG_FC = coefTest(mdl, HG_FC);
disp(' ')
disp(['Significance of Group vs. Factor A interaction: ', ...
num2str(pG_FA, 3)]);
disp(['Significance of Group vs. Factor B interaction: ', ...
num2str(pG_FB, 3)]);
disp(['Significance of Group vs. Factor C interaction: ', ...
num2str(pG_FC, 3)]);
disp(' * Here''s where I would like to test for differences among JKL within Group C')
% Uncomment to check above against MATLAB anova
% disp(' ')
% disp(anova(mdl))
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/972465/image.jpeg)
4 commentaires
Scott MacKenzie
le 24 Avr 2022
Sorry, but I don't think I can be of any help here. I don't have experience doing pairwise comparisons in the manner you are describing. But, good luck.
Réponse acceptée
James Akula
le 25 Avr 2022
1 commentaire
MM
le 9 Sep 2023
Hi James,
I thought your approach could work for my problem.
I have 3 factors with each having multiple levels. Factor color has 3 levels, order has 4 levels and condition has 5 levels. I want to look at the main and interaction effect on the outcome variable, individual_correction_error. I have subjects as the random effect.
I have excluded condition_3 as there was no data there.
I would like to use coeftest to do post hoc analysis on the interaction effects (condn*order). How do I go about building the contrast matrix for sepcifically testing the interaction between condition and order?
This is the output of fitlme on the model ('individual_correction_error ~ color*condition*order+(1|subject)')
lm_full_condition =
Linear mixed-effects model fit by ML
Model information:
Number of observations 2286
Fixed effects coefficients 48
Random effects coefficients 45
Covariance parameters 2
Formula:
individual_correction_error ~ 1 + color*condition + color*order + condition*order + color:condition:order + (1 | subject)
Model fit statistics:
AIC BIC LogLikelihood Deviance
12030 12316 -5964.9 11930
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue Lower Upper
{'(Intercept)' } -0.15933 0.74743 -0.21317 2238 0.83121 -1.6251 1.3064
{'color_blue' } -0.1484 0.098592 -1.5052 2238 0.13241 -0.34174 0.044939
{'color_green' } 0.13717 0.098592 1.3912 2238 0.16429 -0.056176 0.33051
{'condition_4' } -2.8125 0.15171 -18.538 2238 1.7385e-71 -3.11 -2.515
{'condition_2' } -0.13007 0.11213 -1.16 2238 0.24619 -0.34996 0.089826
{'condition_1' } 3.0123 0.16372 18.398 2238 1.6508e-70 2.6912 3.3333
{'order_[0 1 2 4 3]' } 1.1165 0.8071 1.3834 2238 0.16669 -0.46624 2.6993
{'order_[0 2 3 1 4]' } 0.14848 1.3039 0.11387 2238 0.90935 -2.4085 2.7055
{'order_[0 3 4 2 1]' } -2.3417 1.2672 -1.8478 2238 0.064757 -4.8267 0.14343
{'color_blue:condition_4' } 0.27389 0.18359 1.4918 2238 0.13588 -0.086139 0.63392
{'color_green:condition_4' } -0.34359 0.18359 -1.8715 2238 0.061407 -0.70362 0.016437
{'color_blue:condition_2' } -0.12677 0.15211 -0.83338 2238 0.40472 -0.42507 0.17153
{'color_green:condition_2' } 0.075611 0.15211 0.49707 2238 0.61919 -0.22269 0.37391
{'color_blue:condition_1' } -0.25319 0.19023 -1.331 2238 0.18332 -0.62623 0.11985
{'color_green:condition_1' } 0.14719 0.19023 0.77377 2238 0.43915 -0.22585 0.52023
{'color_blue:order_[0 1 2 4 3]' } 0.30762 0.17161 1.7925 2238 0.073182 -0.028914 0.64416
{'color_green:order_[0 1 2 4 3]' } -0.30302 0.17161 -1.7657 2238 0.077581 -0.63956 0.033519
{'color_blue:order_[0 2 3 1 4]' } 0.10408 0.17381 0.59878 2238 0.54938 -0.23677 0.44492
{'color_green:order_[0 2 3 1 4]' } -0.013282 0.17381 -0.076417 2238 0.93909 -0.35413 0.32757
{'color_blue:order_[0 3 4 2 1]' } -0.33717 0.1713 -1.9682 2238 0.049165 -0.6731 -0.0012345
{'color_green:order_[0 3 4 2 1]' } 0.26726 0.1713 1.5601 2238 0.11887 -0.068678 0.60319
{'condition_4:order_[0 1 2 4 3]' } 0.11191 0.2537 0.44112 2238 0.65917 -0.3856 0.60943
{'condition_2:order_[0 1 2 4 3]' } 0.64807 0.20125 3.2203 2238 0.001299 0.25342 1.0427
{'condition_1:order_[0 1 2 4 3]' } -0.38955 0.2858 -1.363 2238 0.17301 -0.95002 0.17091
{'condition_4:order_[0 2 3 1 4]' } 0.78129 0.29921 2.6112 2238 0.0090828 0.19454 1.368
{'condition_2:order_[0 2 3 1 4]' } -0.72327 0.1926 -3.7553 2238 0.00017757 -1.101 -0.34557
{'condition_1:order_[0 2 3 1 4]' } -0.50564 0.25323 -1.9968 2238 0.045969 -1.0022 -0.0090558
{'condition_4:order_[0 3 4 2 1]' } -1.2315 0.24498 -5.0268 2238 5.3815e-07 -1.7119 -0.75105
{'condition_2:order_[0 3 4 2 1]' } -1.1713 0.19561 -5.988 2238 2.4683e-09 -1.5549 -0.78773
{'condition_1:order_[0 3 4 2 1]' } 2.3455 0.29693 7.8989 2238 4.3708e-15 1.7632 2.9278
{'color_blue:condition_4:order_[0 1 2 4 3]' } 0.12839 0.3146 0.4081 2238 0.68324 -0.48854 0.74532
{'color_green:condition_4:order_[0 1 2 4 3]'} 0.14677 0.3146 0.46654 2238 0.64088 -0.47016 0.7637
{'color_blue:condition_2:order_[0 1 2 4 3]' } -0.10329 0.2706 -0.3817 2238 0.70272 -0.63394 0.42737
{'color_green:condition_2:order_[0 1 2 4 3]'} 0.14738 0.2706 0.54464 2238 0.58606 -0.38328 0.67804
{'color_blue:condition_1:order_[0 1 2 4 3]' } 0.17675 0.3325 0.53158 2238 0.59507 -0.47529 0.82879
{'color_green:condition_1:order_[0 1 2 4 3]'} -0.22094 0.3325 -0.66448 2238 0.50645 -0.87298 0.4311
{'color_blue:condition_4:order_[0 2 3 1 4]' } -0.033804 0.34999 -0.096585 2238 0.92306 -0.72015 0.65254
{'color_green:condition_4:order_[0 2 3 1 4]'} -0.20315 0.34999 -0.58045 2238 0.56167 -0.8895 0.48319
{'color_blue:condition_2:order_[0 2 3 1 4]' } 0.2268 0.26636 0.85149 2238 0.39459 -0.29554 0.74914
{'color_green:condition_2:order_[0 2 3 1 4]'} -0.13989 0.26636 -0.52518 2238 0.59951 -0.66223 0.38245
{'color_blue:condition_1:order_[0 2 3 1 4]' } 0.30113 0.30137 0.9992 2238 0.3178 -0.28987 0.89213
{'color_green:condition_1:order_[0 2 3 1 4]'} -0.089843 0.30137 -0.29811 2238 0.76565 -0.68084 0.50116
{'color_blue:condition_4:order_[0 3 4 2 1]' } 0.36974 0.30391 1.2166 2238 0.22388 -0.22623 0.96571
{'color_green:condition_4:order_[0 3 4 2 1]'} -0.021389 0.30391 -0.07038 2238 0.9439 -0.61736 0.57458
{'color_blue:condition_2:order_[0 3 4 2 1]' } -0.13751 0.26 -0.52888 2238 0.59694 -0.64738 0.37236
{'color_green:condition_2:order_[0 3 4 2 1]'} -0.11538 0.26 -0.44375 2238 0.65726 -0.62524 0.39449
{'color_blue:condition_1:order_[0 3 4 2 1]' } -0.73679 0.35229 -2.0914 2238 0.036604 -1.4276 -0.045935
{'color_green:condition_1:order_[0 3 4 2 1]'} 0.48379 0.35229 1.3732 2238 0.16981 -0.20707 1.1746
Random effects covariance parameters (95% CIs):
Group: subject (45 Levels)
Name1 Name2 Type Estimate Lower Upper
{'(Intercept)'} {'(Intercept)'} {'std'} 4.9857 4.0418 6.1499
Group: Error
Name Estimate Lower Upper
{'Res Std'} 3.1361 3.0456 3.2293
Plus de réponses (0)
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!