[R] How to use lmer function and multicomp package?

Giovanni Bacaro bacaro at unisi.it
Sun Jun 4 11:39:12 CEST 2006

Dear list members,
First of all thank you for your helpful advices.
After your answeres to my firt mail I studied a lot (R-News n°5) and I
tried to perform my analysis:

First, to fit a GLM with a nested design I decided to use the function
"lmer" in package "lme4"
as suggested by Spencer Graves and Filippo Piro.
I remember to you that my data were:

land use classes, 3 levels (fixed factor) = cla (R variable)
plot number, 98 levels each with 4 replicates (random factor within "cla")
= plotti (R variable)

number of species, totally 392 counts, response variable = sp (R variable)

Now, I started my analysis as follow (after the creation of the data.frame

mod1<-lmer(sp~cla+(1|cla:plotti), data=bacaro, family=poisson(link=log))

> summary(mod1) #sunto del modello

Generalized linear mixed model fit using PQL
Formula: sp ~ cla + (1 | cla:plotti)
          Data: bacaro
 Family: poisson(log link)

      AIC      BIC    logLik deviance
 451.2908 467.1759 -221.6454 443.2908

Random effects:
 Groups     Name        Variance Std.Dev.
 cla:plotti (Intercept) 0.60496  0.77779
number of obs: 392, groups: cla:plotti, 98

Estimated scale (compare to 1)  0.6709309

Fixed effects:
            Estimate Std. Error z value  Pr(>|z|)
(Intercept)  2.06406    0.14606 14.1316 < 2.2e-16 ***
cla2        -0.59173    0.17695 -3.3440 0.0008257 ***
cla3        -0.74230    0.83244 -0.8917 0.3725454
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
     (Intr) cla2
cla2 -0.825
cla3 -0.175  0.145

> anova(mod1,test="Chsqr")
Analysis of Variance Table
    Df Sum Sq Mean Sq
cla  2 11.352   5.676

Now, my questions are:
1) is the mod1 well specified? Have I said to R that "plotti" is my random
factor and that "plotti" is nested inside "cla" ("cla" as grouping
factor)? (can be  "lmer(sp~cla+(plotti|cla:plotti), data=bacaro,
family=poisson(link=log)" an alternative solution?)
2) Why if I try "lmer(sp~cla+(1|cla:plotti),..." or
"lmer(sp~cla+(1|plotti:cla),...." I obtain the same results?
3) why the anova summary don't say if differences in classes are
significance (or not significance)?
4) I'd like to perform a post-hoc test with the package "multicomp" but
the lmer function give me a lmer object (and this kind of object is not
read by the "multicomp" package). How could I perform my analysis in a
different way?

Thank you a lot for your help!

>I'd like to perform a glm analysis with a hierarchically nested design.
In particular,
>I have one fixed factor ("Land Use Classes") with three levels and a
random factor ("quadrat") nested within Land Use Classes with different
levels per classes (class artificial = 1 quadrat; class crops = 67
quadrats; and class seminatural = 30 quadrats).
>I have four replicates per each quadrats (response variable = species
richness per plot)

>Here some question about:
>1) could I analize these data using the class "artificial" (i.e. I have
only 1 level)?

>2) using R I'd like perfor a glm analysis considering my response
variable (count of species) with a Poisson distribution. How can I
develop my model considering the nested nature of my design? I'm sorry
but I don't know the right package to use.

>3) for the Anova analysis I'd like use a post-hoc comparison between
pairwise classes. What is the right procedure to do this? Is this
analysis performed in R?

Dr. Giovanni Bacaro
Università degli Studi di Siena
Dipartimento di Scienze Ambientali "G. Sarfatti"
Via P.A. Mattioli 4 53100 Siena
tel. 0577 235408
email: bacaro a unisi.it

More information about the R-help mailing list