Dear Mixed-models group,

I am working with a dataset containing fixed and nested random effects. I
have
one fixed effect that I am most interested in and the others were collected
to
see if they also might have an influence. I apologize for the novel but
hopefully discussion of this will help others in the future who are as
intimidated as I was/am.
 
The data consist of total counts of anuran individuals from a particular
species
of epiphytic phytotelm plant found in tree canopies at two sites. The site
difference (if any) is my main interest. At each site all trees with
suitable
#ąs of this epiphyte species for sampling were located within a
predetermined
size area. 16 trees were randomly selected from those available at each site
and
5 epiphytes were then randomly sampled for all anurans within them. So, 2
sites
-> 16 trees (at each site) -> 5 epiphytes (in each tree), which equals 80
samples from each site for a total of 160. I also collected 5 other
variables
like size of epiphyte, height in tree, etc. I had originally planned to do a
nested anova analysis and a regression for the other variables until I had
the
data in hand ­ there were a lot of zeroąs. I began by fitting the count data
to
expected values for the Poisson and Negative binomial distributions
according to
Crawley 2007. The neg. binomial was the better fit, in the poisson
comparison my
data had a slightly greater # of zeroąs and a few values of 5 and 7 counts
per
epiphyte. With a mean = 0.59375 anurans/epiphyte and variance = 1.22386 It
appears I also have overdispersion.

I believe I need to run a GLMM analysis with the negative binomial for
nested
random effect (epiphyte nested in tree). I discovered the R package glmmADMB
which looks like it might work but there is very little documentation
available
and I have now tried to install it on my mac in OSX10.6.4 and Windows 7
32-bit
and it will not run the sample model on their site. It will start but not
run
anything. I tried one work around I found but I could not get that to fix it
(I
may have not done it right=complicated). I have read that it can only handle
one
random effect so I am not sure if it will work for my nested random effect.

I decided to proceed with lme4 and use the poisson and quasipoisson so I
could
at least get myself a little familiar with mixed models before posting here
as
what to do next. I began with poisson since it gives p-values and I was able
to
find some good examples on how to work my way through the model selection
process. I began with the following syntax:

m1po<-lmer(count~treat+treedbh+treehgt+numepi+elevepi+hgtepi+leafepi+
(1|tree/epi),family=poisson,data=ecpad2)
 
I included the summary of this model and the final model at the end of this
email. I worked my way through the model reduction and in the end only the
treat
(site) factor was significant.
 
I moved on to the quasipoisson in lme4. Having initially read how this has
an
additional parameter for adjusting for overdispersion I thought this was it!
Then I read more and saw that it appears to have issues in lme4 but is still
better than glmmPQL. Correct? I started with the full model but I am
uncertain
how to begin the reduction process. I thought I read somewhere to start with
the
t-values closest to zero and remove those factors one by one by conducting
an
anova model comparison for each subsequent model (syntax and results at end
of
email). That was the procedure I went with and once again only treat was
significant. Is this correct? I am thinking/hoping that maybe this quasi-
GLMM
will suffice as it addresses the overdispersion issue and my data is not too
far
off from the poisson, but I just keep coming across conflicting mailing list
posts that talk about problems with it but some say to still go with it over
PQL. If quasi- GLMM works can I report the p-value given for the model
comparison? Or should I calculate degrees of freedom and use the t-value
from
the simplest model to generate a p-value? I know some donąt like the
infamous
p-value but I would prefer to avoid the inquisition from manuscript
reviewers.
 
If it is recommended that I should use glmmADMB I will try to get it running
but
certainly will need help with syntax and interpretation.
 
I sincerely appreciate any and all assistance.
 
As one poster said, łI am a mere (student) ecologist, trying to get the
truth
out of his data, and not a statistician, so forgive me my lack of
statistical
background and possible conceptual misunderstandings˛ ­ my sentiments
exactly. 
 
Best,
 
Shawn McCracken
 
PhD Student
Texas State University ­ San Marcos
 
 
 
 
 
 
Full Poisson model
 
> summary(m1po)
Generalized linear mixed model fit by the Laplace approximation
Formula: count ~ treat + treedbh + treehgt + numepi + elevepi + hgtepi +
leafepi + (1 | tree/epi)
   Data: ecpad2 
   AIC   BIC logLik deviance
 231.6 262.3 -105.8    211.6
Random effects:
 Groups    Name        Variance Std.Dev.
 epi:tree (Intercept) 0.472608 0.68747
 tree      (Intercept) 0.053143 0.23053
Number of obs: 160, groups: epi:tree, 80; tree, 16
 
Fixed effects:
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.338606   1.556223  -2.145  0.03193 *
treat      -0.296865  0.271361  -1.094  0.27396
treedbh     -0.008924   0.003160  -2.824  0.00474 **
treehgt      0.076279   0.036195   2.107  0.03508 *
numepi      0.001981   0.003742   0.529  0.59651
elevepi    -0.014964   0.032105  -0.466  0.64116
hgtepi      0.002821   0.010349   0.273  0.78517
leafepi     0.021787   0.016733   1.302  0.19289
---
Signif. codes:  0 Ś***ą 0.001 Ś**ą 0.01 Ś*ą 0.05 Ś.ą 0.1 Ś ą 1
 
Correlation of Fixed Effects:
        (Intr) treat  tredbh trehgt numbrm elvbrm hgtbrm
