[R-sig-ME] Plotting results of interaction from GLMM using model summary coefficients

ONKELINX, Thierry Thierry.ONKELINX at inbo.be
Wed May 16 15:06:38 CEST 2012


Dear Roslyn,

It seems to me like you are confusing coefficients and predictions. How did you make the plots? Did you back-transform the data of the glmer with poisson? Note that the poisson family uses the lok-link by default.

A reproducible example would ben ice.

Best regards,

Thierry

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
+ 32 2 525 02 51
+ 32 54 43 61 85
Thierry.Onkelinx op inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey


-----Oorspronkelijk bericht-----
Van: r-sig-mixed-models-bounces op r-project.org [mailto:r-sig-mixed-models-bounces op r-project.org] Namens Roslyn Anderson
Verzonden: zaterdag 12 mei 2012 19:58
Aan: r-sig-mixed-models op r-project.org
Onderwerp: [R-sig-ME] Plotting results of interaction from GLMM using model summary coefficients

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]]

* * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * *
Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document.
The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document.



More information about the R-sig-mixed-models mailing list