[R-sig-ME] Restrictions on the class of GLMMs

[John Maindonald]

It occurred to me that I should try the calculation that I was after with
the development version (lme4a).  This version does allow the calculation 
(with one level of the grouping factor transect for the random effects
the same as the number of observations) to proceed. So I am hoping that 
lme4a is the way that it will be.  [Below, I compare the result with that from
glmmPQL from MASS.]

For this particular example, there are particular issues with the "Bank" habitat.
The count is zero, and the relevant parameter (on a log scale) has to 
represent 0 as exp(-C), where C is a suitably large number whose estimate 
depends on the convergence criteria.  So it is probably sensible to leave "Bank" 
out of the main calculation.  One has:

[lme4a_0.999375-44]
[Dependencies create problems for loading DAAG into R-devel.
So I saved the image under R-2.10.0 into moths.RData, and then
loaded that image.]
> library(lme4a)
> moths$transect <- 1:41  # Each row is from a different transect 
> moths$habitat <- relevel(moths$habitat, ref="Lowerside") 
> (A.glmer <-  glmer(A~habitat+log(meters)+(1|transect),  
+                    family=poisson, data=subset(moths,subset=habitat!="Bank"))) 
Generalized linear mixed model fit by the Laplace approximation 
Formula: A ~ habitat + log(meters) + (1 | transect) 
  Data: subset(moths, subset = habitat != "Bank") 
AIC BIC logLik deviance
208 223    -95      190
Random effects:
Groups   Name        Variance Std.Dev.
transect (Intercept) 0.229    0.478   
Number of obs: 40, groups: transect, 40

Fixed effects:
                Estimate Std. Error z value Pr(>|z|)
(Intercept)        1.0876     0.3963    2.74  0.00607
habitatDisturbed  -1.2326     0.4699   -2.62  0.00872
habitatNEsoak     -0.8210     0.4402   -1.87  0.06216
habitatNWsoak      1.5166     0.3915    3.87  0.00011
habitatSEsoak      0.0515     0.3505    0.15  0.88321
habitatSWsoak      0.2435     0.4543    0.54  0.59188
habitatUpperside  -0.1669     0.5366   -0.31  0.75570
log(meters)        0.1506     0.1374    1.10  0.27293
. . . .

The same (?) analysis is possible under glmmPQL from MASS.  
(This relies on iterated calls to lme(), from nlme.) The estimates are 
very similar, but the SEs are noticeably different:

> A.glmmPQL <-  glmmPQL(A~habitat+log(meters), random=~1|transect,  
+                    family=poisson, data=subset(moths,subset=habitat!="Bank"))
Random effects:
Formula: ~1 | transect
       (Intercept) Residual
StdDev:       0.448     1.04

Variance function:
Structure: fixed weights
Formula: ~invwt 
Fixed effects: A ~ habitat + log(meters) 
                 Value Std.Error DF t-value p-value
(Intercept)       1.110     0.437 32    2.54  0.0162
habitatDisturbed -1.240     0.528 32   -2.35  0.0253
habitatNEsoak    -0.821     0.488 32   -1.68  0.1027
habitatNWsoak     1.514     0.422 32    3.58  0.0011
habitatSEsoak     0.053     0.385 32    0.14  0.8917
habitatSWsoak     0.240     0.497 32    0.48  0.6326
habitatUpperside -0.166     0.590 32   -0.28  0.7795
log(meters)       0.147     0.151 32    0.97  0.3399

John Maindonald             email: john.maindonald at anu.edu.au
phone : +61 2 (6125)3473    fax  : +61 2(6125)5549
Centre for Mathematics & Its Applications, Room 1194,
John Dedman Mathematical Sciences Building (Building 27)
Australian National University, Canberra ACT 0200.
http://www.maths.anu.edu.au/~johnm

On 07/12/2009, at 8:10 PM, John Maindonald wrote:

> Reasonably recently, it used to be possible, with glmer, to
> associate one random term with each observation.  With the
> current version of glmer(), I find that this is not allowed.  In the
> case I was using it for, this allowed me to fit a random between
> observations variance that was additive on the scale of the
> linear predictor, rather than as with the dispersion estimate
> fudge which estimates a constant multiplier for the theoretical
> variance.  I am wondering what the reason may be for disallowing
> this; does it unduly complicate code somewhere or other?
> 
> I am using lme4_0.999375-32;    Matrix_0.999375-32
> 
> Here is what I had been doing:
> 
> library(DAAG)
> moths$transect <- 1:41  # Each row is from a different transect
> moths$habitat <- relevel(moths$habitat, ref="Lowerside")
> (A.glmer <-  glmer(A~habitat+log(meters)+(1|transect), 
>                               family=poisson, data=moths))
> 
> 
> Generalized linear mixed model fit by the Laplace approximation 
> Formula: A ~ habitat + log(meters) + (1 | transect) 
> Data: moths 
> AIC BIC logLik deviance
> 95 112  -37.5       75
> Random effects:
> Groups   Name        Variance Std.Dev.
> transect (Intercept) 0.234    0.483   
> Number of obs: 41, groups: transect, 41
> 
> ...
> 
> Thanks
> 
> John Maindonald             email: john.maindonald at anu.edu.au
> phone : +61 2 (6125)3473    fax  : +61 2(6125)5549
> Centre for Mathematics & Its Applications, Room 1194,
> John Dedman Mathematical Sciences Building (Building 27)
> Australian National University, Canberra ACT 0200.
> http://www.maths.anu.edu.au/~johnm
> 
> 
> John Maindonald             email: john.maindonald at anu.edu.au
> phone : +61 2 (6125)3473    fax  : +61 2(6125)5549
> Centre for Mathematics & Its Applications, Room 1194,
> John Dedman Mathematical Sciences Building (Building 27)
> Australian National University, Canberra ACT 0200.
> http://www.maths.anu.edu.au/~johnm
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models