Dear List Members,

I have a query regarding the plotting of results from a GLMM model.

I am working with both linear mixed effects models (LMMs) and generalised
linear mixed models (GLMMs), mainly using the nlme and lme4 packages in R
(2.14.2).

The study involved looking at the effects of *altitude* and the *percentage
of traditional sheep on a farm* on various measures of plant diversity.
Here I will concentrate on the response variable *plant species
richness*(which ranges from 4 – 26). It is count data so I have used a
Poisson
distribution.

The experimental design had a nested structure, so I chose to use LMMs and
GLMMs. Model syntax for plant species richness is:

*SR <- glmer (plant.spp.rich ~ log10(alt) * TRADEWE.cat +
(1|site/farm/habitat), family = poisson, data=veg)*

Fixed effects are:
1) log10(alt) - a log10 transformation on altitude (m) (range: 1.56 – 2.78)
and
2) TRADEWE.cat - the percentage of traditional sheep on a farm split into
four categories (1 = 0%, 2 = 50%, 3 = 70%, 4 = 100%). Groups are of unequal
size, therefore data is unbalanced.

Random effects are:
habitats (categories 1-12), which are nested within farms (categories
1-12), which are nested within sites (categories 1-4).

After rigorous model simplification and validation I am happy with the
model summary output and the interaction between log10(alt) and
TRADEWE.cat1 is significantly different to the interaction between
log10(alt) and TRADEWE.cat4. Therefore I wanted to plot this result with
plant species richness on the y-axis, log10(alt) on the x-axis and each of
the four traditional sheep categories as regression lines on the plot.

See summary below from the GLMM model and then the output of a linear
regression model containing the same variables:

> *library(lme4)*

> *SR<-glmer(plant.spp.rich ~ log10(alt) * TRADEWE.cat +
(1|site/farm/habitat), family=poisson, data=veg)*

> *summary(SR)*

Generalized linear mixed model fit by the Laplace approximation
Formula: plant.spp.rich ~ log10(alt) * TRADEWE.cat + (1 |
site/farm/habitat)
Data: veg
AIC   BIC logLik deviance
146.3 178.9 -62.13    124.3
Random effects:
Groups                       Name        Variance Std.Dev.
habitat:(farm:site)    (Intercept)  0.028298 0.16822
farm:site                    (Intercept)  0.000000 0.00000
site                            (Intercept)   0.000000 0.00000
Number of obs: 144, groups: habitat:(farm:site), 48; farm:site, 12; site, 4

Fixed effects:
                                                Estimate Std. Error z value
Pr(>|z|)
(Intercept)                                2.0134     0.5203   3.870
0.000109 ***
log10(alt)                                  0.1717     0.2365   0.726
0.467844
TRADEWE.cat2                     -1.3003     1.7218  -0.755 0.450117
TRADEWE.cat3                     0.5645     1.5340   0.368 0.712872
TRADEWE.cat4                     1.6082     0.6260   2.569 0.010200 *
log10(alt):TRADEWE.cat2   0.5322     0.6907   0.771 0.440985
log10(alt):TRADEWE.cat3   -0.2772     0.7047  -0.393 0.694059
log10(alt):TRADEWE.cat4   -0.6733     0.2814  -2.393 0.016720 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
               (Intr) lg10() TRADEWE.2 TRADEWE.3 TRADEWE.4 l10():TRADEWE.2
l10():TRADEWE.3
log10(alt)      -0.991

TRADEWE.ct2     -0.302  0.299

TRADEWE.ct3     -0.339  0.336  0.102

TRADEWE.ct4     -0.831  0.824  0.251     0.282

l10():TRADEWE.2  0.339 -0.342 -0.997    -0.115    -0.282

l10():TRADEWE.3  0.333 -0.336 -0.100    -0.996    -0.276     0.115

l10():TRADEWE.4  0.833 -0.841 -0.252    -0.282    -0.991     0.288
  0.282


> *SRa<-lm(plant.spp.rich ~ log10(alt) * TRADEWE.cat, data=veg)*

> *summary(SRa)*

Call:
lm(formula = plant.spp.rich ~ log10(alt) * TRADEWE.cat, data = veg)

Residuals:
  Min     1Q Median     3Q    Max
-7.955 -2.091 -0.221  1.943 12.404

Coefficients:
                                                Estimate Std. Error t value
Pr(>|t|)
(Intercept)                                 7.092      4.421   1.604
 0.11096
log10(alt)                                   1.837      2.017   0.911
 0.36397
TRADEWE.cat2                      -17.959     15.078  -1.191  0.23570
TRADEWE.cat3                       6.066     12.795   0.474  0.63618
TRADEWE.cat4                      19.437      5.365   3.623  0.00041 ***
log10(alt):TRADEWE.cat2     7.318      6.058   1.208  0.22913
log10(alt):TRADEWE.cat3    -3.025      5.875  -0.515  0.60747
log10(alt):TRADEWE.cat4    -8.101      2.409  -3.363  0.00100 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.54 on 136 degrees of freedom
Multiple R-squared: 0.1819,     Adjusted R-squared: 0.1398
F-statistic: 4.319 on 7 and 136 DF,  p-value: 0.0002399

My question is why, when using a linear regression model, my coefficients
(intercept and slope) for each traditional sheep category fit nice
regression lines within the data points on the plot but when using an LMM
or a GLMM the coefficients are much lower (and therefore the regression
lines are also much lower) than the actual data points on the plot? Is this
related to the random effects?

I was just wondering if anyone could suggest an alternative method I might
use to plot this result using the coefficients from the GLMM output?

Many, many thanks in advance for any help, advice or guidance.

Best Wishes,

Roz

	[[alternative HTML version deleted]]

