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