We use individual child data from the MAL-ED study, whose protocols, methodology and aggregate growth results have been presented previously ^{8, 17} . MAL-ED was initiated as a multi-country cohort study located across eight low- and middle-income sites with historically high incidence of diarrhoeal disease and undernutrition, with a research focus on investigating determinants of development in children from birth through early years. Ethical approval was obtained from the National Institute of Medical Research for Tanzania (NIMR/HQ/R.8a/Vol.IX/858 and NIMR/HQ/R.8c/Vol.II/1034) and University of Venda Ethics Committee for South Afirca (SMNS/09/MBY/004). Approval was additionally given by the institutional review board of the University of Virginia, USA and all methods used in this study followed the relevant guidelines and regulations. Informed consent was taken from parents of all children prior to enrolment.

In our case studies we use data from Haydom, Tanzania ¹⁹ (n = 224) and Venda, South Africa ²⁰ (n = 236). Our focus here is on anthropometric data from ages 0–24 months where each child included in these analyses had between 20 and 25 monthly observations (within a window of ±14 days), with 83% of children in Haydom and 86% in Venda having at least 24 observations. Weight and length, collected by trained fieldworkers and with minimal measurement error ⁸ , were converted to age and sex standardised z-scores using the WHO growth standards ²¹ . Here we focus on weight-for-length data, reflective of the relative weight of a child given their stature and therefore a child’s current nutritional status ²² , and one of the growth z-scores recommended by the WHO for diagnosing acute malnutrition ²³ . Figure 1 shows all trajectories for weight-for-length z-score (ZWfL), along with the site-specific cross-sectional means.

High within- and between-child variability is the predominant feature of the raw trajectory data, which holds from birth through 24 months and for both sites. Three trajectories are highlighted in each site, and these illustrate the dynamic complexity of each child’s growth.

The standard introductory text for SDEs is by Øksendal ⁵ , which contains a detailed mathematical exposition of SDEs. We focus on application and SDE models are introduced through examples with technical details largely omitted. We begin with a well-studied special case SDE model which we fit to data from three individual children (separately) and compare results with a linear regression model. We then introduce a more general SDE model and fit this to data from all children from the Haydom site.

Example 1 – Ornstein-Uhlenbeck (OU) process for a single trajectory . Our growth outcome of interest is ZWfL, where the first observation in the growth trajectory for a single child is X ₀ when the child is aged t ₀. We now develop a model to estimate likely values of X ₁, ZWfL at a later age, t ₁. In the simplest linear context, we could formulate an expression for the mean of X ₁ conditional on the previous observation, X ₀, and the difference in age between t ₁ and t ₀. Assuming X ₁ is approximately normally distributed with mean μ and variance σ ², N(μ,σ ²), this would give a model of the form X 1 | X 0 ∼ N ( X 0 + b ( t 1 − t 0 ) , σ 2 ) , ( 1 )

which models ZWfL at a subsequent age as a linear function of the current value and the elapsed age. In this model, growth velocity – rate of change per unit time (age) – is described by b, a parameter to be estimated from the observed data. This is a fully specified growth model between any two time points, although too simple to be practically useful; firstly, the growth velocity, b, is assumed constant, and secondly, the variance σ ², is constant and independent of the duration of the elapsed age between t ₁ and t ₀, whereas it might reasonably be expected that two time points closer together in age may be more similar than those further apart.

Consider the same example as above but now where we have X 1 | X 0 ∼ N ( e ( t 0 − t 1 ) α ( X 0 − β ) + β , e − 2 t 1 α ( e 2 t 1 α − e 2 t 0 α ) σ 2 2 α ) . ( 2 )

This model now has three parameters, α, β and σ, and ( 2) is the transitional probability density function (or slice density) for the OU stochastic process. The OU model is well studied and has applications in mathematical finance ²⁴ and theoretical biology ²⁵ . The distribution in ( 2) is a more flexible model compared to ( 1), with non-linear mean and variance terms, both of which depend on the elapsed time between t ₁ and t ₀. The OU process is typically defined in differential form by the stochastic differential equation dX t = α ( β − X t ) dt + σdW t , ( 3 )

where α(β − X _t) is called the drift, σ is the diffusion and W _t is a Wiener (Brownian motion) process. The drift can be thought of as the slow-moving trend in growth velocity, while the diffusion is the continual perturbation of the system giving rise to volatility in velocity and therefore growth. This SDE has correlated movements through time, for example, in an OU process that commenced at time t ₀ the covariance between any two points in time, t and s, is (e ^−(s+t)α(−e ^{2t ₀α} + e ^{2α Min[s,t]})σ ²)/(2α). Equation (3) provides the interpretation of the model parameters as components of the rate of change of growth, where ( 2) can be derived from ( 3) by solving the Fokker-Planck equation (see later).

