Including ordinal independent variables in a linear mixed effects model (using the lme4 package in R)

Question

I am analyzing data from cohort of 500 calves investigating the impact of disease on growth.

My outcome variables are normally distributed, continuous data. I am using hierarchical models with calf nested within farms and testing for the longer term impacts of disease.

The problem I am having is with how to include disease data. I have variables for the number of weeks a calf had disease and the total score over a validated threshold for diagnosis

As I am inexperienced in uploading images, here are the tabulated results of the data above:

---

Disease Duration (weeks) 0   1    2   3   4   5   6 
Frequency               266 128  50  33   8   5   2

---

Total Score  0   1   2   3    4   5   6   7   8   9  10  13  14  15 
Frequency   266  88  51  30  20  13   2   6   4   5   3   1   2   1 

---

Obviously, this data is far from normal. But there a lot of levels to use a dummy coded categorical variable, and I think an ordinal scale better represents the data. What do you think it the best way to include this data as an independent variable in my LME models? (n.b. I don't include both in the same model just one or the other)

The models do return results without convergence errors or other warnings when I include these variables but it doesn't feel like very good practice and I am unsure of what sort of transformation I could do to make this data better (e.g. log transformation leaves the data looking very odd and plots of the raw data make it look like a linear relationship is the most likely)

Here is an example of what I would like to improve:

(adj_w_63 - calf weight, weeks_brd - weeks with disease (as described above), rid - a normally distributed continuous variable, milksolids_total - a normally distributed continuous variable)

library(lme4)
model1<-lmer(adj_w_63 ~ weeks_brd + rid + milksolids_total + (1|farm_ac),
 data=comp)
summary(model1)

Linear mixed model fit by REML ['lmerMod']
Formula: adj_w_63 ~ weeks_brd + rid + milksolids_total + (1 | farm_ac)
   Data: comp

REML criterion at convergence: 3247

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.5180 -0.5525 -0.0458  0.5945  6.1674 

Random effects:
 Groups   Name        Variance Std.Dev.
 farm_ac  (Intercept) 30.10    5.487   
 Residual             83.37    9.131   
Number of obs: 443, groups:  farm_ac, 11

Fixed effects:
                 Estimate Std. Error t value
(Intercept)      68.06279    3.30996  20.563
weeks_brd        -1.00200    0.42089  -2.381
rid               0.11010    0.04981   2.210
milksolids_total  0.19904    0.07679   2.592

Correlation of Fixed Effects:
            (Intr) wks_br rid   
weeks_brd   -0.174              
rid         -0.285  0.141       
mlkslds_ttl -0.795  0.038 -0.016

Thank you so much for your help.

Answer 1

(This answer applies to [generalized] linear models generally, not just mixed models.)

This answer on SO discusses the interpretation of linear models with ordinal independent (predictor) variables. Here are two reasonable approaches, not clear which is best:

• treat the score as numeric. Advantages: simple, parsimonious (only takes one parameter). Disadvantages: assumes that the degree of change between each successive pair of scores is identical.
• convert the score to an ordered factor; in R, this automatically (by default) uses orthogonal polynomial contrasts. This will use the same number of parameters as treating the score as an unordered factor (and will give the same overall predictions, goodness-of-fit, etc.), but will give more interpretable parameters in terms of linear, quadratic, cubic ... terms. You may be able to reduce the number of terms (equivalent to using a lower-order orthogonal polynomial). Advantages: makes no assumptions about the size of differences. Disadvantages: less parsimonious.

Depending on what your goal is (prediction, hypothesis testing, etc.), you may be able to find some level of intermediate complexity (by regularizing/penalizing, or more crudely by reverting to lower-order orthog. polynomials, or by fitting a generalized additive model using a spline function of the scores), but the two options above are the simplest.