Linear mixed effects model on physiological parameter as function of time

Question

I have 1000 individuals with measurements of plasmatic HDL levels over the last 15 years. The number of measurements per individual ranges from 1 to 15, and were not taken in the same dates. I would like to compare the HDL levels (quantitative trait) over time in the individuals and make a ranking of those individuals whose HDL levels are consistently the highest over time. I want to adjust for sex and Tanner (categorical variable related to puberty development, it can be 1,2,3,4 or 5).

I have a data.frame (data) with 5 columns: ID (individual ID), HDL, Age (at the time of measurement), Sex and Tanner. This is a data frame with a toy example:

ID  HDL   Age   Sex  Tanner
1   50    12.3  M    2
1   52.1  15.4  M    4
1   55.3  17.1  M    5
2   45    12.1  M    1
2   46.3  13.1  M    1
3   60    14.3  F    3
3   55    16.2  F    5
... ...   ...   ...  ...

This is the model I figured out:

library(lme4)
model = lmer(HDL ~ Age + Sex + (1|Tanner) + (1|ID), data)

Age and Sex are my fixed variables and I expect to adjust for variation between Tanner. Individuals of the same Tanner should have more similar HDL values for reasons other than Age and Sex.

After running the regression, I extract the intercept and slope from each individual using coef(data)$ID. Then, I multiply each intercept by the corresponding slope and get the final value I will use for the ranking.

Is my reasoning correct?

Answer 1

You are on the right track, but you have to consider carefully the hidden assumptions in your model and whether you need to use a more complicated model to achieve your goals.

Treating ID as a random effect in this study with 1000 individuals makes a lot of sense. You have modeled this as (1|ID), however, so you are only allowing for intercepts to differ among individuals. That would seem to be consistent with your interest in "ranking ... individuals whose HDL levels are consistently the highest over time."

You say, however, that you wish to "extract the intercept and slope from each individual." (Emphasis added.) If that is the case, then you need to specify the slope (presumably for Age here) that you wish to allow to differ among individuals. Analyzing slopes might pose difficulties, however, with only 1 observation for some individuals.

Treating Tanner as a random effect, with only 5 levels, is not so obviously correct. If this is an ordered categorical variable for which you expect an ordered relation to HDL, you probably should model it as an ordinal predictor variable instead.

Specifying Age as a linear predictor makes sense if you think that HDL changes strictly linearly with age, other things being equal. I suppose that is possible, but it seems a bit unlikely. Often it makes more sense to model age in a more flexible way, for example with splines.

Finally, although your use of coef() to extract random effects from the model might work OK in this case, interpreting its output in more complicated models with nesting can be troublesome, as for example in this question. With the R lme4 package, the ranef() function might provide more intuitive results. It might make sense to look at both functions and make sure that you understand the types of results that they return.