[R-lang] mixed logit models

Wed Apr 8 15:47:07 PDT 2009

Dear Maria,

On 08-04-2009 23:06, r-lang-request at ling.ucsd.edu wrote:
> The DV of my expt  is a target response (poresp)  which can be either
> PO (coded as 1) or DO (coded as 0). So PO is a success and DO is a
> failure. The predictors are prime and nounrep (nounrepetition), each
> with 2 levels. I have effect coded the predictors (-.5 and +.5
> respectively for each of the two levels of the predictors, following
> Barr, 2008): with this coding, which I believe is a way of centering,
> the intercept should correspond to the grand mean and the coefficients
> to the differences between the means etc., like in an ANOVA. Please
> correct me if I am wrong. I tried this coding in the past where the
> outcome was a continuous variable and indeed when I checked the means
> and differences in my data , I found a correspondence between the ANOVA
> measures and the regression intercept and coefficients; however, now I
> am dealing  with binomially distributed data, where the model parameters
> are in log odds space, so it is a bit more complicated, because I have
> to do transformations). 
> 
> The problem is that when I check the output of the mixed  logit
> regression (using lmer and family = binomial ) and try to make sense of
> the coefficients by checking, for example, whether the intercept
> corresponds to the grand mean, in log odds, of course, things do not
> match. I checked and re-checked, and I do not know where I am going
> wrong. 
> 
> ==========
> 
>> head(verbdiff)
>      X subject cond item myscore poresp nounrep prime primec nounrepc
> 1 1069       1    2   19       p      1       0     1    0.5     -0.5
> 2 1070       1    3    6       p      1       1     0   -0.5      0.5
> 3 1071       1    3    5       p      1       1     0   -0.5      0.5
> 4 1072       1    2   12       p      1       0     1    0.5     -0.5
> 5 1073       1    4   16       p      1       1     1    0.5      0.5
> 6 1074       1    3   14       p      1       1     0   -0.5      0.5
> 
> #primec and nounrepc are effect coded
> 
> 
> 
> COUNTS OF 0s and 1s:
> 
>> xtabs (~ poresp +prime +nounrep, data=verbdiff)
> 
>      nounrep = 0
> 
>       prime
> poresp   0   1
>      0  89  58
>      1 207 235
> 
> , , nounrep = 1
> 
>       prime
> poresp   0   1
>      0  91  64
>      1 202 228
> 
> 
> MODEL:
> 
>> verbdiff.lmer = lmer(poresp ~ primec *nounrepc + (1|subject) +
> (1|item), data= verbdiff, family="binomial")
> 
>> print (verbdiff.lmer)
> 
> 
> Generalized linear mixed model fit using Laplace 
> Formula: poresp ~ primec * nounrepc + (1 | subject) + (1 | item) 
>    Data: verbdiff 
>  Family: binomial(logit link)
>   AIC  BIC logLik deviance
>  1092 1122   -540     1080
> Random effects:
>  Groups  Name        Variance Std.Dev.
>  item    (Intercept) 1.2143   1.1020  
>  subject (Intercept) 1.8470   1.3590  
> number of obs: 1174, groups: item, 40; subject, 32
> 
> Estimated scale (compare to  1 )  0.917418 
> 
> Fixed effects:
>                 Estimate Std. Error z value Pr(>|z|)    
> (Intercept)       1.6605     0.3121   5.320 1.04e-07 ***
> primec            0.7854     0.1645   4.773 1.81e-06 ***
> nounrepc         -0.1054     0.1612  -0.653    0.513    
> primec:nounrepc  -0.2138     0.3224  -0.663    0.507    
> ---
> Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.?
> 0.1 ? ? 1 
> 
> Correlation of Fixed Effects:
>             (Intr) primec nonrpc
> primec       0.053              
> nounrepc    -0.009 -0.023       
> primc:nnrpc -0.012 -0.028  0.128
> 
> =============
> 
> QUESTIONS:
> 
> 1. THE ESTIMATED INTERCEPT IN THE LOGIT MODEL IS 1.66 IN LOG ODDS,
> WHICH should correspond to grand mean, i.e. likelihood of successes (=PO
> responses coded as "1") independent of the predictors.
> 
> THERE WERE OVERALL 872 SUCCESSES AND 302 FAILURES IN THE EXPT, SO ODDS
> SHOULD BE 872/302=2.88 or (in probability space) .74/.26 = 2.85;
> THIS SHOULD GIVE A LOG OF ODDS OF APPROX 1.05, BUT THE INTERCEPT
> PREDICTED BY THE MODEL  IS MUCH HIGHER (1.66)

Indeed:
 > logit(872/(872+302))
[1] 1.060362

This may be because with your coding, you are estimating the 
intercept if primec==0 & nounrepc==0. This means that if prime and 
DV are related (e.g. because there are more hits if prime==1), or 
likewise if noun and DV are related, then the estimated intercept is 
no longer the grand mean, but rather the estimated value of your DV 
if all predictors would be zero.

> 
> SAME THING WHEN I USE THE CORRESPONDING DUMMY CODED (0-1) V
> ARIABLES
> -INTERCEPT IS MUCH HIGHER THAN WHAT MY DATA SUGGEST 
> 

In my experience, models are far easier to interpret if you do NOT 
centralize the binary dummy predictors. The estimated intercept then 
equals the mean logit for the particular cell where all dummies are 
zero. So, my advice is to use your dummy coded variables and not the 
centered dummies.

The main point I guess is that the intercept does not reflect the 
grand mean, in your model, but it reflects the estimated DV if all 
predictors would be zero, which is different from the grand mean IF 
predictors and DV are related.

Did you try this model...
 >> verbdiff.m0.lmer = lmer(poresp ~ 1+(1|subject) +
 > (1|item), data= verbdiff, family="binomial")

It should give you ONLY the intercept which should then be about 1.06.

More background and details are available at
   http://www.hugoquene.nl/mixedeffects/
and references in the articles mentioned there.

Hope this helps, with kind regards, Hugo Quene

-- 
Dr Hugo Quené | assoc prof Phonetics | Utrecht inst of Linguistics 
OTS | Utrecht University | Trans 10 | 3512 JK Utrecht | The 
Netherlands | T +31 30 253 6070 | F +31 30 253 6000 | H.Quene at uu.nl 
| www.hugoquene.nl | www.hum.uu.nl