[R] glmmPQL and Convergence

[Douglas Bates]

On 8/19/05, rab45+ at pitt.edu <rab45+ at pitt.edu> wrote:
> I fit the following model using glmmPQL from MASS:
> 
> fit.glmmPQL <-
> glmmPQL(ifelse(class=="Disease",1,0)~age+x1+x2,random=~1|subject,family=binomial)
> summary(fit.glmmPQL)
> 
> The response is paired (pairing denoted by subject), although some
> subjects only have one response. Also, there is a perfect positive
> correlation between the paired responses. x1 and x2 can and do differ
> within each pair. Here is the output:
> 
> > summary(fit.glmmPQL)
> Linear mixed-effects model fit by maximum likelihood
>  Data: fernando
>        AIC      BIC    logLik
>   30.51277 49.25655 -9.256384
> 
> Random effects:
>  Formula: ~1 | subject
>         (Intercept)     Residual
> StdDev:    8.284993 4.113725e-09

Notice the value of the residual standard deviation.  It's far too
small (it should be approximately 1 for a binomial-response model fit
by IRLS).  You have perfect prediction in your model and surprisingly
that is a problem in these models.

> 
> Variance function:
>  Structure: fixed weights
>  Formula: ~invwt
> Fixed effects: ifelse(class == "Disease", 1, 0) ~ age + x1 + x2
>                   Value Std.Error  DF   t-value p-value
> (Intercept)   -35.01862 2.4414559 123     -14.3       0
> age             0.59026 0.0441817 123      13.4       0
> x1              1.39317 0.0000014  41 1000507.2       0
> x2              0.93695 0.0000010  41  915150.3       0
>  Correlation:
>               (Intr) age        x2
> age           -0.952
> x1             0.000  0.000
> x2             0.000  0.000 -0.057
> 
> Standardized Within-Group Residuals:
>           Min            Q1           Med            Q3           Max
> -2.939213e+00 -2.509951e-07 -1.169248e-07  2.999710e-06  3.825035e+00
> 
> Number of Observations: 168
> Number of Groups: 125
> 
> 
> The t-values are huge and the se's are correspondingly tiny. The model
> does a great job of discriminating between disease and no disease. But I
> have a feeling there is something wrong here. Is there something wrong
> with the type of model I'm trying to fit? If it weren't for the pairing I
> would just have used glm. Any insights would be appreciated.
> 
> 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
>