Nonlinear mixedeffects model
Part of a series on 
Regression analysis 

Models 
Estimation 
Background 

Nonlinear mixedeffects models constitute a class of statistical models generalizing linear mixedeffects models. Like linear mixedeffects models, they are particularly useful in settings where there are multiple measurements within the same statistical units or when there are dependencies between measurements on related statistical units. Nonlinear mixedeffects models are applied in many fields including medicine, public health, pharmacology, and ecology.^{[1]}^{[2]}
Definition
While any statistical model containing both fixed effects and random effects is an example of a nonlinear mixedeffects model, the most commonly used models are members of the class of nonlinear mixedeffects models for repeated measures^{[1]}
 [math]\displaystyle{ {y}_{ij} = f(\phi_{ij},{v}_{ij}) + \epsilon_{ij},\quad i =1,\ldots, M, \, j = 1,\ldots, n_i }[/math]
where
 [math]\displaystyle{ M }[/math] is the number of groups/subjects,
 [math]\displaystyle{ n_i }[/math] is the number of observations for the [math]\displaystyle{ i }[/math]th group/subject,
 [math]\displaystyle{ f }[/math] is a realvalued differentiable function of a groupspecific parameter vector [math]\displaystyle{ \theta_{ij} }[/math] and a covariate vector [math]\displaystyle{ v_{ij} }[/math],
 [math]\displaystyle{ \phi_{ij} }[/math] is modeled as a linear mixedeffects model [math]\displaystyle{ \phi_{ij}= \boldsymbol{A}_{ij}\beta + \boldsymbol{B}_{ij} \boldsymbol{b}_{i}, }[/math] where [math]\displaystyle{ \beta }[/math] is a vector of fixed effects and [math]\displaystyle{ \boldsymbol{b}_{i} }[/math] is a vector of random effects associated with group [math]\displaystyle{ i }[/math], and
 [math]\displaystyle{ \epsilon_{ij} }[/math] is a random variable describing additive noise.
Estimation
When the model is only nonlinear in fixed effects and the random effects are Gaussian, maximumlikelihood estimation can be done using nonlinear least squares methods, although asymptotic properties of estimators and test statistics may differ from the conventional general linear model. In the more general setting, there exist several methods for doing maximumlikelihood estimation or maximum a posteriori estimation in certain classes of nonlinear mixedeffects models – typically under the assumption of normally distributed random variables. A popular approach is the LindstromBates algorithm^{[3]} which relies on iteratively optimizing a nonlinear problem, locally linearizing the model around this optimum and then employing conventional methods from linear mixedeffects models to do maximum likelihood estimation. Stochastic approximation of the expectationmaximization algorithm gives an alternative approach for doing maximumlikelihood estimation.^{[4]}
Applications
Example: Disease progression modeling
Nonlinear mixedeffects models have been used for modeling progression of disease.^{[5]} In progressive disease, the temporal patterns of progression on outcome variables may follow a nonlinear temporal shape that is similar between patients. However, the stage of disease of an individual may not be known or only partially known from what can be measured. Therefore, a latent time variable that describe individual disease stage (i.e. where the patient is along the nonlinear mean curve) can be included in the model.
Example: Modeling cognitive decline in Alzheimer's disease
Alzheimer's disease is characterized by a progressive cognitive deterioration. However, patients may differ widely in cognitive ability and reserve, so cognitive testing at a single time point can often only be used to coarsely group individuals in different stages of disease. Now suppose we have a set of longitudinal cognitive data [math]\displaystyle{ (y_{i1}, \ldots, y_{in_i}) }[/math] from [math]\displaystyle{ i=1,\ldots,M }[/math] individuals that are each categorized as having either normal cognition (CN), mild cognitive impairment (MCI) or dementia (DEM) at the baseline visit (time [math]\displaystyle{ t_{i1} =0 }[/math] corresponding to measurement [math]\displaystyle{ y_{i1} }[/math]). These longitudinal trajectories can be modeled using a nonlinear mixed effects model that allows differences in disease state based on baseline categorization:
 [math]\displaystyle{ {y}_{ij} = f_{\tilde\beta}(t_{ij} + A^{MCI}_i \beta^{MCI} + A^{DEM}_i \beta^{DEM} + b_i) + \epsilon_{ij},\quad i =1,\ldots, M, \, j = 1,\ldots, n_i }[/math]
