Dear List Members,
I am working with generalised linear mixed models (GLMMs), mainly using the
lme4 package in R (2.9.2 and 2.12.1).
The study involves looking at the effects of varying levels of grazing
intensity (grazing states) on the biodiversity (birds, beetles and plants)
of hill sheep farms. Here I will concentrate on the birds. The data is of a
nested nature so I chose to use linear mixed effects models. None of the
response variables are normally distributed so I could:
i) transform the response variables and run a linear
mixed effects model or
ii) use a generalised linear mixed model
Random effects are: *transects*, which are nested within *farms*, which are
nested within *sites*
Fixed effects are: *grazing state* and *altitude*
Response variable 1 is *bird density* (continuous, range: 1.61-36.05)
Error structure: GAMMA DISTRIBUTION
Response variable 2 is *bird diversity* (Simpson’s Index) (continuous,
range: 1-15.07)
Error structure: GAMMA DISTRIBUTION
Response variable 3 is *bird species richness* (counts, range: 1-19)
Error structure: POISSON OR NEGATIVE BINOMIAL DISTRIBUTION
Response variable 4 is *bird evenness* (proportions, range: 0.52-1)
Error structure: BINOMIAL DISTRIBUTION
The birds were surveyed over two years (and three farms from 2007 were
repeated in 2008):
*2007*
grazing.state = (1,2,5)
alt = altitude (continuous in metres)
site = area containing 3 farms (1,2,3,4)
farm = farm number (1,2,3,4,5,6,7,8,9,10,11,12)
t.id = transect identification number (1 – 10)
*2008*
grazing.state = (1,2,3,5)
alt = altitude (continuous in metres)
site = area containing 3 farms (1,2,3 (repeated from 2007), 5,6,7)
farm = farm number (3, 4, 8 (repeated from 2007), 13, 14, 15, 16, 17, 18,
19, 20, 21)
t.id = transect identification number (1 – 10)
Examples of model syntax:
*library(nlme)*
a) *linear regression model*
> with(avian,{model1 <- lm(*bird.density* ~ grazing.state*alt)})
b) *generalised linear model*
> with(avian,{model2<-glm(*bird.density* ~ grazing.state * alt, family =
Gamma, data = avian)})
c) *generalised linear model (poisson distribution, as bird species richness
is count data)*
> with(avian,{model3 <- glm(*bird.sp.rich* ~ grazing.state*alt, family =
poisson, data = avian)})
*library(MASS)*
d) *generalised linear model (negative binomial distribution)*
> with(avian,{model4 <- glm.nb(*bird.sp.rich* ~ grazing.state * alt ,
link="log",data= avian)})
e) *linear mixed effects model*
> with(avian,{model5<-lme(*bird.density* ~ grazing.state * alt, random = ~
1|site/farm/t.id, data = avian)})
*library(lme4)*
As density is continuous, I am using a Gamma distribution. However I can’t
seem to get the GLMM to work in lme4:
f) *generalised linear mixed model*
> with(avian,{model6a<-glmer(*bird.density* ~ grazing.state * alt +
(1|site/farm/t.id), family = Gamma, link = log, data = avian)})
Error in function (fr, FL, glmFit, start, nAGQ, verbose) :
Number of levels of a grouping factor for the random effects
must be less than n, the number of observations
Removing transect id gets rid of that error message, although I’m not too
confident about the random error structure now. I then get a new error
message:
> with(avian,{model6b<-glmer(*bird.density* ~ grazing.state * alt +
(1|site/farm), family = Gamma, link = log, data = avian)})
Error in mer_finalize(ans) :
mu[i] must be positive: mu = -0.0211949, i = 100
> avian$bird.density
[1] 17.14 23.17 9.89 13.70 4.82 6.50 3.81 3.81 7.84 9.46 7.24 9.89
[13] 4.02 10.84 14.37 12.46 6.23 7.24 4.62 3.81 4.62 4.02 8.65
19.12
[25] 14.51 7.39 1.61 3.81 9.46 6.00 9.65 12.67 3.81 6.21 4.02
6.21
[37] 3.44 3.81 3.81 6.00 21.90 11.34 10.04 15.69 23.96 12.87 11.63
9.27
[49] 12.07 9.44 11.63 19.95 10.68 17.72 22.56 7.84 8.03 7.02 5.69
7.82
[61] 35.26 22.56 5.42 5.26 4.10 2.19 2.19 2.19 5.05 3.81 20.96
19.70
[73] 9.66 6.87 15.67 3.81 4.84 4.03 2.15 9.02 23.16 16.31 16.35
27.56
[85] 18.88 8.01 3.81 5.63 5.42 6.00 36.05 20.36 13.25 11.50 16.33
6.46
[97] 4.02 2.40 7.04 3.81 34.43 33.96 35.55 20.94 25.20 7.82 6.46
5.42
[109] 5.42 5.55 27.03 26.74 25.35 29.16 20.55 6.00 8.78 12.87 5.20
8.94
There are no negative values here?
I get similar error messages when I replace bird.density with the response
variable bird.diversity (also continuous).
Bird species richness (count data) does work:
> with(avian,{model6c<-glmer(*bird.sp.rich* ~ grazing.state * alt +
(1|site/farm), family = poisson, link = log, data = avian)})
I understand how to use the quasipoisson command to look at the dispersion
parameter in a glm but I wasn’t sure how to do this in a glmm. I found this
syntax on an earlier thread from the r-sig-mixed-models list:
> (phi<-lme4:::sigma(model))
This worked but it gave me a figure of 1, which I was a little suspicious
of, however, plotting the residuals shows no patterns so perhaps this is
okay?
Bird evenness (proportions) doesn’t work either:
> with(avian,{model6d<-glmer(*bird.even* ~ grazing.state * alt +
(1|site/farm), family = binomial, link = logit, data = avian)})
Error in glmer(bird.even ~ grazing.state * alt + (1 | site/farm), family =
binomial, :
object 'logit' not found
> with(avian,{model6e<-glmer(*bird.even* ~ grazing.state * alt +
(1|site/farm), family = binomial, data = avian)})
Warning messages:
1: In eval(expr, envir, enclos) :
non-integer #successes in a binomial glm!
2: In mer_finalize(ans) : singular convergence (7)
I also tried link = quasibinomial, which works in R 2.9.2 but not in R
2.12.1.
Going back to bird density. It seems to works in *MASS* but a message
appears saying iteration 1? It still gives me a model summary though.
> with(avian,{model7<-glmmPQL(*bird.density* ~ grazing.state * alt, random =
~ 1|site/farm/t.id, family = Gamma, data = avian)})
iteration 1
Bird diversity, species richness and evenness give larger error messages.
I have also run the models using *glmmML* (although species richness is the
only one which works as it only allows poisson and binomial distributions).
> with(avian,{model8<-glmmML(*bird.sp.rich* ~ grazing.state * alt, cluster =
site, family = poisson, data = avian)})
I tried *glmmADMB* too (in R 2.9.2) but couldn’t get it to work at all.
So, my questions are:
Q1) Do you think I am having so many problems fitting glmm models to my data
because the response variables are not counts or proper proportions? I had
originally thought that Gamma would be able to deal with continuous data. Do
you think that given all my problems, it may be more sensible to transform
the response variables and run linear mixed effects models instead?
Q2) My random effects structure is nested (1|site/farm/t.id), as opposed to
crossed (1|site) + (1|farm) + (1|t.id). However warning messages appear in
most of the glmm unless I remove t.id, and in some cases farm too. I had
read about adding a random effect with a distinct level for each of my 120
observations, however this doesn’t seem to work either. Have I made a
fundamental error in my random effects structure?
Q3) I also have a fifth response variable which I would like to look at. It
is bird beta diversity (Sorensen’s Index). Data are in the form of
proportions (range: 0.18-0.53). However, the data show a clear
*bimodal*distribution. I understand that transformations are not
really an option but
I have read that *finite mixture models* may be a good alternative and was
just wondering if anyone may be able to suggest some R syntax I might use to
tackle this?
Many, many thanks in advance for any help, advice or guidance.
Best Wishes,
Roz
--
Roslyn Anderson
Postgraduate Researcher
Environmental Research Institute, Lee Road, University College Cork, Cork,
Ireland
and
Department of Zoology, Ecology, and Plant Sciences, University College Cork,
Distillery Fields, North Mall, Cork
E-mail: r.anderson@umail.ucc.ie Tel: +353 (0)21 490 1944 Fax: +353 (0)21 490
1932
[[alternative HTML version deleted]]