Dear mixed models list,


I have inherited a data set from a field experiment examining the effects
of fire exclusion on grassland vegetation.  The design has uneven
replication at one level and this has caused me some minor grief - it
seemed such a simple problem initially -  I would *really* appreciate some
feedback on the random effects structures I have tried.


There are 10 sites (which I am treating as random) and these occur on one
of two soils (fixed effect with 2 levels - (1) alluvial , 6 sites or
(2)aeolian derived sand, 4 sites).


 Within each site the grasslands were burnt on a biennial to triennial
basis with the exception of an area of approximately 1ha.  Within the burnt
area there were two quadrats (100m2 quadrats) and a single quadrat (again
100m2) was placed within the unburnt section (so three quadrats in each
site,  but 2 in burnt).


Within each quadrat there were 15 random plots (0.25m2) within which data
on vegetation was collected.  I only have data for the 10th year – so no
temporal component to the data.

I have been analysing the data using glmer and lmer in lme4.  The example I
provide is for the proportional cover of the dominant grass species *Themeda
triandra*.  *T. triandra*  cover was measured as a proportion (out of 25)
in each small plot.


My fixed effects model structure was simply  Them_tria/n~ Burnt+ soilgrp +
Burnt: soilgrp

Where n = 25.


I have tried a lot of alternative random effects structures including:

a)  (1 | sitef/Quadf)+( 1 | obs)  # observation level added to account for
overdispersion in this case

And

b)  (Burnt | sitef) + (1|obs)

Or

c)  (1 | sitef)+( 1 | obs)

In (a) I was hoping to account for the nested structure ie quadrats within
sites but because there is only one unburnt  quadrat in each site it seems
to me that I am possibly allocating variation that should be attributed to
fixed effects to that stratum of the random effects.  Indeed the random
effects model typically gives an AIC similar to that of the full model
(with both random and fixed effects), and while the fit of the data is
excellent, errors for fixed effects suggest this is a very conservative
approach, but is currently the one I prefer.

The second choice (b) allows the slopes associated with burning to vary
within each site –but model AIC is always  substantially higher than for
(a) and does not quite tackle the experimental structure.  I would probably
be happier with this if there were only a single burnt quadrat in each
site.

The last choice (c) seems just plain wrong and does not account for the
design structure – standard errors for fixed effects are relatively small
and generate confidence intervals that seem too tight given the data .  Also
this generates, not surprisingly, the largest AIC and the poorest fit to
the data.

I would be really appreciative of any suggestions for an alternative random
effects structure that might provide a better approach than (a) but still
capture the data structure.  Or even just comments on the approaches I have
suggested.

Below are model outputs for each of the alternatives (having fitted the
full fixed effects structure only), using the full dataset.  I have also
attached a subset of the data and code at the bottom of the message  (that
only includes  4 sites and 5 plots per quadrat).





fm.a <- glmer(Them_tria/n~ Burnt+ soilgrp + Burnt: soilgrp

+ (1 | sitef/Quadf)+ (1|obs), data, family= binomial(link="logit"), weight=
as.numeric(n))



> fm.a

Generalized linear mixed model fit by the Laplace approximation

Formula: Them_tria/n ~ Burnt + soilgrp + Burnt:soilgrp + (1 |
sitef/Quadf) +
(1 | obs)

   Data: data

  AIC  BIC logLik deviance

 1089 1118 -537.4     1075

Random effects:

 Groups      Name        Variance Std.Dev.

 obs         (Intercept)  8.2616  2.8743

 Quadf:sitef (Intercept) 17.3863  4.1697

 sitef       (Intercept)  0.0000  0.0000

Number of obs: 450, groups: obs, 450; Quadf:sitef, 30; sitef, 10



Fixed effects:

                  Estimate Std. Error z value Pr(>|z|)

(Intercept)          1.826      1.775   1.029   0.3034

BurntYes             4.996      2.179   2.293   0.0218 *

soilgrp2            -3.066      2.771  -1.107   0.2685

BurntYes:soilgrp2   -3.713      3.402  -1.091   0.2751

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



Correlation of Fixed Effects:

            (Intr) BrntYs slgrp2

BurntYes    -0.815

soilgrp2    -0.640  0.522

BrntYs:slg2  0.522 -0.640 -0.814





fm.b <- glmer(Them_tria/n~ Burnt+ soilgrp + Burnt: soilgrp + (Burnt | sitef)

 + (1|obs), data, family= binomial(link="logit"), weight= as.numeric(n))



> fm.b

Generalized linear mixed model fit by the Laplace approximation

Formula: Them_tria/n ~ Burnt + soilgrp + Burnt:soilgrp + (Burnt | sitef) +
(1 | obs)

   Data: data

  AIC  BIC logLik deviance

 1171 1204 -577.4     1155

Random effects:

 Groups Name        Variance Std.Dev. Corr

 obs    (Intercept) 16.891   4.1098

 sitef  (Intercept) 63.067   7.9414

        BurntYes    58.109   7.6229   -0.998

Number of obs: 450, groups: obs, 450; sitef, 10



Fixed effects:

                  Estimate Std. Error z value Pr(>|z|)

(Intercept)        5.46246    3.33825   1.636    0.102

BurntYes           1.94663    3.24980   0.599    0.549

soilgrp2          -7.09645    5.22212  -1.359    0.174

BurntYes:soilgrp2 -0.05792    5.06001  -0.011    0.991



Correlation of Fixed Effects:

            (Intr) BrntYs slgrp2

BurntYes    -0.987

soilgrp2    -0.639  0.631

BrntYs:slg2  0.634 -0.642 -0.991





fm.c <- glmer(Them_tria/n~ Burnt+ soilgrp + Burnt: soilgrp +

(1 | sitef) + (1|obs), data, family= binomial(link="logit"), weight=
as.numeric(n))





