Solved – Odds ratio confidence intervals and p-values suggest different conclusions in a binary logistic mixed-effects model (glmer)

confidence intervallme4-nlmemixed modelodds-ratior

I am running a generalized linear mixed-effects model in R using the glmer function of lme4. The outcome variable is trial-level accuracy in a task (incorrect trials are 0, correct trials are 1), and I have 3, trial-level predictors:

A, a binary categorical variable
B, another binary categorical variable
C, a categorical variable with 3 levels

(All variables, including accuracy are factors.)

Here is the model I'm running:

    acc_model <- glmer(accuracy ~ 1 + A*B*C + (1|subid), data = acc_data, family = binomial,
                   control = glmerControl(optimizer ="bobyqa"))
    #subid simply refers to subject IDs as trials are nested within participants
    #I use bobyqa because the model failed to converge with the default optimizer

Then I use the Anova() function from the "car" package to see if the terms are significant.

library(car)
Anova(acc_model)

This has worked nicely for me in the past (anova() does not give p-values for glmer models), and the results seem sensible in this case too:

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: acc
                              Chisq Df Pr(>Chisq)    
A                           10.0468  1   0.001526 ** 
B                           18.6358  1  1.582e-05 ***
C                            6.7366  2   0.034449 *  
A:B                          4.5702  1   0.032532 *  
A:C                          0.0142  2   0.992915    
B:C                          0.6303  2   0.729685    
A:B:C                        2.8599  2   0.239319

The main effects of A, B, C are sig, and so is the A:B two-way interaction.

However, when I calculate odds ratios and corresponding 95% confidence intervals for these effects, the CIs ALL include 1.

This is the code I use to get ORs and CIs:

#CI and OR for variable A
exp(summary(acc_model)$coefficients["A1",1] + qnorm(c(0.025,0.5,0.975)) * summary(acc_model)$coefficients["A1",2])

What can I do here? Is it okay to interpret the Anova() output? I suspect the significance of main effects A, B, and C in the Anova() output could be due to that being a Type II ANOVA (i.e., interactions including a given variable are not controlled for when the main effect of that variable is being evaluated, if I get that correclty), however, the CI for the OR of the A:B interaction ALSO contains 1, and that is NOT involved in any significant higher order interactions.

NOTE: when I use lmerTest, and get summary(acc_model), neither of these variables is significant. I would like to avoid using lmerTest, however, because that does not give an omnibus p-value for categorical variables with more than 2 levels (e.g., Variable C here, that would have 2 separate lines for the Reference_Category – Level1 and the Reference_Category – Level2 contrasts), so anova() and Anova() are neater.

Best Answer

The ANOVA is Type II sum of squares and the method you used below is Type III sum of squares. In general, they don't agree. Models for independent data give the same results. For mixed models, they differ. Type II and Type III sum of squares handle the random effect differently. The Type III method assumes, (incorrectly in most cases) that if the "true" data generating fixed effect were different, such as if it were at the lower end of the 95% confidence interval, that the random effect would be unchanged. This is well and good provided the random effects and fixed effects aren't correlated. But they usually are, and the default empirical variance estimates do not suffice to summarize variability.

If you use the confint method in R, there are three confint procedures for mixed models (merMod objects in the nlme package). Wald (Type III, wrong), profile, and bootstrap. Try profile and bootstrap and report back the findings. Note they still may not agree with Type II SS, but the reasons are more refined.

https://www.rdocumentation.org/packages/lme4/versions/1.1-21/topics/confint.merMod

Related Solutions

Solved – Mixed, repeated measure model specification and results interpretation using LMER in R

This mainly pertains to your first question.

You want to specify a "Dx interaction with any of the above variables", and seem to think this is done with $Dx*var$. But it is not, $Dx*var$ expands to $Dx + var + Dx:var$. The $Dx:var$ specifies the interaction.

So let's change your model to

Correct ~ Dx + No_of_Stim + Trial_Type + Probe_Loc + Dist +
          Dx:No_of_Stim + Dx:Trial_Type + Dx:Probe_Loc + Dx:Dist +
          (1 | Trial) + (1 | SubID)

(I'm not sure what lmer made out of your model specification.) So this fits a model with the fixed predictors you asked for, plus it allows for a individual intercept for each Trial, and additionally for each SubID.

Because I want to also examine a continuous variable, the total Cartesian Distance > between stimuli per trial (divided by number of stimuli to control for varying > numbers), I have opted to use a mixed linear model with repeated measures.

I don't get that. So you would like to include a continuous predictor - why does that make it a repeated measures design? Do you measure the same subject on the same tasks more than once? On differing tasks more than once?

3) If I attempt to reduce the model, I should favor the model with the lower REML, yes?

