# crazy odds ratios with mixed-effects logistic regression

I just run two mixed-effects logistic regressions. I narrowed down the best models using a stepwise manual backwards method and anova comparisons

The two regressions measure two different kinds of linguistic productivity: morphological (morph_num: 0 or 1) and syntactical (synt_num, 0 or 1).

The main predictor are : tvl_scaled: a re-scaled version of participants' vocabulary score(it was necessary for the model to be computationally feasible).

Wo: either “vo” (0) or “vs” (1)

Verb: “test”(1) vs “control” (0).

Design

Each participant got a trial with one test and one control verb, each verb presented with either wo. That is either “control+vs and test+vo” or “control+vo and test+vs”. Hence verb is within participants.

As previously suggested to me here I run models using nAGQ=10

The two results are presented below. I shall refer to the first output as morph and to the second as synt.

MORPH

Generalized linear mixed model fit by maximum likelihood (Adaptive
Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
Family: binomial  ( logit )
Formula:
morph_num ~ tvl_scaled + verb + wo + tvl_scaled:verb + (1 | participants)
Data: opz

AIC      BIC   logLik deviance df.resid
92.3    109.1    -40.2     80.3      114

Scaled residuals:
Min      1Q  Median      3Q     Max
-2.8857  0.0796  0.1717  0.2693  1.5567

Random effects:
Groups       Name        Variance Std.Dev.
participants (Intercept) 2.325    1.525
Number of obs: 120, groups:  participants, 67

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)             1.738      2.649   0.656   0.5117
tvl_scaled              3.321      2.667   1.245   0.2131
verb[T.t]              -9.235      4.588  -2.013   0.0441 *
wo[T.vs]               -2.179      1.068  -2.040   0.0413 *
tvl_scaled:verb[T.t]    9.887      5.422   1.824   0.0682 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
(Intr) tvl_sc vr[T.] w[T.v]
tvl_scaled  -0.755
verb[T.t]   -0.721  0.392
wo[T.vs]    -0.595  0.044  0.551
tvl_sc:[T.]  0.694 -0.456 -0.981 -0.483


SYNT

Generalized linear mixed model fit by maximum likelihood (Adaptive
Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
Family: binomial  ( logit )
Formula: synt_num ~ tvl_scaled + verb + wo + tvl_scaled:wo + (1 | participants)
Data: opz

AIC      BIC   logLik deviance df.resid
115.7    132.5    -51.9    103.7      114

Scaled residuals:
Min      1Q  Median      3Q     Max
-2.8500 -0.3247  0.1830  0.4142  1.4667

Random effects:
Groups       Name        Variance Std.Dev.
participants (Intercept) 1.522    1.234
Number of obs: 120, groups:  participants, 67

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)          -10.344      4.146  -2.495  0.01260 *
tvl_scaled            16.958      6.286   2.698  0.00698 **
verb[T.t]             -2.647      1.075  -2.462  0.01380 *
wo[T.vs]              10.361      4.643   2.231  0.02565 *
tvl_scaled:wo[T.vs]  -13.569      6.051  -2.242  0.02494 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
(Intr) tvl_sc vr[T.] w[T.v]
tvl_scaled  -0.981
verb[T.t]    0.600 -0.720
wo[T.vs]    -0.930  0.928 -0.622
tvl_sc:[T.]  0.928 -0.938  0.631 -0.992


PROBLEM: estimates are ridiculous and odds ratios get up to 9,795,948,000,000,000,000, which is nonsense.

I tried to do a bit of diagnosis, but I am quite confused.

COOKS DISTANCE

I calculate the cooks distance for both morph and synt (I report example for morph only) and check whether any case had a cooks distance higher than 0.7 (instead of 1) with the following codes

influence.morph=influence(morph4, obs = TRUE)

cooks.morph=cooks.distance(influence.morph, sort=TRUE)

check.morph.cooks = cooks.morph > .7


none of the participants had a cooks distance higher than 0.7 (that is, they’re ok)

sum (check.morph.cooks)
0


DFBETA VALUES

I calculated dfbeta for each cell (120 participants for 4 predictor + intercept = 600 cells) with

