Everything I know about glmms is from the internet, and after extensive searching, I haven't come across a good clearcut guide for how to visualize your data in a way that is relevant to hypotheses and fits the model results. Essentially, if you just plot means for various groups, that doesn't account for random effects. I am seeking confirmation that using the lsmeans package makes sense but also, what specifically I should plot from lsmeans.
Data available: https://docs.google.com/spreadsheets/d/1qhq2qhwriotkRQUFYP1a_2oqBN_TQ37JhqOVFBis7Ds/edit?usp=sharing
Specifics of my data:
I studied the effect of four pollination treatments on the production of a specific type of flowers (CHFlwr). Size (BSD) and height (hgt) of plants may affect floral production so I included those in my models. There were max ~15 flowers per plant so I used Poisson errors, then followed with Tukey's post hoc using multcomp.
CHfoct<-glmer(POL$CHFlwr~Trtmt+(1|Hgt)+(1|BSD), family=poisson, data=POL)
multchoct<-glht(CHfoct, linfct = mcp(Trtmt = "Tukey" ))
Diagnostics plots seem normal (qqnorm plots/resid v fitted, etc). Results as follow:
> CHfoct<-glmer(POL$CHFlwr~Trtmt+(1|Hgt)+(1|BSD), family=poisson, data=POL)
> summary(CHfoct)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: poisson ( log )
Formula: POL$CHFlwr ~ Trtmt + (1 | Hgt) + (1 | BSD)
Data: POL
AIC BIC logLik deviance df.resid
954.6 972.4 -471.3 942.6 138
Scaled residuals:
Min 1Q Median 3Q Max
-2.5876 -0.5811 -0.0711 0.3527 3.2776
Random effects:
Groups Name Variance Std.Dev.
Hgt (Intercept) 0.3415 0.5844
BSD (Intercept) 0.1367 0.3698
Number of obs: 144, groups: Hgt, 103; BSD, 15
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.9717 0.1368 14.412 <2e-16 ***
TrtmtBagged 0.2486 0.1264 1.966 0.0492 *
TrtmtOx_Supplemented -0.1848 0.1250 -1.478 0.1393
TrtmtSelf_Supplemented -0.1406 0.1323 -1.063 0.2878
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) TrtmtB TrtO_S
TrtmtBagged -0.329
TrtmtOx_Spp -0.202 0.308
TrtmtSlf_Sp -0.253 0.425 0.275
> Anova(CHfoct)
Analysis of Deviance Table (Type II Wald chisquare tests)
Response: POL$CHFlwr
Chisq Df Pr(>Chisq)
Trtmt 11.747 3 0.008302 **
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> multchoct<-glht(CHfoct, linfct = mcp(Trtmt = "Tukey" ))
> summary(multchoct)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: glmer(formula = POL$CHFlwr ~ Trtmt + (1 | Hgt) + (1 | BSD), data = POL,
family = poisson)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
Bagged - Ambient_Pollinated == 0 0.24863 0.12644 1.966 0.1983
Ox_Supplemented - Ambient_Pollinated == 0 -0.18477 0.12499 -1.478 0.4473
Self_Supplemented - Ambient_Pollinated == -0.14059 0.13226 -1.063 0.7093
Ox_Supplemented - Bagged == 0 -0.43340 0.14790 -2.930 0.0179*
Self_Supplemented - Bagged == 0 -0.38922 0.13882 -2.804 0.0257*
Self_Supplemented - Ox_Supplemented == 0 0.04418 0.15502 0.285 0.9918
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
In the lsmeans package, I was able to get lsmeans for treatment but not a backtransformed value (because I used Poisson errors in my original model). Would it make sense to manually backtransform and simply graph those values with SE in another program?
In short, am I taking the appropriate direction and any hints on resources? I would much appreciate advice or direction to a tutorial/book/website for novices that could help.
