Fitting Bayesian Stochastic Differential Equation Models with Mixed Effects through a Filtering Approach

Abstract

Recent advances in technology contribute to a fast-growing number of studies utilizing intensive longitudinal data, and call for more flexible methods to address the demands that come with them. One issue that arises from collecting longitudinal data from multiple units in time is nested data, where the variability observed in such data is a mixture of within-unit changes and between-unit differences. This article aims to provide a model-fitting approach that simultaneously models the within-unit changes with differential equation models and accounts for between-unit differences with mixed effects. This approach combines a variant of the Kalman filter, the continuous-discrete extended Kalman filter (CDEKF), and the Markov chain Monte Carlo method often employed in the Bayesian framework through the platform Stan. At the same time, it utilizes Stan’s functionality of numerical solvers for the implementation of CDEKF. For an empirical illustration, we applied this method in the context of differential equation models to an empirical dataset to explore the physiological dynamics and co-regulation between couples.

Keywords:

• Mixed models
• dynamic modeling
• Bayesian analysis
• differential equation models
• Kalman filter

Article information

Conflict of interest disclosures: Each author signed a form for disclosure of potential conflicts of interest.

Ethical principles: The authors affirm having followed professional ethical guidelines in preparing this work. These guidelines include obtaining informed consent from human participants, maintaining ethical treatment and respect for the rights of human or animal participants, and ensuring the privacy of participants and their data, such as ensuring that individual participants cannot be identified in reported results or from publicly available original or archival data.

Funding: This work was supported by Grant U24AA027684 and R01HD076994 from the National Institute of Health (NIH), as well as Grant IGE-1806874 and BCS-05-27766 from the National Science Foundation (NSF).

Role of the funders/sponsors: None of the funders or sponsors of this research had any role in the design and conduct of the study; collection, management, analysis, and interpretation of data; preparation, review, or approval of the manuscript; or decision to submit the manuscript for publication.

Acknowledgments: The authors would like to thank Drs. Murali Haran and Erika Lunkenheimer for their comments on prior versions of this manuscript. The ideas and opinions expressed herein are those of the authors alone, and endorsement by the authors’ institutions and the funding agency is not intended and should not be inferred.

Notes

1 In principle, the estimation method described here can be implemented using other Bayesian software as well. Stan was selected from the available Bayesian platforms because of the relative ease to write user-defined functions, and the readily available built-in tools for solving differential equations.

2 Code to replicate this example can be found at https://github.com/mchdfs/MBR-MLSDE-Stan-Supplementary

3 Using high-density intervals yielded the same conclusions in this illustration.

4 In our implementation, we did not constrain the priors for the self-regulation parameters to only take on non-negative values, based on our definition of self-regulation and previous exploration on the parameter space demonstrated by Figure 6. However, we do realize the problem-specific need for researchers to impose a stricter assumption on self-regulation by constraining the relevant parameters. Sample scripts with constrained and unconstrained priors are both provided on https://github.com/mchdfs/MBR-MLSDE-Stan-Supplementary

5 The default for the starting value of an unconstrained parameter is a draw from U(-2, 2). In our case, standard deviation parameters have a lower bound of zero, thus starting values for those drew from U(0, 2).

6 The 26 positive $b_{1i}s$ were not necessarily from the same dyads as the 27 positive $a_{2i}s.$

7 This function requires the starting time index to be a real number and the end time index a real-numbered array, as the ODE solver is often used to solve for the entire sequence instead of the next time index. Therefore, we specified ts as an array with dimensions $T\times 1$ to get around these requirements.

8 https://github.com/mchdfs/MBR-MLSDE-Stan-Supplementary