How can I run post hoc tests for a mixed-effects model using fitlme?

Hi,

I have a model with several independent categorical variables. Using fitlme and anova, I find significant effects. However, I do not understand how to run post hoc tests, such as which level was sig. different within a variable, or which part of the interaction was significant.

I see that fitlm has a multcompare function. I would want to do something similar. Or at the very least, run contrast tests. I know coefTest is supposed to be able to test contrasts, however, I'm really unclear on how to set up a contrast that would correspond to one of the above tests. I'd really appreciate it if someone could walk me through a contrast for different levels (let's say out of 3 levels) for a particular variable, or for a 2x2 interaction.

Below is code for one of the models that I ran, as an example. All variables are binary except OrdValence, which has 3 levels. I'm particularly interested in an OrdValence by LikingVSWanting interaction.

Thanks in advance!! Galit

----

maxModelStr = ['DecisionNum ~ Sex * LikingVSWanting * ArmAction * OrdValence + (1|Subject:Sex) + (1|Subject:LikingVSWanting) +'... '(1|Subject:ArmAction) + (1|Subject:OrdValence) + (1|Subject:LikingVSWanting:Sex) + (1|Subject:LikingVSWanting:ArmAction) +'... '(1|Subject:LikingVSWanting:OrdValence) + (1|Subject:Sex:ArmAction) + (1|Subject:Sex:OrdValence) + (1|Subject:ArmAction:OrdValence) +'... '(1|Subject:Sex:LikingVSWanting:ArmAction) + (1|Subject:Sex:LikingVSWanting:OrdValence) + (1|Subject:LikingVSWanting:ArmAction:OrdValence)']; maxModel = fitlme(raw_data, maxModelStr,'FitMethod', 'ML');

anova_max = anova(maxModel,'DFMethod', 'Satterthwaite');

[Alexis]

Given how popular GLME models are, it would be really great if TMW addressed the poor documentation on this topic. I have cobbled together examples from here, Stack Overflow, etc., and it still isn't really clear how to use coefTest. Most of the related questions on MATLAB Answers remain unresolved.

Answer 1

[Pär Nyström]

Dear Galit,

I am resurrecting this question, because I think many people want to know the answer, and because I have to give up on MATLAB and statistics if this is not solved. So let's hope that this will lead to good discussion... While I am not completely confident that this is the correct approach, and I want others to justify it, I think that it could be a step closer to what you want. I'm fitting a model with fitlme on the form "relconst ~ group*site + (1|subject)", which is much less complex than your model, where the group factor has 3 levels, and site 2 levels. I have pasted the stats output below.

I'm using coefTest to do post-hoc comparisons on the significant group effect, and by using different contrasts I can recreate the p-values found using the "anova" function and almost all post hoc tests in JASP (free stats software) with bonferroni correction, except for the HR-ASD vs HR-no-ASD which is weirdly very different from JASP. However, SPSS will give other p-values for the post-hoc tests if there are more than one categorical factor, such as "site", suggesting that there are some differing modelling parameters between statistics software...

Anyway, I used the code below to do the works:

        mdl1 = fitlme(tbl, mdlstring1, 'DummyVarCoding', 'effects', 'FitMethod', 'REML');         disp(mdl1);
        fprintf('ANOVA:s\n');
        disp(anova(mdl1));
        for y = 1:numel(mdl1.CoefficientNames)
            fprintf('Contrast position %i: %s\n', y, char(mdl1.CoefficientNames{y}));
        end
        bonfCorrection = 1;
        [pVal F df1 df2] = coefTest(mdl1, [0  0  1 0 0 0]);   pVal = min(1,pVal*bonfCorrection);           fprintf('Planned comparison HR-ASD vs HR-no-ASD,  F(%i, %i)=%0.2f,\t p=%0.7f (contrast 0  0  1 0 0 0)\n', df1, df2, F, pVal);
        [pVal F df1 df2] = coefTest(mdl1, [0  1  0 0 0 0]);   pVal = min(1,pVal*bonfCorrection);           fprintf('Planned comparison HR-ASD vs TD,         F(%i, %i)=%0.2f,\t p=%0.7f (contrast 0  1  0 0 0 0)\n', df1, df2, F, pVal); 
        [pVal F df1 df2] = coefTest(mdl1, [0  1 -1 0 0 0]);   pVal = min(1,pVal*bonfCorrection);           fprintf('Planned comparison TD vs HR-no-ASD,      F(%i, %i)=%0.2f,\t p=%0.7f (contrast 0 -1  1 0 0 0)\n', df1, df2, F, pVal);
        [pVal F df1 df2] = coefTest(mdl1, [0  0  0 0 0 1; 0 0 0 0 1 0]); fprintf('INTERACTION,                             F(%i, %i)=%0.2f,\t p=%0.7f (contrast 0 0 0 0 0 1; 0 0 0 0 1 0)\n', df1, df2, F, pVal);
        [pVal F df1 df2] = coefTest(mdl1, [0  1  0 0 0 0; 0 0 1 0 0 0]); fprintf('MAIN GROUP,                              F(%i, %i)=%0.2f,\t p=%0.7f (contrast 0  1 0 0 0 0; 0 0 1 0 0 0)\n', df1, df2, F, pVal);
        [pVal F df1 df2] = coefTest(mdl1, [0  0  0 1 0 0]);              fprintf('MAIN SITE,                               F(%i, %i)=%0.2f,\t p=%0.7f (contrast 0 0 0 1 0 0)\n', df1, df2, F, pVal);

This will generate the output below, and you can see that the coefTest on these contrasts correspond with the anova (for the main effects and interaction) and the post-hoc tests correspond to the output from the linear model. If you want bonferroni correction for the post-hoc tests you can change the bonfCorrection parameter from 1 to 3 (actually, the number of levels in your factor, which is 3 in my case).

Linear mixed-effects model fit by REML

Model information:
    Number of observations             187
    Fixed effects coefficients           6
    Random effects coefficients          0
    Covariance parameters                1

Formula:
    medianRelconst10 ~ 1 + groupname*site

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    42.759    65.148    -14.379          28.759  

Fixed effects coefficients (95% CIs):
    Name                                   Estimate     SE          tStat     DF     pValue         Lower        Upper    
     '(Intercept)'                             1.0969    0.022045     49.756    181    1.6013e-107       1.0534       1.1404
     'groupname_TD'                          -0.10793    0.031778    -3.3963    181     0.00083963     -0.17063    -0.045225
     'groupname_HR-no-ASD'                   0.014524    0.027115    0.53564    181        0.59286    -0.038978     0.068026
     'site_EASE'                              0.11274    0.022045     5.1143    181     7.9731e-07     0.069246      0.15624
     'groupname_TD:site_EASE'               -0.056058    0.031778     -1.764    181       0.079412     -0.11876    0.0066453
     'groupname_HR-no-ASD:site_EASE'         0.027101    0.027115    0.99947    181         0.3189    -0.026401     0.080603

Random effects covariance parameters (95% CIs):
Group: Error
    Name             Estimate    Lower    Upper  
    'Res Std'        0.24277     0.219    0.26911

ANOVA:s

    ANOVA marginal tests: DFMethod = 'Residual'

    Term                    FStat     DF1    DF2    pValue     
    '(Intercept)'           2475.7    1      181    1.6013e-107
    'groupname'             5.9674    2      181      0.0030921
    'site'                  26.157    1      181     7.9731e-07
    'groupname:site'        1.6512    2      181         0.1947

Contrast position 1: (Intercept)
Contrast position 2: groupname_TD
Contrast position 3: groupname_HR-no-ASD
Contrast position 4: site_EASE
Contrast position 5: groupname_TD:site_EASE
Contrast position 6: groupname_HR-no-ASD:site_EASE

Planned comparison HR-ASD vs HR-no-ASD,  F(1, 181)=0.29,   p=0.5928635 (contrast 0  0  1 0 0 0)
Planned comparison HR-ASD vs TD,         F(1, 181)=11.53,   p=0.0008396 (contrast 0  1  0 0 0 0)
Planned comparison TD vs HR-no-ASD,      F(1, 181)=6.46,   p=0.0118413 (contrast 0 -1  1 0 0 0)
INTERACTION,                             F(2, 181)=1.65,   p=0.1946977 (contrast 0 0 0 0 0 1; 0 0 0 0 1 0)
MAIN GROUP,                              F(2, 181)=5.97,   p=0.0030921 (contrast 0  1 0 0 0 0; 0 0 1 0 0 0)
MAIN SITE,                               F(1, 181)=26.16,   p=0.0000008 (contrast 0 0 0 1 0 0)

I hope this provides a concrete example on how post-hoc tests can be done, but beware the fact that the HR-ASD vs HR-no-ASD test was not replicated in another statistics software, and I don't know whether MATLAB or JASP is wrong. Or if this should be done differently.

Btw, did you solve this yourself? If so, I would be happy to know how you did!

EDIT: I just realized that is is only the post-hoc contrast [0 -1 1 0 0 0] that always gives the right p-value is you reorder the levels :( So I guess that we need to refit the model with different level orders. In the code below I commented out all but one order, and ran the code three times... The correct contrast coding should probably be done with 'DummyVarCoding', 'full', but I don't get that to work. So I will suggest going for the reordering of levels...:

% tbl.groupname = reorderlevels(nominal(tbl.groupname), {'TD', 'HR-ASD', 'HR-no-ASD'}); % tbl.groupname = reorderlevels(nominal(tbl.groupname), {'HR-ASD', 'HR-no-ASD', 'TD'}); tbl.groupname = reorderlevels(nominal(tbl.groupname), {'HR-no-ASD', 'TD', 'HR-ASD'});

[David Lewis]

I think your contrast for HR-ASD vs HR-no-ASD should have been [0 1 2 0 0 0]. One of your compared levels was the reference level (HR-ASD). So if you're using effects dummy coding (which you are), then the level (HR-no-ASD) that you're comparing against the reference gets a 2 in the contrast, while other uncompared levels (TD in this case) of the same treatment factor gets a 1. By the same logic, your contrast for HR-ASD vs. TD should have been [0 2 1 0 0 0]. I'm really novice at this, so take me with a huge grain of salt. But this is how I read the description of how they set up the "supplier A vs. supplier B" contrast in this mathworks example: https://www.mathworks.com/help/stats/generalizedlinearmixedmodel.coeftest.html#namevaluepairarguments

Answer 2

[John Hartman]

There is no need to fit multiple models for post-hoc tests involving reference levels of predictor variables, just define the contrasts carefully. You can do this using coefTest but it isn't explained well enough in the documentation for generalized linear mixed effect models (at least for complicated cases). It's usually easier to consider contrasts of the estimated marginal means, which you can do with this, https://www.mathworks.com/matlabcentral/fileexchange/71970-emmeans

[David Lewis]

I downloaded your emmeans package, John. Thanks for writing it. Is there any way to get a full matrix of pairwise comparisons? In your demo.m file, for example the H0_usa just gives a single p-value output.

[John Hartman]

Hi David,

In the demo, H0_usa is jointly testing all of those different hypotheses (i.e. testing "is at least one of these contrasts significant?") which is why you only get a single p-value.

To get the full set of pairwise comparisons, you would have to loop through each contrast on its own. Here's some code to do that, assuming you're interested in pairs of the estimated marginal means.

numberofbetas = 5;
fullpairs = nchoosek(1:numberofbetas,2);
for ii=1:size(fullpairs,1)
    L = zeros(1,numberofbetas);
    L(fullpairs(ii,1)) = 1;
    L(fullpairs(ii,2)) = -1;
    
    H0 = contrasts_wald(mdl,emm,L);
end

If you are testing the full set of pairwise comparisons, it sounds like you might not have specified your hypotheses before your analysis, and therefore you probably want to do some sort of multiple comparison correction, e.g. Bonferroni, Tukey, Holm, etc.

Best,

Jack

[Alexis]

Dear Jack, I'd like to use emmeans(), but I'm not totally clear on how to run post-hoc tests on more complex interactions (e.g., 3-way).

I'm using dummy variable coding ('effects') in my model -- am still a bit confused by ML's documentation for coefTest() on this -- so if you can provide a relevent example, I'd really appreciate that.

[John Hartman]

Hi Aldema, Did you check out the demo.m script in the package? The purpose of the emmeans package is essentially not to use coefTest because, as you see, it's very confusing and poorly documented.

[John Hartman]

The data here sort of messes up how this works but this code should sort of explain how to think about both approaches.

%% For 2-way interactions
twowaycombos = permn(1:3,2);
datatable.indep_var1 = categorical( twowaycombos(:,1), 1:3, {'Ham' 'Egg' 'Cheese'} );
datatable.indep_var2 = categorical( twowaycombos(:,2), 1:3, {'Red' 'Blue' 'Green'} );
datatable.dep_var = randperm(9)';
datatable = struct2table(datatable);
mdl = fitglme( datatable,'dep_var~indep_var1*indep_var2', 'DummyVarCoding', 'effects');

% using MATLAB method, test difference between RedHam and BlueEgg
% to get these vectors you match to the betas from mdl
RedHam = [1 2 1 2 1 2 1 1 1]; % for effects coding
% RedHam = [1 1 0 1 0 1 0 0 0]; % for reference coding
BlueEgg = [1 1 2 1 2 1 1 1 2]; % for effects coding
% BlueEgg = [1 0 1 0 1 0 0 0 1]; % for reference coding
L = RedHam-BlueEgg; % contrast
H = coefTest(mdl,L); 

% using emmeans
% to get these vectors you match to the rows from emm.table
emm = emmeans( mdl, 'effects' );
L = [1 0 0 0 -1 0 0 0 0]; 
% L = [ (ismember(emm.table.indep_var1,'Ham') & ismember(emm.table.indep_var2,'Red'))... 
%     -(ismember(emm.table.indep_var1,'Egg') & ismember(emm.table.indep_var2,'Blue')) ]; % another way to get same contrast
H = contrasts_wald(mdl,emm,L);

Hope that helps!

[Alexis]

Thanks so much! I will give this a try and provide feedback at your File Exchange entry.

[Adam Ranson]

Hi John,

The original question above was regarding contrasts in linear mixed effects models from fitlme() - does your emmeans function work with these? It appears to be designed for generalized linear mixed effects models fit using fitglme()?

Thanks a lot

Adam