>fm.c

Generalized linear mixed model fit by the Laplace approximation

Formula: Them_tria/n ~ Burnt + soilgrp + Burnt:soilgrp + (1 | sitef) +      (1
| obs)

   Data: data

  AIC  BIC logLik deviance

 1265 1289 -626.4     1253

Random effects:

 Groups Name        Variance Std.Dev.

 obs    (Intercept) 20.7804  4.5586

 sitef  (Intercept)  4.4972  2.1207

Number of obs: 450, groups: obs, 450; sitef, 10



Fixed effects:

                  Estimate Std. Error z value Pr(>|z|)

(Intercept)         1.9293     1.0164   1.898   0.0577 .

BurntYes            6.5087     0.7763   8.385  < 2e-16 ***

soilgrp2           -3.4173     1.5988  -2.137   0.0326 *

BurntYes:soilgrp2  -4.8413     1.0917  -4.435 9.22e-06 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



Correlation of Fixed Effects:

            (Intr) BrntYs slgrp2

BurntYes    -0.355

soilgrp2    -0.636  0.226

BrntYs:slg2  0.252 -0.711 -0.389



#Input data (subset) and code


structure(list(Site = c(1L, 1L, 1L, 2L, 2L, 2L, 9L, 9L, 9L, 10L,

10L, 10L, 1L, 1L, 1L, 2L, 2L, 2L, 9L, 9L, 9L, 10L, 10L, 10L,

1L, 1L, 1L, 2L, 2L, 2L, 9L, 9L, 9L, 10L, 10L, 10L, 1L, 1L, 1L,

2L, 2L, 2L, 9L, 9L, 9L, 10L, 10L, 10L, 1L, 1L, 1L, 2L, 2L, 2L,

9L, 9L, 9L, 10L, 10L), Quadrat = c(1.1, 1.2, 1.3, 2.1, 2.2, 2.3,

9.1, 9.2, 9.3, 10.1, 10.2, 10.3, 1.1, 1.2, 1.3, 2.1, 2.2, 2.3,

9.1, 9.2, 9.3, 10.1, 10.2, 10.3, 1.1, 1.2, 1.3, 2.1, 2.2, 2.3,

9.1, 9.2, 9.3, 10.1, 10.2, 10.3, 1.1, 1.2, 1.3, 2.1, 2.2, 2.3,

9.1, 9.2, 9.3, 10.1, 10.2, 10.3, 1.1, 1.2, 1.3, 2.1, 2.2, 2.3,

9.1, 9.2, 9.3, 10.1, 10.2), Plot = c(1L, 1L, 1L, 1L, 1L, 1L,

1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,

2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L,

4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L,

5L, 5L, 5L, 5L, 5L), obs = c(1L, 16L, 31L, 46L, 61L, 76L, 361L,

376L, 391L, 406L, 421L, 436L, 2L, 17L, 32L, 47L, 62L, 77L, 362L,

377L, 392L, 407L, 422L, 437L, 3L, 18L, 33L, 48L, 63L, 78L, 363L,

378L, 393L, 408L, 423L, 438L, 4L, 19L, 34L, 49L, 64L, 79L, 364L,

379L, 394L, 409L, 424L, 439L, 5L, 20L, 35L, 50L, 65L, 80L, 365L,

380L, 395L, 410L, 425L), Burnt = structure(c(2L, 2L, 1L, 2L,

2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L,

1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L,

2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L,

2L, 1L, 2L, 2L, 1L, 2L, 2L), .Label = c("No", "Yes"), class = "factor"),

    soilgrp = c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,

    1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,

    1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,

    2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,

    2L, 2L), Them_tria = c(25L, 25L, 25L, 25L, 25L, 8L, 25L,

    17L, 0L, 0L, 24L, 5L, 25L, 24L, 1L, 25L, 25L, 0L, 25L, 3L,

    0L, 0L, 25L, 4L, 25L, 25L, 11L, 22L, 0L, 7L, 25L, 0L, 0L,

    0L, 25L, 3L, 25L, 25L, 2L, 11L, 25L, 25L, 25L, 4L, 0L, 0L,

    25L, 0L, 25L, 23L, 0L, 25L, 25L, 14L, 25L, 10L, 5L, 0L, 25L

    ), n = c(25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L,

    25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L,

    25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L,

    25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L,

    25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L,

    25L)), .Names = c("Site", "Quadrat", "Plot", "obs", "Burnt",

"soilgrp", "Them_tria", "n"), class = "data.frame", row.names = c(NA,

-59L))

#create factors

data$soilgrp<-as.factor(data$soilgrp)

data$sitef<-as.factor(data$Site)

data$Quadf<-as.factor(data$Quadrat)



#random effects model (a)

fm.a <- glmer(Them_tria/n~ Burnt+ soilgrp + Burnt: soilgrp

+ (1 | sitef/Quadf)+ (1|obs), data, family= binomial(link="logit"), weight=
as.numeric(n))



#random effects model (b)

fm.b <- glmer(Them_tria/n~ Burnt+ soilgrp + Burnt: soilgrp + (Burnt | sitef)

 + (1|obs), data, family= binomial(link="logit"), weight= as.numeric(n))



#random effects model (c)

fm.c <- glmer(Them_tria/n~ Burnt+ soilgrp + Burnt: soilgrp +

(1 | sitef) + (1|obs), data, family= binomial(link="logit"), weight=
as.numeric(n))





Any comments on the choices above or alternatives would be highly
appreciated .


Kind regards,


Josh



Josh Dorrough
Visiting Research Scientist
CSIRO Ecosystem Sciences

	[[alternative HTML version deleted]]