There are various ways to compare two models. The simplest is to set REML=FALSE for both, and compare them via a likelihood ratio test by anova(mod1, mod2). Another way is to use information criteria such as the AIC (Akaike information criterion) or BIC (Bayesian information criterion). Those try to balance model fit with model complexity, and a lower number indicates a "better" model. I don't think it's a good idea to use the REML criterion as a model comparison metric, afaik this is, at least in general, not valid (perhaps someone more knowledgeable can provide an explanation for that).

Update after ADDENDUM

Remove $(1 | Trial)$ from the random effects specification (unless you intend to model temporal correlation between trial runs; but even then you would need something other than the running trial index number as a grouping variable). It's pointless to include this as a grouping variable, since you only have one observation per level:

Correct ~ Dx + No_of_Stim + Trial_Type + Probe_Loc + Dist +
          Dx:No_of_Stim + Dx:Trial_Type + Dx:Probe_Loc + Dx:Dist +
          (1 | SubID)

This has the minimal random effects structure justified by your experiment: It fits a model with the fixed effects predictors, and adds an additional "random" (read: individual) intercept per individual subject, effectively shifting the regression line determined by the fixed effects predictors up and down by a subject-specific amount. The fact that subjects were repeatedly measured is thereby accounted for.

You can build upon this model significantly. The maximal random effects structure would allow for correlated random intercepts and slopes of all within-subject predictors, grouped by subject. Barr et al., (2013) recommend, for confirmatory hypothesis testing of fixed effects, to include for all tested fixed effects correlated random intercepts and slopes per subject.

Presuming you intend to test all your within-subjects variables number of stimuli, trial type, probe location, distance, and omit any interactions between them, this would be

Correct ~ Dx + No_of_Stim + Trial_Type + Probe_Loc + Dist +
          Dx:No_of_Stim + Dx:Trial_Type + Dx:Probe_Loc + Dx:Dist +
          (No_of_Stim + Trial_Type + Probe_Loc + Dist | SubID)

where intercepts and slopes for No_of_Stim, Trial_Type, Probe_Loc, and Dist are allowed to correlate. A simpler model, imposing independence between intercept and slope, would be

Correct ~ Dx + No_of_Stim + Trial_Type + Probe_Loc + Dist +
          Dx:No_of_Stim + Dx:Trial_Type + Dx:Probe_Loc + Dx:Dist +
          (1 | SubID) + (0 + No_of_Stim | SubID) + (0 + Trial_Type | SubID) +
          (0 + Probe_Loc | SubID) + (0 + Dist | SubID))

Relevant google terms for finding "the best" model are model building, top-down strategy, bottom-up strategy, but there is a lot of controversy regarding the best strategy.

Regarding interpretation of the coefficients in the output

By default, R uses treatment contrasts (see $help(contr.treatment)$ and http://talklab.psy.gla.ac.uk/tvw/catpred/). This means that one level of a factor is chosen as the reference factor, and the coefficients indicate the contrast to this reference level. If you are instead asking about the interpretation of the coefficients of fixed within-subject predictors that are additionally allowed to vary randomly per subject: It's essentially the same as in a simple linear model. They tell you about whether there's any systematic effect across subjects of this predictor on the outcome, after accounting for random variation by subjects.

You need to use a logistic regression model!

The elephant in the room, so to speak, is that you are trying to predict a binary outcome: correct/incorrect. You need to use a logistic mixed regression model (try $glmer$ and look at http://www.ats.ucla.edu/stat/r/dae/melogit.htm).

So, I hope this helped you a bit. Maybe also have a look at http://glmm.wikidot.com/faq and R's lmer cheat-sheet.

(I am by no means an expert on this, and in a state of learning as well, so if there's something to disagree, please do so freely.)

Solved – Different p-values for fixed effects in summary() of glmer() and likelihood ratio test comparison in R

It looks like you are seeing the difference between Wald p-values (based on the curvature of the log-likelihood surface at the maximum likelihood estimate) and likelihood ratio test p-values (based on comparisons between the full and reduced models).

take a look at tpr <- profile(a25,which="beta_"); lattice::xyplot(tpr). You should see that the lines are far from straight (straight lines would indicate a log-quadratic likelihood surface, which is what's assumed by Wald p-values)
compare the results of confint(a25,which="beta_") (likelihood ratio intervals) and confint(a25,which="beta_",method="Wald"); they should be quite different.

LRT CI/p-values are essentially always better than the Wald equivalents (but much slower to compute, which is why Wald p-values are the default in summary()).

Best Answer

Related Solutions

Solved – Mixed, repeated measure model specification and results interpretation using LMER in R

Solved – Different p-values for fixed effects in summary() of glmer() and likelihood ratio test comparison in R

Related Question