Bayesian Multilevel Compositional Data Analysis:
Introduction, Evaluation, and Application

The data come from three studies with similar daily intensive designs and repeated measures: Activity, Coping, Emotions, Stress, and Sleep (ACES. N = 187); Diet, Exercise, Stress, Emotions, Speech, and Sleep (DESTRESS, N = 78); and Stress and Health Study (SHS, N = 96). Study materials are available on the Open Science framework for ACES (https://doi.org/10.17605/OSF.IO/H5497), DESTRESS (https://doi.org/10.17605/OSF.IO/QM63W), and SHS (https://doi.org/10.17605/OSF.IO/TZ48Y). Details of the data collection have been described previously [81, 113]. This data set had the data structure found in typical applications of multilevel analysis in psychological research (i.e., daily observations of movement behaviours are nested within individuals). For the purposes of this illustration, we used the data of 345 healthy adults, from whom we have repeated measurements of sleepiness and five movement behaviours of total sleep time, time awake in bed, moderate-to-vigorous physical activity (MVPA), light physical activity (LPA), and sedentary behaviour (SB). Sleepiness was a single item and self-reported 3-4 times daily, which was averaged to obtain the average daily level of sleepiness. The five behaviours were recorded via an actigraph for 7-12 days and scored using the GGIR R package [107, 106, 108, 109, 92]. These data are available from the corresponding author upon request.

The movement behaviours make up of a 5-part composition ($D$ = 5), which corresponds to a 4-dimensional set of $ilr$ coordinates. Individuals with missing data and zero values of any behaviours were excluded, as missing data and zeros hamper the analysis of compositional data (as the $ilr$ transformation cannot compute 0s). The two sets of between- and within-person $ilr$ coordinates were constructed using a SBP that represents the relative information of compositional parts as follows:

and

The $ilr$ coordinates represent the relative effects of behaviours (increasing some while at decreasing others), accounting for the constrained nature between behaviours within the 24-hour day. Specifically, across the between- and within-person levels, they represent the effects of (1) increasing sleep and time awake in bed while proportionally decreasing MVPA, LPA, and SB, (2) increasing sleep while proportionally decreasing time awake in bed, (3) increasing MVPA while proportionally decreasing LPA and SB, and (4) increasing LPA while proportionally decreasing SB.

We considered a Bayesian compositional MLM with a varying-intercept. The predictors are the 5-part movement composition, expressed as a total of 8 between- and within-person $ilr$ coordinates. The outcome of the model is next-day sleepiness. A varying-intercept by participants was included to account for non-independence. The model was fitted with weakly informative priors, 4 chains, and 4 cores, with 3000 iterations including 500 warmups (total of 10000 post-warmup draws), using CmdStanR [102] as back-end. We used weakly-informative priors, which play a minimal role in the computation of the posterior distribution, and maximise the influence of the data. For the population-level effects, student’s t distribution were used for the constant (i.e., fixed) intercept, and flat priors (improper priors over the reals) were used for the constant parameters of the predictors. Group-level effects also have their standard deviation parameters (i.e., varying intercept and residual), which were specified using student’s t distribution. The priors for the standard deviation parameters are restricted to be non-negative and have a half student-t prior with 3 degrees of freedom and a scale parameter that depends on the standard deviation of the outcome. These priors are only weakly informative, but provide some regularisation to improve convergence and sampling efficiency [67]. Prior information is given in Table 2.

The Bayesian compositional substitution MLM was then conducted for both between- and within-person levels, examining the predicted change in sleepiness associated with the pairwise reallocation from 1 to 30 minutes between the composition of 24-hour movement behaviours.

Significance of individual parameters was assessed using the Bayesian 95% posterior credible interval, with 95% credible intervals (CIs) not containing 0 providing evidence for the probability that the true estimate would lie within the interval. All analyses were performed in R v4.3.1 [97], using package multilevelcoda v1.1.0 (model estimation, workflow outlined in Figure 2), brms [66], future [[, parallel processing,]]future, and ggplot2 [[, results visualisation,]]ggplot. All analysis code is available at: https://github.com/florale/multilevelcoda-sim.

Results from the MLM predicting next-day sleepiness from a 5-part composition are presented in Table 3, supporting the effects of all within-person $ilr$ coordinates (indicated by 95% CIs not containing 0s), but not any between-person $ilr$ coordinates (indicated by 95% CIs containing 0s). This demonstrated the relationships between movement behaviours and next-day sleepiness occurred only at within-person level, but not between-person level. Overall, the $1^{st}$ within-person $ilr$ coordinate (longer time spent on sleep behaviours than usual [sleep and time awake in bed], relative to wake behaviours [MVPA, LPA, and SB]) predicted -0.59 [95% CI -0.70, -0.49] lower next-day sleepiness. The $2^{nd}$ within-person $ilr$ coordinate (longer sleep than usual, relative to spending time staying awake in bed), also predicted lower -0.44 [95% CI -0.55, -0.34] next-day sleepiness. Similarly, the $3^{rd}$ within-person $ilr$ coordinate (higher-than-usual MVPA, relative to LPA and SB) and the $4^{th}$ within $ilr$ coordinate (higher-than-usual LPA relative to SB), predicted lower sleepiness (-0.27 [95% CI -0.39, -0.16] and -0.20 [95% CI -0.35, -0.06], respectively).

Consistent with the main MLM, reallocation of time between movement behaviours predicted changes in sleepiness at the within-person level, but not the between-person level. Individuals slept longer-than-usual at the expense of any behaviours, except MVPA, at within level, experienced lower level of next-day sleepiness. However, when individuals sacrificed their sleep on a given day is for any other behaviours (i.e., including MVPA), they experienced a higher level of sleepiness. Additionally, individuals who spent longer time in LPA at the expense of time awake in bed on a given day, also experienced a higher level of sleepiness the next day, and vice versa. Results of the substitution model for 30-minute reallocation are in Table 4. For brevity, we presented only the significant results for the reallocation from 1 to 30 minutes of Sleep and Awake in bed, respectively, in Figure 3.

We created a range of simulation conditions including different values of the number of clusters ($J$), cluster size ($I$), the number of compositional parts ($D$), and the magnitude of sample variability (assessed by the varying-intercept variance $\sigma^{2}_{u}$ and residual variance $\sigma^{2}_{\varepsilon}$). The values for the number of clusters and cluster sizes were informed by our review of a systematic review and meta-analyses on daily sleep and physical activity [60]. Given the different number of compositional parts used in existing studies, we constructed models with different numbers of possible compositional parts and assessed their performances using different sets of ground truth values. Finally, we examined the influences of sample variability, including varying-intercept variance ($\sigma^{2}_{u}$) and residual variance ($\sigma^{2}_{\varepsilon}$) on the estimation of our models. Table 5 summarises the factors and their levels considered in this simulation study. The combination of these factors resulted in a total of 240 scenarios. For each cell of the simulation design, 2000 replications were generated ($n_{sim}=2000$), resulting in $4[I]\times 4[J]\times 3[D]\times 5[\sigma]\times 2000=$ 480 000 data sets to be analysed.

In the following, we described the simulation procedure to generate data sets resembling the data structure used in real data study. The varying-intercept $u_{0j}$ was generated from $\text{Normal}(0,\sigma^{2}_{u})$. The design matrices of the predictors, the between-person $ilr$ ($z^{(b)}$) and within-person $ilr$ ($z^{(w)}$) corresponding to 5-part composition (total sleep time, time in bed awake, MVPA, LPA, and SB) were generated, respectively, as follows:

and

with values of the means and covariances informed by the data set used in the real data study. Compositional data were then generated by inverse-transforming the 4-dimension $ilr$ coordinates. At this step, when necessary, the 4-part and 3-part compositions were created by collapsing variables. The 4-part composition was obtained by collapsing total sleep time and wake in bed to a single variable named sleep. The 3-part composition was created by collapsing MVPA and LPA to a single variable named physical activity. These compositions were transformed again to $ilr$ coordinates for model estimation.

The outcome vector $\boldsymbol{y}$ was generated from a normal distribution:

with the values for the constant parameters set to be close to those found in the real data study.

The primary estimands of the simulation study are the parameters of the Bayesian MLM models, including the constant parameter estimates: the intercept ($\gamma_{0}$), the between-person and within-person $ilr$ coordinates ($\beta$s), and the varying parameters: varying-intercept ($\sigma_{u}$) and residual error ($\sigma_{\varepsilon}$). For the Bayesian compositional substitution MLM, estimation of predicted change in outcome at between-level ($\Delta{\hat{y}^{(b)}_{ij}}$) and within-level ($\Delta{\hat{y}^{(w)}_{ij}}$) were evaluated for all possible pairwise substitution between compositional parts, totalling to $2\times D\times(D-1)$ parameters.

Model performance of 2000 replications across 240 conditions was evaluated using the following criteria.

1. 1.

   Quality of the MCMC-based sampling procedure of the Bayesian compositional MLM. We considered the proportion of replication that sufficiently converged [[, $\hat{R}<1.05$, ]]vehtari2021 and had no divergent transition. Effective sample size (ESS) was investigated both at the bulk of the distribution (e.g., for the mean or median) and in the tails (e.g., for posterior interval estimates and inferences about extreme quatiles). Any parameters with ESS $>$ 400 indicated sampling inefficiency and required further diagnostics [110].

2. 2.

   Quality of model performance was evaluated in terms of accuracy in parameter estimates and inference, using three performance measures: bias, coverage, and bias-eliminated coverage [93]. Monte Carlo standard errors were used to calculate 95% confidence interval.

Using package multilevelcoda, each simulated data set was fitted in Bayesian MLM with a varying-intercept to predict next-day sleepiness from the $D$-part behaviour composition, expressed as a total of $2(D-1)$ between- and within-person $ilr$ coordinates. The Bayesian substitution MLM was then conducted, and the model performance in estimating the predicted change in outcome for 30-minute reallocation was evaluated. The simulation study results were summarised using package rsimsum [73] and visualised using package ggplot2 [112]. Reproducible material for this study is available at: https://github.com/florale/multilevelcoda-sim.

Divergences were observed in 1312 replications (0.00%), predominantly with small number of clusters (73.6% $J$: 30), small cluster size (90.5% $I$: 3), and large residual variation (97.6% $\sigma^{2}_{\varepsilon}:1.5$). An additional 17 (0.00%) replications had $\hat{R}>1.05$, demonstrating convergence issues. These replications were excluded for the evaluation of parameter estimates and inference.

In contrast, low bulk ESS were observed as sample size increases with large between-person heterogeneity and small within-person heterogeneity. Particularly, 27 651 replications (0.06%) had bulk $\text{ESS}<400$ for some parameters, predominantly with large number of clusters (51.1% $J$: 1200), large cluster size (70.1% $I$: 14), and small residual variation (95.6% $\sigma^{2}_{\varepsilon}:0.5$). The low ESS values under these conditions represent a technical difficulty posed by the MCMC sampling methods, wherein small variation in the sample (i.e., $\sigma^{2}_{\varepsilon}$) cause the sampler to produce higher within-chain correlation [64]. Additionally, the default non-centered parameterisation [[, i.e., separation of population parameters and hyperparameters in the prior, ]]papaspiliopoulos2007 used in our model estimation procedure can be less efficient for large data sets and strong likelihood (i.e., small sample variability), compared to centered parameterisation [64]. Therefore, we conducted a case study [82] into a replication generated using a 3-part composition, 1200 clusters and cluster size of 14, with large varying intercept variation ($\sigma^{2}_{u}=1.5$) and residual variation ($\sigma^{2}_{\varepsilon}=0.5$), wherein model produced low bulk ESS values for 4 out of 7 parameters. Posterior predictive distributions were checked and two methods to improve the MCMC sampling were tested: centered-parameterisation and increased iterations. We found no evidence of non-convergence (e.g., poor mixing of chains or funnel degeneracy in the posterior). Both reparameterisation or increasing iterations and warm-ups improved ESS, with centered parameterisation showing substantial gain of ESS for the same number of iterations. A sensitivity analysis comparing the model performance with and without the replications with low ESS revealed that ESS did not have an influence on the quality of parameter estimates and inference. Replications with low ESS were, therefore, kept in the subsequent evaluation of model performance.

Across the simulated conditions, both Bayesian compositional MLMs and Bayesian compositional substitution MLMs yielded negligible biases in the estimation of all parameters, with a mean of 0.00 and a range from -0.09 to 0.05 and mean of 0.00 and range from -0.03 to 0.04, respectively. Both models had coverage and bias-eliminated coverage close to the nominal value, with means of 0.95 and ranges from 0.93 to 0.97.

As there was no impact of simulation conditions on the model performance, for brevity, we presented the results for individual parameters estimated using composition with 4 parts ($D$ = 4) and a medium level of modelled variance ($\sigma^{2}_{u}=1$ and $\sigma^{2}_{\varepsilon}=1$) under different conditions of the number of clusters ($J$) and cluster size ($I$). Full results can be accessed via the shiny app by locally running in R using package multilevelcoda. Both models performed well consistently across the number of clusters and cluster size, as demonstrated in Figure 4, 5, 6, and 7.