[R-sig-ME] Plotting results of interaction from GLMM using model summary coefficients
Roslyn Anderson
r.anderson at umail.ucc.ie
Wed May 16 16:36:16 CEST 2012
Dear Thierry,
Many thanks for your email. I think I may be confusing coefficients and
predictions. Please see below the code that I used to plot the figure and
the figure as an attachment. My maximum value for the response variable,
plant species richness, is 26:
> veg$plant.spp.rich
[1] 15 14 11 10 11 11 11 15 12 15 14 12 14 16 14 5 4 5 17 12 7 12 11 9
8
[26] 7 10 12 13 12 7 12 11 20 15 15 13 8 9 18 16 15 15 13 17 14 9 11
23 26
[51] 16 5 9 7 6 12 10 13 12 11 11 11 9 10 8 8 11 14 14 12 10 9 12
15 13
[76] 11 9 7 8 14 13 9 11 11 12 18 11 16 13 10 13 14 15 10 9 9 8 7 7
11
[101] 10 9 7 9 11 12 13 22 11 10 10 16 15 13 9 11 9 9 9 10 9 5 7
11 11
[126] 11 13 24 18 14 11 6 18 21 19 11 10 10 14 13 16 12 8 11
> tapply(veg$plant.spp.rich, veg$TRADEWE.cat, max)
1 2 3 4
24 20 14 26
However, the intercepts and slopes that I have calculated for each
TRADEWE.cat (1-4) are all quite a bit below this.
> 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
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
l10():TRADEWE.3
log10(alt)
TRADEWE.ct2
TRADEWE.ct3
TRADEWE.ct4
l10():TRADEWE.2
l10():TRADEWE.3
l10():TRADEWE.4 0.282
# for TRADEWE.cat = 1
plant.spp.rich = 2.0134 + 0.1717*log10(alt) #intercept = 2.0134, slope =
0.1717
# for TRADEWE.cat = 2
plant.spp.rich = (2.0134 – 1.3003) + (0.1717 + 0.5322)*log10(alt)
#intercept = 0.7131, slope = 0.7039
# for TRADEWE.cat = 3
plant.spp.rich = (2.0134 + 0.5645) + (0.1717 – 0.2772)*log10(alt)
#intercept = 2.5779, slope = -0.1055
# for TRADEWE.cat = 4
plant.spp.rich = (2.0134 + 1.6082) + (0.1717 – 0.6733)*log10(alt)
#intercept = 3.6216, slope = -0.5016
> par(mfrow=c(2,2), mar=c(2,3,2,1), oma=c(3,1.5,0,1))
> plot(plant.spp.rich[TRADEWE.cat=="1"]~log10(alt)[TRADEWE.cat=="1"],
pch=1, xlim=c(1.0, 4.0), ylim=c(0,30), xlab="", ylab="", cex.axis=0.6,
main="0% traditional ewes", cex.main = 0.75, data=veg)
> abline(2.01, 0.17, lty=1, lwd=1)
> plot(plant.spp.rich[TRADEWE.cat=="2"]~log10(alt)[TRADEWE.cat=="2"],pch=2,
xlim=c(1.0, 4.0), ylim=c(0,30),xlab="", ylab="", cex.axis=0.6,main="50%
traditional ewes",cex.main = 0.75,data=veg)
> abline(0.71, 0.70, lty=2, lwd=2)
> plot(plant.spp.rich[TRADEWE.cat=="3"]~log10(alt)[TRADEWE.cat=="3"],pch=3,
xlim=c(1.0, 4.0), ylim=c(0,30), xlab="", ylab="", cex.axis=0.6,main="70%
traditional ewes",cex.main = 0.75, data=veg)
> abline(2.58, -0.12, lty=3, lwd=3)
> mtext(side=2, "plant species richness", cex=0.5, at =30, font=2, line=3)
> plot(plant.spp.rich[TRADEWE.cat=="4"]~log10(alt)[TRADEWE.cat=="4"],pch=4,
xlim=c(1.0, 4.0), ylim=c(0,30), xlab="", ylab="", cex.axis=0.6,main="100%
traditional ewes",cex.main = 0.75, data=veg)
> abline(3.62, -0.50, lty=4, lwd=4)
> mtext(side=1, expression("log"[10]*" altitude (m)"), cex=0.5, at =0,
font=2, line=3)
I haven’t back-transformed the data of the glmer with poisson. I apologise
for my basic knowledge of statistical procedures in R but what would be the
best way to go about doing this? Could I then use this to create regression
lines which run through the real data points?
Please let me know if you would like me to provide any further information.
Many thanks again for all your advice and guidance,
Best Wishes,
Roz
On 16 May 2012 14:06, ONKELINX, Thierry <Thierry.ONKELINX at inbo.be> wrote:
> 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 at 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 at r-project.org [mailto:
> r-sig-mixed-models-bounces at r-project.org] Namens Roslyn Anderson
> Verzonden: zaterdag 12 mei 2012 19:58
> Aan: r-sig-mixed-models at 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