Seeking feedback on the interpretation and reporting of a glmm output

Question

This might be a very beginner level question. I have fitted a model with glmmTMB R package, considering the interactive effect of mean wind speed and average daily rain as well as an interactive effect of mean temperature and mean relative humidity on plant disease severity. Here is my fitted model

mod1 <-
  glmmTMB(disease_severity ~ mean_rh * mean_temp + mean_ws * avg_daily_rain + (1|year),
          family = beta_family(), data = dat_seasonal)
summary(mod1)

Here is the model output.

Formula:          disease_severity ~ mean_rh * mean_temp + mean_ws * avg_daily_rain +  
    (1 | year)
Data: dat_seasonal
 AIC      BIC   logLik deviance df.resid 

-41.9    -27.4     29.9    -59.9       28
Random effects:
Conditional model:
 Groups Name        Variance Std.Dev.
 year   (Intercept) 0.7874   0.8874

Number of obs: 37, groups:  year, 11
Dispersion parameter for beta family (): 7.03
Conditional model:
                        Estimate Std. Error z value Pr(>|z|)

(Intercept)            -336.9044   151.8930  -2.218   0.0266 *
mean_rh                   4.3731     1.9177   2.280   0.0226 *
mean_temp                22.8305     9.9872   2.286   0.0223 *
mean_ws                 -19.8944    11.1741  -1.780   0.0750 .
avg_daily_rain           -3.7751     2.5388  -1.487   0.1370

mean_rh:mean_temp        -0.2854     0.1258  -2.269   0.0232 *
mean_ws:avg_daily_rain    5.0808     2.3155   2.194   0.0282 *

Signif. codes:  0 ‘*’ 0.001 ‘’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

For reporting,

I would like to report my model results. Can I say that the predictors mean_rh and mean_temp had a significant positive effect individual effect (effect on their own). The predictors mean_ws and avg_daily_rain had no significant effect on their own? There was a signficant negative interaction effect of mean_rh & mean_temp as well as a significant positive interaction effect of mean_ws and avg_daily_rain.

For plotting,

I just need to plot interactive terms of my model. Plotting individual effects alone would be misleading. For example, if I'm using ggeffects package.

Figure 1 <- plot(ggpredict(mod1, c("mean_rh", "mean_temp"))) +
  theme_pubclean() +
  labs(x = "Mean relative humidity", title = "")
Figure 2 <- plot(ggpredict(mod1, c("mean_ws", "avg_daily_rain"))) +
  theme_pubclean() + 
  labs(x = "Mean wind speed", title = "")

In summary, it's okay to report individual as well as interactive effects (main & interaction effect). But it is NOT okay to plot main effects alone in my case. Is that right?

Edit: Thank you very much for the feedback @EdM. I am editing the question based on the test you suggested, so all can benefit. I ran car::Anova(mod1) on my model and I got the following output.

Analysis of Deviance Table (Type II Wald chisquare tests)
Response: disease_severity
                        Chisq Df Pr(>Chisq)

mean_rh                0.0535  1   0.817057

mean_temp              0.2366  1   0.626650

mean_ws                1.1419  1   0.285252

avg_daily_rain         9.8801  1   0.001671 **
mean_rh:mean_temp      5.1498  1   0.023249 * 
mean_ws:avg_daily_rain 4.8148  1   0.028217 *

Signif. codes:  0 ‘*’ 0.001 ‘’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

This output suggest significant interaction effect of the interactive terms as shown in the original model above. However, it doesn't give a positive or a negative sign, but I assume there is a significant negative interaction effect of mean_rh:mean_temp and a significant positive effect of mean_ws:avg_daily_rain as shown in the original model above? So if need to report this, I'd report that there was a positive main effect of avg_daily_rain and significant interaction effect of mean_rh:mean_temp & mean_ws:avg_daily_rain. Is that right?

Regarding overfitting, I fitted the same model (but excluded random effect of year) using beta regression as shown below.

