# Joint Bayesian weight and height postnatal growth model to study the effects of maternal smoking during pregnancy

Sophie Carles, Marie-Aline Charles, Barbara Heude, Ismaël Ahmed, Jérémie Botton

► To cite this version:

Sophie Carles, Marie-Aline Charles, Barbara Heude, Ismaël Ahmed, Jérémie Botton. Joint Bayesian weight and height postnatal growth model to study the effects of maternal smoking during pregnancy. Statistics in Medicine, 2017, 36 (25), pp.3990-4006. 10.1002/sim.7407. hal-03382825

HAL Id: hal-03382825

https://hal.inrae.fr/hal-03382825

Submitted on 26 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Joint Bayesian weight and height postnatal growth model to study the effects of maternal smoking during pregnancy

Sophie Carles<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID | Marie-Aline Charles<fnref n="1" />,<fnref n="2" /> ID | Barbara Heude<fnref n="1" />,<fnref n="2" /> ID | Ismaël Ahmed<fnref n="4" /> ID |
Jérémie Botton<fnref n="1" />,<fnref n="5" /> ID | the EDEN mother-child study group

**Correspondence**

Sophie Carles, Early origins of the child's health and development Team (ORCHAD), UMR1153 Epidemiology and Biostatistics Sorbonne Paris Cité Center (CRESS), INSERM, Paris, France; Paris Descartes University, Paris, France; or Paris-Sud University, Paris, France.
Email: sophie2.carles@gmail.com

Growth models used for describing the dynamics of body weight and height generally consider each trait independently. We proposed modeling height and weight trajectories jointly with a nonlinear heteroscedastic mixed model based on the Jenss-Bayley growth function with correlated individual random effects and using Bayesian inference techniques. Simulations showed that our model provides good estimates of the growth parameters. We illustrated how it can be used to assess the associations between maternal smoking during pregnancy, an early-life factor potentially involved in prenatal programming of obesity, and children's growth from birth to 5 years of age. We used real data from the EDEN study, a large French mother-child cohort study with a high number of height and weight measurements (a total of approximately 30 000 measurements for each of the 2 traits across the 1666 children). Our results supported the existence of a relationship between maternal smoking during pregnancy and growth from birth to 5 years of age. Children from mothers who smoked throughout pregnancy were shown to display a higher body mass index from the first few months of life onwards compared to children from nonsmokers. At 5 years of age, their mean body mass index was 0.21 kg/m<fnref n="2" /> higher than unexposed children. It was mainly explained by the fact that these children tended to be smaller at birth but rapidly exceeded the weight of children from nonsmokers postnatally.

**KEYWORDS**

Bayesian modeling, joint model, postnatal growth trajectories, multivariate, nonlinear modeling

## 1 | INTRODUCTION

Growth of children is used as an indicator of health and development and is influenced by genetic and environmental factors. In epidemiology, growth is studied by different approaches that generally aim to identify critical periods when growth is associated with later health or diseases, or as an outcome to study its determinants. Among them, structural modeling is attractive to provide parameters with meaningful information on key growth features.<fnref n="1" /> Most growth models are generally implemented to assess patterns of height and weight growth separately. However, there is an interest in jointly studying their evolution given the high correlation between these 2 traits. Indeed, multivariate models have proved their ability (1) to study the relationship between longitudinal series over time and (2) to enhance the prediction of a longitudinal outcome when no or few past observations are available by incorporating the information from profiles of other series.<fnref n="2" />,<fnref n="3" /> Appropriately modeling children's growth during the first few years of life creates statistical challenges

<footnote n="1">¹ Early origins of the child's health and development Team (ORCHAD), UMR1153 Epidemiology and Biostatistics Sorbonne Paris Cité Center (CRESS), INSERM, Paris, France</footnote>

<footnote n="2">² Paris Descartes University, Paris, France</footnote>

<footnote n="3">³ Paris-Sud University, Paris, France</footnote>

<footnote n="4">⁴ Biostatistics, Biomathematics, Pharmacoepidemiology and Infectious Diseases (B2PHI), INSERM, UVSQ, Institut Pasteur, Université Paris-Saclay, Villejuif, France</footnote>

<footnote n="5">⁵ Faculty of Pharmacy, Paris-Sud University, Paris-Saclay University, Châtenay-Malabry, France</footnote>
due to the complexity of growth patterns. Repeated collections of height and weight measures are needed to study growth and must account for the nonindependence of measures of a same individual across time. Mixed-effects models provide a flexible and powerful tool to analyze these data.⁴

Growth patterns have already been fitted using nonlinear models, based on the structural Jenss-Bayley equation,⁵ using mixed-effect models,⁶⁻⁸ that adequately fit growth data from birth to 8 years of age and describe individual growth trajectories using 4 subject-specific parameters with clinical interpretations.¹ The likelihood of nonlinear mixed models (NLMMs) does not have a closed-form expression because of the integration of random effects, making inference about parameters somewhat challenging.⁹ Stochastic versions of the Expectation Maximization (EM) algorithms as Monte Carlo EM have been proposed.¹⁰ The “stochastic approximation EM algorithm” (SAEM) is one that is more efficient computationally than Monte Carlo EM because of smoothing on the E step based on recycling simulations from between iterations. This algorithm is now implemented in R¹¹ in the package SAEMIX.¹²

Bayesian approaches to joint NLMMs have been proposed especially for HIV dynamics and have shown to be a powerful way to analyze such complex data. Indeed, specifying prior distributions for the model parameters allows users to gain insight from the results of previous studies but also to handle some parameter identifiability problems. The Bayesian approach also offers a flexible way to fit complex models, especially hierarchical or multilevel ones.¹³ One major interest, inherent to Bayesian reasoning, is that the full posterior distribution of parameters is available. Thus, this approach is suitable for predictions that account for uncertainty in the parameters.¹⁴

Herein, we propose a joint Bayesian model of height and weight growth and illustrate how it can be used to assess the associations between maternal smoking during pregnancy, an early-life factor potentially involved in prenatal programming of obesity, and children's growth from birth to 5 years of age. We used real data from the EDEN study, a large French mother-child cohort study with a high number of height and weight measurements (a total of approximately 30 000 measurements for each of the 2 traits across the 1666 children). We hypothesized that the Bayesian approach would offer computational facilities especially when the number of factors to consider for their association with the growth trajectory is high. Our approach should allow us to simultaneously assess the associations between the determinant of interest and growth in weight and height and infer the association with body mass index (BMI).

In Section 2, we introduce the bivariate model based on the Jenss-Bayley growth function. In Section 3, we apply this model to simulated data to assess its ability to adequately fit growth data by providing accurate parameter estimates. In Section 4, we apply the joint model on growth data resulting from the follow-up of children from the EDEN mother-child cohort study. In that section, we compare the separate modeling of height and weight to the joint modeling, and we apply the joint model to assess the adjusted association between maternal smoking and height, weight and BMI growth. Section 5 provides some elements of discussion.

# 2 | MODEL FOR POSTNATAL GROWTH DURING CHILDHOOD

## 2.1 | A general 3-stage hierarchical model

The joint model consists of the following 2 nonlinear mixed-effects submodels: one for height and the other for weight. Compared to the univariate model, the joint model allows (1) considering height and weight longitudinal patterns as nonindependent profiles, and (2) capturing and characterizing the relationship between them and predicting BMI in different groups of interest. It should allow us to identify whether children in subgroups of interest have atypical weight trajectories considering their height trajectory and describe the repercussions on the BMI. In addition, it should allow us to improve the predictive ability for growth trajectories. Notably, when there are no or few observations for a growth profile, incorporating the information available for the other one should improve the prediction of its theoretical pattern. It could be a way to correct the bias linked to the exclusion of informative data and is a key advantage in cohort studies where there are loss to follow-up children. The model can be displayed as a general 3-stage hierarchical Bayesian model along the same line as in Davidian and Giltinan.¹⁵

### 2.1.1 | Stage 1: intraindividual variability

We start by writing the NLMM under its classical formulation, as follows:

