Pairwise post-hoc testing using coefTest

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))

[Scott MacKenzie]

There's a lot to digest here. I'm wondering why you aren't using an analysis of variance. The description of your study sounds like a factorial experiment and you speak about looking for significant main effects, interaction effects, etc. So, doing an analysis of variance is the norm.

As for pairwise post hoc testing, the usual way to do this is in MATLAB is with the multcompare function. Alternatively, you can use t-tests and apply the Bonferroni correction to maintain the desired family-wise alpha.

[James Akula]

Hi Scott,

Thank you for taking the time to answer. Apologies if I packed too much into my qeustion. It is actually a fairly simple one (see bolded text at end of post). However, first, let me just assert that I want to use an LME model because GLM (in this case ANOVA) doesn't support the designs I am interested in (e.g., mixed longitudinal). As to "main effects" and "interaction effects," I am interested in the LME equivalents. In my example above, if I run ANOVA in on the LME I get:

>> ANOVA = anova(mdl)

ANOVA =

ANOVA marginal tests: DFMethod = 'Residual'

Term FStat DF1 DF2 pValue

{'(Intercept)' } 12.467 1 479 0.00045435

{'Group' } 17.425 2 479 4.9564e-08

{'FactorA' } 1.4629 2 479 0.2326

{'FactorB' } 0.74766 2 479 0.47402

{'FactorC' } 12.581 2 479 4.7285e-06

{'Group:FactorA'} 1.365 4 479 0.24502

{'Group:FactorB'} 1.0619 4 479 0.37481

{'Group:FactorC'} 8.1556 4 479 2.2822e-06

I recognize that ANOVA isn't perfectly appropriate, but it is an "accepted" way to get stats out of an LME model (I could use Satterthwaite approximation for the degrees of freedom, if it helps). As you will note, I get a significant effect of Group (p=0.000000050), a significant effect of Factor C (p=0.0000047) and a significant Group×Factor C interaction (p=0.0000023), which is what I would expect from my data, because I designed it to be that way (I actually didn't intend for there to be an effect of Factor C, but its clear from looking at the plots that it should, in the end, be in there). By the way, I know how to produce the p values for each of the effects shown in the ANOVA table using coefTest.

Notably, I also know how to look at the group effect and, noting that it is significant, perform a set of post-hoc comparisons to determine which groups differ from which other groups; those p values are stored in pAB, pAC, and pBC. As expected, A and B do not differ, but both groups differ from C.

Stated simply, my question is, noting the significant Group×Factor C interaction, how can I test, pairwise, which levels of Factor C differ within Group C. That is, how can I determine, within Group C, values for pJK, pJL, and pKL?

[Scott MacKenzie]

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.

[James Akula]

Thanks. I do beleive I solved it.

Answer 1

[James Akula]

I do believe I answered the question. Adding the following code to the problem seems to give me the answer I want. Uncommenting the symbolic bits leads to the answer.

% To do a pairwise comparison deep in the interaction effect gets tricky.
% Here, we need to figure out what the comparisons should be for a nested
% interaction effect.  However, we can solve for the needed factors using
% the same approach as above:
% syms Group_A_FactorC_J Group_B_FactorC_J Group_C_FactorC_J ...
%      Group_A_FactorC_K Group_B_FactorC_K Group_C_FactorC_K ...
%      Group_A_FactorC_L Group_B_FactorC_L Group_C_FactorC_L
% 
% Group_C_FactorC_J = -Group_A_FactorC_J - Group_B_FactorC_J;
% Group_C_FactorC_K = -Group_A_FactorC_K - Group_B_FactorC_K;
% Group_C_FactorC_L = -Group_C_FactorC_J - Group_C_FactorC_K;
% 
% fGC_FC_JK = Group_C_FactorC_J - Group_C_FactorC_K
% fGC_FC_JL = Group_C_FactorC_J - Group_C_FactorC_L
% fGC_FC_KL = Group_C_FactorC_K - Group_C_FactorC_L

% The answers are:
% fGC_FC_JK = -Group_A_FactorC_J   - Group_B_FactorC_J   + Group_A_FactorC_K   + Group_B_FactorC_K
% fGC_FC_JL = -2*Group_A_FactorC_J - 2*Group_B_FactorC_J - Group_A_FactorC_K   - Group_B_FactorC_K
% fGC_FC_KL = -Group_A_FactorC_J   - Group_B_FactorC_J   - 2*Group_A_FactorC_K - 2*Group_B_FactorC_K

% 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]
HGC_FC_JK = [  0,   0,   0,    0,    0,    0,    0,    0,    0,        0,        0,        0,        0,        0,        0,        0,        0,       -1,       -1,        1,        1];
HGC_FC_JL = [  0,   0,   0,    0,    0,    0,    0,    0,    0,        0,        0,        0,        0,        0,        0,        0,        0,       -2,       -2,       -1,       -1];
HGC_FC_KL = [  0,   0,   0,    0,    0,    0,    0,    0,    0,        0,        0,        0,        0,        0,        0,        0,        0,       -1,       -1,       -2,       -2];

pGC_FC_JK = coefTest(mdl, HGC_FC_JK);
pGC_FC_JL = coefTest(mdl, HGC_FC_JL);
pGC_FC_KL = coefTest(mdl, HGC_FC_KL);

disp([' * Significance of Factor C-J vs. Factor C-K difference within Group C: ', ...
    num2str(pGC_FC_JK, 3)])
disp([' * Significance of Factor C-J vs. Factor C-L difference within Group C: ', ...
    num2str(pGC_FC_JL, 3)])
disp([' * Significance of Factor C-K vs. Factor C-L difference within Group C: ', ...
    num2str(pGC_FC_KL, 3)])

[MM]

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