mod2 <-
  betareg(disease_severity ~ mean_rh * mean_temp + mean_ws * avg_daily_rain,
          data = dat_seasonal)
summary(mod2)

The output of interactive effects is the same as mod1 (above model with random effect)

Call:
betareg(formula = disease_severity ~ mean_rh * mean_temp + mean_ws * avg_daily_rain, 
    data = dat_seasonal)
Standardized weighted residuals 2:
    Min      1Q  Median      3Q     Max 
-1.8117 -0.9641 -0.0679  0.9726  1.8027
Coefficients (mean model with logit link):
                         Estimate Std. Error z value Pr(>|z|)

(Intercept)            -439.66169   87.71840  -5.012 5.38e-07 ***
mean_rh                   5.65518    1.11886   5.054 4.32e-07 ***
mean_temp                29.66389    5.84693   5.073 3.91e-07 ***
mean_ws                 -23.39753    6.14241  -3.809 0.000139 ***
avg_daily_rain           -4.89648    1.40426  -3.487 0.000489 ***
mean_rh:mean_temp        -0.36917    0.07362  -5.015 5.31e-07 ***
mean_ws:avg_daily_rain    6.39018    1.31982   4.842 1.29e-06 ***
Phi coefficients (precision model with identity link):
      Estimate Std. Error z value Pr(>|z|)

(phi)   3.7717     0.8953   4.213 2.52e-05 ***

Signif. codes:  0 '*' 0.001 '' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Type of estimator: ML (maximum likelihood)
Log-likelihood: 24.96 on 8 Df
Pseudo R-squared: 0.5381
Number of iterations: 358 (BFGS) + 56 (Fisher scoring) 

When I ran car::Anova(mod2). I got the following output

Analysis of Deviance Table (Type II tests)
Response: disease_severity
                       Df   Chisq Pr(>Chisq)

mean_rh                 1  0.5807   0.446035

mean_temp               1  3.7047   0.054259 .

mean_ws                 1  9.0713   0.002597 ** 
avg_daily_rain          1 32.5315  1.173e-08 ***
mean_rh:mean_temp       1 25.1495  5.305e-07 ***
mean_ws:avg_daily_rain  1 23.4420  1.287e-06 ***

Signif. codes:  0 ‘*’ 0.001 ‘’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

I'm happy to use the second model if ignoring year is not an issue, and increases the reliability of results. Interaction effects are the same in both models, but wind speed became significant too after running car::anova on the second model. The study was conducted in 11 years.

Answer 1

Can I say that the predictors mean_rh and mean_temp had a significant positive effect individual effect (effect on their own). The predictors mean_ws and avg_daily_rain had no significant effect on their own?

NO.

This is a common misinterpretation of the coefficients reported by linear modeling software when there are interactions. With continuous predictors, an individual coefficient represents the association with outcome when all its interacting predictors are at a value of 0. At least, you need to qualify your claims about "significance" accordingly.

More thoroughly, you need to evaluate each predictor including all of its interactions, plus evaluate the interactions. The default Type II analysis used by the Anova() function of the car package is a good choice for that. This vignette shows that car::Anova() works on this type of model. You should also look at the other methods in that vignette for documenting model adequacy. This answer illustrates some of that for a beta model with glmmTMB.

With interactions it's usually best to display model predictions at useful values of the covariates. I don't use the ggeffects package. Make sure that the ranges of all covariate values used for the displays make sense.

You have an additional issue to deal with: you are at risk of severely over-fitting your data with this model. As a rule of thumb, you need about 15 observations per coefficient (beyond the intercept) to avoid overfitting. It appears that you are trying to estimate 6 such coefficients plus a random-effect variance based on 37 observations. If I understand the software correctly, the coefficients are in a logit (log-odds) scale. Some of the coefficients look very large on that scale, which might be a sign of overfitting (although the magnitudes of coefficients depend on the scales of the predictors and can be misleading with interactions as you have here).