A Bayesian functional approach to test models of life course epidemiology over continuous time

Abstract Background Life course epidemiology examines associations between repeated measures of risk and health outcomes across different phases of life. Empirical research, however, is often based on discrete-time models that assume that sporadic measurement occasions fully capture underlying long-term continuous processes of risk. Methods We propose (i) the functional relevant life course model (fRLM), which treats repeated, discrete measures of risk as unobserved continuous processes, and (ii) a testing procedure to assign probabilities that the data correspond to conceptual models of life course epidemiology (critical period, sensitive period and accumulation models). The performance of the fRLM is evaluated with simulations, and the approach is illustrated with empirical applications relating body mass index (BMI) to mRNA-seq signatures of chronic kidney disease, inflammation and breast cancer. Results Simulations reveal that fRLM identifies the correct life course model with three to five repeated assessments of risk and 400 subjects. The empirical examples reveal that chronic kidney disease reflects a critical period process and inflammation and breast cancer likely reflect sensitive period mechanisms. Conclusions The proposed fRLM treats repeated measures of risk as continuous processes and, under realistic data scenarios, the method provides accurate probabilities that the data correspond to commonly studied models of life course epidemiology. fRLM is implemented with publicly-available software.


Simulation design
We further assess the proposed method, which estimates curves as Gaussian Processes (hereafter, GP), with additional simulations and provide comparisons approaches for estimating the functional curves, including the Principal Analysis by Conditional Estimation (PACE) method 1 , the Bayesian implementation of the Penalized Functional Regression (PFR) 2 , as well as with the RLM 3 , which is the discrete version of the fRLM.The PACE and PFR methods are two mainstream methods in functional regression for longitudinal data 4,5 and we adapt them to estimate the fRLM.The PFR method uses B-splines to model the curves, and the PACE method performs a functional regression on curves estimated with functional PCA.
We consider a standard scenario where ω(t) = cos(3πt + π/4) + 1 is non-sparse.Errors are from a normal distribution with variance σ 2 = 1.Parameters are set to δ = 3, α = 1 with C i = 1.Curves X i (t) were generated as random trigonometric functions, We consider different sample sizes (n ∈ {100, 400, 800}) and again a different number of measurement occasions: a sparse scenario, (N i ∈ {3, 4, 5}), and a moderately sparse scenario, (N i ∈ {6, 7, 8}).We generated random observed time points by t i,j as in the main simulation, but we restricted to the space such that there is at least one t i,j in each bin [0, 1/3], (1/3, 2/3], (2/3, 1], for each i, in order to permit the comparison with the discrete RLM.

Numerical implementation of competing methods
For the PFR, PACE, and GP methods, we use the same prior distribution for the Bayesian functional regression step, that is: For all approaches, we selected L = 9 B-splines basis function in order to balance flexibility of the approximations and computational burden.
For the PFR, we follow their approach and approximate the curves with K = 9 B-splines bases functions.For the PACE method, we used 15 functional Principal Components computed on a regular grid of size 150.
Finally, we provide comparisons with a discrete RLM with three periods.We aim to mimic the common approach in the literature when applying discrete models to longitudinal data with irregularly spaced observations.We divided the time interval [0, 1] into three periods of equal lengths: T 1 = [0, 1/3], T 2 = (1/3, 2/3], and T 3 = (2/3, 1].We average the observations collected on each interval and use them as a time varying exposure, that is we use the model where x t,i is the average of the observations on interval T t .

Results
We performed 100 Monte Carlo runs for each scenario and used 4 chains of 500 iterations with 250 as warmup.For each method, we computed the mean square error for δ and ω when available and report them in Table S1.We also report in Table S2 the mean integrated squared error of X i (t), where the integral is approximated with Riemann sums over a grid of length 150.Results were, as expected, very similar across all n and thus reported only for the largest n.Gaussian Processes are shown to outperform both PFR and the PACE method for both settings (very sparse and moderately sparse).
From the results shown in Table S1, we see that GP and PFR methods outperform the PACE method regarding ω(t).Regarding ω(t), the GP approach appears more advantageous in the 6 − 8 setting and PFR for the 3 − 5 setting.However, regarding δ, our approach outperforms the other models and the PFR method appeared to be the worst in all scenarios.Therefore, modeling curves with Gaussian Process methods appears to be superior for sparse data with irregularly spaced observations.This is likely to stem from the ability of the Gaussian Processes to better estimate the functional curves in sparse designs 6 .

Convergence diagnostics
We provide a summary of the convergence diagnostics for the Hamiltonian Monte Carlo (HMC) implementation for Bayesian estimation in the simulations and the data analysis.This includes quantitative statistics: the Gelman-Rubin statistics 7 (or Rhat), effective sample size (ESS bulk and tail), which indicates the number of independent draws in the sample 8 , and the proportion of number of divergent iterations after warmup (i.e.ratio of divergent iterations to total iteration after warmup).These statistics are reported for δ and ω and are provided in Tables S3, S4 and S5 for Simulation 1, the data analysis and Simulation 2 respectively.
The Rhat is a statistical measure used to determine if the chains are well mixed.Rhat should ideally be close to 1 and Rhat values greater than 1.2 suggest a lack of convergence.In our study, all parameters had Rhat values between 1.00 and 1.15, suggesting convergence for our chains.The effective sample size (ESS) is a measure that considers both the number of iterations and the autocorrelation between iterations.A higher ESS generally suggests better chain mixing, and is especially important in HMC, where autocorrelation can be substantial.ESS was always higher than 100, which is usually considered sufficient.The proportion divergent was always close to 0. For the data analyses, over 8000 total iterations, there was 14, 50, 46 divergent iterations for CKD, Inflammation, and Breast Cancer respectively.In addition, we report the traceplots for δ for each disease signature in Figure S1.Traceplots show that the chains are well mixed, indicating convergence.

Simulation based calibration
To validate our Bayesian model we used Simulation Based Calibration (SBC), which uses the selfrecovering property of Bayesian models 9,10 .For this purpose we used the R package SBC 11 .We used 100 simulations of a simple fRLM model with L = 4 and an intercept.If the algorithm is correctly calibrated, SBC rank statistics should be uniformly distributed.Results are shown in Figure S2 and S3 and suggest uniformly distributed rank statistics.

Figure S2 :
Figure S2: The distribution of SBC ranks.

Table S1 :
Performance metrics (median and MAD) over 100 replications for different estimation methods: GP, PFR, PACE and RLM.

Table S2 :
Summary statistics (median and MAD) of the Mean squared errors of the estimated curves for n = 800.

Table S3 :
Simulation study 1: median of the convergence statistics.

Table S4 :
Data analysis: median of the convergence statistics for each parameter.For the covariates parameters(alpha) we report the median over all coefficient.

Table S5 :
Simulation study 2: median of the convergence statistics.