Numerical discretization-based estimation methods for ordinary differential equation models via penalized spline smoothing with applications in biomedical research

Differential equations are extensively used for modeling dynamics of physical processes in many scientific fields such as engineering, physics, and biomedical sciences. Parameter estimation of differential equation models is a challenging problem because of high computational cost and high-dimensional parameter space. In this article, we propose a novel class of methods for estimating parameters in ordinary differential equation (ODE) models, which is motivated by HIV dynamics modeling. The new methods exploit the form of numerical discretization algorithms for an ODE solver to formulate estimating equations. First, a penalized-spline approach is employed to estimate the state variables and the estimated state variables are then plugged in a discretization formula of an ODE solver to obtain the ODE parameter estimates via a regression approach. We consider three different order of discretization methods, Euler's method, trapezoidal rule, and Runge-Kutta method. A higher-order numerical algorithm reduces numerical error in the approximation of the derivative, which produces a more accurate estimate, but its computational cost is higher. To balance the computational cost and estimation accuracy, we demonstrate, via simulation studies, that the trapezoidal discretization-based estimate is the best and is recommended for practical use. The asymptotic properties for the proposed numerical discretization-based estimators are established. Comparisons between the proposed methods and existing methods show a clear benefit of the proposed methods in regards to the trade-off between computational cost and estimation accuracy. We apply the proposed methods t an HIV study to further illustrate the usefulness of the proposed approaches.     

Dynamic Lung Modeling and Tumor Tracking Using Deformable Image Registration and Geometric Smoothing

A greyscale-based fully automatic deformable image registration algorithm, based on an optical flow method together with geometric smoothing, is developed for dynamic lung modeling and tumor tracking. In our computational processing pipeline, the input data is a set of 4D CT images with 10 phases. The triangle mesh of the lung model is directly extracted from the more stable exhale phase (Phase 5). In addition, we represent the lung surface model in 3D volumetric format by applying a signed distance function and then generate tetrahedral meshes. Our registration algorithm works for both triangle and tetrahedral meshes. In CT images, the intensity value reflects the local tissue density. For each grid point, we calculate the displacement from the static image (Phase 5) to match with the moving image (other phases) by using merely intensity values of the CT images. The optical flow computation is followed by a regularization of the deformation field using geometric smoothing. Lung volume change and the maximum lung tissue movement are used to evaluate the accuracy of the application. Our testing results suggest that the application of deformable registration algorithm is an effective way for delineating and tracking tumor motion in image-guided radiotherapy.     

Estimating a predator-prey dynamical model with the parameter cascades method

Summary .   Ordinary differential equations (ODEs) are widely used in ecology to describe the dynamical behavior of systems of interacting populations. However, systems of ODEs rarely provide quantitative solutions that are close to real field observations or experimental data because natural systems are subject to environmental and demographic noise and ecologists are often uncertain about the correct parameterization. In this article we introduce "parameter cascades" as an improved method to estimate ODE parameters such that the corresponding ODE solutions fit the real data well. This method is based on the modified penalized smoothing with the penalty defined by ODEs and a generalization of profiled estimation, which leads to fast estimation and good precision for ODE parameters from noisy data. This method is applied to a set of ODEs originally developed to describe an experimental predator–prey system that undergoes oscillatory dynamics. The new parameterization considerably improves the fit of the ODE model to the experimental data sets. At the same time, our method reveals that important structural assumptions that underlie the original ODE model are essentially correct. The mathematical formulations of the two nonlinear interaction terms (functional responses) that link the ODEs in the predator–prey model are validated by estimating the functional responses nonparametrically from the real data. We suggest two major applications of "parameter cascades" to ecological modeling: It can be used to estimate parameters when original data are noisy, missing, or when no reliable priori estimates are available; it can help to validate the structural soundness of the mathematical modeling approach.     

Beals smoothing revisited

