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-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 own 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 poissonPoisson 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|Block1 | Block),
family = poisson (link=sqrtlink = sqrt), nAGQ=0nAGQ = 0, data=Xdata = 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 seemseems to work, but don't make much difference.
Does anyone have any ideas for what I should try next?