[R] lmer - Is this reasonable output?

Dimitris Rizopoulos dimitris.rizopoulos at med.kuleuven.be
Thu Jun 29 16:31:13 CEST 2006


Hello Rick,

some comments below:

o in lmer(), 'method = "ML"' and its counterpart "REML" refer to the 
case of linear mixed models; for GLMMs the methods currently available 
are PQL and Laplace. The control argument 'usePQL = FALSE' can be used 
when 'method = Laplace', and in fact specifies not to use PQL as a 
refinement to the initial values.

o SAS proc NLMIXED performs (adaptive) Gaussian Quadrature and is 
conceptually closer to 'method = Laplace' in lmer(); for PQL you have 
to use proc GLIMMIX.

o for GLMMs both lmer() and NLMIXED work with an approximation to the 
observed data likelihood, since this likelihood involves an integral 
with no closed-form solution (except from very specific cases).

o with only 4 clusters I think it'd be difficult to estimate variance 
components; If you want to correct for site, you could just put it as 
an ordinary covariate in your logistic regression.

I hope my comments help.

Best,
Dimitris

----
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/(0)16/336899
Fax: +32/(0)16/337015
Web: http://med.kuleuven.be/biostat/
     http://www.student.kuleuven.be/~m0390867/dimitris.htm


----- Original Message ----- 
From: "Rick Bilonick" <rab45+ at pitt.edu>
To: "R Help" <r-help at stat.math.ethz.ch>
Sent: Thursday, June 29, 2006 3:52 PM
Subject: [R] lmer - Is this reasonable output?


> I'm estimating two models for data with n = 179 with four clusters 
> (21,
> 70, 36, and 52) named siteid. I'm estimating a logistic regression 
> model
> with random intercept and another version with random intercept and
> random slope for one of the independent variables.
>
>
> fit.1 <- lmer(glaucoma~(1|siteid)+x1
> +x2,family=binomial,data=set1,method="ML",
>  control=list(usePQL=FALSE,msVerbose=TRUE))
>
> Generalized linear mixed model fit using PQL
> Formula: glaucoma ~ (1 | siteid) + x1 + x2
>          Data: set1
> Family: binomial(logit link)
>      AIC      BIC    logLik deviance
> 236.7448 249.4944 -114.3724 228.7448
> Random effects:
> Groups Name        Variance Std.Dev.
> siteid (Intercept) 0.05959  0.24411
> number of obs: 179, groups: siteid, 4
>
> Estimated scale (compare to 1)  0.464267
>
> Fixed effects:
>             Estimate Std. Error z value Pr(>|z|)
> (Intercept) -2.213779   0.688158 -3.2170 0.001296 **
> x1           0.609028   0.293250  2.0768 0.037818 *
> x2           0.025027   0.009683  2.5846 0.009749 **
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> Correlation of Fixed Effects:
>       (Intr)     x1
> x1     -0.871
> x2     -0.372 -0.024
>
>> ranef(fit.1)
> An object of class “ranef.lmer”
> [[1]]
>  (Intercept)
> 1 -0.05053772
> 2 -0.21592157
> 3  0.36643051
> 4 -0.04141520
>
>
> fit.2 <- lmer(glaucoma~(x1|siteid)+x1
> +x2,family=binomial,data=set1,method="ML",
>  control=list(usePQL=FALSE,msVerbose=TRUE))
>
> Generalized linear mixed model fit using PQL
> Formula: glaucoma ~ (x1 | siteid) + x1 + x2
>          Data: set1
> Family: binomial(logit link)
>      AIC      BIC    logLik deviance
> 239.3785 258.5029 -113.6893 227.3785
> Random effects:
> Groups Name        Variance Std.Dev. Corr
> siteid (Intercept) 0.059590 0.24411
>        x1          0.013079 0.11436  0.000
> number of obs: 179, groups: siteid, 4
>
> Estimated scale (compare to 1)  0.4599856
>
> Fixed effects:
>              Estimate Std. Error z value Pr(>|z|)
> (Intercept) -2.2137787  0.6911360 -3.2031 0.001360 **
> x1           0.6090279  0.2995553  2.0331 0.042042 *
> x2           0.0250268  0.0097569  2.5650 0.010316 *
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> Correlation of Fixed Effects:
>       (Intr)     x1
> x1     -0.854
> x2     -0.372 -0.023
>
>> ranef(fit.2)
> An object of class “ranef.lmer”
> [[1]]
>  (Intercept)           x1
> 1 -0.05534800  0.009265667
> 2 -0.16991678 -0.052723584
> 3  0.27692026  0.137597965
> 4  0.01648737 -0.062232012
>
> Note that the fixed coefficient estimates are identical for both 
> models
> and they are exactly equal to what glm gives ignoring the sites (but 
> the
> standard errors given by lmer are definitely larger - which would 
> seem
> reasonable). The se's for the fixed factors differ slightly between 
> the
> two models. Note also that the estimated random effect sd for siteid
> intercept is identical for both models.
>
> I ran both models using PROC NLMIXED in SAS. It gives be similar
> estimates but not identical for the fixed effects and random 
> effects.
> The confidence intervals on the random effects for each site are 
> very
> large.
>
> Am I getting these results because I don't really need to fit a 
> random
> effect for siteid? The estimated random effects for slope seem to 
> say
> that siteid matters. When I plot the data for each site with 
> smoothing
> splines it indicates reasonably different patterns between sites.
>
> I don't think these models are that complicated and I have a 
> reasonable
> amount of data. I fit the fixed and random curves to the separate 
> sites
> (along with separate glm estimates for each site). As would be 
> expected,
> the random curves lie between the fixed effect curve and the glm 
> curve.
> But there seems to be a lot of shrinkage toward the fixed effect 
> curve.
> (The fixed effect curve fits the smoothing spline curve for all 
> sites
> combined very closely.)
>
> If I specify method="ML" and use PQL, I get similar fixed effect
> estimates (but not identical to glm). The random intercept sd about
> doubles. I think SAS NLMIXED uses an approximation to the likelihood 
> so
> that may explain some of the differences.
>
> One other thing that seems odd to me is the random intercept sd for 
> site
> id. It equals 0.24411 in both models. If I change x1 to, say x3 (an
> entirely different variable), I still get 0.24411. However, the 
> random
> slope sd does change.
>
> I just want to make sure I'm fitting a reasonable model and getting
> reasonable estimates.
>
> Rick B.
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html 


Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm



More information about the R-help mailing list