# post-hoc tests for a GLMM with polynomial term

I have count data (winter bud production) from a greenhouse experiment in which 48 plant genotypes were subjected to 4 salinity treatments ranging from low to high salinity. For optimal model fitting I converted Treatment to numeric and centered it using scale(), and called the new variable Treatment_Num.

I fit a zero-inflated poisson GLMM:

M8_tur_pois_nb <- glmmTMB(production ~ poly(Treatment_Num,2,raw=TRUE)+Source_Salinity +Rarity+poly(Treatment_Num,2,raw=TRUE)*Max_ramets_num+replanted+(1|Genotype)+(1|Unit/Replicate),
data = turions_produced_1_turion_planting,zi=~1,
family=poisson)

summary(M8_tur_pois_nb)

Family: poisson  ( log )
Formula:
production ~ poly(Treatment_Num, 2, raw = TRUE) + Source_Salinity +
Rarity + poly(Treatment_Num, 2, raw = TRUE) * Max_ramets_num +
replanted + (1 | Genotype) + (1 | Unit/Replicate)
Zero inflation:              ~1
Data: turions_produced_1_turion_planting

AIC      BIC   logLik deviance df.resid
2029.7   2085.8  -1001.8   2003.7      539

Random effects:

Conditional model:
Groups         Name        Variance Std.Dev.
Genotype       (Intercept) 0.08155  0.2856
Replicate:Unit (Intercept) 0.01176  0.1084
Unit           (Intercept) 0.11805  0.3436
Number of obs: 552, groups:  Genotype, 48; Replicate:Unit, 72; Unit, 24

Conditional model:
Estimate Std. Error z value
(Intercept)                                         1.59163    0.14688  10.836
poly(Treatment_Num, 2, raw = TRUE)1                -1.84356    0.17127 -10.764
poly(Treatment_Num, 2, raw = TRUE)2                -1.15687    0.17732  -6.524
Source_SalinityFresh                               -0.39566    0.10262  -3.856
RarityRare                                          0.05603    0.10198   0.549
Max_ramets_num                                      0.13522    0.04049   3.339
replantedYes                                       -0.20325    0.09780  -2.078
poly(Treatment_Num, 2, raw = TRUE)1:Max_ramets_num  0.17057    0.07176   2.377
poly(Treatment_Num, 2, raw = TRUE)2:Max_ramets_num  0.09091    0.05791   1.570
Pr(>|z|)
(Intercept)                                         < 2e-16 ***
poly(Treatment_Num, 2, raw = TRUE)1                 < 2e-16 ***
poly(Treatment_Num, 2, raw = TRUE)2                6.84e-11 ***
Source_SalinityFresh                               0.000115 ***
RarityRare                                         0.582698
Max_ramets_num                                     0.000840 ***
replantedYes                                       0.037696 *
poly(Treatment_Num, 2, raw = TRUE)1:Max_ramets_num 0.017456 *
poly(Treatment_Num, 2, raw = TRUE)2:Max_ramets_num 0.116433
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -4.0789     0.5499  -7.418 1.19e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


I had planned to do post-hoc pairwise comparisons (for example using emmeans) to find out if the counts differ from one treatment to the next. Is this appropriate with the polynomial term, or should I only look at differences among treatments visually, for example by plotting the fitted polynomials using sjPlot? Or both?

• Does it make sense to have Treatment_Num = 1.4, say? If not, I strongly recommend fitting it as a factor, not as a numeric predictor. And if so, then I think your post hoc analysis should focus on trends and curvature, not on comparing treatment means. Jun 24 at 15:54