A Comparison of Two-Stage Approaches for Fitting Nonlinear Ordinary Differential Equation Models with Mixed Effects

Abstract

Several approaches exist for estimating the derivatives of observed data for model exploration purposes, including functional data analysis (FDA; Ramsay & Silverman, 2005 ), generalized local linear approximation (GLLA; Boker, Deboeck, Edler, & Peel, 2010 ), and generalized orthogonal local derivative approximation (GOLD; Deboeck, 2010 ). These derivative estimation procedures can be used in a two-stage process to fit mixed effects ordinary differential equation (ODE) models. While the performance and utility of these routines for estimating linear ODEs have been established, they have not yet been evaluated in the context of nonlinear ODEs with mixed effects. We compared properties of the GLLA and GOLD to an FDA-based two-stage approach denoted herein as functional ordinary differential equation with mixed effects (FODEmixed) in a Monte Carlo (MC) study using a nonlinear coupled oscillators model with mixed effects. Simulation results showed that overall, the FODEmixed outperformed both the GLLA and GOLD across all the embedding dimensions considered, but a novel use of a fourth-order GLLA approach combined with very high embedding dimensions yielded estimation results that almost paralleled those from the FODEmixed. We discuss the strengths and limitations of each approach and demonstrate how output from each stage of FODEmixed may be used to inform empirical modeling of young children's self-regulation.

Keywords: Differential equation; derivatives; dynamic; dynamical systems; functional data analysis.

Figure 1

(A) Over-time trajectories of x(t) from the proportional change model generated using based on a numerical ODE solver; (B) phase plane of levels and first derivatives in the proportional change model. The solid arrow points to one of the equilibrium points of the system at x(t) = x^* = 3, defined as points at which $\frac{dx}{dt}=0$. (C) Over-time trajectories of x(t) from the damped oscillator model generated using based on a numerical ODE solver at different values of ζ and η set to −0.8; (D) phase plane of levels and second derivatives in the damped oscillator model. (E) phase plane of first and second derivatives in the damped oscillator model (F)–(G): component + residuals plot revealing the associations between the residuals and x(t) and dx(t)/dt, respectively, with the linear effect of the other predictor partialled out.

Figure 2

(A)–(B) Over-time trajectories of the predator and prey population sizes from the classical predator-prey model and the corresponding phase plane plot generated with r₁ = 4, r₂ = −3.5,a₁₂ = −2 and a₂₁ = 1.5; (C)–(D) over-time trajectories of the two species’ sizes and the corresponding phase plane plot after adding a quadratic element, $-0.20x_1^2(t)$ to the prey equation in Equation 10; (E) Over-time trajectories of x₁(t) from the van der Pol model generated with η = 1, ζ = 1.0, 3.0, and 5.0; (F) the phase plane of the second and first derivatives; (G) component + residuals plot revealing the association between d²x(t)/dt² and dx(t)/dt after the linear effect of x(t) has been partialled out; (H) plot of the residuals from an expanded model (a model in which x(t), dx(t)/dt and x²(t)dx(t)/dt were included as predictors) against dx(t)/dt.

Figure 3

(A) True and estimated first derivatives from the van der Pol oscillator model with no measurement noise; (B) true and estimated second derivatives with no measurement noise; (C)–(D) true and estimated first and second derivatives when reliability = 0.9. The lowest panels show boxplots of the second derivative estimates from the GOLD, the GLLA and the 4th-order GLLA when d, the embedding dimension, was set to 5. Here, reliability was computed as the ratio between the variance of the true x(t) generated by means of a numerical ODE solver and the sum of the true x(t) and measurement error variance. The same legends were used in plots (A)–(D) but are only shown in plots (A)–(B) due to space limitations.

Figure 4

Results from fitting the coupled nonlinear oscillators model (see Equation 23) across sample size conditions, with (A) the RMSEs of the fixed effects point estimates; (B) the RMSEs of the random effect SD parameters, (C) biases in the SE estimates and (D) proportions of non-convergent replications.

Figure 5

(A) Over-time trajectories of EP and PR from 5 randomly selected children; component + residuals plot revealing the partial associations between d²PR_i(t)/dt² and: (B) demeaned PR_i(t), (C) demeaned EP_i(t), and (D) $\frac{{\mathit{dPR}}_i(t)}{dt}$; (E) the component + residuals plot revealing the partial associations between $\frac{d^2{EP}_i(t)}{{dt}^2}$ and $\frac{{\mathit{dEP}}_i(t)}{dt}$.

Figure 6

(A) Observed data and predicted trajectories for one randomly selected partipant generated using parameter estimates from the final model; (B)–(C) simulated trajectories illustrating the effects of γ_EP,1 and ζ_EP,2, respectively, for a hypothetical child with high (+2SD above the sample average) and low (−2SD below the sample average) effortful control. (D) simulated trajectories generated using the estimated parameters from the final nonlinear model in comparison to trajectories generated by setting the values of γ_EP,3 and γ_{P R,3} to zero; (E–F) discrepancies between the observed data and predicted trajectories for the same participant using the final nonlinear model and the linear comparison model, for EP and PR, respectively. The shaded regions of plots (B)–(D) correspond to regions where the predicted EP trajectories marked with solid triangles show accelerations (positive second derivatives). The unshaded regions correspond to regions with decelerations (negative second derivatives) or second derivative value of zero (points of inflection). EC = effortful control.