where
 [math]\displaystyle{ f_{\tilde\beta} }[/math] is a function that models the mean timeprofile of cognitive decline whose shape is determined by the parameters [math]\displaystyle{ \tilde\beta }[/math],
 [math]\displaystyle{ t_{ij} }[/math] represents observation time (e.g. time since baseline in the study),
 [math]\displaystyle{ A^{MCI}_i }[/math] and [math]\displaystyle{ A^{DEM}_i }[/math] are dummy variables that are 1 if individual [math]\displaystyle{ i }[/math] has MCI or dementia at baseline and 0 otherwise,
 [math]\displaystyle{ \beta^{MCI} }[/math] and [math]\displaystyle{ \beta^{DEM} }[/math] are parameters that model the difference in disease progression of the MCI and dementia groups relative to the cognitively normal,
 [math]\displaystyle{ b_{i} }[/math] is the difference in disease stage of individual [math]\displaystyle{ i }[/math] relative to his/her baseline category, and
 [math]\displaystyle{ \epsilon_{ij} }[/math] is a random variable describing additive noise.
An example of such a model with an exponential mean function fitted to longitudinal measurements of the Alzheimer's Disease Assessment ScaleCognitive Subscale (ADASCog) is shown in the box. As shown, the inclusion of fixed effects of baseline categorization (MCI or dementia relative to normal cognition) and the random effect of individual continuous disease stage [math]\displaystyle{ b_{i} }[/math] aligns the trajectories of cognitive deterioration to reveal a common pattern of cognitive decline.
Example: Growth analysis
Growth phenomena often follow nonlinear patters (e.g. logistic growth, exponential growth, and hyperbolic growth). Factors such as nutrient deficiency may both directly affect the measured outcome (e.g. organisms with lack of nutrients end up smaller), but possibly also timing (e.g. organisms with lack of nutrients grow at a slower pace). If a model fails to account for the differences in timing, the estimated populationlevel curves may smooth out finer details due to lack of synchronization between organisms. Nonlinear mixedeffects models enable simultaneous modeling of individual differences in growth outcomes and timing.
Example: Modeling human height
Models for estimating the mean curves of human height and weight as a function of age and the natural variation around the mean are used to create growth charts. The growth of children can however become desynchronized due to both genetic and environmental factors. For example, age at onset of puberty and its associated height spurt can vary several years between adolescents. Therefore, crosssectional studies may underestimate the magnitude of the pubertal height spurt because age is not synchronized with biological development. The differences in biological development can be modeled using random effects [math]\displaystyle{ \boldsymbol{w}_i }[/math] that describe a mapping of observed age to a latent biological age using a socalled warping function [math]\displaystyle{ v(\cdot, \boldsymbol{w}_i) }[/math]. A simple nonlinear mixedeffects model with this structure is given by
 [math]\displaystyle{ {y}_{ij} = f_{\beta}(v(t_{ij}, \boldsymbol{w}_i)) + \epsilon_{ij},\quad i =1,\ldots, M, \, j = 1,\ldots, n_i }[/math]
where
 [math]\displaystyle{ f_{\beta} }[/math] is a function that represents the height development of a typical child as a function of age. Its shape is determined by the parameters [math]\displaystyle{ \beta }[/math],
 [math]\displaystyle{ t_{ij} }[/math] is the age of child [math]\displaystyle{ i }[/math] corresponding to the height measurement [math]\displaystyle{ y_{ij} }[/math],
 [math]\displaystyle{ v(\cdot, \boldsymbol{w}_i) }[/math] is a warping function that maps age to biological development to synchronize. Its shape is determined by the random effects [math]\displaystyle{ \boldsymbol{w}_i }[/math],
 [math]\displaystyle{ \epsilon_{ij} }[/math] is a random variable describing additive variation (e.g. consistent differences in height between children and measurement noise).