To demonstrate the practical application of an SDE model to real data, we fitted OU models to three different ZWfL trajectories from the Haydom data. We compared these with the fit of a classical linear regression (LR) model, e.g. ZWfL(t) ~ N(a ₀ + a ₁t, σ ²), with both models having three parameters. Figure 2 shows the raw trajectory data, along with fitted values from the OU and LR models (separately for each child). See SI.1 ( Extended data) ¹⁸ for the model fitting code. In the OU model each successive data point (say, ZWfL ₁) depends on the previous point (ZWfL ₀), which means the OU model fit is not a smooth curve, but transitions from point to point. In the LR model, the fitted value at age t ₁ does not depend directly on the previously observed value ZWfL ₀, but rather it is from a globally parameterized smooth curve computed across all observations. The OU model is a special case of a more general SDE formulation, which we present next, along with the key concepts in fitting SDEs to data.

Defining an SDE model – basic concepts . A general formulation for an SDE model (see 7.1.2 in Øksendal ⁵ ) is dX t = μ ( t , X t , θ ) dt + g ( t , X t , θ ) dW t , ( 4 )

where X _t is our growth outcome variable of interest (e.g. ZWfL). Functions μ(t,X _t,θ) and g(t,X _t,θ) generalise the drift and diffusion terms from ( 3), respectively, where θ is a set of parameters to be estimated and, as before, W _t is a Wiener process. The OU model from ( 3) is therefore a special case of ( 4) where μ(t,X _t,θ)=α(β−X _t) and g(t,X _t,θ)=σ. Equation (4) can also be written as a stochastic integral equation X t = X 0 + ∫ 0 t μ ( s , X s , θ ) ds + ∫ 0 t g ( s , X s , θ ) dW s , ( 5 )

which usefully emphasises that these are models of the evolution of a continuous time stochastic process – here growth of a child. Equation (5) says that the value of growth (a random variable) for a child at t time units into the future from the currently observed time ( t=0) is the current value of growth (X ₀) plus the sum over this time interval of the drift and diffusion components of velocity. As the diffusion is integrated with respect to a “white noise” process W _s, then X _t follows a probability distribution at time t. To make this more concrete, if ZWfL _t1 is our growth outcome at time t1, then f(ZWfL _t1) is the probability distribution of ZWfL _t1 conditional on ZWfL _t0, the previously observed value of growth at time t0. This function, f(X _t), which is essential for model fitting, is a solution to the partial differential equation ∂ f ∂ t = ∂ ∂ x ( 1 2 ∂ ( g 2 ( t , X t , θ ) f ) ∂ x − μ ( t , X t , θ ) f ) . ( 6 )

Equation (6) is called the Fokker-Planck or forward-Kolmogorov equation (whose complete specification includes initial and boundary conditions, which have been omitted). The SDE in equation ( 4) defines the terms in Equation (6). Solving Equation (6) gives the expression for f( X _t). If we consider again the OU process and plug in the relevant terms from ( 3) into ( 6) and do the necessary calculus, then we arrive at exactly the normally distributed slice density stated in ( 2).

For model fitting we need to compute the likelihood function given the slice density. If we consider first the likelihood function for trajectory data from a single child, and where we have N observations over time, then the negative log-likelihood for a single child can be written as (see Hurn ¹⁰ ) − log ℒ ( θ ) = − log ​ f 0 ( X 0 | θ ) − ∑ k = 0 Ν − 1 log f ( X k + 1 | X k , θ ) , ( 7 )

where f ₀(X ₀|θ) is the probability density of the growth outcome variable at the first available data point, f(X _k+1|X _k,θ) ≡ f((X _k+1, t _k+1) | (X _k, t _k), θ) is the value of the slice density function for a stochastic process starting at (X _k,t _k) and evolving to (X _k+1,t _k+1). Equation (7) allows us to compute the likelihood function for all trajectories, including, if necessary, covariance structures across children through the inclusion of random effects (see later). One minor remark is how to deal with the first available observation, X ₀, as SDE models are defined in terms of transitions. We follow Schneider ²⁶ and the existing literature in maximum likelihood estimation in SDEs and treat X ₀ as a constant, which is also typical in the time series literature.

In summary, the key steps for working with SDE models focussed on model fitting are: (i) choose a form of μ(t,X _t,θ) and g(t,X _t,θ) in ( 4) and then; (ii) determine the corresponding slice density which satisfies the Fokker-Planck Equation (6) and then; (iii) fit the model to the data using the likelihood function in ( 7). Here we only consider forms of μ(t,X _t,θ) and g(t,X _t,θ) which have known analytical solutions (slice densities) to ( 6), which then makes fitting SDE models to data no different from a standard maximum likelihood problem using standard statistical software. Mathematica software (version 11.3, Wolfram Research Inc.), for example, can be used to compute slice densities for a wide range of SDEs, and a selection of these solutions is provided for reference in SI.2 (see Extended data) ¹⁸ as illustration. While in theory it is possible to fit SDEs that do not have an analytical solution to (6) to data, in practice this is numerically challenging (see the Discussion).

