# Finding GLMMs that fit my count data for multiple datasets

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?

• What is block and how was the data generated? Is block a block/patch of soil with multiple plants that you randomly assigned to receiving the sane nitrogen level (then you may need Moore variation between individual plants)? Are there since theories how the relationship with nitrogen could look like- e.g. could it be non monotonic increase up to a point, then decrease again (then spline could be an idea)? – Björn Apr 5 '19 at 5:09
• 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 number of nodules with nodules, but a decrease at the higher levels. This is the case for 4 of my 5 datasets, there's an increase in nodules across the first two N levels, then nodules decrease. How would I go about applying a spline, then? Thanks! – Niall Millar Apr 5 '19 at 15:45
• Firstly, you may consider (1|Block) + (1|Plant), where plant is a factor with 24*15 =360 levels to capture whether data are overdispersed vs. the Poisson distribution, Secondly, a spline could be done e.g. using + bs(nitrogen,knots = c(1,5,10)) (no idea whether 1, 5 and 10 are sensible spots for knots). Also, you need a lot of levels of a continuous variable for a generic spline to make sense. I guess 24 combinations means something like 12 values vs 2 species? If so, great. If not, with e.g. 3 values vs. 8 species, a spline may not work and you may be better of treating N.level as a factor. – Björn Apr 5 '19 at 16:32
• There are a lot of posts here about count-data, so search this site! Specifically see. And, please add the new information in comments as an edit to the original Q, that way more eople will see it. Can you also include a (link to) the data? – kjetil b halvorsen Apr 5 '19 at 17:43

A couple of points:

• It is better to fit the model using the adaptive Gaussian quadrature with a sufficient number of quadrature points, e.g., 10 or 15. This would provide a better approximation of the log-likelihood of the model. You could also give a try in the GLMMadaptive package that can also fit random slopes with the adaptive Gaussian quadrature.
• From your description, it seems that variable N.Level has a nonlinear effect on the log expected counts. You could account for that using splines or polynomials (the former are prefered). For example, you could first load the splines package, and then define your formula as Nodules ~ ns(N.Level, 3) * Species + (1 | Block).