title: Pairwise post-hoc testing using coefTest (MATLAB)
body: <p>I found <a href="https://www.mathworks.com/matlabcentral/answers/185041-how-can-i-run-post-hoc-tests-for-a-mixed-effects-model-using-fitlme" rel="nofollow noreferrer">this post</a> 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.</p>
  <p>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.</p>
  <p>My mock dataset can be found <a href="https://1drv.ms/u/s!AnVfp9IB0y_GnaRNT5x5y_Z1cgB9PQ?e=gXuC9f" rel="nofollow noreferrer">here</a>.</p>
  <p>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.</p>
  <p>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?</p>
  <p>Thank you!!!</p>
  <pre><code>    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
  
  %     [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)])
  
  %        [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))
  </code></pre>
  <p><a href="https://i.stack.imgur.com/46QS3.jpg" rel="nofollow noreferrer">Some helpful output to visualize expected outcomes</a></p>