Example 2 - OU process for multiple trajectories . In the Haydom data we have 224 trajectories (children) across 0–24 months, and we now fit (non-linear) mixed model variants of the OU process, along with standard linear mixed models to these data. Table 1 gives a summary of different parameterizations and goodness of fit metrics. The modelling code is provided in SI.3a along with model output SI.3b (see Extended data) ¹⁸ .

For the same number of parameters, the OU process gives substantially better Akaike information criterion (AIC) and Bayesian information criterion (BIC) metrics, and fitting these SDE models including random effects is straightforward, requiring only with a few lines of code in SAS's proc nlmixed. These mixed models can also be implemented in the Stan ²⁷ language, with an OU specific example using Stan provided by Goodman (2018) ²⁸ .

Our main results comprise of an illustrative case study where we considered the combined data from Haydom and Venda. The general model formulations considered, and model search process are detailed below. To keep the analysis as clear as possible we considered only one covariate, (in addition to age) in the modelling, a categorical variable indicating site.

Linear mixed models (LMMs) . We considered LMMs where the most general formulation for ZWfL (response) was: a polynomial function of age (continuous) and site (two categories); with interactions between age and site; fixed and random (normally distributed) effects for the age terms in the polynomial; where random effects for the age terms (including intercept) had an unstructured covariance matrix; and within child errors had an AR(1) – autoregressive first order – covariance structure to allow serial dependence between errors. Increasing orders of polynomial (up to fourth order) were examined, guided by AIC and BIC metrics, in addition to trimming terms with high p-values (>0.1). These model formulations can be readily fitted to trajectory data using proc mixed in SAS. Relevant SAS code is provided in SI.4a (see Extended data) ¹⁸ . This can also be achieved using Stan, for example adapting Goodman (2018) ²⁸ .

SDE (non-linear) mixed models . The most general formulation considered for the SDE models was d X t = ( α 2 + α 3 t + α 4 t 2 + α 5 t 3 + α 1 X t ) dt + σ d W t , ( 8 )

where X _t is ZWfL at age t and so the rate-of-change per unit time for ZWfL depends on both the current age of the child and the child’s current value of ZWfL. Specifically, we considered a polynomial (a cubic in equation 8) function of age, plus a linear term (α ₁X _t) in ZWfL. This polynomial formulation of SDE in ( 8) has a closed form of slice density (a normal distribution) if we hold (α ₁, α ₂, α ₃, α ₄, α ₅, σ) constant. Mathematica was used to determine the slice density which is

X 1 | X 0 ∼ N ( 1 α 1 4 e − t 0 α 1 ( e t 1 α 1 x 0 α 1 4 − 6 ( e t 0 α 1 − e t 1 α 1 ) α 5 + α 1 ( − 2 ( e t 0 α 1 − e t 1 α 1 ) α 4 − 6 ( − e t 1 α 1 t 0 + e t 0 α 1 t 1 ) α 5 )   +   α 1 2 ( ( − e t 0 α 1 + e t 1 α 1 ) α 3 +   2 e t 1 α 1 t 0 α 4 + 3 e t 1 α 1 t 0 2 α 5 − e t 0 α 1 t 1 ( 2 α 4 + 3 t 1 α 5 ) ) + α 1 3 ( ( − e t 0 α 1 + e t 1 α 1 ) α 2 + e t 1 α 1 t 0 α 3 + e t 1 α 1 t 0 2 α 4 + e t 1 α 1 t 0 3 α 5 − e t 0 α 1 t 1 ( α 3 + t 1 ( α 4 + t 1 α 5 ) ) ) ) , ( − 1 + e 2 ( − t 0 + t 1 ) α 1 ) σ 2 2 α 1 ) . ( 9 )

Comparing Equation (8) with ( 9) demonstrates how compact the differential form of the model is, but it is the slice density which is required for model fitting and parameter estimation. The most general form of ( 8) we considered additionally allowed each of (α ₁, α ₂, α ₃, α ₄, α ₅, σ) to be included in the model as both a fixed and random (normally distributed) effect, thereby allowing trajectories to be tailored to individual children. An unstructured covariance matrix for the random effects with simplifications down to a diagonal covariance matrix were considered. The model search considered increasing orders of polynomial up to cubic and in keeping with the LMM search was guided by AIC and BIC metrics, in addition to trimming terms with high p-values (>0.1). Relevant SAS code is provided in SI.4a (see Extended data) ¹⁸ . This can also be achieved using Stan, for example, adapting Goodman (2018) ²⁸ .

