Grouping Factor throws model off

One grouping factor of my model appears to pose a problem for model fit. Let me elaborate:

I tried to fit a GLMM to ecological data for a behavioural study on termites. I did 80 experiments split between two castes (soldiers and workers) from 8 colonies across four sampling sites, investigating how aggression between different colonies might be linked to relatedness, their mating system and differences in caste.

I chose a Gamma distribution with a log link, since my response variable is strictly positive, right skewed and bound, this was also suggested by using descdist from the fitdistrplus package. I'm working with glmmTMB instead lme4, but when I tried glmer the same problem occured. The LHS of the model formula looks like this: meanW * Pairing + meanW * Costructure + Caste + (1 | YF). Pairing is coded as [0,1], were 0 indicates tests between nestmates and 1 tests between non-nestmates. meanW is a numeric measurement of relatedness, Costructure a factor of two levels, indicating whether a colony is polygamous or monogamous, and Caste a factor of the two levels, worker and soldier. YF represents the sampling site.

I'm fairly new to these approaches, so if parts of it are jumbled, forgive me, but my reasoning behind the interaction terms is this: relatedness is expected to vary within levels of control, since one would expect members of one colony to be closer related to each other than to non-nestmates. The same is true for MatingSystem, since individuals of a polygamous colony should be less related to each other than members from a monogamous colony. This is the model summary:

   Family: Gamma  ( log )
Formula:          rAi ~ meanW * Pairing + meanW * Costructure + Caste + (1 | YF)
Data: combdat

AIC      BIC   logLik deviance df.resid
-323.6   -302.3    170.8   -341.6       70

Random effects:

Conditional model:
Groups Name        Variance Std.Dev.
YF     (Intercept) 0.01489  0.122
Number of obs: 79, groups:  YF, 4

Dispersion estimate for Gamma family (sigma^2): 0.0392

Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept)        -2.80315    0.43162  -6.494 8.34e-11 ***
meanW               0.87847    0.72938   1.204  0.22843
Pairing1            0.22977    0.21451   1.071  0.28411
Costructure1        1.12301    0.41903   2.680  0.00736 **
Casteworkers        0.33556    0.04566   7.349 1.99e-13 ***
meanW:Pairing1     -2.46882    1.97951  -1.247  0.21233
meanW:Costructure1 -2.12110    0.78427  -2.705  0.00684 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Interpretation aside, the crux for me lies in the DHARMa-diagnostics results: the residuals are not distributed as expected and residuals plotted against predicted values deviate heavily.

However, if I try the same model on a dataset including only workers, therefore omitting the factor Caste entirely, the DHARMa diagnostics look like this, much more in line with what we'd want to see:

Should this be an indication to me that analysing both castes separately is the way to go? I'm not quite sure how to interpret this phenomenon, as I am just beginning to learn about these approaches. I understand that having a random factor of 4 levels is subpar, and I might be overlooking some obvious flaws in how I constructed the model. My main question here is how one would go about fixing this. Looking at my descriptive statistics, nothing but the response value seems to differ between the castes. The data for the soldiers includes one very obvious outlier, but I decided to not remove it, since there were no clear experimental reasons to take it out. I wouldn't think that one grouping factor can make such a big difference, considering that there's almost an equal number of tests for both castes across all other grouping factors.

If I'm being ignorant, excuse me, I have struggled with fitting this model for days now and am not quite sure what is going wrong. I'd be very happy if this turns out to be a very obvious mistake on my part.

A few points:

• It doesn't matter that the distribution of the response is skewed. What matters is the distribution of the residuals.

• It also doesn't necessarily matter that the response is strictly positive. Think about the distribution of height among, say, male students in a university. This is a typical example where we would expect a normal distribution even though height very clearly must be strictly positive and indeed very small values are basically impossible.

• Another problem here is that you only have four groups for the factor FY. This is definitely insufficient in my opinion. You are asking the software to estimate a variance for a normally distributed variable from only 4 observations, so this should be a fixed effect, not random. If you really want to treat it as random then a Bayesian approach would be more appropriate.

I would suggest starting from scratch with a simple linear model:

lm(rAi ~ meanW * Pairing + meanW * Costructure + Caste + YF, data = mydata)


and inspect the usual residual plots. If these are way off then think about transformations you could apply to make the residuals plausibly normal and homoskedastic. If this doesn't lead anywhere, then think about a gamma or beta GLM

• Thanks for the quick response, but this poses more questions to me: if the distribution of the response variable doesn't matter, why choose GLMs at all, and why even utilise tools like fitdistrplus? Or am I misunderstanding your point here? I acknowledged the problem with the random factor in the post, but since it appears to work well in the context of the model without soldiers, I thought that any "negative side-effects" by the low amount of levels wouldn't be too great. I will take the advise an try out the basic approach asap, thank you. May 22, 2021 at 20:42
• You're welcome. We use a GLM so that we can model the response variable via a link function and a distribution so that the variance of an individual response depends on its mean. When we have discrete data such as binary or count data, this is obvious but for continuous data it isn't always. It may be that a GLM fits the data better, in which case, of course we should use it, but it's not as simple as just looking at the distribution of the response. I hadn't come across fitdistrplus before and after googling it I can see what it does. I've never found a use for anything like that. May 22, 2021 at 20:49
• Yes you acknowledged the problem of too few groups, but it is not a matter of it being "sub par". It is much worse than that. But as I said if you adopt a Bayesian approach you could still use random intercepts, though I really don't know why you would want to. As Einstein once said "Everything should be made as simple as possible, but no simpler" and this is probably more true now than it ever was. May 22, 2021 at 20:51
• This clarifies a few things to me, thank you very much again. I will also try the Bayesian approach then. However, just out of interest and for future reference, are the big problems caused by the levels of YF apparent anywhere or just assumed? I would want to include random factors for sites since we cannot rule out any environmental factors within sites which might affect the response variable. May 22, 2021 at 21:00
• You're very welcome :) Think of the problem like this: if you had to estimate the variance in the height of students in a classrom, but you only sampled 4 students, you would not get a very good estimate. Maybe you don't care about this estimate in a mixed model, and that's fine, but if you have convergence issues, that could be the cause, for example. There is really no need to model sites as random for the reason you just said. It is perfectly OK to fit them as fixed. If this is intended to go for publication I would be very surprised if a reviewer was OK with it. May 22, 2021 at 21:12