$$y_{ijk} = f(\phi_{ik}, t_{ijk}) + e_{ijk}, \qquad (1)$$
$k=(H, W)$ denotes the 2 outcomes of interest, ie, height and weight. $y_{ik} = (y_{1ik}, ..., y_{in_{ik}})^T$ denotes a $n_{ik}$ vector of observed measures of height or weight for individual $i$ at time $t_{ijk}$, $i = 1,..., N$, $j = 1,..., n_{ik}$. $f$ is a nonlinear function of a $m$ vector of individual specific parameters $\phi_{ik}$. $e_{ik} = (e_{1ik}, ..., e_{in_{ik}})^T$ is a $n_{ik}$ vector of residual errors of height or weight that are assumed to be independently and normally distributed with a null mean and variance $\sigma_{e_k}^2$ so that

$$e_{ik} = (e_{ijk}) \sim N(0, \sigma_{e_k}^2 I_{n_{ik}}), \qquad (2a)$$

where $I_{n_{ik}}$ is the identity matrix.

To describe residual variances of height and weight, we used general variance functions involving residual variations with age to account for a potential heteroscedasticity, as described by Foulley et al and Duval et al$^{16,17}$:

$$\log(\sigma_{e_{ijk}}^2) = \delta_{0k} + \delta_{1k} t_{ijk} + \delta_{2k} t_{ijk}^2.$$

Applied to height and weight models separately, we selected, among the 3 underlying submodels described below, the one that minimized the fit parameters assessed by the Watanabe Akaike Information Criterion (WAIC) (see supplementary material I and II for details on WAIC, method, and values)

M1: $\log(\sigma_{e_{ijk}}^2) = \delta_{0k}$ (where $\delta_{1k}=\delta_{2k}=0$), ie, homogeneous residual variances

M2: $\log(\sigma_{e_{ijk}}^2) = \delta_{0k} + \delta_{1k}t_{ijk}$ (where $\delta_{2k}=0$), ie, linear adjustment of logvariances according to age

M3: $\log(\sigma_{e_{ijk}}^2) = \delta_{0k} + \delta_{1k}t_{ijk} + \delta_{2k}t_{ijk}^2$, ie, quadratic adjustment of logvariances according to age

On our data, the best improvement of the fit parameters was obtained using the model M2 for height and M3 for weight leading to the following equations:

Height: $\log(\sigma_{e_{ijH}}^2) = \delta_{0H} + \delta_{1H}t_{ijH}$;

Weight: $\log(\sigma_{e_{ijW}}^2) = \delta_{0W} + \delta_{1W}t_{ijW} + \delta_{2W}t_{ijW}^2$.

Formulas 1 and 2a are equivalent in the Bayesian hierarchical models to writing the conditional distribution of the data given the individual parameters $\phi_{ik}$ as

$$y_{ik} | \phi_{ik} \sim N(f_{ik}, \sigma_{e_k}^2 I_{n_{ik}}), \qquad (2b)$$

where $y_{ik}=(y_{ijk})$ for $j = 1,..., n_{ik}$ is a $n_{ik}$ vector and similarly $f_{ik}$ is defined as $f_{ik} = (f(\phi_{ik}, t_{ijk}))$ for $j = 1,..., n_{ik}$.

For the weight and height modeling, $f$ is a modification of the initial nonlinear Jenss-Bayley equation that improves the interpretation of parameters and accounts for mathematical constraints.$^8$ This growth function allows the expression of expectations of height in centimeters (cm) and weight in kilograms (kg) for the individual $i$ at time $t_{ijk}$, as follows:

$$f(\phi_{ik}, t_{ijk}) = \exp(A_{ik}) + \exp(B_{ik})t_{ijk} + \exp(C_{ik})(1 - \exp(-\exp(D_{ik})t_{ijk})) \qquad (3)$$

with $k=(H, W)$ and $\phi_{ik}=(A_{ik}, B_{ik}, C_{ik}, D_{ik})^T$.

## 2.1.2 | Stage 2: interindividual variability

Individual effects are treated as random with the following multinormal distribution:

$$\phi_i \sim N(\mu_i, \Omega), \qquad (4)$$

where $\phi_i=(\phi_{iH} \phi_{iW})^T$ with height $\phi_{iH}$ and weight $\phi_{iW}$ vector components having mean vector $\mu_i = (\mu_{iH} \mu_{iW})^T = (\mu_{iH}^A, ..., \mu_{iH}^D, \mu_{iW}^A, ..., \mu_{iW}^D)^T$ and variance covariance matrix $\Omega = \begin{bmatrix} \Omega_{HH} & \Omega_{HW} \\ \Omega_{WH} & \Omega_{WW} \end{bmatrix}$ of dimension 2 m × 2 m.

In Formula 3, $\exp(A_{iH})$ describes the length at birth of subject $i$. For weight modeling, to bypass the neonatal weight loss not correctly fitted by the model, birth weight measures were excluded and the minimum weight in the first 4 days of
life was selected for each child. Consequently, in the weight model of subject *i*, $\exp(A_{iW})$ represents the extrapolated weight at time 0 without accounting for neonatal weight loss; for any trait *k*, $\exp(B_{ik})$ is the growth velocity of height or weight beyond 2 years of age; $\exp(C_{ik})$ is the spurt of height/weight growth in the first 2 years of life, and $\exp(D_{ik})$ is the trajectory curvature in early life (combination of growth velocity and deceleration rate in the beginning of life). Means $\mu_{ik}$ of the distribution shown in (4) are related to a matrix of explanatory variables $X_i$ of dimension $m \times p$ (intercept and discrete or continuous covariates) via a classical linear model involving a $p$ dimensional vector of parameters $\beta_k$.

$$ \mu_{ik} = X_i \beta_k \qquad (5) $$

In models without covariates, the vector of “fixed” effects $\beta_k$ is restricted to intercept elements $(\beta_{0Ak}, \beta_{0Bk}, \beta_{0Ck}, \beta_{0Dk})^T$ corresponding to the population means of growth parameters. In models with covariates on growth parameters (cf Section 4), $\beta_k$ includes both a vector of $m$ parameters intercepts $\beta_{0Ak}, \beta_{0Bk}, \beta_{0Ck}, \beta_{0Dk}$ and a vector of $mp_{cov}$ regression coefficients $\beta_{1Ak} \cdots \beta_{p_{covAk}}, \cdots, \beta_{1Dk} \cdots \beta_{p_{covDk}}$ pertaining to the $p_{cov}$ covariates of interest.

### 2.1.3 | Stage 3: hyperprior distributions

We specified standard prior distributions usually used for Bayesian hierarchical models¹⁴,¹⁸ and parameter values that take advantage of published studies in the field of growth, including ours.⁸

The prior distributions are described in Table 1 (in rjags notations). Briefly, for the weight submodel, a mean value of 1 for the $A$ parameter specifies the prior expectation of the extrapolation of the average minimal weight after birth to be 2.7 kg and a precision (inverse of variance) of 16 (95% probability in the range [1.70, 4.40]). We defined a mean $B$ parameter value of -5, a mean $C$ parameter of 1.5, and a mean $D$ parameter of -5 with a precision of 16. The same approach was applied for the height submodel. The priors for the height model parameters (see Table 1 for prior distribution) were adapted to provide plausible birth lengths and height trajectories.

In the presence of covariates, to keep the values generating realistic growth trajectories, we assumed the prior distribution for all the coefficients of the $\beta_k$ vector except the intercept terms to be normally distributed with a mean zero and a precision of 16, (ie, a variance of 0.0625 or a standard deviation of 0.25), ie, in the range [-0.49, 0.49]) with a 95% probability for $k = (H, W)$.