Beals smoothing is a multivariate transformation specially designed for species presence/absence community data containing noise and/or a lot of zeros. This transformation replaces the observed values of the target species by predictions of occurrence on the basis of its co-occurrences with the remaining species. In many applications, the transformed values are used as input for multivariate analyses. As Beals smoothing values provide a sense of “probability of occurrence”, they have also been used for inference. However, this transformation can produce spurious results, and it must be used with caution. Here we study the statistical and ecological bases underlying the Beals smoothing function, and the factors that may affect the reliability of transformed values are explored using simulated data sets. Our simulations demonstrate that Beals predictions are unreliable for target species that are not related to the overall ecological structure. Furthermore, the presence of these “random” species may diminish the quality of Beals smoothing values for the remaining species. A statistical test is proposed to determine when observed values can be replaced with Beals smoothing predictions. Two real-data example applications are presented to illustrate the potentially false predictions of Beals smoothing and the necessary checking step performed by the new test. Electronic supplementary material The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

Measuring size and composition of species pools: a comparison of dark diversity estimates

Ecological theory and biodiversity conservation have traditionally relied on the number of species recorded at a site, but it is agreed that site richness represents only a portion of the species that can inhabit particular ecological conditions, that is, the habitat‐specific species pool. Knowledge of the species pool at different sites enables meaningful comparisons of biodiversity and provides insights into processes of biodiversity formation. Empirical studies, however, are limited due to conceptual and methodological difficulties in determining both the size and composition of the absent part of species pools, the so‐called dark diversity. We used >50,000 vegetation plots from 18 types of habitats throughout the Czech Republic, most of which served as a training dataset and 1083 as a subset of test sites. These data were used to compare predicted results from three quantitative methods with those of previously published expert estimates based on species habitat preferences: (1) species co‐occurrence based on Beals' smoothing approach; (2) species ecological requirements, with envelopes around community mean Ellenberg values; and (3) species distribution models, using species environmental niches modeled by Biomod software. Dark diversity estimates were compared at both plot and habitat levels, and each method was applied in different configurations. While there were some differences in the results obtained by different methods, particularly at the plot level, there was a clear convergence, especially at the habitat level. The better convergence at the habitat level reflects less variation in local environmental conditions, whereas variation at the plot level is an effect of each particular method. The co‐occurrence agreed closest the expert estimate, followed by the method based on species ecological requirements. We conclude that several analytical methods can estimate species pools of given habitats. However, the strengths and weaknesses of different methods need attention, especially when dark diversity is estimated at the plot level.     

Semiparametric additive time-varying coefficients model for longitudinal data with censored time origin

Statistical analysis of longitudinal data often involves modeling treatment effects on clinically relevant longitudinal biomarkers since an initial event (the time origin). In some studies including preventive HIV vaccine efficacy trials, some participants have biomarkers measured starting at the time origin, whereas others have biomarkers measured starting later with the time origin unknown. The semiparametric additive time-varying coefficient model is investigated where the effects of some covariates vary nonparametrically with time while the effects of others remain constant. Weighted profile least squares estimators coupled with kernel smoothing are developed. The method uses the expectation maximization approach to deal with the censored time origin. The Kaplan–Meier estimator and other failure time regression models such as the Cox model can be utilized to estimate the distribution and the conditional distribution of left censored event time related to the censored time origin. Asymptotic properties of the parametric and nonparametric estimators and consistent asymptotic variance estimators are derived. A two-stage estimation procedure for choosing weight is proposed to improve estimation efficiency. Numerical simulations are conducted to examine finite sample properties of the proposed estimators. The simulation results show that the theory and methods work well. The efficiency gain of the two-stage estimation procedure depends on the distribution of the longitudinal error processes. The method is applied to analyze data from the Merck 023/HVTN 502 Step HIV vaccine study.     

Dynamic logistic state space prediction model for clinical decision making

Prediction modeling for clinical decision making is of great importance and needed to be updated frequently with the changes of patient population and clinical practice. Existing methods are either done in an ad hoc fashion, such as model recalibration or focus on studying the relationship between predictors and outcome and less so for the purpose of prediction. In this article, we propose a dynamic logistic state space model to continuously update the parameters whenever new information becomes available. The proposed model allows for both time-varying and time-invariant coefficients. The varying coefficients are modeled using smoothing splines to account for their smooth trends over time. The smoothing parameters are objectively chosen by maximum likelihood. The model is updated using batch data accumulated at prespecified time intervals, which allows for better approximation of the underlying binomial density function. In the simulation, we show that the new model has significantly higher prediction accuracy compared to existing methods. We apply the method to predict 1 year survival after lung transplantation using the United Network for Organ Sharing data.