Uncertainty computation at finite distance in nonlinear mixed effects models - a new method evaluated on simulations and applied to the evolution of clinical status of patients hospitalised for COVID-19
Mélanie Guhl (1), Julie Bertrand (1), Jérémie Guedj (1), France Mentré (1,2), Emmanuelle Comets(1,3)
(1) Université Paris Cité and Université Sorbonne Paris Nord, Inserm, IAME, F-75018 Paris, France (2) AP-HP Hôpital Bichat, Département d’Epidémiologie Biostatistiques et Recherche Clinique, Paris, France (3) Univ Rennes, Inserm, EHESP, Irset - UMR S 1085, F-35000 Rennes, France
Number of characters: 4429/4500
Objectives:
The standard error (SE) of the maximum likelihood estimate (MLE) of the population parameter vector θ in nonlinear mixed effect models (NLMEM) is usually estimated as the inverse of the Fisher Information Matrix (FIM). However, at finite distance, the FIM can underestimate the SE of NLMEM [1,2]. As the limit distributions of the MLE and the maximum a posterior estimator in a Bayesian framework are equivalent (Bernstein-von Mises theorem), the standard deviation of the posterior distribution, obtained in Stan via the HMC algorithm, has been shown to be a proxy for the SE [3,4]. Here, we develop a similar method using the Metropolis Hastings (MH) algorithm and implement it in the saemix R package [5]. We assess it with different simulation sets and a real dataset.
Methods:
We use the MH algorithm, already embedded in the SAEM algorithm for other purposes, to draw from the posterior distribution of θ. The parameters to be calibrated are the length of the chain, the prior distribution and the kernel distribution. For now, the Bayesian step uses the frequentist estimations as parameters of the proposal kernel. We obtain a frequentist and a Bayesian estimate of the SE as respectively the inverse of the FIM and the standard deviation of the chain sampled.
Our first set of simulations used the one-compartment theophylline pharmacokinetic (PK) model with linear absorption and elimination, and a proportional error model. The number of patients varied from N=150 with n=10 sampling points per patient, to N=12 and n=3. The evaluation of the method was based on the MH acceptation rates, the boxplots of SE (the target being the empirical SE obtained with SAEM) and the 95% coverage rates computed over 1000 simulated datasets. We compared our method to the FIM (Asympt), the HMC algorithm in Stan (Post) and the sampling importance resampling (SIR) [6].
We then assessed the method in more challenging settings with larger inter-individual variability and correlations between the random effects.
We also applied this method on the evolution of the NEWS-2 [7], a composite clinical score ranging from 0 (best) to 20 (worst), in patients hospitalised for Covid19, using real data from the Discovery trial [8], a European clinical trial aiming to evaluate antiviral drugs for the treatment of COVID-19. We modelled the evolution of NEWS-2 and tested for a treatment effect between the standard of care (SoC, N=408 patients) and SoC + remdesivir (N=402 patients) arms.
Results:
On theophylline simulations with N=150 and n=10, the SE obtained from all methods were close to the target. Coverage rates (CR) were appropriate. These results validate the Bernstein von-Mises theorem supporting our approach.
On N=12 and n=3, Asympt, SIR and MH underestimated the SE. Coverage rates were below the prediction interval. MH gave improved results when inflating the variance of the kernel. In this scenario, the SE were overestimated by Post, giving coverage rates over or under the target. Further work is needed to investigate suitable priors.
In both cases, the acceptation rates of MH were between 15% and 40% which seems reasonable given the number of dimensions we are working with. They dropped however with increased variability or strong correlations, and acceptation ratios were very sensitive in this scenario.
These first results show that, as expected, the FIM underestimates the SE at finite distance. In this case, the Bayesian paradigm seems promising but we need to further investigate how to calibrate our implemented method on challenging settings, using the acceptation rate as a tool to diagnose the optimal inflation on the kernel variance.
On Discovery data, no significant treatment effect of remdesivir on the evolution of the NEWS2 score was shown. The MH method gave lower SE than all other methods and acceptation rates were critically low, which could be explained by the very high inter-individual variabilities in the model selected (100 to 200%). There were also strong correlations between random effects (one at -0.98).
Conclusions:
Post was the most stable method in both simulated and real data settings, confirming previous studies [3,4] and the validity of the Bayesian paradigm. Our algorithm still requires calibration in more challenging settings, where acceptation rates became critically low. Perspectives for improvement include block sampling or conditional univariate kernels.
References:
[1] Bertrand J, Comets E, Lafont CM, Chenel M, Mentré F, Pharmacogenetics and population pharmacokinetics: impact of the design on three tests using the SAEM algorithm, Biometrics 2011.
[2] Dubois A, Lavielle M, Gsteiger S, Pigeolet E, Mentré F, Model-based analyses of bioequivalence crossover trials using the stochastic approximation expectation maximisation algorithm, Stat in Med 2012.
[3] Loingeville F, Bertrand J, Nguyen T, Sharan S, Feng K et al., New model-based bioequivalence statistical approaches for pharmacokinetic studies with sparse sampling, AAPS J 2020.
[4] Guhl M, Mercier F, Hofmann C, Sharan S, Donnelly M et al., Impact of model misspecification on model-based tests in PK studies with parallel design: real case and simulation studies, J Pharmacokinet Pharmacodyn 2022.
[5] Comets E, Lavenu A and Lavielle M, Parameter estimation in nonlinear mixed effect models using saemix, an R implementation of the SAEM algorithm, J Stat Soft 2017.
[6] Dosne AG, Bergstrand M, Harling K, Karlsson MO, Improving the estimation of parameter uncertainty distributions in nonlinear mixed effects models using sampling importance resampling, J Pharmacokin Pharmacodyn 2016.
[7] Smith GB, Prytherch DR, Meredith P, Schmidt PE and Featherstone PI, The ability of the National Early Warning Score (NEWS) to discriminate patients at risk of early cardiac arrest, unanticipated intensive care unit admission, and death, Resuscitation 2013.
[8] Ader F, Bouscambert-Duchamp M, Hites M, Peiffer-Smadja N, Poissy J et al., Remdesivir plus standard of care versus standard of care alone for the treatment of patients admitted to hospital with COVID-19 (DisCoVeRy): a phase 3, randomised, controlled, open-label trial, Lancet Infect Dis 2022.