For the simulation study, prior normal distributions with a mean of zero and precision of 0.1 were specified on the parameters of the residual variance function $\delta_{0H}$ and $\delta_{1H}$ for height and $\delta_{0W}$, $\delta_{1W}$ and $\delta_{2W}$ for weight. For the application, we used precision parameter of 16. The precision matrix $\Gamma = \Omega^{-1}$ of the 8 individual parameters ($A, B, C, D$ for height and

<table><caption>TABLE 1 Prior distributions used for hyperparameters involved in the Bayesian model (rjags notation)</caption><thead><tr><th>Parameters</th><th>Weight (k = W)</th><th>Height (k = H)</th></tr></thead><tbody><tr><td>β<sub>0Ak</sub></td><td>dnorm(1; 16)</td><td>dnorm(4; 16)</td></tr><tr><td>β<sub>0Bk</sub></td><td>dnorm(-5; 16)</td><td>dnorm(-4; 16)</td></tr><tr><td>β<sub>0Ck</sub></td><td>dnorm(1.5; 16)</td><td>dnorm(3; 16)</td></tr><tr><td>β<sub>0Dk</sub></td><td>dnorm(-5; 16)</td><td>dnorm(-5; 16)</td></tr><tr><td>β<sub>1(A-D)k</sub> to β<sub>p<sub>cov(A-D)k</sub></sub></td><td>dnorm(0; 16)</td><td>dnorm(0; 16)</td></tr><tr><td>Ω<sup>-1</sup></td><td colspan="2">dwish(R<sup>-1</sup>, 9)</td></tr><tr><td>Residual variance</td><td colspan="2"></td></tr><tr><td>Simulation study</td><td colspan="2"></td></tr><tr><td>δ<sub>0k</sub></td><td>dnorm(0; 0.1)</td><td>dnorm(0; 0.1)</td></tr><tr><td>δ<sub>1k</sub></td><td>dnorm(0; 0.1)</td><td>dnorm(0; 0.1)</td></tr><tr><td>δ<sub>2W</sub></td><td>dnorm(0; 0.1)</td><td>NA</td></tr><tr><td>Application</td><td colspan="2"></td></tr><tr><td>δ<sub>0k</sub></td><td>dnorm(0; 16)</td><td>dnorm(0; 16)</td></tr><tr><td>δ<sub>1k</sub></td><td>dnorm(0; 16)</td><td>dnorm(0; 16)</td></tr><tr><td>δ<sub>2W</sub></td><td>dnorm(0; 16)</td><td>NA</td></tr></tbody></table>
weight) was supposed to follow a Wishart distribution $W(R, v)$, a convenient distribution for positive-definite matrices due to conjugacy properties. $\Gamma$ is thus a priori centered on $vR$ and has $v$ degrees of freedom.¹⁹ Here, $R = \Gamma_0/v$, where $\Gamma_0 = \Omega_0^{-1}$ is the prior mean of the precision matrix. Setting the prior on the precision $\Gamma$ is equivalent to assuming that the covariance $\Omega$ has an inverse Wishart distribution $IW(R^{-1}, v)$ with $E(\Omega) = R^{-1}/(v - 2m - 1)$. Notably, in Winbugs and rjags, the Wishart distribution is specified by $R^{-1} = v\Omega_0$ instead of $R$. The parameter $v$ should be equal or higher than the dimension of the precision matrix for the prior being proper; low values indicate less informative prior information.²⁰ Here, $v$ was set equal to the number of parameters (2 traits with m components each) plus one ie, $v=2m+1$ as advocated by Gelman and Hill.²¹ Therefore, we used, according to rjags notations, $\Gamma = \Omega^{-1} \sim W(R^{-1}, 9)$ with $R^{-1}$ chosen from previous estimations of the variance-covariance between the height and weight model parameters⁷ (see supplementary material III for detailed values).

## 2.2 | Bayesian implementation

### 2.2.1 | Algorithm and software

The full conditional posterior distribution of all model parameters was calculated using Markov Chain Monte Carlo (MCMC) methods within the rjags package (version 3-15) in R software (version 3.1.0).²² In this package, MCMC methods are simulated using the Gibbs sampler. Within this iterative process, after choosing arbitrary starting values, each parameter of the model is sampled for each iteration t from the conditional distribution given all the others parameters at iteration t-1. In other words, for each iteration s, each component $\theta_j$ of the whole vector of parameters $\theta$ will be sampled conditional on the most recent values of the other components of $\theta$.²⁰,²³

### 2.2.2 | Computational aspects

Models were run for 60 000 iterations (with 20 000 burn-in iterations not used for posterior inference) using 3 chains initialized with different starting values. Chains were thinned by selecting every 20th iteration such that the autocorrelation was sufficiently low. Inferences were made on a total of 1500 samples (500 samples per chain).

### 2.2.3 | Convergence assessment

Convergence assessments were based on visual inspections of the simulation trace plot for each parameter and on the Brooks-Gelman-Rubin criteria (ie, criterion lower than 1.05).²³⁻²⁵ This criterion compares between and within chain variability using the ratio of the total to the within chain variations. A criterion close to 1 indicates a good convergence.

## 2.3 | Model selection and evaluation

In the application (see Section 4), to test the effects of smoking covariates on height and weight trajectories, models with and without those covariates were compared using the Deviance Information Criterion (DIC)²⁶ and the WAIC²⁷ (see supplementary material I for details on both criterions).

To validate the selected model, posterior predictive checks derived by Gelman et al were performed.²³ At each simulation (retained iteration) ($s=1, ..., S$), simulated replicated data $y^{\text{rep } s}$ were generated under the model from the predictive distribution of $\theta^s$ at the corresponding simulation. At each simulation, using a Chi-squared discrepancy measure, we assessed, for height and weight separately, the test quantity for replicated datasets $T(y_k^{\text{rep } s}, \theta^s)$ and for observed data $T(y_k, \theta^s)$. We calculated the posterior predictive $P$ value for each trait $k$ ($p_{B_k}$) that corresponds to the proportion of simulations $S$ for, which

$$T(y_k^{\text{rep } s}, \theta^s) \geq T(y_k, \theta^s),$$

$$\text{where } T(y_k, \theta^s) = \sum_{i=1}^{N} \sum_{j=1}^{n} \frac{(y_{ijk} - E(y_{ijk} | \theta^s))^2}{\text{Var}(y_{ijk} | \theta^s)} \quad \text{and} \quad T(y_k^{\text{rep } s}, \theta^s) = \sum_{i=1}^{N} \sum_{j=1}^{n} \frac{(y_{ijk}^{\text{rep } s} - E(y_{ijk} | \theta^s))^2}{\text{Var}(y_{ijk} | \theta^s)}.$$

To assess the sensitivity of posterior inferences to prior distribution, a sensitivity analysis using noninformative priors was performed (see supplementary material IV).
# 3 | SIMULATION STUDY

## 3.1 | Simulated growth data

We performed a simulation study to assess and illustrate the ability of the model (see supplementary material III for model syntax) to adequately estimate the key growth features with respect to classical properties of estimators (bias, standard, and mean square error). Datasets were simulated mimicking height and weight growth of typical children. A total of 200 datasets were generated each containing the anthropometric data of 500 children with 10 measurements at time points ranging from day 5 to day 1825. Specifically, data were generated at days 5, 30.44 (1 month), 91.32 (3 months), 182.64 (6 months), 273.96 (9 months), 365.25 (1 year), 730.5 (2 years), 1095.75 (3 years), 1461 (4 years), and 1826.25 (5 years). Population parameters values were guided by the results obtained from a previous study through a frequentist SAEM algorithm.⁷

1. Each individual parameter was sampled from a multinormal distribution with the mean parameter values as described in Table 2 (see true values). The 8-dimensional variance covariance matrix between the individual height and weight parameters was empirically chosen according to previous research⁷ (see values of $\Omega_0$ in supplementary material III).

2. Heights and weights were then predicted at each time point (days 5-1825) for each child in the simulated dataset by applying the Jenssen-Bayley equation. We considered these values as the true height and weight measurements free from measurement error.

