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.

enter image description here

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)?

New contributor
Dangerbuee is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.

1 Answer 1


Consider ordinal regression. You have data that are ordered but, from your description, it doesn't seem that the difference between scores of 1 and 2 is the same as the difference between scores of, say, between 4 and 5 or between 8 and 9. Frank Harrell recommends this as an approach when you have this type of data, even when there are so many ordered outcome levels that you might think of the outcomes as continuous. Chapter 13 of his Regression Modeling Strategies provides much detail.

The R GLMMadaptive package can fit mixed models with ordinal outcomes. This page shows how to proceed.

  • 2
    $\begingroup$ The ordinal package can also handle this kind of model. $\endgroup$
    – Ben Bolker
    May 26 at 1:30

Your Answer

Dangerbuee is a new contributor. Be nice, and check out our Code of Conduct.

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.