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
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.