<table><caption>**TABLE 2** True values, bias, and probability coverage obtained from the Bayesian joint model of height and weight on 200 simulated datasets</caption><thead><tr><th>Parameter</th><th>True Value</th><th>Bias</th><th>Relative Bias</th><th>95% Coverage Probability</th><th>RMSE</th></tr></thead><tbody><tr><td>A<sub>W</sub></td><td>1.00000</td><td>-0.00063</td><td>-0.00063</td><td>95.5</td><td>0.00892</td></tr><tr><td>B<sub>W</sub></td><td>-5.00000</td><td>0.00028</td><td>0.00006</td><td>93.5</td><td>0.00395</td></tr><tr><td>C<sub>W</sub></td><td>1.50000</td><td>-0.00005</td><td>-0.00003</td><td>94.5</td><td>0.00069</td></tr><tr><td>D<sub>W</sub></td><td>-5.00000</td><td>0.00089</td><td>0.00018</td><td>92.5</td><td>0.01261</td></tr><tr><td>A<sub>H</sub></td><td>4.00000</td><td>-0.00002</td><td>-0.00001</td><td>95.5</td><td>0.00030</td></tr><tr><td>B<sub>H</sub></td><td>-4.00000</td><td>-0.00023</td><td>-0.00006</td><td>95.5</td><td>0.00338</td></tr><tr><td>C<sub>H</sub></td><td>3.00000</td><td>-0.00005</td><td>-0.00002</td><td>95.0</td><td>0.00081</td></tr><tr><td>D<sub>H</sub></td><td>-5.00000</td><td>-0.00036</td><td>-0.00007</td><td>95.5</td><td>0.00504</td></tr><tr><td>δ<sub>0W</sub></td><td>-3.00000000</td><td>0.01541804</td><td>-0.00513935</td><td>94.0</td><td>0.21804405</td></tr><tr><td>δ<sub>1W</sub></td><td>0.00200000</td><td>-0.00003799</td><td>-0.01899257</td><td>95.5</td><td>0.00053719</td></tr><tr><td>δ<sub>2W</sub></td><td>-0.00000050</td><td>0.00000002</td><td>-0.04195978</td><td>96.0</td><td>0.00000030</td></tr><tr><td>δ<sub>0H</sub></td><td>0.40000000</td><td>0.00312113</td><td>0.00780282</td><td>94.5</td><td>0.04413940</td></tr><tr><td>δ<sub>1H</sub></td><td>0.00040000</td><td>0.00000213</td><td>0.00533664</td><td>97.5</td><td>0.00003019</td></tr><tr><td>A<sub>W</sub>A<sub>W</sub></td><td>0.01444</td><td>0.00011</td><td>0.00795</td><td>95.5</td><td>0.00162</td></tr><tr><td>B<sub>W</sub>B<sub>W</sub></td><td>0.04000</td><td>0.00079</td><td>0.01985</td><td>95.5</td><td>0.01123</td></tr><tr><td>C<sub>W</sub>C<sub>W</sub></td><td>0.05444</td><td>0.00082</td><td>0.01508</td><td>94.0</td><td>0.01161</td></tr><tr><td>D<sub>W</sub>D<sub>W</sub></td><td>0.09556</td><td>0.00105</td><td>0.01094</td><td>95.5</td><td>0.01478</td></tr><tr><td>A<sub>H</sub>A<sub>H</sub></td><td>0.00111</td><td>0.00001</td><td>0.00983</td><td>94.0</td><td>0.00015</td></tr><tr><td>B<sub>H</sub>B<sub>H</sub></td><td>0.01111</td><td>0.00032</td><td>0.02925</td><td>94.0</td><td>0.00460</td></tr><tr><td>C<sub>H</sub>C<sub>H</sub></td><td>0.01778</td><td>0.00044</td><td>0.02455</td><td>95.5</td><td>0.00617</td></tr><tr><td>D<sub>H</sub>D<sub>H</sub></td><td>0.06778</td><td>0.00152</td><td>0.02245</td><td>97.0</td><td>0.02152</td></tr></tbody></table>

$(A - D)_k(A - D)_k$ refers to the variance of the 4 parameters describing the trait $k$ corresponding to the diagonal values of the variance-covariance matrix $\Omega_{kk}$, where $k = H$ for height and $k = W$ for weight.
3. Consequently, we introduced some measurement error to generate data that would be closer to those reported in practice. For each individual in the simulated dataset, the observed data at given time points were obtained by sampling one value from a normal distribution. The mean was the predicted value (2) and the standard error varied with age according to the selected residual variance function (see Section 2.1). To generate standard error values for height ($\sigma_{e_{ijH}}$) at each time, the parameters $\delta_{0H}$ and $\delta_{1H}$ describing the residual variances were fixed to 0.4 and 0.0004, respectively. To calculate standard error for weight ($\sigma_{e_{ijW}}$), $\delta_{0W}$, $\delta_{1W}$, and $\delta_{2W}$ were fixed to -3, 0.002, and $-5 \times 10^{-7}$, respectively. These values correspond to a residual standard deviation increasing from 5 days to 5 years from 1.22 to 1.76 cm for height and from 224 to 602 g for weight.

We ran the Bayesian model on the 200 simulated datasets. We fixed a burn-in period of 10 000 iterations for 3 chains, ran the chains for an additional 20 000 iterations, and collected samples at intervals of 20 iterations.

The bias (ie, difference between the true value of the parameters (step 1) and the mean estimate obtained), the relative bias (ie, ratio of the bias over the true value), the RMSE (ie, root mean square error), and the probability coverage (ie, the proportion of 95% credible interval containing the true value) were additionally computed. The expectation of the posterior distribution was used to compute the estimates for each parameter.

## 3.2 | Simulation results

Estimates obtained for the 8 main parameters (A, B, C, and D for height and weight parameters) were quite satisfactory from a frequentist point of view. Absolute and relative biases were closed to 0, and the coverage probabilities ranged from 92.5 to 95.5% (Table 2 summarizes the results for the main parameters and their variance; results for the covariance between parameters are described in supplementary material V). Most of the other parameters (ie, variances and covariances) were also adequately estimated by the model.

Figure 1 illustrates the simulated and fitted (mean trajectories and 95% credible interval of the trajectories) weight growth trajectories for 3 randomly selected children and provides a visual validation of the model adequacy. Height growth trajectories were also well fitted by the model (see supplementary material VI).

<figure><img src="image_2.png" /><figcaption>**FIGURE 1** Weight measures (simulated data) of randomly selected individuals and datasets (open circles), mean curves fitted by the joint Bayesian model and 95% credible intervals</figcaption></figure>
# 4 | APPLICATION TO GROWTH DATA

## 4.1 | The EDEN study dataset

We applied the joint Bayesian model to the data from the EDEN study, a birth cohort designed to assess pre and early postnatal determinants of child health and development.²⁸ Within the frame of the study, between 2003 and 2006, 2002 pregnant women were recruited in the maternity wards of 2 university hospitals in the French cities of Nancy and Poitiers. Women were recruited prior to 24 weeks of gestation defined by the date of the last menstrual period. The child's birth weight, daily weights on days 1 to 5, and birth length were extracted from the obstetric records. Length or height and weight postnatal measurements were collected from health booklet via mailed questionnaires and during the study's clinical visits at 1, 3, and 5 years.

## 4.2 | Comparison between univariate and multivariate models

Our aim was to assess the better predictive ability of the joint model compared to models of height and weight separately. We worked on a subsample of children from the EDEN study. More precisely, we considered boys only with at least 2 weight and/or 2 height measures: 942 boys had at least 2 weight measures (16 751 measures as a whole) and 917 had at least 2 height measures (15 567 measures as a whole). We performed 2 types of evaluations:

First, we performed a comparison of the intraindividual variability of the height and weight predictions obtained with the joint model of height and weight and with separate height and weight models. All these models were fitted within the Bayesian framework. Three MCMC chains were launched for 20 000 iterations to compute the posterior distributions (after a burning phase of 20 000 iterations). To limit autocorrelation, each chain was thinned by selecting every 20th iteration. Individual height, weight, and BMI and the corresponding instantaneous velocities (details on calculations methods are given in previous studies⁷,⁸) were predicted at 3 months, 1, 2, 3, 4, and 5 years. From the posterior distributions of the individual predictions, we get the mean intraindividual variability and computed the mean over all subjects (Table 3).

The intraindividual variability was systematically lower with the joint model for every growth characteristic considered in this work. In addition, we could compute directly the predictions for BMI.

Secondly, in the presence of missing values, we studied whether the precision of the predictions obtained for the height growth was better when adding the available information on the weight growth and reciprocally. We randomly selected 20% of the available height measures from our sample that we considered as missing values. It represented 3113 measures among the 15 567 available ones. All weight measures were kept. The joint model of height and weight was used to predict the removed height measures. A height-only growth modeling was used to predict the same information. We assessed the differences between the observed measures of height that were removed and the mean predictions obtained with the joint model on the one hand and with the univariate model on the other hand (Table 4).