There exists several methods and software packages for fitting such models. The socalled SITAR model^{[7]} can fit such models using warping functions that are affine transformations of time (i.e. additive shifts in biological age and differences in rate of maturation), while the socalled pavpop model^{[6]} can fit models with smoothlyvarying warping functions. An example of the latter is shown in the box.
Example: Population Pharmacokinetic/pharmacodynamic modeling
PK/PD models for describing exposureresponse relationships such as the Emax model can be formulated as nonlinear mixedeffects models.^{[8]} The mixedmodel approach allows modeling of both population level and individual differences in effects that have a nonlinear effect on the observed outcomes, for example the rate at which a compound is being metabolized or distributed in the body.
Example: COVID19 epidemiological modeling
The platform of the nonlinear mixed effect models can be used to describe infection trajectories of subjects and understand some common features shared across the subjects. In epidemiological problems, subjects can be countries, states, or counties, etc. This can be particularly useful in estimating a future trend of the epidemic in an early stage of pendemic where nearly little information is known regarding the disease.^{[9]}
Example: Prediction of oil production curve of shale oil wells at a new location with latent kriging
The eventual success of petroleum development projects relies on a large degree of well construction costs. As for unconventional oil and gas reservoirs, because of very low permeability, and a flow mechanism very different from that of conventional reservoirs, estimates for the well construction cost often contain high levels of uncertainty, and oil companies need to make heavy investment in the drilling and completion phase of the wells. The overall recent commercial success rate of horizontal wells in the United States is known to be 65%, which implies that only 2 out of 3 drilled wells will be commercially successful. For this reason, one of the crucial tasks of petroleum engineers is to quantify the uncertainty associated with oil or gas production from shale reservoirs, and further, to predict an approximated production behavior of a new well at a new location given specific completion data before actual drilling takes place to save a large degree of well construction costs.
The platform of the nonlinear mixed effect models can be extended to consider the spatial association by incorporating the geostatistical processes such as Gaussian process on the second stage of the model as follows:^{[10]}
 [math]\displaystyle{ {y}_{it} = \mu(t;\theta_{1i},\theta_{2i},\theta_{3i}) + \epsilon_{it},\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad i =1,\ldots, N, \, t = 1,\ldots, T_i, }[/math]
[math]\displaystyle{ \theta_{li}=\theta_{l}(s_i) = \alpha_l + \sum_{j=1}^{p}\beta_{lj}x_j + \epsilon_{l}(s_i) + \eta_{l}(s_i), \quad \epsilon_{l}(\cdot) \sim GWN(\sigma_l^2), \quad\quad l=1,2,3, }[/math] [math]\displaystyle{ \eta_{l}(\cdot) \sim GP(0,K_{\gamma_{l}}(\cdot, \cdot)),\quad K_{\gamma_{l}}(s_i,s_j) = \gamma_l^2 \exp (e^{\rho_l} \ s_i  s_j \^2), \quad\quad\quad l=1,2,3, }[/math] [math]\displaystyle{ \beta_{lj}\lambda_{lj},\tau_l,\sigma_l \sim N(0,\sigma_l^2 \tau_l^2 \lambda_{lj}^2 ),\quad \sigma,\lambda_{lj},\tau_l,\sigma_l\sim C^{+}(0,1), \quad\quad\quad\quad\quad\quad\quad l=1,2,3,\, j =1,\cdots, p, }[/math] [math]\displaystyle{ \alpha_l \sim \pi(\alpha)\propto 1, \quad \sigma_l^2 \sim \pi(\sigma^2) \propto 1/\sigma^2, \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad l=1,2,3, }[/math]
where
 [math]\displaystyle{ \mu(t;\theta_{1},\theta_{2},\theta_{3}) }[/math] is a function that models the mean timeprofile of logscaled oil production rate whose shape is determined by the parameters [math]\displaystyle{ (\theta_{1},\theta_{2},\theta_{3}) }[/math]. The function is obtained from taking logarithm to the rate decline curve used in decline curve analysis,
 [math]\displaystyle{ x_i = (x_{i1},\cdots,x_{ip})^{\top} }[/math] represents covariates obtained from the completion process of the hydraulic fracturing and horizontal directional drilling for the [math]\displaystyle{ i }[/math]th well,
 [math]\displaystyle{ s_i = (s_{i1},s_{i2})^{\top} }[/math] represents the spatial location (longitude, latitude) of the [math]\displaystyle{ i }[/math]th well,
 [math]\displaystyle{ \epsilon_{l}(\cdot) }[/math] represents the Gaussian white noise with error variance [math]\displaystyle{ \sigma_l^2 }[/math] (also called the nugget effect),
 [math]\displaystyle{ \eta_{l}(\cdot) }[/math] represents the Gaussian process with Gaussian covariance function [math]\displaystyle{ K_{\gamma_{l}}(\cdot, \cdot) }[/math],
 [math]\displaystyle{ \beta }[/math] represents the horseshoe shrinkage prior.
