I am trying to compare the lodging resistance scores of different wheat cultivars in an agronomic trial. Lodging is the phenomenon in which wheat plant can bend and lean closer to the ground as a result of strong wind and rains, for example. My score goes from 1 (wheat plants are lying flat on the ground) to 9 (wheat plants proudly standing upright). My goal is to compare the means of different varieties for this score (4 reps per variety), for which I aim to run post-hoc tests on the fitted model. I do not wish to make predictions with this model.
My data is thus bound from 1 to 9 and is continuous (but quasi integer: I can get .5 or .25 but no more decimals). Another property of my data is that the variance is quite heterogeneous, with a tendency to increase when mean scores get closer to 1. As 1 is the lower bound, variance rapidly decreases when it gets there, but this is rare in my data set.
For now, I have tried fitting a generalized linear model with a Poisson distribution. Indeed, Poisson is bounded on one side (0) and variance=mean, which fits pretty well with my data after the following transformation:
mydata$TransfVar = abs(mydata$lodging_score-9)*100
(The idea here is to convert the 9 into a zero by subtracting the whole data by 9, then "flip" my data with abs(), and multiply it by 100 to get integer values only)
I then fit the model as follows:
model_lodging = glmer(data=my_data,formula = TransfVar ~ Variety + (1|Block), family = poisson())
This seems to work pretty well, as I get the following the following Fitted~Residuals plot, which is much better than when I first tried fitting a linear model.
I also get the following summary:
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: poisson ( log )
Formula: TransfVar ~ Variety + (1 | Block)
Data: my_data
AIC BIC logLik deviance df.resid
2745.6 2774.3 -1357.8 2715.6 35
Scaled residuals:
Min 1Q Median 3Q Max
-17.0608 -5.0263 -0.9449 5.0668 14.3831
Random effects:
Groups Name Variance Std.Dev.
Block (Intercept) 0.01485 0.1219
Number of obs: 50, groups: Block, 4
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.74225 0.09796 38.203 < 2e-16 ***
VarietyCH - PA 2.65454 0.08036 33.034 < 2e-16 ***
VarietyChiddam d'automne a epi blanc 2.16400 0.08098 26.724 < 2e-16 ***
VarietyJ - C 1.57351 0.08427 18.673 < 2e-16 ***
VarietyJ - PA 2.39522 0.08011 29.899 < 2e-16 ***
VarietyJ - PA - RSL - C 1.23198 0.09045 13.620 < 2e-16 ***
VarietyJ - RSL 1.97698 0.08323 23.753 < 2e-16 ***
VarietyJaphet 2.55357 0.07962 32.073 < 2e-16 ***
VarietyPA - R -0.35157 0.12817 -2.743 0.00609 **
VarietyPA - RSL 1.76694 0.08298 21.293 < 2e-16 ***
VarietyPA - RSL - S - R 0.62584 0.10038 6.235 4.52e-10 ***
VarietyPrince Albert 2.56268 0.07959 32.198 < 2e-16 ***
VarietyRouge de St-Laud 0.25783 0.10211 2.525 0.01157 *
VarietyRSL - S -0.30029 0.12610 -2.381 0.01725 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
However, this feels quite unorthodox. My data is not count data, is not integer, and has an upper bound, unlike Poisson.
Another issue is that some of my varieties did not lodge at all and have all 9 as lodging scores (so no variance), which glmer does not like. From my understanding, a zero-inflated model would not apply in this case. For now my solution has been to remove these groups (hence the problem is not apparent in the above summary).
Are there other distribution families that could better match the type of data I have here, or is my approach valid for the intended use (mean comparisons between groups, not predictions)?
