As of January 1970 to December 2021, PubMed.gov (https://pubmed.ncbi.nlm.nih.gov/, accessed on 23 February 2022) searched of “nonlinear mixed-effect models” yielded 6288 publications. Among the published articles, nearly 94% of works used frequentist approaches (5929 articles), while only 6% of works adopted Bayesian approaches (359 articles).Figure 1 displays a bar plot based on the published articles, categorized by the frequentist and Bayesian approaches over time. In the panel, it is observed that, until the late 1980s, the Bayesian research was nearly dormant, but since the early 1990s, Bayesian works begin to re-emerge, and the gap between frequentist and Bayesian works becomes gradually narrower as time evolves.

The dormancy of the Bayesian approaches until the late 1980s is mainly due to the time-consuming nature of the calculations based on a sampling scheme, which was previously impossible by the limitation of computing power. Fortunately, a breakthrough in computer processors (for, e.g., Am386 in 1991, Pentium Processor in 1993, etc.) took place in the early 1990s, driving the computing revolution to solve computationally intense problems, and this has given the statistical community the ability to solve statistical questions by using Bayesian methods. This timeline is also aligned with the widespread Markov chain Monte Carlo (MCMC) sampling techniques in the Bayesian community [32,33]. Since then, the Bayesian community has been gradually gaining the momentum to leverage the rapidly growing developments of computing power, and now, assorted Bayesian software packages (e.g.,JAGS [34],BUGS [35], and Stan [17]) are available for researchers to answer scientific questions arising from industrial and academic research.

To understand the rise of the Bayesian approaches, we first want to understand what will be some of the advantages of using Bayesian methods over frequentist methods in the context of nonlinear mixed effects models. As the primary focus of this review paper is to provide the readers with some insight on methodologies and practical implementation of using Bayesian approaches, our comparison and exposition below are described from an operational viewpoint.Table 1 summarizes the modeling strategies of using the frequentist and Bayesian approaches for the nonlinear mixed effects models. Broadly speaking, the usual estimation method of the frequentist computation is optimization, while that of the Bayesian computation is sampling. Normally, it is known that the former is much faster than the latter. This is not surprising because a sampling scheme, by its nature, needs to explore a wide range of the parameter space, whereas the optimization only needs to find the best point estimate, which is often described by the maximum likelihood estimate. In many practical problems, widely used frequentist optimization algorithms are the first-order approximation [36], Laplace approximation [37], and stochastic approximation of expectation-maximization algorithm [38]. They will be briefly discussed inSection 5.4. As for the Bayesian sampling algorithms, combinations of Gibbs sampler [39], Metropolis-Hastings algorithm [40], Hamiltonian Monte Carlo [41], and No-U-Turn sampler [42] are popularly used, among many others [43,44,45]. We explain these in detail inSection 6.

The stark difference between using frequentist and Bayesian approaches may be the procedure of describing an uncertainty underlying the parameter estimation for the nonlinear mixed effects models. Here, the parameter which is of primary interest is the population-level parameters (also called fixed effects), typical values for the individual-level parameters. In many cases, frequentist$95\%$ confidence intervals for the parameters of the models are constructed by assuming that asymptotic normality of maximum likelihood estimator holds in a finite sample study, which is actually the most accurate in large sample scenario [51]. Most frequentist software packages, such asNONMEM [13,47],Monolix [48], and nlmixr [18], by default, may print out a$95\%$ confidence interval of the form, “Estimate ± 1.96 × Standard Error”, or some transformation of the lower and upper bounds, if necessary, such that theStandard Error is calculated by using (observed) Fisher information matrix [52,53,54]. Using such a scheme in small sample studies is highly likely to overlook the gap between the reality of the data and the idealistic asymptotic situation.

In contrast, as for the Bayesian approaches, the large-sample theory is not needed for the uncertainty quantification, and the procedure to obtain$95\%$ posterior credible intervals is a lot easier than obtaining$95\%$ confidence intervals (See Chapter 4 of [55]). Furthermore, Bayesian credible intervals based on percentiles of posterior samples allow for a strongly skewed distribution, wherein frequentist confidence intervals (based on large-sample theory) may induce a non-negligible approximating error due to the deviation from the asymptotic normality. Along with that, Bayesian methods are highly appreciated when researchers wish to incorporate prior knowledge from previous studies into the model so that posterior inference provides the researchers with an updated view on the problem, possibly with a more accurate estimation. Using prior information would be particularly useful in small-sample contexts [56,57]. For example, in medical device clinical trials, some opportunities and challenges in developing a new medical device are: (i) there is often a great deal of prior information for a medical device; (ii) a medical device evolves in relatively small increments from previous generations to a new generation; (iii) there are only a few numbers of patients for the trials; and (iv) companies need to make a rational decision promptly to reduce cost. In those settings, Bayesian methods have been demonstrated to be suitable, and their proper use is guided by Food and Drug Administration [58,59].

We outline the first two steps in the Bayesian workflow of using Bayesian nonlinear mixed effect models described inFigure 2. The panel includes some mathematical notations that are consistently used throughout the paper. These notations will be clearly understood later. The aim of our explanation at this point is to provide readers with a blueprinted plan to implement Bayesian modeling strategies for the nonlinear mixed effects models. We assume that readers are familiar with basic concepts and generic workflow in Bayesian statistics; see [55,60,61] for those basic concepts and refer to the review paper by [62] and references therein for detailed concepts and general terminologies used in workflow, such as prior and posterior predictive checks, and prior elicitation, etc.

The first step of the Bayesian research cycle is (a) standard research cycle [63,64]. Some early activities at this step involve reviewing literature, defining a problem, and specifying a research question and a hypothesis. After that, researchers specify which analytic strategy would be taken to solve the research question and suggest possible model types, followed by data collection. The data type arising in this process may include a response variable and some covariates that are grouped longitudinally, which then formulates repeated measures data of a population of interest. Furthermore, if there appears to be some nonlinear temporal tendency at each subject, then a possible model type for the analysis is a nonlinear mixed effects model [15,65].

The second step of the Bayesian research cycle is (b) Bayesian-specific workflow. Logically, the first thing to do at this step is to determine prior distributions (see Step (b)–(i) inFigure 2). The selection of priors is often viewed as one of the most crucial choices that a researcher makes when implementing a Bayesian model, as it can have a substantial impact on the final results [66]. As exemplified earlier in the context of Bayesian medical device trials, using a prior in small sample studies may improve the estimation accuracy, but an unthoughtful choice of priors would lead to a significant bias in estimation. Prior elicitation effort would require Bayesian expertise to formulate domain expert’s knowledge in a probabilistic form [67]. Strategies for prior elicitation include asking domain experts to provide suitable values for the hyperparameters of the prior [68,69]. After prior is specified, one can check the appropriateness of the priors through prior predictive checking process [70]. For almost all practical problems, prior distribution of Bayesian nonlinear mixed effect models can be hierarchically represented as follow: (1) a prior for the parameters used in likelihood, often called ‘population-level model’ in the literature of mixed effects modeling; and (2) a prior for the parameters used in the population-level model and for the parameters describing the residual errors used in likelihood. It is important to note that the former type of prior distribution (that is, (1)) is also a requirement to implement frequentist approaches for the nonlinear mixed effects model, as a name of ‘distribution for random effects’. Essentially, the defining factor of the Bayesian framework is the latter type of prior distribution (that is, (2)), which is fixed in the frequentist framework, as a name of ‘fixed effects’. Some prior options of the latter type will be discussed inSection 7.

The second task is to determine the likelihood function (see Step (b)–(ii) in the panel). At this time, the raw dataset collected in (a) standard research cycle should be cleaned and preprocessed. Before embarking on more serious statistical modeling, it is a common practice to get some insight about the research question via exploratory data analysis and have a discussion with domain experts such as clinical pharmacologists, clinicians, physicians, engineers, etc. To some extent, eventually, all these efforts are to determine a nonlinear function (denoted asf in this paper) that best describes the temporal profiles of all subjects. This nonlinear function is aknown function because it should be specified by researchers. In other words, the branch of the nonlinear mixed effects models belongs to parametric statistics. However, one technical challenge is that, in many problems, such a nonlinear function is represented as a solution of a differential equation system [71,72], and therefore there is no guarantee that we can conveniently work with a closed-form expression of the nonlinear function. For example, if researchers wish to work with nonlinear differential equations [73,74], then some approximation via differential equation solver [75,76] may be needed to calculate the nonlinear function. As such, most software packages dedicated to implementing a nonlinear mixed effect model, or, more generally, a Bayesian hierarchical model, are equipped with several built-in differential equation solvers [17,47,77]. For instance, visit the website (https://mc-stan.org/docs/2_29/stan-users-guide/ode-solver.html, accessed on 20 February 2022) to see some functionality supported in Stan [17].

Finally, the likelihood is combined with the prior to form the posterior distribution (see Step (b)–(iii) in the panel). Given the important roles that the prior and the likelihood have in determining the posterior, this step must be conducted with care. The implementational challenge at this step is to construct an efficient MCMC sampling algorithm. The basic idea behind MCMC here is the construction of a sampler that simulates a Markov chain that is converging to the posterior distribution. One can use software packages if prior distributions to be implemented in Bayesian models exist in the list of prior options available in the packages. Otherwise, professional programmers and Bayesian statisticians are needed to make codes manually; this review paper will be useful for that purpose. Another activity important at this step is to compare multiple models with different priors and nonlinear functions, specified in Step (b)–(i) and (ii), and select the best model out of them. This topic is broadly called the model selection [78], which will be discussed inSection 8.

To exemplify circumstances for which the nonlinear mixed effects model is a suitable modeling framework, we review challenges from several diverse applications.Table 2 summarize four real-life problems that will be illustrated in the next subsections.

Studies of the pharmacokinetics of drugs help us learn about the variability in drug disposition in a population [92].Figure 3 shows theophylline concentration in the plasma as a function of time after oral administration of the same amount of anti-asthmatic theophylline for 12 subjects (the data considered here are courtesy of Dr. Robert A. Upton of the University of California, San Francisco.) As seen in the panel, concentration trajectories have a similar functional shape for all individuals. However,$C_{max}$ and$t_{max}$ (peak concentration and time when it is achieved), absorption, and elimination phases are substantially different across subjects. Clinical pharmacologists believe that these differences are attributable to between-subject variation in the underlying pharmacokinetic processes, explained by Absorption, Distribution, Metabolism, and Excretion (ADME), understanding of which is crucial in a new drug development in the pharmaceutical industry.

In pharmacokinetics analysis, often abbreviated by ‘PK analysis’, it is routine to use compartmental modeling to describe the amount of drug in the body by dividing the whole body into one or more compartments [93]. For theophylline, a one-compartment model is normally used, which assumes that the entire body acts like a single, uniform compartment; see page 30 from [94] for a detailed explanation about the model: where$C\left(t\right)$ is drug concentration at timet for a single subject following oral doseD at$t=0$. Here,F is the bioavailability which expresses the proportion of a drug that gains access to the systemic circulation.$k_a$ is the absorption rate constant describing how quickly the drug is absorbed from the gut into the systemic circulation.V is the volume of the central compartment.$Cl$ is the clearance rate representing the volume of plasma from which the drug is eliminated per unit time. Eventually, the pharmacokinetic processes for a given subject is summarized by the 4-dimensional vector with ‘PK parameters’$(F,k_a,V,Cl)$. Obviously, it is the modeler’s discretion to proceed with a more complex PK model such as a three compartment models with nonlinear clearance to fit the data, but in this case, over-parameterization should be carefully examined [95].

Typically, the dataset collected in a drug development program includes demographic and clinical covariates obtained from each subject, for, e.g., body weight, height, age, sex, creatinine clearance, albumin, etc.; and furthermore, one can also involve genetic information in an individual’s response to drugs. Most covariates are measured at baseline, before assigning the drug, while some covariates can be measured at every sampling time. One of the crucial goals of PK analysis is to illustrate the effect of such covariates on the PK parameters [96]. The causal relationship inferred by the covariate analysis can be used to support physicians in making the necessary judgments about the medicines that they prescribe, tailored to individual patients [97].

In a PK report for a new drug application to government authorities like U.S. Food and Drug Administration (FDA) or European Medicines Agency (EMA), the PK parameters are summarized by mean or median, and very importantly, estimates of parameter precision. Estimates of parameter precision can provide valuable information regarding the adequacy of the data to support those parameters [98]. Parameter uncertainty can be estimated through several methods, including bootstrap procedures [99], log-likelihood profiling [100], or using the asymptotic standard errors of parameter estimates, and recently, Bayesian approaches draw a lot of attention from the pharmaceutical industry [101]. Particularly, Bayesian approach for the population PK analysis can be very useful when there is prior knowledge about PK parameters learned from preclinical studies, published works, etc., and one wants to incorporate them into the prior specification for PK parameters [28].

The US shale boom—a product of technological advances in horizontal drilling and hydraulic fracturing that unlocked new stores of energy—has greatly benefited the growth in the US economy. Horizontal drilling is a directional drilling technology where a well is drilled parallel to the reservoir bedding plane [102]. Well productivity of a horizontal well is known to often be 3 to 5 times greater than that of a vertical well [103,104], but it also costs 1.5 to 2.5 times more than a vertical well [105]. Therefore, the eventual success of the drilling project of unconventional shale wells relies on a large degree of well construction costs [106]. Because of very low permeability, and a flow mechanism very different from that of conventional reservoirs, estimates for the shale well construction cost often contain high levels of uncertainty. For this reason, one of the crucial tasks of petroleum engineers is to quantify the uncertainty associated with the process of oil or gas production to reduce the extra initial risk for the projects.

Figure 4 shows monthly production rate trajectories of 360 shale oil wells completed in the Eagle Ford Shale of South Texas, studied by [31]. The declining pattern manifested in the trajectories is commonly observed in almost all oil production rate time series data following well completion. Here, the completion is terminology in petroleum engineering, meaning the process of transforming a well ready for the initial production [107]. Decline curve analysis (DCA), introduced by [83] around 100 years ago, is one of the most popularly utilized methods for petroleum engineers. Its purpose is to (i) theorize a curve describing the declining pattern, (ii) analyze the declining production rates, (iii) characterize the well-productivity, and (iv) forecast the future performance of oil and gas wells. Particularly, estimation and uncertainty quantification of estimated ultimate recovery (EUR) (here, EUR is a special jargon defined as an approximated quantity of oil from a well that is potentially recoverable by the end of its producing life [108]) is the utmost important task and a starting point in the decision-making process for future drilling projects. In addition, the oil and gas companies comply with financial regulations for EUR outlined by the U.S. Securities and Exchange Commission: seehttps://www.sec.gov/rules/final/2008/33-8995.pdf (accessed on 20 February 2022), for the regulations.

Most curves used in DCA are derived from solving certain differential equations that describe a hidden dynamic from production rate trajectory [109,110,111,112,113,114]. See [84,85,115,116] for an overview of such curves. Ref. [31] studied Arps’ hyperbolic, stretched exponentiated decline, Duong, and Weibull curves to fit the trajectories shown in theFigure 4. Particularly, the Duong model was developed for unconventional reservoirs with very low permeability: where$P\left(t\right)$ is the production rate at timet for a single well following completion.$q_1$ is the initial rate coefficient, and m anda are additional model parameters. We note that the parameters,$q_1$,m anda, have their own meanings in terms of well-productivity: see [114] for the interpretation. That being said, the well-productivity for a given well is summarized by the 3-dimensional parameter vector,$(q_1,m,a)$. In a modeling perspective, the variation of the well-productivity across different wells is attributable to the different values for$(q_1,m,a)$. To explain this variability, one can regress the values$(q_1,m,a)$ on the well-design parameters such as true vertical depth, measure depth, etc. The causal relationship inferred by the covariate analysis will be used in a future drilling project. Geological information of wells can also be incorporated to make a spatial prediction for the EUR at a new location, as researched by [31].

Macroeconomists, financial economists, and market participants all attempt to build good models of the ‘yield curve’ [117]. The yield curve on a given day is a curve showing the interest rates across different maturity spans (one month, one year, five years, etc.) for a similar debt contract at a particular date. It determines the interest rate pattern (i.e., cost of borrowing), which can be used to calculate a bond’s price [118].Figure 5 shows daily treasury par yield curve rates spanning from 3 to 13 January 2022, with maturities up to 30 years. The data source is from the U.S. Department of the treasury (https://www.treasury.gov/resource-center/data-chart-center/interest-rates/Pages/TextView.aspx?data=yield, accessed on 20 February 2022). As seen from the panel, the shape of the yield curve displays a slightly delayed humped shape. Economists believe that such a shape of the yield curve has an important implication on the economic growth [119].

The Nelson–Siegel model [86] is a very popular model in the literature to fit the term structure: where$Y\left(\tau \right)$ denotes the (zero-coupon) yield evaluated at$\tau$, and $\tau$ denotes the time to maturity. The model parameters have a specific financial meaning:${\beta }_0$,${\beta }_1$, and ${\beta }_2$ are related long-term, short-term, and mid-term effects on the interest rate, respectively, and $\lambda$ is referred to as a decay factor [87]. Each of the yield curves is summarized by the 4-dimensional parameter$({\beta }_0,{\beta }_1,{\beta }_2,\lambda )$, and it is known that the model can capture a wide range of possible shapes of the yield curve [86,87,120,121]. Therefore, the Nelson–Siegel model is extensively used by central banks and monetary policymakers [122]. For example, The Federal Reserve updates estimates of$({\beta }_0,{\beta }_1,{\beta }_2,\lambda )$ once per week: visit the website (https://www.federalreserve.gov/data/yield-curve-tables/feds200628_1.html, accessed on 20 February 2022). In recent years, there has been a great deal of interest in the uncertainty quantification of the Nelson–Siegel parameters over time, and their relationship with macroeconomic variables such as inflation and real activity, etc, in financial applications: refer to  [121,123,124,125] for some of those works.

Novel coronavirus disease 2019 (COVID-19) is a big threat to global health. The rapid spread of the virus has created a pandemic, and countries all over the world are struggling with a surge in COVID-19 infected cases.Figure 6 displays the daily infection trajectories describing the cumulative numbers of infected cases for 40 countries, spanning from 22 January to 14 May 2020, studied by [30]. The data source is from COVID-19 Data Repository by the Center for Systems Science and Engineering at Johns Hopkins University (https://coronavirus.jhu.edu/map.html, accessed on 20 February 2022). Refer to Table S.1 in [30] for the list of 40 countries. The time frame of the authors’ research was the early stage of the pandemic when there was no drug or other therapeutics approved by the US FDA.

In general, during an early phase of a pandemic, information regarding the disease is very limited and scattered, even if it exists. In spite of that, it is crucial to predict future cases of infection or death. In such a situation, one consideration is to use data integration (also called ‘borrowing information’), combining data from diverse sources and eliciting useful information with a unified view of them. Additionally, it is very important to find risk factors relevant to the disease. Reliable and early risk assessment of a developing infectious disease outbreak allow policymakers to make swift and well-informed decisions that would be needed to ensure epidemic control. Quantifying uncertainty about the final epidemic size is also very important.

Richards growth curve [126], so-called the generalized logistic curve [127], is a popularly used growth curve for population studies in situations where growth is not symmetrical about the point of inflection [128,129]. There are variant reparamerized forms of the Richards curve in the literature [130,131,132,133], and one of the frequently used form is where$I\left(t\right)$ is the cumulative number of infected cases at timet. Here, epidemiological meanings of the parameters,a,b, and c, are the final epidemic size, infection rate, and lag phase of the trajectory, respectively. The parameter$\xi$ is the shape parameter, and there seems no clear epidemiological meaning [134]. Each infection trajectory inFigure 6 can be characterized by the 4-dimensional parameters$(a,b,c,\xi )$ if the Richards curve is used. Due to its flexibility originating from the shape parameter$\xi$, Richards curve has been widely used in epidemiology for real-time prediction of outbreak of diseases, possibly at an early phase of the pandemic when there is no second wave. Examples include SARS [135,136], dengue fever [137,138], pandemic influenza H1N1 [139], and COVID-19 outbreak [30,140].

In the previous subsections, we presented a range of examples in which nonlinear mixed effects models can be exploited. They have their own challenges to solve the problems that are representative of issues many researchers have to deal with in other areas: for example, (1) how to describe a possible nonlinear clearance with a limited number of patients; (2) how to handle an enormously large number of shale oil wells and make a spatial prediction of EUR at a new location; (3) how to describe the dynamic of the financial parameters over time; and (4) how to integrate data from different sources to produce more accurate forecast on the epidemic size.

An emerging issue accompanied by these problems, requested from researchers, government agencies, domain experts, etc., is how to quantify the uncertainty associated with parameter estimation and prediction. Although the traditional nonlinear mixed effects models, based on the maximum likelihood method, can provide confidence intervals and statistical tests, calculations of those generally involve approximations that are most accurate for large sample sizes, as discussed inSection 2.1. On the other hand, in the Bayesian approach—in which the prior automatically imposes the parameter constraints—inferences about parameter values based on the posterior distribution usually require integration rather than maximization, and no further approximation is involved. For that reason, the Bayesian approach is often suggested as a viable alternative to the frequentist approach to solving the problems.

We now formulate these problems as a statistical problem. First, we summarize common features of the dataset for the analysis.

The subject of sampling units considered in the statistical analysis is quite comprehensive. We have seen that it can be a patient, a shale oil well, a particular date, and a country. As the unique identifier, we assign the indexi to each individual. By denotingN as the number of individuals (i.e., the sample size), the indexi will take an integer from 1 toN. The sample sizeN available for the data analysis substantially varies across different industrial problems as well as subfields within the same industry. For example, the number of shale oil wells on Eagle Ford Shale Play can be as large as 6000 [31]. As for the pharmaceutical industry, in phase I cancer clinical trials, the number of cancer patientsN may be strictly confined to 25 [141], but for phase III trials for non-oncology drug studies,N can be as large as 2000 [142].

Here the term ‘time’ is meant in the broadest sense. It can be a calendar time, a nominal time, a time after some event (e.g., the time after dose fromFigure 3 and the time after well completion fromFigure 4), or a time to some event (e.g., the time to maturity fromFigure 5). Essentially, time can be defined as a physical quantity that can be indexed with consecutive integers to produce a temporal record. Another important characteristic of the time is that each subject may have different time points where observations are measured. In this article, we use$t_{ij}$ to represent the time point, where the integer$j=1,2,\cdots ,M_i$ indexes the time point from the earliest to the last observations. Thus,$M_i$ represents the number of repeated observations for thei-th individual. When$M_i$ is relatively small (or large), we say the repeated measures are sparsely observed (or densely observed). For example, the theophylline and yield curve data shown inFigure 3 andFigure 5 are sparse data, while the oil production and COVID-19 data shown inFigure 4 andFigure 6 are dense data.

As for the repeated measures,$y_{ij}$ denotes the continuous response of thei-th subject at the time point$t_{ij}$. We assume that$y_{ij}$ has been already pre-processed so that it is ready to be used for statistical modeling. For most applications, it may be necessary first to transform the data into some new representation before training the model. For example, as seen fromFigure 4, oil productions vary substantially across different wells. For that reason, the authors [31] take a logarithm on the productions to derive the response$y_{ij}$, followed by appropriate statistical modeling on the log-scale. To some extent, data pre-processing may enhance the performance of the model.

Suppose that researchers collectedP number of covariates at the baseline from each subjecti ($i=1,\cdots ,N$). Here, the baseline refers to the time point$t_{i1}$ (or possibly right before the time point$t_{i1}$), where at the first response$y_{i1}$ has not been observed yet. Let$x_{ib}$ denote theb-th covariate of thei-th subject ($b=1,\cdots ,P$). In general, there are two types of covariates: time-invariant and time-varying covariates. This article mainly concerns the former type. As similar toN, the number of covariatesP substantially varies across industries and specific problems. For instance, in pharmacogenetics analysis, the number of protein-coding genesP would be around 20,000 [143]. In the oil and gas industry, if we consider most of the covariates obtained from the well completion procedure,P could be at least 100 [31].

In conclusion, the dataset for the statistical analysis can be represented by the collection of theN triplets${\left\{({\mathbf{y}}_i,{\mathbf{t}}_i,{\mathbf{x}}_i)\right\}}_{i=1}^N$. Here, for each subjecti ($i=1,\cdots ,N$), we formulated two$M_i$-dimensional vectors${\mathbf{y}}_i={(y_{i1},\cdots ,y_{ij},\cdots ,y_{i{M}_i})}^{\top }$ and${\mathbf{t}}_i={(t_{i1},\cdots ,t_{ij},\cdots ,t_{i{M}_i})}^{\top }$, and aP-dimensional vector${\mathbf{x}}_i={(x_{i1},\cdots ,x_{ib},\cdots ,x_{iP})}^{\top }$.

The data structures described up to this point are commonly encountered in longitudinal data studies [7]. Essentially, the feature of dataset motivating the use of nonlinear mixed effects models is that, for each subjecti, the response vector${\mathbf{y}}_i$ displays some nonlinear tendency over time${\mathbf{t}}_i$, as seen inFigure 3,Figure 4,Figure 5 andFigure 6. To explain this nonlinearity, a researcher needs to theorize some nonlinear function, denoted asf, such as one compartment, Duong, Nelson-Siegel, and Richards models, depending on the contexts. The construction of such functions relies on human modelers’ abstraction of data into a suitable dynamical system, which is often represented by a differential equation. Such a differential equation has a finite number of parameters that control the dynamic of the solution of the system, the understanding of which is vital for causal inference for the nature of the system by associating with covariates${\mathbf{x}}_i$.

Figure 7 displays a pictorial description about how a PK modeler would see the theophylline concentration trajectory from the modeling perspective, where she theorized that one compartment model (1) would be suitable to describe the trajectories${\mathbf{y}}_i$ over time${\mathbf{t}}_i$ for each subjecti ($i=1,\cdots ,N$). Then the 10-dimensional vector${\mathbf{y}}_i$ is summarized by a 4-dimensional PK parameter vector$(F_i,k_{ai},V_i,C{l}_i)$; the dimension reduction is intrinsically embedded in this process. As each of the parameters$F_i$,$k_{ai}$,$V_i$, and $C{l}_i$ has an important clinical meaning, it is very natural to ask how they are related withP covariates${\mathbf{x}}_i$ to induce a causal relationship. For the purpose of modeling, it may be necessary to transform the original PK parameters$(F_i,k_{ai},V_i,C{l}_i)\in [0,1]\times {(0,\infty )}^3$ to model parameters$({\theta }_{1i},{\theta }_{2i},{\theta }_{3i},{\theta }_{4i})\in {\mathbb{R}}^K$ ($K=4$) so that elements${\theta }_{li}$ ($l=1,\cdots ,K$) are supported on the real number by taking transformations${\theta }_{1i}=log\{F_i/(1-F_i)\}$,${\theta }_{2i}=log{k}_{ai}$,${\theta }_{3i}=log{V}_i$, and ${\theta }_{4i}=logC{l}_i$. As these transformations were taken only for the modeling purpose, interpretations on the PK parameter for the PK report should be carried out after transforming back to the original scale.

Assume that we have dataset for a statistical analysis${\left\{({\mathbf{y}}_i,{\mathbf{t}}_i,{\mathbf{x}}_i)\right\}}_{i=1}^N$ fromN subjects, as explained inSection 3.6. We consider a basic version of the model here. Extensions are discussed inSection 7.3 and Section 9. The usual Bayesian nonlinear hierarchical model may then be written as a three-stage hierarchical model as follows:

Directed asymmetric graphical (DAG) model representation of the basic model (2)–(4) is depicted inFigure 8. Following the grammar of the graphical model (Chapter 8 of [144]), the circled variables indicate stochastic variables, while the observed ones are additionally colored in gray. Non-stochastic quantities are uncircled. The arrows indicate the conditional dependency between the variables.

We will often wish to write the hierarchy (2)–(4) for thei-th individual’s entire response vector and represent it with an equivalent vector-form. This turns out to be useful to develop relevant computational algorithms. We first introduce a$K\times N$ dimensional matrix frequently used throughout this article:

The matrix (5) is referred to asmodel matrix because it comprises of scalar model parameters${\left\{{\theta }_{li}\right\}}_{l=1,i=1}^{K,N}$ from all subjects. Indeed, most of the computational techniques either via frequentist or Bayesian setting in the literature have been developed to overcome an obstacle of a nonlinear association of the model matrix$\mathrm{\Theta }$ into the mean functionf.

In (5), the subject indexi is stacked column-wisely, while model parameter indexl is stacked row-wisely, different from the conventional way adopted in most statistics. The column indexing for the subjects (i.e., stacking individual-based vector horizontally) shown in (5) is often adopted in modern computation theory of deep learning [145], and one of the main advantages of using this indexing is that it may give some pedagogical insights on the use of vectorization toward the entries${\left\{{\theta }_{li}\right\}}_{l=1,i=1}^{K,N}$ to exploit parallel computations, stochastic updating, etc., in optimization or sampling techniques.

The model matrix$\mathrm{\Theta }$ (5) can be re-expressed as$\mathrm{\Theta }=[{\mathbf{\theta }}^1\cdots {\mathbf{\theta }}^i\cdots {\mathbf{\theta }}^N]\in {\mathbb{R}}^{K\times N}$, obtained by stacking the individual model parameter vector in Stage 1 (2). Alternatively, we can represent the matrix with$\mathrm{\Theta }={[{\mathbf{\theta }}_1\cdots {\mathbf{\theta }}_l\cdots {\mathbf{\theta }}_K]}^{\top }\in {\mathbb{R}}^{K\times N}$ by defining aN-dimensional vector corresponding thel-th model parameter across all subjects${\mathbf{\theta }}_l={({\theta }_{l1},{\theta }_{l2},\cdots ,{\theta }_{lN})}^{\top }\in {\mathbb{R}}^N$ ($l=1,\cdots ,K$). Former and latter indexing method are referred to asi-indexing andl-indexing, respectively.

We are now in a position to re-write the hierarchy (2)–(4) using the vector notations:

To summarize, we derived two equivalent vectorized formulations representing the basic model (2)–(4) according to how the model matrix$\mathrm{\Theta }$ (5) is vectorized:

Figure 9 displays the DAG representations of the two vector forms of the basic model. In vector-form (a),K latent nodes${\left\{{\mathbf{\theta }}_l\right\}}_{l=1}^K$ are fully connected toward theN response vectors${\left\{{\mathbf{y}}_i\right\}}_{i=1}^N$. On the other hand, in vector-form (b),N latent nodes${\left\{{\mathbf{\theta }}^i\right\}}_{i=1}^N$ are bijectively connected to theN response vectors${\left\{{\mathbf{y}}_i\right\}}_{i=1}^N$ for each subjecti. These two ways of looking at the framework of the Bayesian nonlinear hierarchical models complement each other for a more proper understanding of the modeling framework and provide modelers with a statistical insight. For example, vector-form (a) is useful to understand some mathematics underlying the linear regression usingP regressors, while vector-form (b) makes it easy to comprehend the role of the covariance matrix$\mathrm{\Omega }$.

In this section, we investigate a likelihood function based on the basic model (2)–(4). As illustrated in Bayesian workflow inSection 2.2, the likelihood theory is fundamental of Bayesian inference (seeFigure 2). Therefore, it is worth spending time to re-study the formulation of the likelihood function. Here, one caveat is that, due to the hierarchical nature of the nonlinear mixed effects model, a notion of the likelihood function depends on what part of the model specification is considered to be part of the likelihood, and what is not. Most papers directly consider the marginal likelihood that will be discussed inSection 5.4. In this paper, before marching there, we study other two formulations of the likelihood inSection 5.2 andSection 5.3 to get some pedagogical insights. We will briefly discuss popularly used frequentist computing strategies inSection 5.4.

As in most of the statistical models, a natural starting point for inference is maximum likelihood estimation. We start with considering only Stage 1 from the basic model inSection 4.1 and ignore Stages 2 and 3 for now. Then the likelihood function for thei-th subject is

Therefore, the likelihood function based on theN subjects${\mathbf{y}}_{1:N}={\left\{{\mathbf{y}}_i\right\}}_{i=1}^N$ is

Now, we maximize the likelihood (10) with respect to the model matrix$\mathrm{\Theta }$ (5) given${\sigma }^2$ fixed: where${\parallel ·\parallel }_2$ is the Euclidean norm.

The estimator$\widehat{\mathrm{\Theta }}=[{\widehat{\mathbf{\theta }}}^1\cdots {\widehat{\mathbf{\theta }}}^i\cdots {\widehat{\mathbf{\theta }}}^N]\in {\mathbb{R}}^{K\times N}$ of the model matrix$\mathrm{\Theta }$ (5) can be obtained by various optimization techniques such as Newton-Raphson method or Gradient descent method [146]. Noting from the summation acrossi in (11),N estimators${\left\{{\widehat{\mathbf{\theta }}}^i\right\}}_{i=1}^N$ are independent. We can obtain an estimator of the variance${\sigma }^2$ by plugging$\widehat{\mathrm{\Theta }}$ into the likelihood (10), and then maximize with respect to the${\sigma }^2$. To investigate a denoised temporal tendency for the trajectory${\mathbf{y}}_i$, we can simply plug${\widehat{\mathbf{\theta }}}^i$ into the function$f(t_{ij};{\mathbf{\theta }}^i)$$(j=1,\cdots ,M_i)$. To see a future pattern, we can extrapolate the function by extending the time index beyond the last time point$t_{i{M}_i}$. Eventually, the illustrated approach is based on traditional least squares estimation.

Unfortunately, there are three major drawbacks in this approach. First, it forfeits the opportunity to use ‘information borrowing’ [30] to improve a predictive accuracy due to the ignorance of Stage 2. What happens in Stage 2 (3) is to borrow strength acrossN individuals to produce a better estimator for$\mathrm{\Theta }$ (5) than an estimator simply based on individual data. A similar issue can be found in the Clemente problem from [11] where the James–Stein estimator [147] predicts better than an individual hitter-based estimator. Another example applied to epidemic data can be found in [30]. Second, it is not well-aligned with the generic motivation to use the mixed effects models, whose primary purpose is to understand “typical” values for the model parameters inf, representing whole subjects, which should be addressed by making an inference about the parameters$\mathbf{\alpha }$,$\mathit{B}$, and $\mathrm{\Omega }$. Third, it only produces point estimates for the parameters, failing to describe the underlying uncertainty.

A remedy of the first two drawbacks is the consideration of Stages 1 (2) and 2 (3) hierarchically in a single model, leading to a frequentist version of nonlinear mixed effects models, which will be discussed inSection 5.3 andSection 5.4. To describe relevant uncertainty within the frequentist framework one may have to resort to bootstrap methods [99] or use a large-sample theory. Uncertainty quantification in frequentist analysis often needs to be done step-wisely. Another solution resolving all the three drawbacks at once is to incorporate Stage 1 (2), 2 (3), and 3 (4) into a single model in a fully Bayesian way, resulting in a Bayesian version of nonlinear mixed effects models, which is the main topic in this paper.

A likehood function based on vector-form (a) is derived. More specifically, we examine the frequentist setting where the assumptions of Stage 1—(6) and Stage 2—(7) are considered, while the parameters introduced in Stage 3—(9) are fixed (i.e., no prior assumptions).

The individual model on Stage 1 (6) yields a conditional density$p\left({\mathbf{y}}_i\right|{\mathbf{\theta }}^i,{\sigma }^2)={\mathcal{N}}_{M_i}\left({\mathbf{y}}_i\right|{\mathit{f}}_i({\mathbf{t}}_i,{\mathbf{\theta }}^i),{\sigma }^{2}I)$ for each subject$i=1,\cdots ,N$. Under the population assumption on Stage 2 (7), we have the density$p\left({\mathbf{\theta }}_l\right|{\alpha }_l,{\mathbf{\beta }}_l,{\omega }_l^2)={\mathcal{N}}_N\left({\mathbf{\theta }}_l\right|\mathbf{1}{\alpha }_l+\mathbf{X}{\mathbf{\beta }}_l,{\omega }_l^2\mathbf{I})$ for each model parameter index$l=1,\cdots ,K$. The joint density of$({\mathbf{y}}_{1:N},{\mathbf{\theta }}_{1:K})$ given parameters${\sigma }^2,\mathbf{\alpha },\mathit{B}$ and$\mathrm{\Omega }$ is a product-form distribution: where${\mathbf{y}}_{1:N}=\{{\mathbf{y}}_1,{\mathbf{y}}_2,\cdots ,{\mathbf{y}}_N\}$ and${\mathbf{\theta }}_{1:K}=\{{\mathbf{\theta }}_1,{\mathbf{\theta }}_2,\cdots ,{\mathbf{\theta }}_K\}$. Now, the next step is to integrate out the latent model parameters${\mathbf{\theta }}_{1:K}$ (or equivalently, the model matrix$\mathrm{\Theta }$ (5)) from the density above to get a likelihood for the$({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$:

In most cases, the integral (12) is not tractable due to the non-linearity of the function${\mathit{f}}_i({\mathbf{t}}_i,{\mathbf{\theta }}^i)$ with respect to the${\mathbf{\theta }}^i$. Although it may be possible to use numerical techniques for the evaluation of the integral (12), this might require enormous computational effort, which is not really appreciated in the literature due to the high-dimensionality of the integral involving the$KN$ dimensional model parameters${\mathbf{\theta }}_{1:K}$.

A likehood function based on vector-form (b) adoptingi-indexing is derived here. Similar toSection 5.3, we preserve the assumption of Stage 1—(6) and Stage 2${}^{\prime }$—(8), but work with fixing the parameters in Stage 3—(9). In these specifications, for each index$i=1,\cdots ,N$, the individual model on Stage 1 (6) and population model on Stage 2 (8) lead to densities$p\left({\mathbf{y}}_i\right|{\mathbf{\theta }}^i,{\sigma }^2)={\mathcal{N}}_{M_i}\left({\mathbf{y}}_i\right|{\mathit{f}}_i({\mathbf{t}}_i,{\mathbf{\theta }}^i),{\sigma }^{2}I)$ and$p\left({\mathbf{\theta }}^i\right|\mathbf{\alpha },\mathit{B},\mathrm{\Omega })={\mathcal{N}}_K\left({\mathbf{\theta }}^i\right|\mathbf{\alpha }+\mathit{B}{\mathbf{x}}_i,\mathrm{\Omega })$, respectively. Thus, the joint density of$({\mathbf{y}}_i,{\mathbf{\theta }}^i)$ given parameters${\sigma }^2,\mathbf{\alpha },\mathit{B}$ and$\mathrm{\Omega }$ is

Given the parameters$({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$, the ordered pairs in the collection${\left\{({\mathbf{y}}_i,{\mathbf{\theta }}^i)\right\}}_{i=1}^N$ are conditionally independent across individuals. Therefore, a likelihood for the$({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$ is based on the marginal density of${\mathbf{y}}_{1:N}=\{{\mathbf{y}}_1,{\mathbf{y}}_2,\cdots ,{\mathbf{y}}_N\}$: where the last equality is derived by using the change of variable (8). The last expression (14) is a standard mathematical formulation that many frequentist computing strategies are constructed with: see Equation (3.2) from [14]. Essentially, the part that makes the MLE computation complicated is the mean vector${\mathit{f}}_i({\mathbf{t}}_i,\mathbf{\alpha }+\mathit{B}{\mathbf{x}}_i+{\mathbf{\eta }}^i)\in {\mathbb{R}}^{M_i}$.

As the model parameter${\mathbf{\theta }}^i$ in (13) (or similarly,${\mathbf{\eta }}^i$ in (14) which is often called random effect in the frequentist framework) participates to the functionf in a non-linear fashion, the integral generally cannot be obtained in a closed-form. Benefiting from a conditional independence [148], dimensionality of theN integrals (13) is much lower than that of the integral (12) based on vector-form (a). Analytically, the likelihood functions of the basic model (2)–(4) based on vector-form (a) (12) and vector-form (b) (13) may be equivalent. That being said, minimization of the two functions with respect to the parameters$({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$ yields the same solution,$(\widehat{{\sigma }^2},\widehat{\mathbf{\alpha }},\widehat{\mathit{B}},\widehat{\mathrm{\Omega }})$, so-called maximum likelihood estimators (MLE).

We shall briefly discuss MLE computations. One approach would be to perform a multivariate numerical integration (e.g., Gauss-Hermite quadrature [149]) to each of theN integrals (13), and then obtain the MLE by maximizing the product of theN numerical integrals with respect to the parameters$({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$ [150]. This approach turns out to be computationally expensive and may have poor converge properties due to the following two reasons [151]. First, the numerical integration necessitates increasingly expensive iterative procedures within a MLE algorithm as the correlation of the model parameters (or equivalently, random effects) increases. Second, convergence property may be highly deteriorated when the number of model parametersK is large (i.e., high-dimensional integral) and the number of sampling times$M_i$ is small (i.e., sparse data) due to the ‘curse of dimensionality’ [152].

A class of common approaches for the MLE computations is based on analytical approximation to each of theN integrals (14) [13,36,153,154,155], and some of them have been successfully adopted to industrial software likeNONMEM [47,101] andSAS [46]. Here, we illustrate a key idea of the first-order method attributed to [36]. Let us define a mapping${\mathbf{g}}_i\left({\mathbf{\eta }}^i\right)={\mathit{f}}_i({\mathbf{t}}_i,\mathbf{\alpha }+\mathit{B}{\mathbf{x}}_i+{\mathbf{\eta }}^i):A\subset {\mathbb{R}}^K\to {\mathbb{R}}^{M_i}$ for each subjecti, whereA is an open set with$\mathbf{0}\in A$. Suppose that${\mathbf{g}}_i\left({\mathbf{\eta }}^i\right)$ is smooth on the setA: then, by Taylor’s theorem (page 375 of [156]), we have the best linear approximation of the mapping${\mathbf{g}}_i\left({\mathbf{\eta }}^i\right)$ at the origin$\mathbf{0}$ given by${\mathbf{g}}_i\left({\mathbf{\eta }}^i\right)\approx {\mathbf{g}}_i\left(\mathbf{0}\right)+\mathbf{D}{\mathbf{g}}_i\left(\mathbf{0}\right){\mathbf{\eta }}^i$, where$\mathbf{D}{\mathbf{g}}_i\left(\mathbf{0}\right)\in {\mathbb{R}}^{M_i\times K}$ is the Jacobian matrix of${\mathbf{g}}_i\left({\mathbf{\eta }}^i\right)$ at$\mathbf{0}$. Now, we shall replace the function${\mathbf{g}}_i\left({\mathbf{\eta }}^i\right)={\mathbf{f}}_i({\mathbf{t}}_i,\mathbf{\alpha }+\mathit{B}{\mathbf{x}}_i+{\mathbf{\eta }}^i)$ in integral (14) with the resulting approximation${\mathbf{g}}_i\left(\mathbf{0}\right)+\mathbf{D}{\mathbf{g}}_i\left(\mathbf{0}\right){\mathbf{\eta }}^i$ for eachi ($i=1,\cdots ,N$), leading to a closed-form expression

To summarize, a linearization was used to convert the nonlinear mixed effects model to a linear mixed effects model, in some sense, equivalent to the Lindley–Smith form [157]. This enables us to integrate out the random vector${\mathbf{\eta }}^i$ from theN integrals (15), deriving a marginal likelihood (15) to approximate the exact marginal likelihood (13). MLE$(\widehat{{\sigma }^2},\widehat{\mathbf{\alpha }},\widehat{\mathit{B}},\widehat{\mathrm{\Omega }})$ can be obtained by jointly maximizing$\tilde{\mathcal{L}}({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega }|{\mathbf{y}}_{1:N})$ (15) assuming the approximation is exact.

Another way to compute the MLE is through the use of expectation-maximization (EM) algorithm [158]. Borrowing terms from EM updating process [159],${\mathbf{y}}_i$,${\mathbf{\theta }}^i$,$({\mathbf{y}}_i,{\mathbf{\theta }}^i)$,$({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$ and$p({\mathbf{y}}_i,{\mathbf{\theta }}^i|{\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$ (i.e., the integrand in (13)) can be viewed as observable incomplete data, missing data, complete data, unknown parameters, and density of complete data, respectively, for thei-th subject. The goal is to maximize the exact marginal likelihood$\mathcal{L}({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega }|{\mathbf{y}}_{1:N})$ (13) by iterating E-step and M-step, leading to the MLE$(\widehat{{\sigma }^2},\widehat{\mathbf{\alpha }},\widehat{\mathit{B}},\widehat{\mathrm{\Omega }})$. The E-step computes a conditional expected log-likelihood of$({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$ based on the hierarchy (2)–(3), followed by the M-step that maximizes the function with respect to$({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$. The nonlinearity associated with the model matrix$\mathrm{\Theta }$ (5) makes the E-step intractable. As a remedy, variant versions of the EM algorithm are proposed; see [160,161,162,163] for a technical detail applied to a hierarchy similar to the basic model. Among them, the scheme of stochastic approximation EM algorithm proposed by [38], splitting the E-step into two steps, namely a simulation step and a stochastic approximation step, is widely used in many applications for its numerical stability, fast computation, and theoretical soundness [164,165], which has been successfully deployed as industrial software includingMonolix [48] as well as open source software such asR packagenlmixr [18].

We briefly overview two contrasting workflows of Bayesian and frequentist approaches for nonlinear mixed effects models before moving to a technical detail. Both settings allow the randomness in the model matrix$\mathrm{\Theta }$ (5), but then, they diverge when it comes to how parameters$({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$ are treated. Bayesians treat$({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$ as random, while frequentists regard it as fixed. To conceptualize a subtlety arising from this difference, let us recap frequentist computing strategies discussed inSection 5.4. There, the model matrix$\mathrm{\Theta }$ was eventually integrated out from the joint density of$({\mathbf{y}}_{1:N},\mathrm{\Theta })$, either approximately or exactly, to derive a marginal likelihood of$({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$ from which the MLE$(\widehat{{\sigma }^2},\widehat{\mathbf{\alpha }},\widehat{\mathit{B}},\widehat{\mathrm{\Omega }})$ is computed via various optimization methods. After that, frequentists apply standard Bayesian formulas, such as posterior density, posterior mean, and so on, to estimate$\mathrm{\Theta }$ [7].

In contrast, to drive the Bayesian engine, one would need an appropriate prior$\pi ({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$. After that, the entire collection of parameters$(\mathrm{\Theta },{\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$ will be updated through the Bayes’ theorem post observing the data${\mathbf{y}}_{1:N}$, leading to the posterior density$\pi (\mathrm{\Theta },{\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega }|{\mathbf{y}}_{1:N})$ [166] (SeeFigure 2). The essence of the Bayesian viewpoint is that there is no logical distinction between$\mathrm{\Theta }$ and$({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$, which are associated with the random and fixed effects, respectively, from the frequentist perspective. In Bayesian framework, both$\mathrm{\Theta }$ and$({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$ are random quantities. It is important to point out that the likelihood principle is naturally incorporated in the Bayes’ theorem [167]. Clearly, modern data complications such as enormous volume, large dimensionality, and multi-level structures may necessitate a sophistication on the prior specifications.

We are now in a position to describe the Bayesian analysis for the basic model (2)–(4), that assumed independence, a priori, for each parameter blocks${\sigma }^2$,$\mathbf{\alpha },$$\mathit{B}$, and $\mathrm{\Omega }$,$\pi ({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })=\pi \left({\sigma }^2\right)·\pi \left(\mathbf{\alpha }\right)·\pi \left(\mathit{B}\right)·\pi \left(\mathrm{\Omega }\right)$. (Our logic below can be generalized to a more complex prior setting.) As discussed inSection 4, it is the discretion of the modeler of how she would treat the model matrix$\mathrm{\Theta }$ (5) withl-indexing${\mathbf{\theta }}_{1:K}={\left\{{\mathbf{\theta }}_l\right\}}_{l=1}^K$ ori-indexing${\mathbf{\theta }}^{1:N}={\left\{{\mathbf{\theta }}^i\right\}}_{i=1}^N$, leading to vector forms (a) and (b), respectively. For the sake of readability, we illustrate the Bayesian inference by using the vector-form (a), but we will sometimes use the vector-form (b) when this seems more understandable.

A central task in the application of the Bayesian nonlinear mixed effects models is to evaluate the posterior density, or indeed to compute expectation with respect to the density: where the last equation can be detailed as follows

From a Bayesian perspective, all inferential problems regarding the parameter$(\mathrm{\Theta },{\sigma }^2,$$\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$ may be addressed in terms of the posterior distribution$\pi (\mathrm{\Theta },{\sigma }^2,$$\mathbf{\alpha },\mathit{B},\mathrm{\Omega }|{\mathbf{y}}_{1:N})$ (16). Unfortunately, for almost all problems, the distribution is intractable. In such situations, we need to resort to approximation techniques, and these fall broadly into two classes, according to whether they rely on stochastic [41,45,168,169] or deterministic [170,171,172,173] approximations. See [174,175] for review papers of these techniques. In this article, we mainly focused on the stochastic approximation. The basic idea behind the methodology is to construct a Markov chain whose stationary distribution is the posterior distribution$\pi (\mathrm{\Theta },{\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega }|{\mathbf{y}}_{1:N})$ (16).

We resort to MCMC technique [174] to sample from the full joint density$\pi (\mathrm{\Theta },{\sigma }^2,\mathbf{\alpha },\mathit{B},$$\mathrm{\Omega }|{\mathbf{y}}_{1:N})$ (16). Among many MCMC techniques, we use the Gibbs sampling algorithm [39,168] to exploit the conditional independence [148] induced by the hierarchical formulation. A generic Gibbs sampler would cycle in turn through each of the conditional distributions for the parameter blocks$\mathrm{\Theta },{\sigma }^2,\mathbf{\alpha },\mathit{B}$, and $\mathrm{\Omega }$ as follows:

Step 1. Sample$\mathrm{\Theta }$ from its full conditional distribution

Step 2. Sample${\sigma }^2$ from its full conditional distribution

Step 3. Sample$\mathbf{\alpha }$ from its full conditional distribution

Step 4. Sample$\mathit{B}$ from its full conditional distribution

Step 5. Sample$\mathrm{\Omega }$ from its full conditional distribution

Sampling the model matrix$\mathrm{\Theta }$ (5) atStep 1 is independent of the choice of the priors, which we discuss shortly. On the other hand, sampling the parameters$({\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega })$ atSteps 2,3,4, and 5 depends on the prior choices for the parameters in Stage 3 (9); we discuss this topic inSection 7.

One of the most computer-intensive steps to implement the Gibbs sampler inSection 6.2 isStep 1 to sample the model matrix$\mathrm{\Theta }\in {\mathbb{R}}^{K\times N}$ (5), or equivalently its entries${\left\{{\theta }_{li}\right\}}_{l=1,i=1}^{K,N}$, from the full conditional distribution$\pi \left(\mathrm{\Theta }\right|{\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega },{\mathbf{y}}_{1:N})$ (17). Clearly, the nonlinear participation of the model parameters to the functionf makes the conditional distribution intractable, hence, non-conjugate sampling is unavoidable, which may suffer from a slow convergence. At the same times, due to the Markovian nature of the Gibbs algorithm, it is difficult to parallelize the whole steps of the Gibbs sampler, which creates difficulties in slower languages likeR [176]. Nevertheless, the increasing number of parallel cores that are available at a very low price drives more and more interest in ‘parallel sampling algorithms’ that can benefit from the available parallel processing units on computers [177,178].

We suggest a framework of parallel computations to efficiently update the model matrix$\mathrm{\Theta }\in {\mathbb{R}}^{K\times N}$. This framework can be particularly appreciated under the setting of Bayesian nonlinear mixed effects models when the number of subjectsN is a lot larger than the number of model parametersK ($N\gg K$).

The first version of parallel sampling algorithms is based on scalar updating. For the derivation, we start with analyzing the full conditional posterior distribution of a single element${\theta }_{li}$ ($l=1,\cdots ,K;i=1,\cdots ,N$): where the notation${\theta }_{-li}$ represents the all the entries of$\mathrm{\Theta }$ except for${\theta }_{li}$, that is,${\theta }_{-li}={\left\{{\theta }_{li}\right\}}_{l=1,i=1}^{K,N}-\left\{{\theta }_{li}\right\}$. Here, we used a conventional notation in Bayesian computation: ‘$\pi \left({\theta }_{li}\right|-)$’ indicates the density$\pi \left({\theta }_{li}\right|{\theta }_{-li},{\sigma }^2,\mathbf{\alpha },\mathit{B},\mathrm{\Omega },{\mathbf{y}}_{1:N})$, where the notation ‘−’ in$\pi \left({\theta }_{li}\right|-)$ implies all the parameters except for the${\theta }_{li}$ in the basic model (2)–(4) along withN observations.