[R-sig-ME] Restrictions on the class of GLMMs
John Maindonald
john.maindonald at anu.edu.au
Tue Dec 8 23:19:31 CET 2009
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
More information about the R-sig-mixed-models
mailing list