Looking for a properly working post hoc test for my experimental analysis. Any suggestions?

Question

I posted my question earlier, but I didn't get an answer after I edited my question, so I'm posting it again.

I had to change from lmer to glmmTMB because of normality problems. The normality looks fine now, but the current post hoc test (emmeans - simple contrast) is not working somehow. Any suggestions what else I can use? The reviewers suggested especially the emmeans package, the Tukey comparison wasn't answering my question unfortunately so I decided to use the contrast. It would be great, but not working with my current models. What else might be good in my case?

I would like to check wheter there were significant differences between the treatment and control plots (TREATED variable) within each monitoring year (TIME.SINCE variable).

Here you can find the whole script and data set: https://www.dropbox.com/scl/fo/a54dku8dxwre5wmq8sx2k/h?rlkey=s6rw85iglxl68r1jfea6y2zpg&dl=0

The post hoc script:

library(emmeans)
#time factor!  
#post hoc test is not working with continous variables!!!!!!!! !!!!!! !!!!!
dt2$time.since <- as.factor(dt2$TIME.SINCE)
model2_4_ <- glmmTMB(sqrt(PW_c)~TREATED*time.since+(1|SITE/CODE),data=dt2)
summary(model2_4_)
AIC(model2_4_)
Anova (model2_4_)
#Response: sqrt(PW)
#Chisq Df Pr(>Chisq)
#TREATED 0.9165 1 0.3384012
#time.since 48.6197 11 1.107e-06 ***
#TREATED:time.since 33.8686 10 0.0001944 ***
#simple emmeans contrast
exp2_PW.emm <- emmeans(model2_4_, ~ TREATED * time.since, data=dt2)
contrast(exp2_PW.emm, "consec", simple = "each", combine = TRUE, adjust = "mvt")`

The contrast results:

> contrast(exp2_PW.emm, "consec", simple = "each", combine = TRUE, adjust = "mvt")
 time.since TREATED contrast                    estimate    SE  df t.ratio p.value
 1          .       TREATED1 - TREATED0           nonEst    NA  NA      NA      NA
 2          .       TREATED1 - TREATED0          -0.7130 0.415 254  -1.717      NA
 3          .       TREATED1 - TREATED0          -0.2078 0.323 254  -0.643      NA
 4          .       TREATED1 - TREATED0          -0.4262 0.254 254  -1.680      NA
 5          .       TREATED1 - TREATED0          -0.8181 0.218 254  -3.751      NA
 6          .       TREATED1 - TREATED0          -0.7933 0.211 254  -3.751      NA
 7          .       TREATED1 - TREATED0          -0.4185 0.207 254  -2.022      NA
 8          .       TREATED1 - TREATED0          -0.3727 0.190 254  -1.957      NA
 9          .       TREATED1 - TREATED0          -0.1256 0.164 254  -0.766      NA
 11         .       TREATED1 - TREATED0           0.0197 0.157 254   0.126      NA
 13         .       TREATED1 - TREATED0          -0.0538 0.135 254  -0.398      NA
 21         .       TREATED1 - TREATED0           0.1862 0.131 254   1.419      NA
 .          0       time.since2 - time.since1     0.4418 0.388 254   1.139      NA
 .          0       time.since3 - time.since2     0.3699 0.341 254   1.085      NA
 .          0       time.since4 - time.since3     0.1594 0.246 254   0.647      NA
 .          0       time.since5 - time.since4     0.2679 0.207 254   1.292      NA
 .          0       time.since6 - time.since5     0.1431 0.207 254   0.690      NA
 .          0       time.since7 - time.since6    -0.5211 0.202 254  -2.582      NA
 .          0       time.since8 - time.since7     0.0807 0.191 254   0.423      NA
 .          0       time.since9 - time.since8     0.0373 0.170 254   0.220      NA
 .          0       time.since11 - time.since9   -0.4131 0.143 254  -2.883      NA
 .          0       time.since13 - time.since11   0.4008 0.130 254   3.093      NA
 .          0       time.since21 - time.since13  -0.3651 0.133 254  -2.754      NA
 .          1       time.since2 - time.since1     nonEst    NA  NA      NA      NA
 .          1       time.since3 - time.since2     0.8752 0.358 254   2.447      NA
 .          1       time.since4 - time.since3    -0.0590 0.287 254  -0.206      NA
 .          1       time.since5 - time.since4    -0.1240 0.218 254  -0.568      NA
 .          1       time.since6 - time.since5     0.1679 0.167 254   1.006      NA
 .          1       time.since7 - time.since6    -0.1463 0.165 254  -0.889      NA
 .          1       time.since8 - time.since7     0.1266 0.157 254   0.808      NA
 .          1       time.since9 - time.since8     0.2844 0.135 254   2.099      NA
 .          1       time.since11 - time.since9   -0.2678 0.127 254  -2.106      NA
 .          1       time.since13 - time.since11   0.3272 0.118 254   2.781      NA
 .          1       time.since21 - time.since13  -0.1250 0.109 254  -1.145      NA

Note: contrasts are still on the sqrt scale 
P value adjustment: mvt method for 34 tests 

Summary of the problematic glmmTMB model:

> summary(modEXP2_4<-glmmTMB(sqrt(PW_c)~TREATED*TIME.SINCE+(1|SITE/CODE),data=dt2)) 
 Family: gaussian  ( identity )
Formula:          sqrt(PW_c) ~ TREATED * TIME.SINCE + (1 | SITE/CODE)
Data: dt2

     AIC      BIC   logLik deviance df.resid 
   418.3    443.7   -202.1    404.3      273 

Random effects:

Conditional model:
 Groups    Name        Variance Std.Dev.
 CODE:SITE (Intercept) 0.06559  0.2561  
 SITE      (Intercept) 0.02366  0.1538  
 Residual              0.20246  0.4500  
Number of obs: 280, groups:  CODE:SITE, 76; SITE, 3

Dispersion estimate for gaussian family (sigma^2): 0.202 

Conditional model:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)          1.376753   0.166070   8.290  < 2e-16 ***
TREATED1            -0.567529   0.163699  -3.467 0.000526 ***
TIME.SINCE          -0.018044   0.008721  -2.069 0.038551 *  
TREATED1:TIME.SINCE  0.037464   0.010879   3.444 0.000574 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Warning messages:
1: In sqrt(PW_c) : NaNs produced
2: In sqrt(PW_c) : NaNs produced

Thanks for your help.

I tried other post hoc test (e.g. Tukey pair wise), but it compares everything with everything in bulk, and that doesn't suit me.

Answer 1

You centered your outcome and then took the square root. However, centering PW results in many negative values, and taking the sqrt() produces a lot of NA values. That might be a reason for your problem here:

d <- read.csv(
  "c:/Users/Daniel/Downloads/data_2.csv",
  sep = ";",
  dec = ",",
  na.string = "na"
)
d$PW_c <- as.vector(scale(d$PW))
head(sqrt(d$PW_c))
#> [1] NaN NaN NaN NaN NaN NaN
#> Warning message:
#> In sqrt(d$PW_c) : NaNs produced

I'm not sure what the purpose of taking the sqrt is. Your data is highly right-skewed:

Taking the sqrt makes the skewness even worse:

What is the natural scale/range of your outcome? Maybe you could think about a Beta or Gamma regression.

Regarding the contrasts with numeric predictors, this should work:

library(glmmTMB)
library(ggeffects)
d <- read.csv("c:/Users/Daniel/Downloads/data_2.csv", sep = ";", dec = ",", na.string = "na")

d$PW_c <- as.vector(scale(d$PW))
d$ts <- d$TIME.SINCE

m <- glmmTMB(sqrt(PW_c) ~ TREATED * ts + (1 | SITE / CODE), data = d)
hypothesis_test(m, c("TREATED", "ts"), vcov = vcov(m), test = NULL)
hypothesis_test(m, c("TREATED", "ts"), vcov = vcov(m))