Creating nonlinear centiles/reference charts with longitudinal data

Question

I have a continuous dependent variable measured for several thousand subjects. The dependent variable has been measured one or more times per subject at nonuniform points in time. I am looking to generate nonlinear centile or reference charts. The structure of the data naturally direct me towards a mixed model approach.

GAMLSS::gamlss and GAMLSS:centiles appears to be used quite frequently to develop nonlinear centile or reference charts. The GAMLSS::gamlss also interfaces with NLME::lme to specify random effects in the GAMLSS model.

I have tried to produce a mixed effects model with the fixed effect of time and a random intercept and slope over time by subject with the following syntax, although it does not produce the anticipated linear centile plot and I am unsure how to add smoothing terms to this mixed model to produce nonlinear centile plots.

model <- gamlss(DV ~ re(fixed= ~ time, random= ~ 1 + time|subject), 
                    ...)

plot <- centiles(model, xvar=time)  

Is GAMLSS the right approach for this problem? If yes, what is the correct syntax for the mixed model and how are smoothing terms applied to this function? If GAMLSS is not the right approach, is there another recommended approach?

Answer 1

Centile estimation for longitudinal data is surprisingly tricky.

1. If the number of subjects is large and the number of measurements (per subject) is generally small, then:

a) One approach is to ignore the longitudinal aspect of the data and treat it as cross-sectional (i.e. independent) observations. This approach was taken by the World Health Organization when they developed their World Child Growth Standards charts (they had measurements on 8440 children), published in 3 books, 2006, 2007, 2009. The centiles (treating the data as independent observations) can be obtained using gamlss, see Stasinopoulos et al. (2017) Chapter 13 (which focusses on a single continuous explanatory variable), and Rigby et al. (2019) Section 5.9 (which focusses on two continuous explanatory variables), although in gamlss any number and type of explanatory variables are allowed.

b) If you are not happy with ignoring the longitudinal aspect of the data, then another approach is to randomly select one observation from each subject, so the data is then cross-sectional.

c) You could repeat b) m times, and then average the resulting m fitted centiles, or average the m fitted cdf’s and obtain the centiles from that average cdf.

I doubt whether b), and especially c), would give centiles which are much different from a), but it would be a useful check.

2. If the number of subjects is small and the number of measurements (per subject) is large, then a random effects model may be important for obtaining centiles.

However this is difficult as there are 2 problems:

a) How to find a suitable random effects model.

b) How to obtain the centiles of marginal distribution of Y (integrating out the random effects).

Answer 2

In "naturally direct me towards a mixed model approach" you have made a leap of faith. There is nothing natural about the use of random effects to handle longitudinal data except in the case where the time span is short and the multiple responses per subject represent repetitions under identical conditions. That would make the responses possibly qualify as being exchangeable so that their time order doesn't matter. Put another way, a random intercepts model assumes a compound symmetric correlation structure, i.e., the correlation between two responses within subject is independent of the time gap between the two measurements. As opposed to that the data are much more likely to follow a serial correlation pattern such as AR(1) and can be analyzed more simply without the use of random effects. Candidate approaches are generalized least squares and Markov processes. The Markov approach is the most general.

You can easily fit a first-order Markov process using a semiparametric regression model such as the proportional odds model. This has the additional advantage of not assuming a distribution for Y, and modeling its whole distribution. This will allow you to estimate means, exceedance probabilities, and quantiles (centiles). Details and extensive case studies are here.