rpairwiselsmeansglmmtmb

glmmTMB: Pairwise comparison post-hoc tests for factors interaction

In my example:

# Packages
library(glmmTMB)
library(DHARMa)
library(multcomp)
library(lsmeans)
library(car)

# My data set
ds <- read.csv("https://raw.githubusercontent.com/Leprechault/trash/main/temp_ger_ds.csv")
str(ds)
#'data.frame':  140 obs. of  4 variables:
# $ temp       : chr  "constante" "constante" "constante" "constante" ...
# $ generation : chr  "G0" "G0" "G0" "G0" ...
# $ development: int  22 24 22 27 27 24 25 26 27 18 ...

First fit the ziGamma model:

mTCFd <- glmmTMB(development ~ temp * generation, data = ds,
               family = ziGamma(link = "log")) 
Anova(mTCFd,test="Chi")
# Analysis of Deviance Table (Type II Wald chisquare tests)
# Response: development
#                   Chisq Df Pr(>Chisq)    
# temp            198.413  1  < 2.2e-16 ***
# generation       18.347  4   0.001056 ** 
# temp:generation  31.250  4  2.723e-06 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Check adjustment with DHARMa:

plot(s1 <- simulateResiduals(mTCFd))

karma Look like OK!!

Pairwise Comparison Post Hoc Tests:

1) For temp:
lsm.TCFd.temp <- lsmeans(mTCFd, c("temp"))
cld(lsm.TCFd.temp, Letters=letters)
#  temp      lsmean     SE  df lower.CL upper.CL .group
#  constante   3.18 0.0082 129     3.17     3.20  a    
#  flutuante   3.37 0.0131 129     3.34     3.39   b  

2) For generation:
lsm.TCFd.gen <- lsmeans(mTCFd, c("generation"))
cld(lsm.TCFd.gen, Letters=letters)
#  generation lsmean     SE  df lower.CL upper.CL .group
#  G3           3.23 0.0159 129     3.20     3.26  a    
#  G1           3.27 0.0198 129     3.23     3.31  ab   
#  G0           3.27 0.0135 129     3.25     3.30  ab   
#  G4           3.29 0.0217 129     3.25     3.34  ab   
#  G2           3.31 0.0141 129     3.28     3.34   b   

3) For temp:generation interaction:
ds$temp_gen <- paste0(ds$temp,"_",ds$generation)
mTCFd.int <- glmmTMB(development ~ temp_gen, data = ds,
               family = ziGamma(link = "log")) 
lsm.TCFd.temp.gen <- lsmeans(mTCFd.int, c("temp_gen"))
cld(lsm.TCFd.temp.gen, Letters=letters)
#  temp_gen     lsmean     SE  df lower.CL upper.CL .group 
#  constante_G3   3.13 0.0180 129     3.09     3.16  a     
#  constante_G2   3.14 0.0180 129     3.11     3.18  ab    
#  constante_G0   3.19 0.0191 129     3.15     3.23  abc   
#  constante_G1   3.22 0.0180 129     3.18     3.25   bc   
#  constante_G4   3.23 0.0185 129     3.19     3.27    cd  
#  flutuante_G1   3.32 0.0352 129     3.25     3.39    cde 
#  flutuante_G3   3.34 0.0262 129     3.28     3.39      e 
#  flutuante_G0   3.36 0.0191 129     3.32     3.39      e 
#  flutuante_G4   3.36 0.0393 129     3.28     3.44     def
#  flutuante_G2   3.47 0.0218 129     3.43     3.52       f

Ok it works, but I'd like to know if is possible for the pairwise comparison directly with the final model (mTCFd) without a new interaction model adjustment (mTCFd.int).

Please, any help with it?

Solution

  • You don't need to adjust another model but just only specify mTCFd factor's condition (c("temp","generation")).

    lsm.TCFd.temp <- lsmeans(mTCFd, c("temp","generation"))
cld(lsm.TCFd.temp, Letters=letters)
 temp      generation lsmean     SE  df lower.CL upper.CL .group 
 constante G3           3.13 0.0180 129     3.09     3.16  a     
 constante G2           3.14 0.0180 129     3.11     3.18  ab    
 constante G0           3.19 0.0191 129     3.15     3.23  abc   
 constante G1           3.22 0.0180 129     3.18     3.25   bc   
 constante G4           3.23 0.0185 129     3.19     3.27    cd  
 flutuante G1           3.32 0.0352 129     3.25     3.39    cde 
 flutuante G3           3.34 0.0262 129     3.28     3.39      e 
 flutuante G0           3.36 0.0191 129     3.32     3.39      e 
 flutuante G4           3.36 0.0393 129     3.28     3.44     def
 flutuante G2           3.47 0.0218 129     3.43     3.52       f

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 10 estimates 
significance level used: alpha = 0.05 
NOTE: Compact letter displays can be misleading
      because they show NON-findings rather than findings.
      Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.