0

I am currently working on a generalized mixed effects model with 3 fixed effects and one random effect. Just as background, the random effect is the pond number (mesocosm), the response variable is the amount of chlorophyll-a (continuous data) and my fixed effects are two treatments: 1)nutrient enrichment and 2) WLF, both of which are binary yes/no type data. the third fixed effect is depth (0.1 meters and 0.8 meters) Anyway, I've created a glmer where all three fixed effects and the random effect have full interaction as follows:

GHQ<-glmer(CHa.t ~ Nutrients * WLF * depth * (1|mesocosm), family = gaussian(link = "log"), data = microalgae, nAGQ = 100)

A summary of this model shows that there is a significant interaction between nutrient enrichment, wlf and depth on Chlorophyll-a.

My question is how to carry out a post-hoc test when there are interactions between three fixed terms? I have tried glht and ls_means but both have resulted in an error.

here is the relevant data (sorry, I had to copy all of the rows for it to work):

structure(list(mesocosm = structure(c(2L, 2L, 10L, 10L, 15L, 15L, 18L, 18L, 22L, 22L, 23L, 23L, 24L, 24L, 5L, 5L, 6L, 6L, 8L, 8L, 9L, 9L, 14L, 14L, 25L, 25L, 26L, 26L, 1L, 1L, 3L, 3L, 7L, 7L, 11L, 11L, 13L, 13L, 16L, 16L, 20L, 20L, 4L, 4L, 12L, 12L, 17L, 17L, 19L, 19L, 21L, 21L, 27L, 27L, 28L, 28L), .Label = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28"), class = "factor"), Nutrients = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("0", "1"), class = "factor"), WLF = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("0", "1"), class = "factor"), depth = c(0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8, 0.1, 0.8), Chl.a = c(3.301666667, 4.026666667, 2.54, 3.936666667, 0.933333333, 1.67, 1.841666667, 1.968333333, 1.493333333, 3.321666667, 2.78, 2.078333333, 3.633333333, 2.621666667, 2.641666667, 4.991666667, 1.521666667, 2.118333333, 3.061666667, 6.083333333, 2.343333333, 4.406666667, 2.243333333, 3.03, 5.155, 3.723333333, 2.43, 2.06, 8.943333333, 5.346666667, 6.161666667, 4.730833333, 7.336666667, 6.958333333, 8.888333333, 1.898333333, 9.978333333, 7.951666667, 10.09333333, 7.516666667, 10.62666667, 3.043333333, 4.621666667, 7.798333333, 2.79, 7.926666667, 4.403333333, 3.85, 5.581666667, 8.18, 2.795, 5.946666667, 4.093333333, 8.5, 7.356666667, 9.861666667 )), row.names = c(NA, 56L), class = "data.frame")

here are the packages I used:

if(!require(psych)){install.packages("car")}
if(!require(MASS)){install.packages("MASS")}
if(!require(rcompanion)){install.packages("rcompanion")}
if(!require(lme4)){install.packages("lme4")
if(!require(mlmRev)){install.packages("mlmRev")}
if(!require(ggplot2)){install.packages("ggplot2")}
if(!require(nlme)){install.packages("nlme")}
if(!require(agricolae)){install.packages("agricolae")}
if(!require(multcomp)){install.packages("multcomp")}

here is the chunk of code I used for the model:

CHa=microalgae$Chl.a
library("mlmRev")
GHQ<-glmer(CHa ~ Nutrients * WLF * depth * (1|mesocosm), family = gaussian(link = "log"), data = microalgae, nAGQ = 100)
summary(GHQ)
overdisp_fun(GHQ)
plot(GHQ)
qqnorm(resid(GHQ))
qqline(resid(GHQ))

here's my code for "ls_means" and the error that I got:

library(lmerTest)
ls_means(GHQ, pairwise~ depth*Nutrients*WLF, adjust="tukey")

Error in UseMethod("ls_means") : no applicable method for 'ls_means' applied to an object of class "c('glmerMod', 'merMod')"

I actually realise now that I don't really know how to use glht.

Sinah Behsangar
  • 1
  • 1
  • You can probably use the `emmeans` package. Then run something like `emmeans(CHQ, pairwise ~ Nutrients + WLF + depth)`, but note that you quickly get a lot of contrasts. – Axeman Feb 11 '20 at 21:08
  • Are you seeking modeling advice? If so then you should ask such questions over at [stats.se] instead. If you are seeking programming advice then show your code with a [reproducible format](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) and state exactly what errors you are getting. – MrFlick Feb 11 '20 at 21:09
  • Please provide a reproducible example and your code for ghlt and ls_means. – Nakx Feb 12 '20 at 02:13
  • thank you for your replies. I have edited the post to include the data I used and error I got. Axeman, `emmeans(CHQ, pairwise ~ Nutrients + WLF + depth)` also gices an error saying "object 'CHQ' not found" – Sinah Behsangar Feb 12 '20 at 17:40

0 Answers0