The results show that when information on height is lacking at one age, including the information on weight slightly improved the prediction since the variance of the difference was lower (variance = 2.10) compared to including only the other height measures of the individual (variance = 2.36).

This joint modeling seems to have better predictive abilities than the height and weight models separately, representing an important strength of this model.

## 4.3 | Application of the Bayesian joint model to the analysis of the association between maternal smoking during pregnancy and growth

Our aim was to assess whether height, weight, and BMI growth patterns differ between children exposed to different degrees of maternal smoking, namely, exposure during only the first trimester of pregnancy, throughout pregnancy, and unexposed children. We analyzed weight and height records (30 515 and 28 381, respectively) of 1666 children with complete information regarding maternal smoking during pregnancy and potential confounders.

To reduce confounding, we adjusted our analyses for the following variables: maternal education level (high school or less; 2-year university degree; 3-year university degree), maternal prepregnancy BMI (overweight/obese or not), maternal age at delivery (continuous), breastfeeding duration (≥3 months or <3 months), recruitment center, and child's gender. For computational reasons, maternal age was centered on its population mean.
<table><caption>**TABLE 3** Intraindividual variability of the predictions of height, weight, BMI, and their corresponding velocities obtained for the joint and separate modeling of height and weight</caption><thead><tr><th rowspan="2">Variable</th><th colspan="3">Mean of the Intrasubject Variability</th></tr><tr><th>Joint Model of Height and Weight</th><th>Univariate Weight Model</th><th>Univariate Height Model</th></tr></thead><tbody><tr><td>Weight, kg</td><td></td><td></td><td></td></tr><tr><td>3 months</td><td>0.0365</td><td>0.0405</td><td></td></tr><tr><td>1 year</td><td>0.0557</td><td>0.0563</td><td></td></tr><tr><td>2 years</td><td>0.1089</td><td>0.1127</td><td></td></tr><tr><td>3 years</td><td>0.1891</td><td>0.1943</td><td></td></tr><tr><td>4 years</td><td>0.3389</td><td>0.3460</td><td></td></tr><tr><td>5 years</td><td>0.5736</td><td>0.5848</td><td></td></tr><tr><td>Height, cm</td><td></td><td></td><td></td></tr><tr><td>3 months</td><td>0.2297</td><td></td><td>0.2819</td></tr><tr><td>1 year</td><td>0.2795</td><td></td><td>0.3068</td></tr><tr><td>2 years</td><td>0.5767</td><td></td><td>0.6449</td></tr><tr><td>3 years</td><td>0.8184</td><td></td><td>0.8953</td></tr><tr><td>4 years</td><td>1.1940</td><td></td><td>1.2951</td></tr><tr><td>5 years</td><td>1.9782</td><td></td><td>2.1741</td></tr><tr><td>IMC, kg/m<sup>2</sup></td><td></td><td></td><td></td></tr><tr><td>3 months</td><td>0.2226</td><td></td><td></td></tr><tr><td>1 year</td><td>0.1259</td><td></td><td></td></tr><tr><td>2 years</td><td>0.1271</td><td></td><td></td></tr><tr><td>3 years</td><td>0.1310</td><td></td><td></td></tr><tr><td>4 years</td><td>0.1568</td><td></td><td></td></tr><tr><td>5 years</td><td>0.2050</td><td></td><td></td></tr><tr><td>Weight growth velocity, kg/month</td><td></td><td></td><td></td></tr><tr><td>3 months</td><td>0.0016</td><td>0.0021</td><td></td></tr><tr><td>1 year</td><td>0.0005</td><td>0.0006</td><td></td></tr><tr><td>2 years</td><td>0.0002</td><td>0.0002</td><td></td></tr><tr><td>3 years</td><td>0.0003</td><td>0.0003</td><td></td></tr><tr><td>4 years</td><td>0.0003</td><td>0.0003</td><td></td></tr><tr><td>5 years</td><td>0.0003</td><td>0.0003</td><td></td></tr><tr><td>Height growth velocity, cm/month</td><td></td><td></td><td></td></tr><tr><td>3 months</td><td>0.0199</td><td></td><td>0.0243</td></tr><tr><td>1 year</td><td>0.0045</td><td></td><td>0.0053</td></tr><tr><td>2 years</td><td>0.0012</td><td></td><td>0.0013</td></tr><tr><td>3 years</td><td>0.0012</td><td></td><td>0.0013</td></tr><tr><td>4 years</td><td>0.0016</td><td></td><td>0.0019</td></tr><tr><td>5 years</td><td>0.0018</td><td></td><td>0.0022</td></tr><tr><td>IMC growth velocity, kg/m<sup>2</sup>/month</td><td></td><td></td><td></td></tr><tr><td>3 months</td><td>0.0090</td><td></td><td></td></tr><tr><td>1 year</td><td>0.0014</td><td></td><td></td></tr><tr><td>2 years</td><td>0.0003</td><td></td><td></td></tr></tbody></table>

(Continues)
<table><caption>**TABLE 3** (Continued)</caption><thead><tr><th rowspan="2">Variable</th><th colspan="3">Mean of the Intrasubject Variability</th></tr><tr><th>Joint Model of Height and Weight</th><th>Univariate Weight Model</th><th>Univariate Height Model</th></tr></thead><tbody><tr><td>3 years</td><td>0.0002</td><td></td><td></td></tr><tr><td>4 years</td><td>0.0001</td><td></td><td></td></tr><tr><td>5 years</td><td>0.0001</td><td></td><td></td></tr></tbody></table>

<table><caption>**TABLE 4** Distribution of the differences (observed-predicted values) for the subsample of 20% randomly suppressed height measures obtained with the joint modeling and univariate modeling of height</caption><thead><tr><th>Model</th><th>Minimum</th><th>First Quartile</th><th>Median</th><th>Mean</th><th>Third Quartile</th><th>Maximum</th><th>Variability</th></tr></thead><tbody><tr><td>Joint</td><td>-11.23</td><td>-0.82</td><td>0.04</td><td>-0.03</td><td>0.84</td><td>7.49</td><td>2.10</td></tr><tr><td>Univariate</td><td>-12.27</td><td>-0.89</td><td>0.01</td><td>-0.06</td><td>0.85</td><td>7.90</td><td>2.36</td></tr></tbody></table>

## 4.4 | Associations between maternal smoking and height and weight growth parameters

We obtained a DIC of 121 013 for the model without the smoking covariates. The DIC value was smaller for the model that included smoking covariates (120 997). This was also true with WAIC (121 646 versus 121 669 for the models with and without smoking covariates). Hence, the model including the smoking covariates was preferred. These differences in DIC/WAIC suggest an association between maternal smoking and growth.

Furthermore, our model with smoking covariates had a global posterior predictive *P* value of .48 for height and .60 for weight. Graphical plots of T(y<sup>rep</sup>, θ) against T(y, θ) also indicated good predictive ability of the model (joint distribution of simulation draws of T(y, θ) and its replication T(y<sup>rep</sup>, θ) are illustrated in supplementary material VII).

Figure 2 describes the collected weight observations for 6 randomly selected children with a mean growth trajectory predicted by the model (ie, mean of the posterior distribution) and its 95% credibility interval. This finding provides a visual indication that the model adequately fits weight data. Similar conclusions were made for height data (figure related to height is given in supplementary material VI).

<figure><img src="image_3.png" /><figcaption>**FIGURE 2** Weight measures (open circles), mean curves fitted by the joint Bayesian model and 95% credible intervals of randomly selected children from the EDEN study</figcaption></figure>
Table 5 describes the association between maternal smoking and the main growth parameters of height and weight. Several credibility intervals of regression coefficients associated with maternal smoking on growth parameters ($A_H$, $A_W$, $B_W$, and $D_H$) did not include zero suggesting an association between maternal smoking and growth. Women who smoked throughout pregnancy delivered children with lower birth lengths ($A_H$ parameter regression coefficient = -0.015, 95% credibility interval: -0.021, -0.008), compared with nonsmoking women. The $A_W$ parameters were also lower for these mother-child pairs. Maternal smoking throughout pregnancy was also associated with higher $B_W$ parameters in children, indicating a higher growth velocity of weight from 2 years of age onwards. Maternal smoking exclusively during the first trimester was associated with a higher $D_H$ parameter, indicating a higher curvature degree of the height growth curve in the first months of life.

