[R-sig-ME] odds ratios from mixed logistic models

[Beate Glaser]

Dear lmer group,

I have been fitting a very simple mixed logistic regression with lmer
The outcome is educational impairment (edu.cat: 0,1) in students and the 
predictors are age and a base-line mental health measure (mh.s). I used a 
random intercept model, as none of my predictors, apart from age, is 
time-varying. From cross-sectional analysis I know that there is a 
significant change in the effect of base-line mental health on educational 
impairment across time, and this effect is included as a fixed-effect 
interaction term.
I wrote a little function to extract the OR for the mental health measure 
(increase in 1 point) given different ages of the student. I derived the SE 
of this OR as pooled SE using the SEs for the fixed effects "age" and "age 
: mental health" interaction, as well as their covariance. I am not sure 
however if this is too simple to obtain this SE, or if I need to do 
simulations (e.g. in WinBugs, unless someone knows a sneaky way to 
incorporate these combined ORs variables into the mcmcsampler within lmer) 
to obtain more accurate estimates. I would be very grateful for any 
comments and suggestions.

Kind regards,

Beate

This is the model with age centered at the mean (22.4 years):

>lmer.mh1 <- lmer(smfq.cat ~ age.c*( iq.8.s) + (0 + female|aln) + (0 + 
male|aln), data=psy.st.l, family=binomial(link=logit))

Generalized linear mixed model fit by the Laplace approximation
Formula: edu.cat ~ age.c * (mh.s) + (0 + female | aln) + (0 + male | aln)
   Data: psy.st.l
  AIC  BIC logLik deviance
 8399 8442  -4193     8387
Random effects:
 Groups Name   Variance Std.Dev.
 aln    female 3.49     1.87
 aln    male   2.18     1.48
Number of obs: 10149, groups: aln, 3438

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)   -0.9671     0.3113   -3.11   0.0019 **
age.c         -0.7776     0.1692   -4.60  4.3e-06 ***
mh.s          -0.1435     0.0292   -4.91  8.9e-07 ***
age.c:mh.s     0.0955     0.0161    5.93  3.1e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) age.c  mh.s
age.c       -0.132
mh.s        -0.989  0.145
age.c:mh.s   0.143 -0.989 -0.159

#OR function to obtain ORs at different ages

orme<-function(agev,mod) {
  fe<-fixef(mod)

#Derive combined SE
  se.fe<-se.fixef(mod)
  se.mh<-se.fe[names(se.fe)=="mh.s"]
  se.ia<-se.fe[names(sefe)=="age.c:mh.s"]
  vcmat<-vcov(lmer.mh1)
  cov.mh.ia<-vcmat[4,3]
  se.or.ia<-sqrt(se.mh^2 + se.ia^2 + 2*cov.mh.ia)

#Derive OR + 95%CI
  mh.dta1<-0
  mh.dta2<-1
  or<-rep(NA,length(agev))
  ci<-matrix(nrow=2,ncol=length(agev))

for(i in 1:length(agev)) {
    dta1<-c(1,agev[i],mh.dta1,agev[i]*mh.dta1)
    dta2<-c(1,agev[i],mh.dta2,agev[i]*mh.dta2)
    logodds1<-dta1%*%fe
    logodds2<-dta2%*%fe
    or[i]<-exp(logodds2-logodds1)
    ci[1:2,i]<- or[i] + c(-1.96,1.96)*se.or.ia
    }
  return(list(or=or,ci=ci))
  }

#Mean age
print(age.m, digits=6)
22.3942

#OR for age20, age23, age24
  orme(c(-2.4,0.6,1.6),lmer.mh1)

$or
[1] 0.69 0.92 1.01

$ci
     [,1] [,2] [,3]
[1,] 0.63 0.86 0.95
[2,] 0.75 0.98 1.07



R version 2.6.2 (2008-02-08)
i386-pc-mingw32

locale:
LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United 
Kingdom.1252;LC_MONETARY=English_United 
Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] arm_1.1-8         R2WinBUGS_2.1-8   coda_0.13-2       lme4_0.999375-18 
Matrix_0.999375-9 xtable_1.5-2      car_1.2-8         mice_1.16
 [9] nnet_7.2-42       MASS_7.2-42       gdata_2.4.2       lattice_0.17-10 
conf.design_0.0-4 foreign_0.8-26

loaded via a namespace (and not attached):
[1] grid_2.6.2   gtools_2.4.0 tools_2.6.2


----------------------
Beate Glaser
ALSPAC
MRC CAiTE
Dept Social Medicine
Oakfield House
15-23 Oakfield Grove
Bristol BS8 2BN
UK

++44 117 33 10101