3
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ 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" $\endgroup$ – kjetil b halvorsen Dec 16 '16 at 18:41
  • $\begingroup$ 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. $\endgroup$ – Pietre Dec 16 '16 at 22:46
  • 1
    $\begingroup$ 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. $\endgroup$ – kjetil b halvorsen Dec 16 '16 at 22:49
  • $\begingroup$ 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? $\endgroup$ – Pietre Dec 16 '16 at 22:56
  • $\begingroup$ 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. $\endgroup$ – kjetil b halvorsen Dec 16 '16 at 23:00
3
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Pietre Dec 16 '16 at 21:22
  • $\begingroup$ @Pietre Split into quantiles, not rigorous but a useful diagnostic. $\endgroup$ – AdamO Dec 16 '16 at 21:24
  • $\begingroup$ 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? $\endgroup$ – Pietre Dec 16 '16 at 22:42
0
$\begingroup$

(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?

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.