The Gaussian process regressions used on the latent level (the second stage) eventually produce kriging predictors for the curve parameters [math]\displaystyle{ (\theta_{1i},\theta_{2i},\theta_{3i}), (i=1,\cdots,N), }[/math] that dictate the shape of the mean curve [math]\displaystyle{ \mu(t;\theta_{1},\theta_{2},\theta_{3}) }[/math] on the date level (the first level). As the kriging techniques have been employed in the latent level, this technique is called latent kriging. The right panels show the prediction results of the latent kriging method applied to the two test wells in the Eagle Ford Shale Reservoir of South Texas.
Bayesian nonlinear mixedeffects model
The framework of Bayesian hierarchical modeling is frequently used in diverse applications. Particularly, Bayesian nonlinear mixedeffects models have recently received significant attention. A basic version of the Bayesian nonlinear mixedeffects models is represented as the following threestage:
Stage 1: IndividualLevel Model
[math]\displaystyle{ {y}_{ij} = f(t_{ij};\theta_{1i},\theta_{2i},\ldots,\theta_{li},\ldots,\theta_{Ki} ) + \epsilon_{ij},\quad \epsilon_{ij} \sim N(0, \sigma^2), \quad i =1,\ldots, N, \, j = 1,\ldots, M_i. }[/math]
Stage 2: Population Model
[math]\displaystyle{ \theta_{li}= \alpha_l + \sum_{b=1}^{P}\beta_{lb}x_{ib} + \eta_{li}, \quad \eta_{li} \sim N(0, \omega_l^2), \quad i =1,\ldots, N, \, l=1,\ldots, K. }[/math]
Stage 3: Prior
[math]\displaystyle{ \sigma^2 \sim \pi(\sigma^2),\quad \alpha_l \sim \pi(\alpha_l), \quad (\beta_{l1},\ldots,\beta_{lb},\ldots,\beta_{lP}) \sim \pi(\beta_{l1},\ldots,\beta_{lb},\ldots,\beta_{lP}), \quad \omega_l^2 \sim \pi(\omega_l^2), \quad l=1,\ldots, K. }[/math]
Here, [math]\displaystyle{ y_{ij} }[/math] denotes the continuous response of the [math]\displaystyle{ i }[/math]th subject at the time point [math]\displaystyle{ t_{ij} }[/math], and [math]\displaystyle{ x_{ib} }[/math] is the [math]\displaystyle{ b }[/math]th covariate of the [math]\displaystyle{ i }[/math]th subject. Parameters involved in the model are written in Greek letters. [math]\displaystyle{ f(t ; \theta_{1},\ldots,\theta_{K}) }[/math] is a known function parameterized by the [math]\displaystyle{ K }[/math]dimensional vector [math]\displaystyle{ (\theta_{1},\ldots,\theta_{K}) }[/math]. Typically, [math]\displaystyle{ f }[/math] is a `nonlinear' function and describes the temporal trajectory of individuals. In the model, [math]\displaystyle{ \epsilon_{ij} }[/math] and [math]\displaystyle{ \eta_{li} }[/math] describe withinindividual variability and betweenindividual variability, respectively. If Stage 3: Prior is not considered, then the model reduces to a frequentist nonlinear mixedeffect model.
A central task in the application of the Bayesian nonlinear mixedeffect models is to evaluate the posterior density:
[math]\displaystyle{ \pi(\{\theta_{li}\}_{i=1,l=1}^{N,K},\sigma^2, \{\alpha_l\}_{l=1}^K, \{\beta_{lb}\}_{l=1,b=1}^{K,P},\{\omega_l\}_{l=1}^K  \{y_{ij}\}_{i=1,j=1}^{N,M_i}) }[/math]
[math]\displaystyle{ \propto \pi(\{y_{ij}\}_{i=1,j=1}^{N,M_i}, \{\theta_{li}\}_{i=1,l=1}^{N,K},\sigma^2, \{\alpha_l\}_{l=1}^K, \{\beta_{lb}\}_{l=1,b=1}^{K,P},\{\omega_l\}_{l=1}^K) }[/math]
[math]\displaystyle{ = \underbrace{\pi(\{y_{ij}\}_{i=1,j=1}^{N,M_i} \{\theta_{li}\}_{i=1,l=1}^{N,K},\sigma^2)}_{Stage 1: IndividualLevel Model} \times \underbrace{\pi(\{\theta_{li}\}_{i=1,l=1}^{N,K}\{\alpha_l\}_{l=1}^K, \{\beta_{lb}\}_{l=1,b=1}^{K,P},\{\omega_l\}_{l=1}^K)}_{Stage 2: Population Model} \times \underbrace{p(\sigma^2, \{\alpha_l\}_{l=1}^K, \{\beta_{lb}\}_{l=1,b=1}^{K,P},\{\omega_l\}_{l=1}^K)}_{Stage 3: Prior} }[/math]
The panel on the right displays Bayesian research cycle using Bayesian nonlinear mixedeffects model.^{[12]} A research cycle using the Bayesian nonlinear mixedeffects model comprises two steps: (a) standard research cycle and (b) Bayesianspecific workflow. Standard research cycle involves literature review, defining a problem and specifying the research question and hypothesis. Bayesianspecific workflow comprises three substeps: (b)–(i) formalizing prior distributions based on background knowledge and prior elicitation; (b)–(ii) determining the likelihood function based on a nonlinear function [math]\displaystyle{ f }[/math]; and (b)–(iii) making a posterior inference. The resulting posterior inference can be used to start a new research cycle.
See also
 Mixed model
 Fixed effects model
 Generalized linear mixed model
 Linear regression
 Mixeddesign analysis of variance
 Multilevel model
 Random effects model
 Repeated measures design
References
 ↑ ^{1.0} ^{1.1} Pinheiro, J; Bates, DM (2006). Mixedeffects models in S and SPLUS. Statistics and Computing. New York: Springer Science & Business Media. doi:10.1007/b98882. ISBN 0387989579.
 ↑ Bolker, BM (2008). Ecological models and data in R. Princeton University Press. https://ms.mcmaster.ca/~bolker/emdbook/.
 ↑ Lindstrom, MJ; Bates, DM (1990). "Nonlinear mixed effects models for repeated measures data". Biometrics 46 (3): 673–687. doi:10.2307/2532087. PMID 2242409.
 ↑ Kuhn, E; Lavielle, M (2005). "Maximum likelihood estimation in nonlinear mixed effects models". Computational Statistics & Data Analysis 49 (4): 1020–1038. doi:10.1016/j.csda.2004.07.002.
 ↑ ^{5.0} ^{5.1} Raket, LL (2020). "Statistical disease progression modeling in Alzheimer's disease". Frontiers in Big Data 3: 24. doi:10.3389/fdata.2020.00024. PMID 33693397.
 ↑ ^{6.0} ^{6.1} "A nonlinear mixedeffects model for simultaneous smoothing and registration of functional data". Pattern Recognition Letters 38: 1–7. 2014. doi:10.1016/j.patrec.2013.10.018.
 ↑ "SITAR—a useful instrument for growth curve analysis". International Journal of Epidemiology 39 (6): 1558–66. 2010. doi:10.1093/ije/dyq115. PMID 20647267.
 ↑ Jonsson, EN; Karlsson, MO; Wade, JR (2000). "Nonlinearity detection: advantages of nonlinear mixedeffects modeling". AAPS PharmSci 2 (3): E32. doi:10.1208/ps020332. PMID 11741248.
 ↑ Lee, Se Yoon; Lei, Bowen; Mallick, Bani (2020). "Estimation of COVID19 spread curves integrating global data and borrowing information". PLOS ONE 15 (7): e0236860. doi:10.1371/journal.pone.0236860. PMID 32726361.
 ↑ Lee, Se Yoon; Mallick, Bani (2021). "Bayesian Hierarchical Modeling: Application Towards Production Results in the Eagle Ford Shale of South Texas". Sankhya B 84: 1–43. doi:10.1007/s13571020002458.
 ↑ Lee, Se Yoon (2022). "Bayesian Nonlinear Models for Repeated Measurement Data: An Overview, Implementation, and Applications". Mathematics 10 (6): 898. doi:10.3390/math10060898.
 ↑ Lee, Se Yoon (2022). "Bayesian Nonlinear Models for Repeated Measurement Data: An Overview, Implementation, and Applications". Mathematics 10 (6): 898. doi:10.3390/math10060898.
Original source: https://en.wikipedia.org/wiki/Nonlinear mixedeffects model.
Read more 