df.beta.morph = dfbetas(influence.morph, sort=FALSE, to.sort=NA, abs=FALSE)


I determined a cutoff point of 2/srt(n) (Belsley, Kuh, and Welsch). Hence 2/( (sqr(120) ) = .18. Turns out that there are cells whose dfbeta is higher than 0.18.

check.df.beta.morph = df.beta.morph > .18

sum (check.df.beta.morph)
25
Sum (check.df.beta.synt)
32


Could those dfbeta values be the problem?

Now, since verb is the between-participant condition, I can run regular logistic regressions on data on the control and test verb (morph+control, morph+test, synt+control and synt+test), using the non-scaled vocabulary and wo as predictors. When I run them, I get odds ratios (and estimates, of course) that are realistic, whose CIs are ok.

If I check at the correlation tabs of my mixed-effects models, there are variables that are highly correlated with the verb condition (it makes no sense to me, though). I know that multicollinearity can affect B’s estimates. Do I get these results because of multicollinearity? This would explain why I get normal results with the separate logistic regressions (verb is not factored in and therefore there is no correlation between predictors).

Could anyone help me understand what’s going on. Please note that I am clearly not a statistician, so if you could avoid jargon (as far as possible), it’d be greatly appreciated.

ps. It doesn't seem that changing optimizer makes much difference

• Does the same phenomenon occur with the usual glm() logistic regresion? (I woulkd guess that). Then search this site for "logistic regresion complete separation" and " Hauck-Donner" – kjetil b halvorsen Dec 16 '16 at 18:41
• yes, but if I run two separate regressions for the two different verb condition, results are ok. I believe because the control nearly perfectly predicts 1 in the DV. – Pietre Dec 16 '16 at 22:46
• Separation is not really a problem, necessarily, the parameter estimates explodes to infinity but that is because some estimated probabilities becomes 0/1. The problem is with inferencve on the parameters, if you only needs predictions there is no problem. – kjetil b halvorsen Dec 16 '16 at 22:49
• thank you for the observation. Wanted I to use such a point in order to justify my results, what sources could I look up and quote? – Pietre Dec 16 '16 at 22:56
• Brian Ripley}s book on Neural Networks says something like within computer science community (today probably would have said Machine learners) separation was seen as a good property, not as a problem, because it permits very good discrimination ... Then the problem of doing inference is maybe secondary. Dont remember page numbers, but you will find it. – kjetil b halvorsen Dec 16 '16 at 23:00

I highly suggest the issue is one of complete separation. In high-dimensional analyses it is harder and harder to visualize such a phenomenon in hyperspace. Basically, a 2 dimensional analogue is fitting a logit curve to the sequence 0, 0, 0, 1, 1, 1. The OR is the slope of the S-shaped curve fit to those data, of course the best fitting line is actually a Heaviside step function with an OR of infinity. R's internal fitter does not actually converge to infinity, so it stops when the information matrix is numerically singular.

Just creating a frequency table will allow you to see which marginal cells are 0-count, and you will identify which if any of the factors need to be dropped.

• Hi, how can I create a table if one of my predictor is continuous? Secondly, if running separate analyses for control and test verb gives me normal estimates, does that mean that the complete separation is in the verb condition? – Pietre Dec 16 '16 at 21:22
• @Pietre Split into quantiles, not rigorous but a useful diagnostic. – AdamO Dec 16 '16 at 21:24
• I have a categorical version of the continuous variable based on results of regression trees (group a,b, c). there are a few 0s and 1s. is there a test to check for perfect separation? – Pietre Dec 16 '16 at 22:42

(This question is better-suited to cross-validated than stack overflow, because it is a stats question and not a programming question.)

First, your right-hand-side variables are very highly correlated. This will lead the coefficient estimates to blow up.

https://en.wikipedia.org/wiki/Variance_inflation_factor

The wikipedia article discusses it in terms of OLS, but the principle holds in linear models generally.

Basically, you can't really interpret your coefficients in isolation from one another when they are so deeply correlated. Especially with such a small sample size.

Second, look at crosstabs before fitting a model. Do the distributions of the RHS variables share support across both binary conditions? I.e.: is there complete separation?