[R] lmer - Is this reasonable output?

[Rick Bilonick]

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.