R help - Jun 2006 - How to use lmer function and multicomp package?

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
"bacaro"):


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!
Giovanni


>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

Comments below:
> 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)? 
Yes, you have plotti as a random factor nested in cla

(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?
Because this call is communative, so either way does not matter
> 3) why the anova summary don't say if differences in classes 
> are significance (or not significance)?
See a recent post by Doug Bates at
http://finzi.psych.upenn.edu/R/Rhelp02a/archive/76742.html
> 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!
> Giovanni
> 
> 
> 
> >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 at unisi.it
> 
> ______________________________________________
> 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
>