Estimation methods for nonlinear mixed-effects modelling have considerably improved over the last decades. Nowadays several algorithms implemented in different softwares are used. The present study aimed at comparing their performance for dose-response models.
Eight scenarios were considered using a sigmoid Emax model, with varying sigmoidicity factors and residual error models. 100 simulated datasets for each scenario were generated. 100 individuals with observations at 4 doses constituted the rich design and at 2 doses for the sparse design. Nine parametric approaches for maximum likelihood estimation were studied: FOCE in NONMEM and R, LAPLACE in NONMEM and SAS, adaptive Gaussian quadrature (AGQ) in SAS, and SAEM in NONMEM and MONOLIX (both SAEM approaches with default and modified settings). All approaches started first from initial estimates set to the true values, and second using altered values. Results were examined through relative root mean squared error (RRMSE) of the estimates.
With true initial conditions, full completion rate was obtained with all approaches except FOCE in R. Runtimes were shortest with FOCE and LAPLACE, and longest with AGQ. Under the rich design with true initial conditions, all approaches performed well except FOCE in R. When starting from altered initial conditions, AGQ, and then FOCE in NONMEM, LAPLACE in SAS, and SAEM in NONMEM and MONOLIX with tuned settings, consistently displayed lower RRMSE than the other approaches.
For standard dose-response models analyzed through mixed-effects models, differences could be identified in the performance of estimation methods available in current software.
MAXIMUM LIKELIHOOD ESTIMATION; FOCE; LAPLACE; ADAPTIVE GAUSSIAN QUADRATURE; SAEM
Estimation methods for nonlinear mixed-effects modelling have considerably improved over the last decades. Nowadays, several algorithms implemented in different software are used. The present study aimed at comparing their performance for dose–response models. Eight scenarios were considered using a sigmoid Emax model, with varying sigmoidicity and residual error models. One hundred simulated datasets for each scenario were generated. One hundred individuals with observations at four doses constituted the rich design and at two doses, the sparse design. Nine parametric approaches for maximum likelihood estimation were studied: first-order conditional estimation (FOCE) in NONMEM and R, LAPLACE in NONMEM and SAS, adaptive Gaussian quadrature (AGQ) in SAS, and stochastic approximation expectation maximization (SAEM) in NONMEM and MONOLIX (both SAEM approaches with default and modified settings). All approaches started first from initial estimates set to the true values and second, using altered values. Results were examined through relative root mean squared error (RRMSE) of the estimates. With true initial conditions, full completion rate was obtained with all approaches except FOCE in R. Runtimes were shortest with FOCE and LAPLACE and longest with AGQ. Under the rich design, all approaches performed well except FOCE in R. When starting from altered initial conditions, AGQ, and then FOCE in NONMEM, LAPLACE in SAS, and SAEM in NONMEM and MONOLIX with tuned settings, consistently displayed lower RRMSE than the other approaches. For standard dose–response models analyzed through mixed-effects models, differences were identified in the performance of estimation methods available in current software, giving material to modellers to identify suitable approaches based on an accuracy-versus-runtime trade-off.
Electronic supplementary material
The online version of this article (doi:10.1208/s12248-012-9349-2) contains supplementary material, which is available to authorized users.
adaptive Gaussian quadrature; FOCE; LAPLACE; maximum likelihood estimation; SAEM
Shrinkage of empirical Bayes estimates (EBEs) of posterior individual parameters in mixed-effects models has been shown to obscure the apparent correlations among random effects and relationships between random effects and covariates. Empirical quantification equations have been widely used for population pharmacokinetic/pharmacodynamic models. The objectives of this manuscript were (1) to compare the empirical equations with theoretically derived equations, (2) to investigate and confirm the influencing factor on shrinkage, and (3) to evaluate the impact of shrinkage on estimation errors of EBEs using Monte Carlo simulations. A mathematical derivation was first provided for the shrinkage in nonlinear mixed effects model. Using a linear mixed model, the simulation results demonstrated that the shrinkage estimated from the empirical equations matched those based on the theoretically derived equations. Simulations with a two-compartment pharmacokinetic model verified that shrinkage has a reversed relationship with the relative ratio of interindividual variability to residual variability. Fewer numbers of observations per subject were associated with higher amount of shrinkage, consistent with findings from previous research. The influence of sampling times appeared to be larger when fewer PK samples were collected for each individual. As expected, sample size has very limited impact on shrinkage of the PK parameters of the two-compartment model. Assessment of estimation error suggested an average 1:1 relationship between shrinkage and median estimation error of EBEs.
Electronic supplementary material
The online version of this article (doi:10.1208/s12248-012-9407-9) contains supplementary material, which is available to authorized users.
empirical Bayes estimate; mixed-effects models; population PK/PD; shrinkage
Missing covariate data is a common problem in nonlinear mixed effects modelling of clinical data. The aim of this study was to implement and compare methods for handling missing covariate data in nonlinear mixed effects modelling under different missing data mechanisms. Simulations generated data for 200 individuals with a 50% difference in clearance between males and females. Three different types of missing data mechanisms were simulated and information about sex was missing for 50% of the individuals. Six methods for handling the missing covariate were compared in a stochastic simulations and estimations study where 200 data sets were simulated. The methods were compared according to bias and precision of parameter estimates. Multiple imputation based on weight and response, full maximum likelihood modelling using information on weight and full maximum likelihood modelling where the proportion of males among the individuals lacking information about sex was estimated (EST) gave precise and unbiased estimates in the presence of missing data when data were missing completely at random or missing at random. When data were missing not at random, the only method resulting in low bias and high parameter precision was EST.
categorical data; covariates; missing data; NONMEM
Multiple imputation (MI) is an approach widely used in statistical analysis of incomplete data. However, its application to missing data problems in nonlinear mixed-effects modelling is limited. The objective was to implement a four-step MI method for handling missing covariate data in NONMEM and to evaluate the method’s sensitivity to η-shrinkage. Four steps were needed; (1) estimation of empirical Bayes estimates (EBEs) using a base model without the partly missing covariate, (2) a regression model for the covariate values given the EBEs from subjects with covariate information, (3) imputation of covariates using the regression model and (4) estimation of the population model. Steps (3) and (4) were repeated several times. The procedure was automated in PsN and is now available as the mimp functionality (http://psn.sourceforge.net/). The method’s sensitivity to shrinkage in EBEs was evaluated in a simulation study where the covariate was missing according to a missing at random type of missing data mechanism. The η-shrinkage was increased in steps from 4.5 to 54%. Two hundred datasets were simulated and analysed for each scenario. When shrinkage was low the MI method gave unbiased and precise estimates of all population parameters. With increased shrinkage the estimates became less precise but remained unbiased.
covariates; missing data; multiple imputation; NONMEM
Efficient power calculation methods have previously been suggested for Wald test-based inference in mixed-effects models but the only available alternative for Likelihood ratio test-based hypothesis testing has been to perform computer-intensive multiple simulations and re-estimations. The proposed Monte Carlo Mapped Power (MCMP) method is based on the use of the difference in individual objective function values (ΔiOFV) derived from a large dataset simulated from a full model and subsequently re-estimated with the full and reduced models. The ΔiOFV is sampled and summed (∑ΔiOFVs) for each study at each sample size of interest to study, and the percentage of ∑ΔiOFVs greater than the significance criterion is taken as the power. The power versus sample size relationship established via the MCMP method was compared to traditional assessment of model-based power for six different pharmacokinetic and pharmacodynamic models and designs. In each case, 1,000 simulated datasets were analysed with the full and reduced models. There was concordance in power between the traditional and MCMP methods such that for 90% power, the difference in required sample size was in most investigated cases less than 10%. The MCMP method was able to provide relevant power information for a representative pharmacometric model at less than 1% of the run-time of an SSE. The suggested MCMP method provides a fast and accurate prediction of the power and sample size relationship.
likelihood ratio test; NONMEM; pharmacometrics; power; sample size
Absorption models used in the estimation of pharmacokinetic drug characteristics from plasma concentration data are generally empirical and simple, utilizing no prior information on gastro-intestinal (GI) transit patterns. Our aim was to develop and evaluate an estimation strategy based on a mechanism-based model for drug absorption, which takes into account the tablet movement through the GI transit. This work is an extension of a previous model utilizing tablet movement characteristics derived from magnetic marker monitoring (MMM) and pharmacokinetic data. The new approach, which replaces MMM data with a GI transit model, was evaluated in data sets where MMM data were available (felodipine) or not available (diclofenac). Pharmacokinetic profiles in both datasets were well described by the model according to goodness-of-fit plots. Visual predictive checks showed the model to give superior simulation properties compared with a standard empirical approach (first-order absorption rate + lag-time). This model represents a step towards an integrated mechanism-based NLME model, where the use of physiological knowledge and in vitro–in vivo correlation helps fully characterize PK and generate hypotheses for new formulations or specific populations.
Electronic supplementary material
The online version of this article (doi:10.1208/s12248-012-9324-y) contains supplementary material, which is available to authorized users.
absorption; model; non-linear mixed effect; semi-mechanistic
The development of covariate models within the population modeling program like NONMEM is generally a time-consuming and non-trivial task. In this study, a fast procedure to approximate the change in objective function values of covariate–parameter models is presented and evaluated. The proposed method is a first-order conditional estimation (FOCE)-based linear approximation of the influence of covariates on the model predictions. Simulated and real datasets were used to compare this method with the conventional nonlinear mixed effect model using both first-order (FO) and FOCE approximations. The methods were mainly assessed in terms of difference in objective function values (ΔOFV) between base and covariate models. The FOCE linearization was superior to the FO linearization and showed a high degree of concordance with corresponding nonlinear models in ΔOFV. The linear and nonlinear FOCE models provided similar coefficient estimates and identified the same covariate–parameter relations as statistically significant or non-significant for the real and simulated datasets. The time required to fit tesaglitazar and docetaxel datasets with 4 and 15 parameter–covariate relations using the linearization method was 5.1 and 0.5 min compared with 152 and 34 h, respectively, with the nonlinear models. The FOCE linearization method allows for a fast estimation of covariate–parameter relations models with good concordance with the nonlinear models. This allows a more efficient model building and may allow the utilization of model building techniques that would otherwise be too time-consuming.
conditional estimation; covariate model building; NONMEM; population PK/PD
Mixed-effect Markov chain models have been recently proposed to characterize the time course of transition probabilities between sleep stages in insomniac patients. The most recent one, based on multinomial logistic functions, was used as a base to develop a final model combining the strengths of the existing ones. This final model was validated on placebo data applying also new diagnostic methods and then used for the inclusion of potential age, gender, and BMI effects. Internal validation was performed through simplified posterior predictive check (sPPC), visual predictive check (VPC) for categorical data, and new visual methods based on stochastic simulation and estimation and called visual estimation check (VEC). External validation mainly relied on the evaluation of the objective function value and sPPC. Covariate effects were identified through stepwise covariate modeling within NONMEM VI. New model features were introduced in the model, providing significant sPPC improvements. Outcomes from VPC, VEC, and external validation were generally very good. Age, gender, and BMI were found to be statistically significant covariates, but their inclusion did not improve substantially the model’s predictive performance. In summary, an improved model for sleep internal architecture has been developed and suitably validated in insomniac patients treated with placebo. Thereafter, covariate effects have been included into the final model.
Electronic supplementary material
The online version of this article (doi:10.1208/s12248-011-9287-4) contains supplementary material, which is available to authorized users.
covariates; Markov; multinomial; sleep; validation
Informative diagnostic tools are vital to the development of useful mixed-effects models. The Visual Predictive Check (VPC) is a popular tool for evaluating the performance of population PK and PKPD models. Ideally, a VPC will diagnose both the fixed and random effects in a mixed-effects model. In many cases, this can be done by comparing different percentiles of the observed data to percentiles of simulated data, generally grouped together within bins of an independent variable. However, the diagnostic value of a VPC can be hampered by binning across a large variability in dose and/or influential covariates. VPCs can also be misleading if applied to data following adaptive designs such as dose adjustments. The prediction-corrected VPC (pcVPC) offers a solution to these problems while retaining the visual interpretation of the traditional VPC. In a pcVPC, the variability coming from binning across independent variables is removed by normalizing the observed and simulated dependent variable based on the typical population prediction for the median independent variable in the bin. The principal benefit with the pcVPC has been explored by application to both simulated and real examples of PK and PKPD models. The investigated examples demonstrate that pcVPCs have an enhanced ability to diagnose model misspecification especially with respect to random effects models in a range of situations. The pcVPC was in contrast to traditional VPCs shown to be readily applicable to data from studies with a priori and/or a posteriori dose adaptations.
Electronic supplementary material
The online version of this article (doi:10.1208/s12248-011-9255-z) contains supplementary material, which is available to authorized users.
mixed-effects modeling; model diagnostics; pcVPC; prediction correction; VPC
It is not uncommon that the outcome measurements, symptoms or side effects, of a clinical trial belong to the family of event type data, e.g., bleeding episodes or emesis events. Event data is often low in information content and the mixed-effects modeling software NONMEM has previously been shown to perform poorly with low information ordered categorical data. The aim of this investigation was to assess the performance of the Laplace method, the stochastic approximation expectation–maximization (SAEM) method, and the importance sampling method when modeling repeated time-to-event data. The Laplace method already existed, whereas the two latter methods have recently become available in NONMEM 7. A stochastic simulation and estimation study was performed to assess the performance of the three estimation methods when applied to a repeated time-to-event model with a constant hazard associated with an exponential interindividual variability. Various conditions were investigated, ranging from rare to frequent events and from low to high interindividual variability. The method performance was assessed by parameter bias and precision. Due to the lack of information content under conditions where very few events were observed, all three methods exhibit parameter bias and imprecision, however most pronounced by the Laplace method. The performance of the SAEM and importance sampling were generally higher than Laplace when the frequency of individuals with events was less than 43%, while at frequencies above that all methods were equal in performance.
importance sampling; Laplace; mixed-effects modeling; NONMEM; repeated time-to-event; SAEM
This article demonstrates techniques for describing and predicting disease progression in acute stroke by modeling scores measured using clinical assessment scales, accommodating dropout as an additional source of information. Scores assessed using the National Institutes of Health Stroke Scale and the Barthel Index in acute stroke patients were used to model the time course of disease progression. Simultaneous continuous and probabilistic models for describing the nature and magnitude of score changes were developed, and used to model the trajectory of disease progression using scale scores. The models described the observed data well, and exhibited good simulation properties. Applications include longitudinal analysis of stroke scale data, clinical trial simulation, and prognostic forecasting. Based upon experience in other areas, it is likely that application of this modeling methodology will enable reductions in the number of patients needed to carry out clinical studies of treatments for acute stroke.
Electronic supplementary material
The online version of this article (doi:10.1208/s12248-010-9230-0) contains supplementary material, which is available to authorized users.
Barthel index; disease progression; NIH stroke scale; NONMEM; stroke
A nonparametric population method with support points from the empirical Bayes estimates (EBE) has recently been introduced (default method). However, EBE distribution may, with sparse and small datasets, not provide a suitable range of support points. This study aims to develop a method based on a prior parametric analysis capable of providing a nonparametric grid with adequate support points range. A new method extends the nonparametric grid with additional support points generated by simulation from the parametric distribution, hence the name extended-grid method. The joint probability density function is estimated at the extended grid. The performance of the new method was evaluated and compared to the default method via Monte Carlo simulations using simple IV bolus model and sparse (200 subject, two samples per subject) or small (30 subjects, three samples per subjects) datasets and two scenarios based on real case studies. Parameter distributions estimated by the default and the extended-grid method were compared to the true distributions; bias and precision were assessed at different percentiles. With small datasets, the bias was similar between methods (<10%); however, precision was markedly improved with the new method (by 43%). With sparse datasets, both bias (from 5.9% to 3%) and precision (by 60%) were improved. For simulated scenarios based on real study designs, extended-grid predictions were in a good agreement with true values. A new approach to obtain support points for the nonparametric method has been developed, and it displayed good estimation properties. The extended-grid method is automated, using the program PsN, for implementation into the NONMEM.
empirical Bayes estimates; extended grid method; NONMEM; nonparametric method; parameter distribution
Empirical Bayes (“post hoc”) estimates (EBEs) of ηs provide modelers with diagnostics: the EBEs themselves, individual prediction (IPRED), and residual errors (individual weighted residual (IWRES)). When data are uninformative at the individual level, the EBE distribution will shrink towards zero (η-shrinkage, quantified as 1-SD(ηEBE)/ω), IPREDs towards the corresponding observations, and IWRES towards zero (ε-shrinkage, quantified as 1-SD(IWRES)). These diagnostics are widely used in pharmacokinetic (PK) pharmacodynamic (PD) modeling; we investigate here their usefulness in the presence of shrinkage. Datasets were simulated from a range of PK PD models, EBEs estimated in non-linear mixed effects modeling based on the true or a misspecified model, and desired diagnostics evaluated both qualitatively and quantitatively. Identified consequences of η-shrinkage on EBE-based model diagnostics include non-normal and/or asymmetric distribution of EBEs with their mean values (“ETABAR”) significantly different from zero, even for a correctly specified model; EBE–EBE correlations and covariate relationships may be masked, falsely induced, or the shape of the true relationship distorted. Consequences of ε-shrinkage included low power of IPRED and IWRES to diagnose structural and residual error model misspecification, respectively. EBE-based diagnostics should be interpreted with caution whenever substantial η- or ε-shrinkage exists (usually greater than 20% to 30%). Reporting the magnitude of η- and ε-shrinkage will facilitate the informed use and interpretation of EBE-based diagnostics.
empirical Bayes estimate; model building; model evaluation; NONMEM; shrinkage
The purpose of this study is to investigate the impact of observations below the limit of quantification (BQL) occurring in three distinctly different ways and assess the best method for prevention of bias in parameter estimates and for illustrating model fit using visual predictive checks (VPCs). Three typical ways in which BQL can occur in a model was investigated with simulations from three different models and different levels of the limit of quantification (LOQ). Model A was used to represent a case with BQL observations in an absorption phase of a PK model whereas model B represented a case with BQL observations in the elimination phase. The third model, C, an indirect response model illustrated a case where the variable of interest in some cases decreases below the LOQ before returning towards baseline. Different approaches for handling of BQL data were compared with estimation of the full dataset for 100 simulated datasets following models A, B, and C. An improved standard for VPCs was suggested to better evaluate simulation properties both for data above and below LOQ. Omission of BQL data was associated with substantial bias in parameter estimates for all tested models even for seemingly small amounts of censored data. Best performance was seen when the likelihood of being below LOQ was incorporated into the model. In the tested examples this method generated overall unbiased parameter estimates. Results following substitution of BQL observations with LOQ/2 were in some cases shown to introduce bias and were always suboptimal to the best method. The new standard VPCs was found to identify model misfit more clearly than VPCs of data above LOQ only.
likelihood; limit-of-quantification; NONMEM; visual-predictive-check
In nonlinear mixed effects modeling using NONMEM, mixture models can be used for multimodal distributions of parameters. The fraction of individuals belonging to each of the subpopulations can be estimated, and the most probable subpopulation for each patient is output (MIXESTk). The objective function value (OFV) that is minimized is the sum of the OFVs for each patient (OFVi), which in turn is the sum across the k subpopulations (OFVi,k). The OFVi,k values can be used together with the total probability in the population of belonging to subpopulation k to calculate the individual probability of belonging to the subpopulation (IPk). Our objective was to explore the information gained by using IPk instead of or in addition to MIXESTk in the analysis of mixture models. Two real data sets described previously by mixture models as well as simulations were used to explore the use of IPk and the precision of individual parameter values based on IPk and MIXESTk. For both real data-based mixture models, a substantial fraction (11% and 26%) of the patients had IPk values not close to 0 or 1 (IPk between 0.25 and 0.75). Simulations of eight different scenarios showed that individual parameter estimates based on MIXEST were less precise than those based on IPk, as the root mean squared error was reduced for IPk in all scenarios. A probability estimate such as IPk provides more detailed information about each individual than the discrete MIXESTk. Individual parameter estimates based on IPk should be preferable whenever individual parameter estimates are to be used as study output or for simulations.
mixture; NONMEM; pharmacometrics; population modeling; subpopulation
A significant bias in parameters, estimated with the proportional odds model using the software NONMEM, has been reported. Typically, this bias occurs with ordered categorical data, when most of the observations are found at one extreme of the possible outcomes. The aim of this study was to assess, through simulations, the performance of the Back-Step Method (BSM), a novel approach for obtaining unbiased estimates when the standard approach provides biased estimates. BSM is an iterative method involving sequential simulation-estimation steps. BSM was compared with the standard approach in the analysis of a 4-category ordered variable using the Laplacian method in NONMEM. The bias in parameter estimates and the accuracy of model predictions were determined for the 2 methods on 3 conditions: (1) a nonskewed distribution of the response with low interindividual variability (IIV), (2) a skewed distribution with low IIV, and (3) a skewed distribution with high IIV. An increase in bias with increasing skewness and IIV was shown in parameters estimated using the standard approach in NON-MEM. BSM performed without appreciable bias in the estimates under the 3 conditions, and the model predictions were in good agreement with the original data. Each BSM estimation represents a random sample of the population; hence, repeating the BSM estimation reduces the imprecision of the parameter estimates. The BSM is an accurate estimation method when the standard modeling approach in NONMEM gives biased estimates.
ordered categorical; proportional odds model; bias in parameter estimates; NONMEM; Laplacian; pharmacodynamics