More moderate associations were observed for some parameters (Figure 3). The posterior distribution of the regression coefficients for both maternal smoking exclusively during the first trimester and maternal smoking throughout pregnancy shifted toward positive values for height $B_H$. The regression coefficient for the height $D_H$ parameter also shifted toward positive values for maternal smoking throughout pregnancy. The posterior distribution of the regression coefficient for $C_H$ shifted toward negative values for maternal smoking during the first trimester exclusively. This finding provides evidence of a greater height growth velocity after 2 years of age in exposed children from both groups (for $B_H$), a greater height spurt in the first 2 years of life in children exposed in early pregnancy (for $C_H$), and a greater curvature of height in children exposed throughout pregnancy (for $D_H$).

To quantify the differences in height, weight, and BMI between the children exposed during only the first trimester or children exposed throughout pregnancy and the reference category, we simulated height, weight, and BMI at day 5, 15, 25, 30, and then at each month from 30 days to 5 years for fictitious children in each smoking group (10 000 children per group), setting fixed values for covariates (male gender, breastfeeding less than 3 months, birth in Poitiers, and having mothers with a normal BMI, a high school level or less, and an age of 30 years at delivery). The mean height, weight and BMI values at each time point were computed among subgroups of children. The differences in the mean values between children whose mothers smoked exclusively during the first trimester or throughout pregnancy and the reference category (ie, children whose mothers did not smoke) were then computed. Table 6 and Figure 4 illustrate the mean posterior estimates of these differences and their 95% credibility intervals.

Maternal smoking during only the first trimester was not associated with any difference in BMI, height, or weight at the selected ages. Children whose mothers smoked throughout pregnancy were born shorter (mean difference at 5 days: -0.71 cm, 95% credibility interval: -1.03, -0.38) and lighter. These babies had a 0.05 kg/m² lower BMI at 5 days (95% credibility interval: -0.23, 0.12). On average, those children reached the weight and height trajectories of the reference category between 4 and 6 months and approximately 4 years later, respectively. Children exposed to maternal smoking throughout pregnancy tended to weigh more than children of nonsmoking mothers at the end of the study period: At 4 and 5 years of age, they weighed 231 and 301 g more (credibility intervals: -36, 505 and -24, 639, respectively). Children whose mothers smoked throughout pregnancy tended to be smaller and heavier than those from nonsmokers over the

<table><caption>TABLE 5 Mean differencesª in height and weight model parameters between each maternal smoking category and the reference category (no maternal smoking) and 95% credibility intervals, 1666 children 0 to 5 years old from the EDEN study</caption><thead><tr><th rowspan="3">Outcomes</th><th colspan="2">A</th><th colspan="2">B</th><th colspan="2">C</th><th colspan="2">D</th></tr><tr><th colspan="2">Birth length<br>Extrapolation of Birth<br>Weight</th><th colspan="2">Growth Velocity<br>Childhood</th><th colspan="2">Spurt of Growth<br>First Months</th><th colspan="2">Curvature Degree<br>First Months</th></tr><tr><th>β<sub>A</sub></th><th>(CI 95%)</th><th>β<sub>B</sub></th><th>(CI 95%)</th><th>β<sub>C</sub></th><th>(CI 95%)</th><th>β<sub>D</sub></th><th>(CI 95%)</th></tr></thead><tbody><tr><td colspan="9"><strong>Height</strong></td></tr><tr><td>Exclusively the first trimester</td><td>-0.004</td><td>(-0.012, 0.005)</td><td>0.027</td><td>(-0.002, 0.057)</td><td>-0.021</td><td>(-0.057, 0.015)</td><td>0.072</td><td>(0.004, 0.136)</td></tr><tr><td>Throughout pregnancy</td><td>-0.015</td><td>(-0.021, -0.008)</td><td>0.021</td><td>(-0.002, 0.043)</td><td>0.005</td><td>(-0.023, 0.033)</td><td>0.039</td><td>(-0.012, 0.092)</td></tr><tr><td colspan="9"><strong>Weight</strong></td></tr><tr><td>Exclusively the first trimester</td><td>0.006</td><td>(-0.028, 0.040)</td><td>-0.001</td><td>(-0.050, 0.048)</td><td>0.014</td><td>(-0.042, 0.067)</td><td>-0.014</td><td>(-0.097, 0.067)</td></tr><tr><td>Throughout pregnancy</td><td>-0.037</td><td>(-0.061, -0.012)</td><td>0.034</td><td>(0.000, 0.068)</td><td>0.009</td><td>(-0.032, 0.052)</td><td>0.022</td><td>(-0.035, 0.080)</td></tr></tbody></table>

Abbreviations: CI, credibility interval.

<footnote>ªAdjusted for maternal education level, BMI, and age at delivery, breastfeeding duration, recruitment center, and child's gender.</footnote>
<figure><img src="image_4.png" /><figcaption>**FIGURE 3** Posterior density of regression coefficients associated with maternal smoking during pregnancy and selected height and weight parameters</figcaption></figure>

age range. Consequently, maternal smoking throughout pregnancy was positively associated with BMI in those children. Between 570 and 600 days (19–20 months), the BMI of children whose mothers smoked throughout pregnancy already surpassed nonsmokers, reaching 0.21 kg/m² 1825 days after birth, ie, at 60 months (credibility interval: 0.03, 0.40).

Posterior inferences regarding the measures of association between maternal smoking and growth were not modified using noninformative priors in the sensitivity analysis (see supplementary material IV).

**5 | DISCUSSION**

