[R-sig-ME] odds ratios from mixed logistic models
Beate Glaser
b.glaser at bristol.ac.uk
Fri Sep 26 21:59:06 CEST 2008
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
More information about the R-sig-mixed-models
mailing list