treat    -0.409    
treedbh   0.211  0.001
treehgt  -0.619  0.309 -0.409
numepi   0.116 -0.080 -0.261  0.005
elevepi -0.088 -0.294  0.060 -0.560 -0.210
hgtepi  -0.612  0.408 -0.021 0.138 -0.175  0.076
leafepi -0.197  0.161 -0.029  0.116  0.040 -0.175 -0.239
 
Final Poisson model
 
> m9po<-lmer(count~treat+(1|tree/epi),family=poisson,data=ecpad2)
> summary(m9po)
Generalized linear mixed model fit by the Laplace approximation
Formula: count ~ treat + (1 | tree/epi)
   Data: ecpad2 
   AIC   BIC logLik deviance
 231.4 243.7 -111.7    223.4
Random effects:
 Groups    Name        Variance Std.Dev.
 epi:tree (Intercept) 0.44136  0.66435
 tree      (Intercept) 0.11113  0.33336
Number of obs: 160, groups: epi:tree, 80; tree, 16
 
Fixed effects:
           Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.5990     0.1797  -3.333 0.000858 ***
treat       -0.4496    0.2181  -2.062 0.039249 *
---
Signif. codes:  0 Ś***ą 0.001 Ś**ą 0.01 Ś*ą 0.05 Ś.ą 0.1 Ś ą 1
 
Correlation of Fixed Effects:
      (Intr)
treat -0.473
 
Full Quasipoisson model
 
>
m1qp<-lmer(count~treat+treehgt+numepi+elevepi+hgtepi+leafepi+(1|tree/epi),fa
mily=quasipoisson,data=ecpad2)
> summary(m1qp)
Generalized linear mixed model fit by the Laplace approximation
Formula: count ~ treat + treehgt + numepi + elevepi + hgtepi + leafepi +
(1
| tree/epi) 
   Data: ecpad2 
   AIC   BIC logLik deviance
 240.9 271.7 -110.4    220.9
Random effects:
 Groups    Name        Variance  Std.Dev.
 epi:tree (Intercept) 0.2233765 0.472627
 tree      (Intercept) 0.0093212 0.096547
 Residual             0.4827734 0.694819
Number of obs: 160, groups: epi:tree, 80; tree, 16
 
Fixed effects:
            Estimate Std. Error t value
(Intercept) -2.341386   0.996047  -2.351
treat      -0.313394  0.182984  -1.713
treehgt      0.035933   0.021961   1.636
numepi     -0.002087   0.002442  -0.854
elevepi    -0.011191   0.021150  -0.529
hgtepi      0.002057   0.006895   0.298
leafepi     0.020682   0.011586   1.785
 
Correlation of Fixed Effects:
        (Intr) treat  trehgt numbrm elvbrm hgtbrm
treat    -0.368    
treehgt  -0.571  0.267
numepi   0.211 -0.123 -0.160
elevepi -0.107 -0.273 -0.593 -0.185
hgtepi  -0.619  0.377  0.102 -0.166  0.091
leafepi -0.228  0.193  0.101  0.004 -0.167 -0.194
 
Final Quasipoisson model
 
> m6qp<-lmer(count~treat+(1|tree/epi),family=quasipoisson,data=ecpad2)
> summary(m6qp)
Generalized linear mixed model fit by the Laplace approximation
Formula: count ~ treat + (1 | tree/epi)
   Data: ecpad2 
   AIC   BIC logLik deviance
 233.4 248.8 -111.7    223.4
Random effects:
 Groups    Name        Variance Std.Dev.
 epi:tree (Intercept) 0.207451 0.45547
 tree      (Intercept) 0.052232 0.22854
 Residual             0.470031 0.68559
Number of obs: 160, groups: epi:tree, 80; tree, 16
 
Fixed effects:
           Estimate Std. Error t value
(Intercept)  -0.5990     0.1232  -4.862
treat       -0.4496    0.1495  -3.007
 
Correlation of Fixed Effects:
      (Intr)
treat -0.473
> anova(m6qp,m5qp)
Data: ecpad2
Models:
m6qp: count ~ treat + (1 | tree/epi)
m5qp: count ~ treat + leafepi + (1 | tree/epi)
     Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)
m6qp  5 233.39 248.77 -111.70
m5qp  6 234.18 252.63 -111.09 1.2126      1     0.2708
> m7qp<-lmer(count~(1|tree/epi),family=quasipoisson,data=ecpad2)
> summary(m7qp)
Generalized linear mixed model fit by the Laplace approximation
Formula: count ~ (1 | tree/epi)
   Data: ecpad2 
   AIC   BIC logLik deviance
 236.1 248.4 -114.0    228.1
Random effects:
 Groups    Name        Variance Std.Dev.
 epi:tree (Intercept) 0.219259 0.46825
 tree      (Intercept) 0.055208 0.23496
 Residual             0.496793 0.70484
Number of obs: 160, groups: epi:tree, 80; tree, 16
 
Fixed effects:
            Estimate Std. Error t value
(Intercept)  -0.7988     0.1116  -7.156
> anova(m7qp,m6qp)
Data: ecpad2
Models:
m7qp: count ~ (1 | tree/epi)
m6qp: count ~ treat + (1 | tree/epi)
     Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)
m7qp  4 236.07 248.37 -114.04
m6qp  5 233.39 248.77 -111.70 4.6807      1    0.03050 *

	[[alternative HTML version deleted]]