We propose a method to jointly model height and weight in children 0 to 5 years of age. The simulation study indicated that the joint modeling yielded unbiased estimates and a good coverage probability, reflecting its ability to adequately fit theoretical growth data.
<table><caption>TABLE 6 Mean differences in body mass index (kg/m²), weight (g), and height (cm) between children exposed to maternal smoking during only the first trimester, children exposed throughout pregnancy, and the reference category (children from nonsmoking mothers during pregnancy)</caption><thead><tr><th rowspan="3">Variable</th><th colspan="4">Maternal Smoking</th></tr><tr><th colspan="2">First Trimester Exclusively</th><th colspan="2">Throughout Pregnancy</th></tr><tr><th>Posterior Coefficient</th><th>(95% Credibility Interval)</th><th>Posterior Coefficient</th><th>(95% Credibility Interval)</th></tr></thead><tbody><tr><td colspan="5">Difference of mean weight, g</td></tr><tr><td>5 days</td><td>19</td><td>(-77, 119)</td><td>-103</td><td>(-176, -29)</td></tr><tr><td>30 days (1 month)</td><td>20</td><td>(-81, 124)</td><td>-72</td><td>(-150, 6)</td></tr><tr><td>90 days (3 months)</td><td>29</td><td>(-100, 160)</td><td>-29</td><td>(-129, 69)</td></tr><tr><td>180 days (6 months)</td><td>45</td><td>(-107, 201)</td><td>2</td><td>(-115, 116)</td></tr><tr><td>270 days (9 months)</td><td>58</td><td>(-114, 231)</td><td>18</td><td>(-112, 146)</td></tr><tr><td>365 days (1 year)</td><td>69</td><td>(-122, 262)</td><td>31</td><td>(-112, 175)</td></tr><tr><td>730 days (2 years)</td><td>88</td><td>(-155, 328)</td><td>90</td><td>(-90, 269)</td></tr><tr><td>1095 days (3 years)</td><td>94</td><td>(-194, 384)</td><td>158</td><td>(-58, 376)</td></tr><tr><td>1470 days (4 years)</td><td>97</td><td>(-255, 449)</td><td>231</td><td>(-36, 505)</td></tr><tr><td>1825 days (5 years)</td><td>98</td><td>(-323, 520)</td><td>301</td><td>(-24, 639)</td></tr><tr><td colspan="5">Difference of mean height, cm</td></tr><tr><td>5 days</td><td>-0.15</td><td>(-0.57, 0.27)</td><td>-0.71</td><td>(-1.03, -0.38)</td></tr><tr><td>30 days (1 month)</td><td>-0.02</td><td>(-0.42, 0.37)</td><td>-0.60</td><td>(-0.90, -0.29)</td></tr><tr><td>90 days (3 months)</td><td>0.16</td><td>(-0.24, 0.55)</td><td>-0.42</td><td>(-0.73, -0.11)</td></tr><tr><td>180 days (6 months)</td><td>0.23</td><td>(-0.2, 0.67)</td><td>-0.29</td><td>(-0.61, 0.04)</td></tr><tr><td>270 days (9 months)</td><td>0.20</td><td>(-0.24, 0.65)</td><td>-0.24</td><td>(-0.56, 0.09)</td></tr><tr><td>365 days (1 year)</td><td>0.14</td><td>(-0.32, 0.60)</td><td>-0.22</td><td>(-0.56, 0.12)</td></tr><tr><td>730 days (2 years)</td><td>-0.01</td><td>(-0.58, 0.57)</td><td>-0.20</td><td>(-0.63, 0.24)</td></tr><tr><td>1095 days (3 years)</td><td>0.04</td><td>(-0.60, 0.69)</td><td>-0.12</td><td>(-0.61, 0.38)</td></tr><tr><td>1470 days (4 years)</td><td>0.19</td><td>(-0.54, 0.94)</td><td>0</td><td>(-0.56, 0.58)</td></tr><tr><td>1825 days (5 years)</td><td>0.37</td><td>(-0.48, 1.24)</td><td>0.13</td><td>(-0.52, 0.82)</td></tr><tr><td colspan="5">Difference of mean BMI, kg/m<sup>2</sup></td></tr><tr><td>5 days</td><td>0.14</td><td>(-0.09, 0.38)</td><td>-0.05</td><td>(-0.23, 0.12)</td></tr><tr><td>30 days (1 month)</td><td>0.08</td><td>(-0.13, 0.29)</td><td>0.07</td><td>(-0.09, 0.24)</td></tr><tr><td>90 days (3 months)</td><td>-0.01</td><td>(-0.25, 0.24)</td><td>0.16</td><td>(-0.04, 0.36)</td></tr><tr><td>180 days (6 months)</td><td>-0.02</td><td>(-0.27, 0.24)</td><td>0.15</td><td>(-0.04, 0.36)</td></tr><tr><td>270 days (9 months)</td><td>0.01</td><td>(-0.24, 0.27)</td><td>0.15</td><td>(-0.05, 0.35)</td></tr><tr><td>365 days (1 year)</td><td>0.06</td><td>(-0.20, 0.31)</td><td>0.15</td><td>(-0.04, 0.35)</td></tr><tr><td>730 days (2 years)</td><td>0.12</td><td>(-0.12, 0.36)</td><td>0.19</td><td>(0.02, 0.37)</td></tr><tr><td>1095 days (3 years)</td><td>0.09</td><td>(-0.13, 0.31)</td><td>0.21</td><td>(0.05, 0.38)</td></tr><tr><td>1470 days (4 years)</td><td>0.03</td><td>(-0.19, 0.25)</td><td>0.22</td><td>(0.05, 0.39)</td></tr><tr><td>1825 days (5 years)</td><td>-0.02</td><td>(-0.25, 0.21)</td><td>0.21</td><td>(0.03, 0.40)</td></tr></tbody></table>

The Bayesian framework offers a flexible technique to predict expected heights, weights, and BMIs of children, allowing the study of its potential determinants, such as sociodemographical profiles or environmental exposures. In our example, the joint model allowed the assessment of growth trajectories, the prediction of the mean differences over time based on maternal smoking and of the credibility intervals, for weight and height but also for BMI.
<figure><img src="image_5.png" /><figcaption>**FIGURE 4** Mean differences and 95% credibility intervals of body mass index (kg/m²), weight (kg) and height (cm) between children exposed to maternal smoking exclusively during the first trimester and the reference category (children from non-smoking mothers during pregnancy) (upper part of the graph); between children exposed throughout pregnancy and the reference category (lower part of the graph) from birth to 1825 days</figcaption></figure>

Our results showed that children from mothers who smoked throughout pregnancy had higher BMIs already in their first months of life compared to children from nonsmokers. At 5 years, they had a 0.21 kg/m² higher BMI. These results support previous studies that found higher BMIs or risk of being overweight in children from mothers who smoked during pregnancy.²⁹,³⁰ The potential mechanisms supporting the association between maternal smoking in pregnancy and the child's growth trajectory in the frame of developmental programming have been discussed elsewhere.⁷,³¹,³² The model provides key information about the dynamics of the association. As described elsewhere, contrarily to a direct BMI modeling approach, the different outcomes (eg, height and weight) do not need to be measured at the same time points in joint modeling. Moreover, the model does not assume that the same number of repeated measurements are available for the different outcomes.³³ Therefore, joint modeling enables the use of all available information on height and weight to predict BMI trajectory.

One of the main strengths compared to the frequentist approach is the fact that we can infer BMI and obtain cred-
ibility intervals for BMI trajectories of each category of maternal smoking during pregnancy in one step, using joint
models of height and weight and determine the differences between the trajectories of the subgroups of interest using
the posterior predictive distribution of the estimated parameters. Indeed, in the frequentist approach, the estimated
height and weight parameters can be used to estimate BMI in subgroups of interest and differences across subgroups,
but the uncertainty in the parameters is ignored. Consequently, the confidence intervals of the differences cannot be
assessed in a single step like the Bayesian approach and require additional simulation studies. In the Bayesian one,
the uncertainty of the parameters is considered in the predictions.

Another strength of our approach, related to structural models and not specific to the Bayesian framework, relies on
the interpretability of the growth parameters. They provide biological information on the individual growth processes
compared to the main growth modeling parameters of the nonstructural techniques such as splines or fractional polynomials. To our knowledge, a single recently published paper used a similar approach. In this study, the model developed by Cole et al³⁴ called the “Superimposition by Translation and Rotation” (SITAR) was generalized in the Bayesian framework to jointly model multiple outcomes of embryonic growth. This multivariate model described individual growth by 3 subject-specific parameters with an obvious meaning.² This model was applied to study prenatal growth but could be applied to postnatal growth data.

We intended to make the model as simple as possible using standard assumptions in NLMM methodology. However, distribution assumptions made about the data model could have been changed or made more complex in several ways. First, the Jenssen-Bayley growth curve function was parameterized according to a modified form (Equation 3) based largely on sign constraints on parameters A, B, and D resulting at the individual level in lognormal versus normal assumptions on these parameters. Instead, we could have considered the original form:

$$f(\boldsymbol{\phi}_{ik, t_{ijk}}) = A_{ik} + B_{ik} t_{ijk} + \exp(C_{ik}) [1 - \exp(D_{ik} t_{ijk})]$$

with the usual normality assumptions on $\boldsymbol{\phi}_{ik} = (A_{ik}, B_{ik}, C_{ik}, D_{ik})^T$.

However, because of convergence issues using this original form despite the use of truncated normal distribution for a priori distribution of the parameters, we used the modified form for the final model.

We assumed that there was no correlation $\sigma_{eHW} = 0$ between the residual components of height $e_{ijH}$ and weight $e_{ijW}$ recorded at the same occasion j and on the same individual i. This is a reasonable assumption as long as the correlation between the weight and height measurements at the same time originates from the subject i or the factors included in the model.

We selected among 2 candidate functions ($\log(\sigma_{eijk}^2) = \delta_{0k} + \delta_{1k}t_{ijk}$ and $\log(\sigma_{eijk}^2) = \delta_{0k} + \delta_{1k}t_{ijk} + \delta_{2k}t_{ijk}^2$) to describe the residual variance of height and weight but did not perform a formal comparison with other possibilities. Consequently, other functions, including functions with nonlinear terms, could have improved the aforementioned model.

The normality assumption of the residual errors could also have been investigated. For instance, the Student t-distribution has been proposed for heavy tailed or contaminated data with outliers or, more generally, the generalized error distribution for any type of kurtosis or skewness.³⁵

