Quasi-Monte Carlo EM Methods for NLME Analysis
R. H. Leary
Pharsight
Problem and Methodology: The classical EM method incorporates an E-step which requires the evaluation of often analytically intractable intergrals. In the context of PK/PD NLME estimation, the integrals of interest are the normalizing factor for the posterior density of the random effects for each subject, as well as the first and second moments of this density.
The Monte-Carlo EM (MCEM) method proposed by Wei and Tanner[1] replaces these integrals with empirical averages using random draws from the distribution of interest. Often this distribution cannot be sampled directly, and techniques such as MCMC (as in the SAEM algorithms in MONOLIX and NM7) and importance sampling (as in the MCPEM algorithm in NM7) are used. The accuracy of the required integrals is controlled by the number of samples N, and typically the theoretical asymptotic error behavior is O(1/sqrt(N)). While often rapid progress can be made in early iterations of MCEM even with low accuracy integrals and small sample sizes, the accuracy requirements increase in the later stage iterations, and very large sample sizes may be required. Thus a significant improvement in efficiency may be obtained in the later stages by lowering the sample size required to achieve the required accuracy.
Here we consider the use of quasi-random (QR) draws using low discrepancy Sobol sequences, as originally proposed for the PK/PD NLME context in the PEM algorithm of Leary et al. [2] The advantage of QR relative to random draws is that the theroretical asymptotic error behavior is now approximately O(1/N), a very significant improvement. A new implementation of PEM based on importance sampling has been developed for the Pharsight Phoenix NLME platform, which utilizes recent enhancements of the QR technique with the scrambling techniques of Owen [3]. This implementation allows the direct comparison of random and QR versions, as well as adaptive error monitoring.
Results: Numerous PK/PD NLME test problems have been compared with random vs. QR importance sampling. Generally the theoretical O(1/N) vs O(1/sqrt(N)) advantage of the QR technique has been confirmed in practice. Improvements in speed of nearly 2 orders of magnitude have been observed for some posterior density integrals where high accuracy is requred, with the sampling requirements being reduced from tens of thousands of points to several hundred. We note that this also confirms a similar QR vs. random draw efficiency improvement obtained by Jank [4] in a MCEM application to a geostatistical generalized linear mixed model.
Conclusions: We have demonstrated that quasirandom integration offers significant practical efficiency advantages relative to the random integration techniques employed in purely stochastic MCEM methods using importance sampling. It is still an open question and an area of active research in the statistical community whether the efficiency advantages of the QR methodology, in conjunction with scrambling, can be extended to MCMC methods such as SAEM.
References:
[1] C. Wei and M. Tanner, J. American Statistical Assoc., 85, 699-704 (1990)
[2] R. H. Leary, R. Jelliffe, A. Schumitzky, and R. E. Port, PAGE 13 (2004), Abstract 491.
[3] A. B. Owen. J. Complexity, 14(4): 466-489 (1998).
[4] W. Jank, Computational Statistics and Data Analysis, 48(4), 685-701 (2005).