I need some help finding a model that fits my data.
I'm using GLMM's to test the effect of nitrogen level (numerical), species (two-level factor), and the interaction between these two variables, on the number of nodules per plant (count data). My random effect is block.
(Edit: Each plant was in its pot, and received one of six N levels. Block is just randomly assigned position in the greenhouse. Each block contained one rep of every treatment combination (24 plants per block, 15 blocks). We do expect there to be an initial increase in the number of nodules with nitrogen, but a decrease at the higher levels.)
I have 5 datasets, and so far I've found that for two of them, a Poisson distribution with a square-root link function in lme4 fits reasonably well, but not perfectly (the residuals vs. fitted plots look odd). For my other datasets, however, this doesn't fit well at all.
Here's the code I'm using:
lme4::glmer(Nodules ~ N.Level + Species + N.Level:Species + (1 | Block), family = poisson(link = sqrt), nAGQ = 0, data = X)
I've tried Poisson distributions and negative binomial distributions, with and without zero-inflation, and other packages like
glmmTMB, but nothing seems to fit. I don't want to transform my nodule count data, but I have tried log transforms for my numerical $N$ level variable, which seems to work, but don't make much difference.
Does anyone have any ideas for what I should try next?