We took advantage of results from previous studies by specifying prior distributions of the mean and covariate effects. Regarding the interindividual variance covariance matrix, we adopted the inverse Wishart distribution. As previously described, this practice is by far the most common because of its conjugate properties and, therefore, simplicity, but having just one parameter as the degree of freedom forces to have the same uncertainty in all elements of the variance covariance matrix. However, the influence of every prior distribution was limited in our study considering the large amount of information in this dataset (both in terms of the total number of measurements on height and weight and in the number of individuals recorded).

In conclusion, because of the complexity of the nonlinear mixed-effects models used for the height and weight modeling, the computational facilities offered by the Bayesian approach add to the above-described advantages¹³ and is one of the key reasons we applied it. We showed that using datasets with rather large number of subjects and measurements by subjects, our model was suitable for growth modeling. The good properties of this joint model could be challenged in less favorable contexts. We assessed the associations between maternal smoking and multiple dimensions of growth, considering several potential confounders in these associations. The model we used could be extended by adding additional covariates of interest. However, the model is already complex, and the computational time is nonnegligible. The ability to fit more complex models could consequently be limited by practical aspects.

## ACKNOWLEDGMENTS

We thank Jean-Louis Foulley for introducing us to Bayesian modeling and providing us advices for the conduction of the analyses.

We thank Lise Giorgis-Allemand and Anne Forhan for data management, the midwife research assistants (Lorraine Douhaud, Sophie Bedel, Brigitte Lortholary, Sophie Gabriel, Muriel Rogeon, and Monique Malinbaum) for data collection, and the data entry operators (Patricia Lavoine, Josiane Sahuquillo, and Ginette Debotte).
REFERENCES

1. Hauspie RC, Cameron N, Molinari L. *Methods in Human Growth Research*. Cambridge University Press; 2004.

2. Willemsen SP, Eilers PHC, Steegers-Theunissen RPM, Lesaffre E. A multivariate Bayesian model for embryonic growth. *Statist. Med.* 2015;34:1351–1365.

3. Rizopoulos D, Ghosh P. A Bayesian semiparametric multivariate joint model for multiple longitudinal outcomes and a time-to-event. *Stat Med.* 2011;30:1366–1380.

4. Pinheiro JC, Bates DM. *Mixed-Effects Models in S and S-PLUS*. Springer Science & Business Media; 2000.

5. Jenss R, Bayley N. A mathematical method for studying the growth of a child. *Hum Biol.* 1937;9:556–563.

6. van Dommelen P, van Buuren S, Zandwijken GRJ, Verkerk PH. Individual growth curve models for assessing evidence-based referral criteria in growth monitoring. *Stat Med.* 2005;24:3663–3674.

7. Carles S, Charles M-A, Forhan A, et al. A novel method to describe early offspring body mass index (BMI) trajectories and to study its determinants. *PLoS One*. 2016;11: e0157766

8. Botton J, Scherdel P, Regnault N, Heude B, Charles M-A. Postnatal weight and height growth modeling and prediction of body mass index as a function of time for the study of growth determinants. *Ann Nutr Metab.* 2014;65:156–166.

9. Samson A, Lavielle M, Mentré F. The SAEM algorithm for group comparison tests in longitudinal data analysis based on non-linear mixed-effects model. *Stat Med.* 2007;26:4860–4875.

10. Levine RA, Casella G. Implementations of the Monte Carlo EM Algorithm. *Journal of Computational and Graphical Statistics*. 2001;10:422–439.

11. R Core Team (2015). R: A language and environment for statistical computing. [Internet]. R Foundation for Statistical Computing, Vienna, Austria; Available from: https://www.R-project.org/

12. Comets E, Lavenu A, Lavielle M. SAEMIX, an R version of the SAEM algorithm for parameter estimation in nonlinear mixed effect models [Internet]. 2012 [cited 2015 Jan 7]. Available from: http://papersjds14.sfds.asso.fr/submission_54.pdf

13. Huang Y, Liu D, Wu H. Hierarchical Bayesian methods for estimation of parameters in a longitudinal HIV dynamic system. *Biometrics*. 2006;62:413–423.

14. De la Cruz-Mesía R, Marshall G. Non-linear random effects models with continuous time autoregressive errors: a Bayesian approach. *Stat Med.* 2006;25:1471–1484.

15. Davidian M, Giltinan DM. *Nonlinear Models for Repeated Measurement Data*. CRC Press; 1995.

16. Foulley JL, San Cristobal M, Gianola D, Im S. Marginal likelihood and Bayesian approaches to the analysis of heterogeneous residual variances in mixed linear Gaussian models. *Computational Statistics & Data Analysis*. 1992;13:291–305.

17. Duval M, Robert-Granié C, Foulley J-L. Estimation of heterogeneous variances in nonlinear mixed models via the SAEM-MCMC algorithm with applications to growth curves in poultry. *Journal de la Société Française de Statistique*. 2009;150:65–83.

18. Donnet S, Foulley J-L, Samson A. Bayesian analysis of growth curves using mixed models defined by stochastic differential equations. *Biometrics*. 2010;66:733–741.

19. Anderson T. *An Introduction to Multivariate Statistical Analysis*. Second ed. Wiley; 1984.

20. Lunn D, Jackson C, Best N, Thomas A, Spiegelhalter D. *The BUGS Book: A Practical Introduction to Bayesian Analysis*. CRC Press; 2012.

21. Gelman A, Hill J. *Data Analysis Using Regression and Multilevel/Hierarchical Models*. Cambridge University Press; 2006.

22. Plummer M, Stukalov A. rjags: Bayesian graphical models using MCMC [Internet]. 2014 [cited 2014 Oct 9]. Available from: http://cran.r-project.org/web/packages/rjags/index.html

23. Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, Rubin DB. *Bayesian Data Analysis*. Third ed. CRC Press; 2013.

24. Gelman A, Rubin DB. Inference from iterative simulation using multiple sequences. *StatistSci*. 1992;7:457–472.

25. Cowles MK, Carlin BP. Markov Chain Monte Carlo convergence diagnostics: a comparative review. *J Am Stat Assoc*. 1996;91:883.

26. Spiegelhalter DJ, Best NG, Carlin BP, Van Der Linde A. Bayesian measures of model complexity and fit. *J R Stat Soc Series B Stat Methodology*. 2002;64:583–639.

27. Watanabe S. Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. *J Mach Learn Res.* 2010;11:3571–3594.

28. Heude B, Forhan A, Slama R, et al. Cohort profile: the EDEN mother-child cohort on the prenatal and early postnatal determinants of child health and development. *Int J Epidemiol*. 2015; dyv151

29. Chen A, Pennell ML, Klebanoff MA, Rogan WJ, Longnecker MP. Maternal smoking during pregnancy in relation to child overweight: follow-up to age 8 years. *Int J Epidemiol*. 2006;35:121–130.

30. Durmus B, Kruithof CJ, Gillman MH, et al. Parental smoking during pregnancy, early growth, and risk of obesity in preschool children: the generation R study. *Am J Clin Nutr*. 2011;94:164–171.
31. DiFranza JR, Aligne CA, Weitzman M. Prenatal and postnatal environmental tobacco smoke exposure and children's health. *Pediatrics*. 2004;113:1007-1015.

32. Toschke AM, Montgomery SM, Pfeiffer U, von Kries R. Early intrauterine exposure to tobacco-inhaled products and obesity. *Am J Epidemiol*. 2003;158:1068-1074.

33. Fieuws S, Verbeke G. Joint modelling of multivariate longitudinal profiles: pitfalls of the random-effects approach. *Statist Med*. 2004;23:3093-3104.

34. Cole TJ, Donaldson MDC, Ben-Shlomo Y. SITAR—a useful instrument for growth curve analysis. *Int J Epidemiol*. 2010;39:1558-1566.

35. Zhang Z. Bayesian growth curve models with the generalized error distribution. *Journal of Applied Statistics*. 2013;40:1779-1795.

## SUPPORTING INFORMATION

Additional Supporting Information may be found online in the supporting information tab for this article.

**How to cite this article:** Carles S, Charles M-A, Heude B, Ahmed I, Botton J, the EDEN mother-child study group. Joint Bayesian weight and height postnatal growth model to study the effects of maternal smoking during pregnancy. *Statistics in Medicine*. 2017;36:3990-4006. https://doi.org/10.1002/sim.7407