The main real-world application area of SDEs is in predictive modelling (e.g. Iversen et al. ²⁶ ). Here we use SDE models to predict future growth given a child’s current age and current ZWfL. Such predictions have two application areas: 1) as part of an individual child’s risk assessment to determine if they require an intervention; and 2) to elucidate structure hidden within the large variability across growth trajectories, which may then offer new insights into drivers of growth at the population level. We use our best fitting SDE model to predict future growth trajectories across a grid of starting points for age and ZWfL, separately for each of the two sites.

Predictions are calculated using a 10-fold random sampling approach, where we draw from all the parameters estimated (fixed and random) in the best fitting SDE model – one set for each child. Which parameter sets are chosen to generate predictions depends on how likely trajectories generated from each set are to have visited each given starting point across an age-ZWfL grid. This adds an important element of “locality” to our predictions, combined with 10-fold sampling to provide an indication of robustness of our predictions. A more detailed description of the prediction algorithm is given below, with full R code provided in SI.5 (see Extended data) ¹⁸ .

For each point across an age-ZWfL grid we compute the likelihood of observing this point for each set of parameters, using the relevant slice density, where the initial starting point is the first age available for each trajectory. Predictions from each grid point progress through increasing ages using the new slice distribution at each next point in time. For example, for predictions in Venda we have n=236 likelihood values for each age-ZWfL grid point. The most likely parameter set is then used for the prediction of the next ZWfL at the next age, with 10-fold sampling used to indicate how robust this prediction is. The 10-fold sampling splits these n=236 parameter sets into 10 random groups, and within each group we choose the parameter set with the highest likelihood value as the one to be used for the next prediction. In summary, from a fixed starting point in an age-ZWfL grid we have a main prediction at each future age up to 24 months, plus 10 additional predictions at each age as a sensitivity analysis (incrementing age in small steps).

Using individual child trajectory data over 0–24 months for ZWfL from all n=460 children from the two sites, the best fitting LMM model was a cubic polynomial with a single interaction term between age squared and site (with no separate term for site), which gave AIC=21913 and BIC=21983. The best fitting SDE model had corresponding values of AIC=21718 and BIC=21772, where this model had random effects in four of the six model parameters, a diagonal covariance matrix and site included in two of the drift parameters and the diffusion parameter. Full modelling descriptions can be found in SI.4a and results, including parameter estimates, can be found in the SI4.b (see Extended data) ¹⁸ . The AIC and BIC metrics suggest the SDE offers an improved fit to the data. Examining residuals and fitted values the fit of each model is visually similar; however, there are notable qualitative and quantitative differences.

Figure 3 shows observed monthly means compared with estimates of the population mean ZWfL from each model. These are quantitatively similar except that the SDE correctly captures the initial shape of the mean, a short rapid increase then decrease during the first six months in the Venda data, whereas the LMM estimates a steady decrease from birth through six months. Figure 3 also shows confidence intervals for the population mean of ZWfL from each model over age; the SDE model has considerably narrower intervals suggesting it has explained more of the variation in the data than the LMM model.

In summary, our results so far suggest that our SDE model is at least as good, and appears superior in some respects, to a reasonable choice of classical LMM.

Figure 4 shows predictions of the most likely future ZWfL trajectory for a child, conditional on current age and ZWfL values. A grid of starting values for age and ZWfL within the ranges observed (from both sites) was used, with separate predictions for each site. There are several particularly striking features in these predictions; most notably, they suggest the presence of longitudinal clustering within each site, as we find what appear to be a small distinct set of “paths” or “streams” hidden inside the large between-child variability in trajectories observed in Figure 1. These clusters also appear to differ quantitatively and qualitatively between sites, with fewer clusters in Haydom and of different shapes to those in Venda. Further, these clusters imply some degree of canalisation of trajectories, for example, the small number of children who are wasted with ZWfL values of ≤ -2 (Venda, n = 16/236; Haydom, n = 2/224) converge on similarly low predicted values (the rapid increase from extreme values [e.g. ≤ -3] shows that predictions return to values where more data are present).

Data from the MAL-ED study are available from https://clinepidb.org/. Guest users can view data and access analysis tools and record pages, but must obtain approval from the data providers to download data. The request may be submitted via a form that pops up when a user logs in with a registered account and clicks “Download data”.

Zenodo: Introducing a drift and diffusion framework for childhood growth research. https://doi.org/10.5281/zenodo.4278002 ¹⁸ .

This project contains the following extended data: -

SI1.pdf (Model code for OU and linear regression models, in SAS)

Data are available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0).

Source code available from: https://github.com/fraseriainlewis/gatesopenresearch_SI

Archived code at time of publication: https://doi.org/10.5281/zenodo.4278002 ¹⁸

License: Creative Commons Attribution 4.0 International license (CC-BY 4.0).