University of Southern California Dissertations and Theses

Bayesian estimation using Markov chain Monte Carlo methods in pharmacokinetic system analysis

(USC Thesis Other)

Bayesian estimation using Markov chain Monte Carlo methods in pharmacokinetic system analysis

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content INFORMATION TO USERS This manuscript has been reproduced from the microfilm master. UMI films the text directly from the original or copy submitted. Thus, some thesis and dissertation copies are in typewriter face, while others may be from any type of computer printer. The quality of this reproduction is dependent upon the quality of the copy submitted. Broken or indistinct print, colored or poor quality illustrations and photographs, print bleedthrough, substandard margins, and improper alignment can adversely affect reproduction. In the unlikely event that the author did not send UMI a complete manuscript and there are missing pages, these will be noted. Also, if unauthorized copyright material had to be removed, a note will indicate the deletion. Oversize materials (e.g., maps, drawings, charts) are reproduced by sectioning the original, beginning at the upper left-hand comer and continuing from left to right in equal sections with small overlaps. Each original is also photographed in one exposure and is included in reduced form at the back of the book. Photographs included in the original manuscript have been reproduced xerographically in this copy. Higher quality 6” x 9” black and white photographic prints are available for any photographs or illustrations appearing in this copy for an additional charge. Contact UMI directly to order. UMI A Bell & Howell Inform ation Company 300 North Zeeb Road, Ann A rbor MI 48106-1346 USA 313/761-4700 800/521-0600 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. NOTE TO USERS The original manuscript received by UMI contains pages with indistinct, light and/or slanted print. Pages were microfilmed as received. This reproduction is the best copy available UMI Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. BAYESIAN ESTIMATION USING MARKOV CHAIN MONTE CARLO METHODS IN PHARMACOKINETIC SYSTEM ANALYSIS Copyright 1998 by Dongwoo Kang A Thesis Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree Master of Science ( Biomedical Engineering) August 1998 Dongwoo Kang Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. UMI Number: 1393173 UMI Microform 1393173 Copyright 1999, by UMI Company. A H rights reserved. This microform edition is protected against unauthorized copying under Title 17, United States Code. UMI 300 North Zeeb Road Ann Arbor, MI 48103 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. This thesis, written by _Dojigwpo Kang__________________________________ under the guidance of his/her Faculty Committee and approved by all its members, has been presented to and accepted by the School of Engineering in partial fidfUlment of the re quirements for the degree of Master o f Science _ Date August 4, 1998 Faculty Committee mm k'cdjoJsT *.— Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. ACKNOWLEDGEMENTS I thank Professor David Z. D’Argenio for his guidance and support that was indispensable for this work. I also thank Dr. John Rodman at St. Jude Children’s Research Hospital in Memphis, Tennessee, for providing me with the patient data and experiment protocol. My sincere gratitude goes to those who helped build the ADAPT II software package, too. On the base of the previous work, I could add RIMH algorithm and test both of the methods. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. iii TABLE OF CONTENTS 1. Introduction 1 1.1 Pharmacokinetics 1 1.2 Modeling Approaches in Pharmacokinetics 2 1.3 A Modeling Framework for Pharmacokinetic Systems 3 1.4 Specific Aims 6 2. Bayesian Estimation 8 2.1 Bayesian Approach vs. Frequentist Approach 8 2.2 Point Estimation and Interval Estimation 10 2.3 Computational Methods 12 2.3.1 Analytic Approximation 12 2.3.2 Numerical Integration 12 2.3.3 Iterated Sampling and Resampling Approaches 15 3. Markov Chain Monte Carlo ( MCMC ) Methods 17 3.1 Monte Carlo Integration 17 3.2 Maikov Chain and Stationary Distribution 18 3.2.1 Basic Concept and Definitions 18 3.2.2 Stationary Distribution 20 3.2.3 Construction of a Markov Chain 23 3.3 MCMC Algorithms 26 3.4 Independence Metropolis-Hastings ( IM H) Sampler 29 3.4.1 Independence and Acceptance Probability 29 3.4.2 IMH Algorithm 30 3.5 Regenerative Independence Metropolis-Hastings ( RIM H) Sampler 31 3.5.1 Rejection Sampling and Acceptance Probability 31 3.5.2 Regenerative Process 33 3.5.3 RIMH Algorithm 36 4. Application with Busulfan 38 4.1 General Description of the Drug 38 4.2 Experiment and Simulation Protocol 39 4.3 Results from IMH Sampler 44 4.3.1 General Information 44 4.3.2 Convergence, Mixing, and Correlation 50 4.3.3 Independent Samples 57 4.4 Results from RIMH Sampler 62 4.4.1 General Information 62 4.4.2 Regeneration and Asymptotic Confidence Interval 66 4.4.3 Different c 69 5. Discussion and Conclusion 72 References 74 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. iv LIST OF TABLES 4.1 Prior mean for parameters from patient data 40 4.2 Prior covariance matrix for parameters from patient data 40 4.3 Prior mean for parameters from population simulation 42 4.4 Prior covariance matrix for parameters from population simulation 42 4.5 Individual parameter data obtained from population simulation 43 4.6 Individual drug concentration data obtained from population simulation to be used in Bayesian estimation using MCMC methods 43 4.7 Posterior mean and standard deviation calculated from 5000 samples with IMH algorithm 46 4.8 Posterior mean and standard deviation for Subj90. Calculated from 5000 samples, 100,000 samples, and 3333 independent samples originally from 100,000 samples with IMH algorithm 58 4.9 Posterior mean and standard deviation calculated from 5000 samples with RIMH algorithm 62 4.10 Asymptotic 95% confidence intervals for V o f Subj90 for various sample lengths 68 4.11 Number o f subchains and average length of them calculated from 5000 samples with RIMH algorithm for Subj90 with different c 70 4.12 Posterior mean and standard deviation calculated from 5000 samples with RIMH algorithm for Subj90 with different c 70 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. V LIST OF FIGURES 1. One compartment model used in this study . 4.1 Concentration plot from population simulation of 1000 4.2 AUC plot from population simulation of 1000 4.3 Markov chains formed by IMH algorithm for Ke 4.4.a Posterior distribution of Ke and Ka for the three individuals with population distribution. Result from IMH sampler 4.4.b Posterior distribution of V and F for the three individuals with population distribution. Result from IMH sampler 4.4.c Posterior distribution o f drug concentration and AUC at t=l 8 hours for the three individuals with population distribution. Result from IMH sampler 4.5 Sequences of Ke started from different initial values for Subj90 4.6 Autocorrelation plots to assess the rate of mixing 4.7 Scatter plots for Subj90 at the lag o f 30 4.8 Selected prior joint densities ( p(V,Ke), p(Ka,Ke), p(Ka,V), and p(F, V) ) from population simulation 4.9 Selected posterior joint densities ( p(V,Ke|z), p(Ka,Ke|z), p(Ka,V|z), and p(F,V|z)) calculated from 5000 samples for Subj90 4. lO.a Posterior distribution of Ke and Ka for Subj90 with population distribution to consider independent samples 4. lO.b Posterior distribution of V and F for Subj90 with population distribution to consider independent samples 4.10.C Posterior distribution o f drug concentration at t=18 hours for Subj90 with population distribution to consider independent samples 4.11 Selected posterior joint densities ( p(V,Ke|z), p(Ka,Ke|z), p(Ka,V|z), and p(F,V|z)) calculated from the independent 3333 samples for Subj90 4.12.a Posterior distribution of Ke and Ka for the three individuals with population distribution. Result from RIMH sampler 4.12.b Posterior distribution o f V and F for the three individuals with population distribution. Result from RIMH sampler 4.12.c Posterior distribution of drug concentration and AUC at t=l 8 hours for the three individuals with population distribution. Result from RIMH sampler 4.13.a SRQ plot for Subj90,40 subchains from 1000 samples 4.13.b SRQ plot for Subj90, 161 subchains from 5000 samples 4.14 Posterior distribution o f Ke and drug concentration at t=l 8 hours calculated from 5000 samples for Subj90 with population distribution to consider different c 5 41 41 45 47 48 49 51 52 54 55 56 59 60 61 61 63 64 65 67 67 71 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. ABSTRACT Two algorithms, independence Metropolis-Hastings ( IMH ) sampler and regenerative independence Metropolis-Hastings ( RIMH ) sampler, have been implemented and applied with the drug busulfan to assess their performance and integrity. One compartment model with four parameters was used to fit the drug concentration data. Both o f them produced clinically meaningful and consistent results. Prior and posterior distributions and joint densities of parameters from the two algorithms were compared with each other and the performance and integrity of the samplers were evaluated. The convergence to an identical sequence o f the IMH sampler was verified with different initial values, and the rate of mixing was evaluated with autocorrelation plots. The regeneration o f the RIMH sampler was confirmed with scaled regeneration quantile plots. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. I 1. INTRODUCTION 1.1 Pharmacokinetics Pharmacokinetics ( PK ) is the study of the quantitative actions o f pharmacological substance and a biological organism. It is defined as the kinetics of absorption, distribution, metabolism, and excretion o f drugs and their corresponding pharmacologic, therapeutic, or toxic response in animals and man ( Riegelman et al. 1973 ). The field o f pharmacokinetics arose as techniques for interpreting curves o f changing blood concentration with time. The word “pharmacokinetics” was first used by F. H. Dost in his book, Der Blutspiegel-Kinetic der Konzentrationablaufe in der Krieslaufflussigkeit, published in 1953. There had been forerunners to Dost who measured blood concentration curves and derived formulae and laws to quantitatively explain them. Dost, however, was the first one to formulate rules applicable to drugs in general, while those workers had applied their knowledge to the solution of one or two problems ( Gladke 1988 ). He regarded the clearance of a drug from the body as a process of concentration dependent elimination, which could be simply and accurately expressed as an exponential function. This concept became a fundamental and basic tool in studying most of the pharmacokinetic systems. Along with pharmacodynamics ( PD ), which is concerned with the time relationships between the concentration of a pharmacological substance and the intensity o f its effect, pharmacokinetics also plays an important role in efficient drug Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 2 development by providing the necessary information to design effective and safe dosage regimens ( Peck 1993 ). 1.2 Modeling Approaches in Pharmacokinetics There are three major modeling approaches in pharmacokinetics: compartment modeling, physiological modeling, and model-independent pharmacokinetics ( Welling 1986 ). In the compartment modeling approach, the body is assumed to be made up of one or more compartments. These compartments may be spatial or chemical in nature. In most cases, the compartment is used to represent a body volume or group of similar tissues or fluids into which a drug distributes. Within each compartment, the drug is assumed to rapidly distribute into a homogeneous fluid volume, so that we can focus on the overall input-output relationships. Drug movement between compartments is based largely on reversible, if between spatial compartments, or irreversible, if between chemical compartments, first-order linear processes. In the physiological modeling approach, pharmacokinetic interpretation is based on known anatomical or physiological values. Unlike the compartment modeling approach, drug movement is based on blood flow rates through particular organs or tissues and experimentally determined blood-tissue concentration ratios. Model-independent pharmacokinetics is a less complex approach based purely on an empirical description of blood or plasma profiles of drug or metabolites, and calculation of useful pharmacokinetic values without invoking a particular model. This approach is Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 3 useful in situations where kinetic parameters such as absorption and elimination rates and clearances are required, but specific distribution characteristics are less important. The approaches developed in this thesis are applicable to both compartment and physiological models of pharmacokinetic system. 1.3 A Modeling Framework for Pharmacokinetic Systems It is assumed that the pharmacokinetic process under study can be described by the following set of equations ( D’Argenio and Schumitzky 1997 ): = f(.x(t),a,r(t),t), *(0) = c. (1.1) at y{t)=h(x{t),a,r{t),t), ( 1.2 ) where jc(/) is a vector of compartment amounts, a is pharmacokinetic parameter vector, r(/) is input vector, and y{t) is output vector which usually represents drug concentration in a compartment. Equations (1 .1 ) and (1.2 ) are referred to as the state and output equations, respectively. Each element of the input vector r(/) is used to represent the dose regimen of a specific compound. One example is the intravenous ( IV ) infusion regimen in which the rates of infusion remain constant between input times: ry(0 <*,-i < t^d tit i = 2,...,nd + 1 , ( 1.3 ) where nd is the number of dose events. Another example is the bolus-type inputs b(t) which produce instantaneous changes in model states: Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. x(dt*) = x(dt-)+b{dt()y i = 1 ,...,nd, ( 1.4 ) where the - and + indicate times immediately before and after bolus administration. Data to be used in fitting the model described by ( 1.1 ) and ( 1.2 ) are generally collected at discrete times and subject to additive output error: z{tj) = h(x(tj),a,tj) + v(fy), j = 1 ,...,m , ( 1.5 ) where h(tj) represents the model output vector at time ( from ( 1.2 ) ) and v(fy) is the vector o f additive output error terms, representing all sources o f uncertainty including measurement error that can be modeled as a random process added to the true model output directly or by transformation. The common method is to model v(f.) as normally distributed with zero mean and variance defined by an appropriate variance model. In general, the variance of the output error can be an arbitrary function of the model output: var{vi(tj )} = gi(y(t),a,tj ,p ), / = 1,...,/, y = l,...,m , (1.6) where / is number o f outputs, and /? represents additional parameters that are unique to the error variance model. For example, the standard deviation o f error can be related to the output either linearly as var M O } = ( o - ^ + crtlopey(t))2, (1 .7 ) or the variance of error, in power function as va/*{v(0} = 0 ’ o V ( 0 - (1.8) In cases o f multiple responses, each output can be defined by a unique variance model which is independent o f the others. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 1.1 One compartment model used in this study. For example, the model in Figure 1.1 is the one that has been used in this study. It is a one compartment linear model with first-order elimination. The model input is via IV infusion and first-order absorption. The model differential equations and the output equations are: ^ P - = -K 'X «) + K'XA<) + K 0 . (1 .9 ) at ± £ i = -K ax2(t), (1.10) “ j p = *•(')’ n . n ) = <112) y2(t) = Fx>(% , (1.13) where Ke is elimination time constant, Ka is absorption time constant, F is fraction of drug absorbed, V is drug administration volume, and x,(r) and x2(f) represents amounts of drug in compartments 1 and 2, respectively. yx (t) is the drug concentration in compartment 1, and y2(t) is the area-under-the-concentration curve ( AUC ), which can be easily obtained by integrating the concentration curve. x3(r) is a temporary variable Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 6 for calculating the AUC. r(t) represents the rate of infusion of drug into compartment 1. Bolus doses are as indicated in ( 1.4). 1.4 Specific Aims Estimating unknown parameters in a model has been a persistently pursued theme in engineering literature. In pharmacokinetics, precise knowledge of the model parameters has significant meaning both in quantitative description of the system under study and in practical clinical value. But with only limited amount of data available and high variability both within and between patients, pharmacokinetic parameter estimation imposes particular difficulty to investigators. In this context, the Bayesian estimation method has a potential advantage in its ability to incorporate data with previously available information in the form of likelihood and prior distribution. The objective of this study is to apply the Bayesian approach in pharmacokinetic parameter estimation framework through Markov chain Monte Carlo methods. In order to accomplish the goal, the following aims were addressed: 1. Regenerative independence Metropolis-Hastings algorithm was implemented and tested along with the independence Metropolis-Hastings algorithm. 2. Both algorithms were applied to the drug busulfan to identify pharmacokinetic parameters and to provide the necessary clinical information in designing safe and effective drug administration regimen. 3. Performance and integrity of the two algorithms were verified. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 7 In the following chapters, basic theories will be summarized and simulated clinical results will be presented to evaluate the usefulness and merit of the Bayesian estimation method in this field. Chapter 1 introduced and summarized pharmacokinetic modeling techniques and presents a specific model to be used in the application. Chapter 2 is an overview of the Bayesian approach in statistical problems. Chapter 3 has been devoted to explaining the basic concept of Markov chain and two algorithms that have been implemented and used in this work. Application to the drug busulfan and the results from two algorithms are presented in chapter 4. Discussion and concluding remarks are given in Chapter 5. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 8 2. BAYESIAN ESTIMATION 2.1 Bayesian Approach vs. Frequentist Approach In the Bayesian approach to parameter estimation, the model parameters are considered to be random variables and the data, fixed values. This is the fundamental difference from the frequentist approach where parameters are treated as fixed ( Carlin and Louis 1996 ). The frequentist evaluates procedures based on repeated sampling, assuming an infinite replication o f the inferential problem and evaluating unknown parameters over this repeated sampling framework. The Bayesian requires a sampling model represented in the form o f likelihood, and a prior distribution on parameters. Unknown parameters are considered random and all inferences are based on their distribution conditional on observed data, which is the posterior distribution. To summarize, the frequentist conditions on parameters and replicates, or integrates, over the data; the Bayesian conditions on the data and integrates over the parameters. The well-known Bayes’ Theorem connects the prior with the posterior through likelihood as follows: p{a\z) p(z\a)p(a) jp(z\a)p(a)da l(a\z)p(a) jl(a \z)p (a )d a ’ Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 9 where a is parameter, z is observation or data, p(a) is the prior distribution, p(a\z) is the posterior distribution, and l(a\z) is the likelihood which depends on the model definition. Thus, the Bayesian method considers the contribution from both the experimental data in the form o f the likelihood and the prior opinion in the form o f the prior distribution. In application to this work, the prior distribution for the parameter is assumed to be lognormally distributed with mean and covariance calculated from previous patient data. Measurement error is modeled as normally distributed with zero mean and variance model of the form ( 1.6 ) as follows: where y(t) is simulated output with chosen parameter vector sampled from the prior. A primary difficulty in the application o f Bayes' Theorem is the need to perform extensive numerical integration. The associated numerical problems increase in severity as the number of parameters in the model increases. Fortunately, this problem has been largely solved in recent years due in part to the development of methodologies that exploit modem computing power. In particular, sampling-based methods for estimating and computing posterior distribution plays a central role in modem Bayesian analysis. (2 .3 ) cr(0 = 0.1 + 0.15y(0, (2 .4 ) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 10 2.2 Point Estimation and Interval Estimation Having specified the prior distribution, we may now use the Bayes’ Theorem to obtain the posterior distribution of the model parameters. Since the posterior distribution summarizes our current state of knowledge arising from both the observed data and our prior opinion, we can use the Bayesian analogues for the frequentist techniques of point estimation and interval estimation to summarize the posterior distribution. In order to obtain a measure of the accuracy of a point estimate a(z) of a , we might A use the posterior variance with respect to a (z ). Writing the posterior mean Ealf(a) = //, we have that £ alr(o r-a(z))2 = ( a - p + f t - a(z))2 (2.5) = Eab[(a - n )2 + 2 (a - //)(// - a(z)) + (// - a(z))2 ] (2.6) = varalI (a) + 2 [£ alI (a) — //][// — a(z)] + (// - a(z))2 (2.7) A = varalI(ar) + (//- a ( z ) ) 2. (2.8) Hence we have shown that the posterior mean Eail(a) = // minimizes the posterior variance with respect to a(z). Furthermore, this minimum value is varalz(a) = Eaiz[a - Eajz(a)]2, the posterior variance of a . This is a possible argument for preferring the posterior mean as a point estimate, and also partially explains why historically only Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1 1 the posterior mean and variance were reported as the summarized results of a Bayesian analysis. The Bayesian analogue of a frequentist confidence interval ( Cl ) is usually referred to as a credible set. More formally, Definition 2.1 A lOOx (l-2 y )% credible set for a is a subset C o f A such that 1 - 2y < P(C\z) = Ja p{a\z)da, (2 .9 ) where integration is replaced by summation for discrete components of a . This definition enables direct probability statements about the likelihood of a falling in C : ‘The probability that a lies in C given the observed data z is at least (1 - 2y) T h i s is in contrast to the usual frequentist interpretation of Cl: ‘If we could recompute C for a large number of datasets collected in the same way as ours, about 100 x (1 - 2 y)% of them would contain the true value of a ' Therefore, for the Bayesian, the credible set provides an actual probability statement, based only on the observed data and the prior opinion. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 12 2 3 Computational Methods The Bayesian approach to statistical modeling requires specifying the likelihood l(a\z) and the prior probability function p(a) , multiplying the two o f them, and normalizing to form the posterior probability function. The technical problem of carrying out a Bayesian analysis therefore reduces to that o f performing or approximating a number of numerical integrals that usually turn out to be high dimensional. Several approaches have developed in response to this computational challenge ( Smith 1991 ). 23.1 Analytic Approximation For large sample sizes, suppressing the dependence on the given data, the posterior distribution of a is often approximately N(a,H), where a is the maximum likelihood estimate and £ is the inverse of the hessian matrix evaluated at a . This approach has the powerful advantage in that practically all posteriors can be treated with normality A A assumption and that a and £ can be easily computed. As such, major problem with this strategy is the justification of approximate normality. 2.3.2 Numerical Integration The numerical integration and marginalization of p{a\z), as given in ( 2.2 ), involves the computation of integrals of the form S,[q(a)] = J q(a)l(a\z)p(a)dar , ( 2.10 ) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1 3 where / is some index set / c {1,...,£} with a ,,...,a t denoting the components of a , and a r denotes the vector of components of a whose subscripts are not elements of / . Thus for example, 1 = 0 and q(a) = 1 corresponds to the normalizing constant in ( 2.2). It is well known that in the case of the k -dimensional integrals, ■S[?(o0] = J q(a)l(a\z)p(a)da, ( 2.11) efficient quadrature formulae are available if the integrand in ( 2.11 ) can be well approximated by a function of the form g(a) = h{a)n(«p,I\a), ( 2.12 ) where n{< pyl\a) denotes a k -dimensional density with known mean < j> and known covariance matrix Z , and h(a) is a polynomial in the components o f a . In such cases, (2.11 ) can be reduced to the form J g(a)da = J h * (£) exp {-(£ +.. .+£ )}d£ (2.13) where h *(£) is a polynomial in £ and given by a set of following transformations: (2.14) ¥ i = a i +YdSijy/j , / = 2,...,£, (2.15) y -1 Si{ = -cov(aiyq/j\z) / var(if/j\z), ( 2.16) (217) where fi(, o] are the means and variances o f the if/(. The class of functions which can be well approximated by polynomial x normal forms is rather rich and covers many of Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 14 the S,[q(a)\ integrands we typically encounter, provided the individual parameter components have support ( - 00, 00). This restriction can be overcome by redefining the likelihood and prior in terms of transformed individual components, for example, by working with the logarithm of variance components and the logits of proportions. The remaining issues are how to perform individual component transformation and how, after that, to approximate the unknown S ~, and cr,. The importance sampling approach to numerical integration is based on the observation that, if / and g are density functions, where G denotes distribution. This method suggests a statistical approach o f generating a sample from distribution and using the average of the values o f the ratio / / g as an J f(x)dx = S[Xx)/g(x)]g(x)dx = j[Xx)/g(x)]dG(x) = Ea[f(X )/g(X )], (2.19) ( 2.20) (2.18) unbiased estimator o f . Since the variance of such an estimator depends critically on the choice of G , it is desired that g should be similar to / . Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 15 23.3 Iterated Sampling and Resampling Approaches Given random variables o f X , Y, and Z , stochastic substitution algorithm estimates density of X by successive stochastic simulation of random variates drawn from the conditional distributions in the following equations: where brackets denote densities and marginalization by integral is represented by forms such as [X] = J [ X\Y]*[Y] . The algorithm proceeds as follows: 1. Draw X (0) from an arbitrary [X]0. 2. Draw Yi0r, Zm from [Y,Z\Xm ]. 3. Draw X (0 Y , Yil) from [X,Y\Zm ]. 4. Draw X ° \ Zw from [X,Z\I*0]. 5. Iterate. After t steps, with m replications o f the process, we obtain , j = l,...,/n . An estimate o f the marginal density [X] is then provided by ( 2.21 ) ( 2.22 ) (2.23 ) (2 .2 4 ) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 16 In application to Bayesian inference, [X> Y,Z] would denote the posterior distribution o f unknown quantities of interest. Gibbs sampling is a markovian updating scheme, and it is a variant o f the Metropolis algorithm which will be discussed in the next chapter. Gibbs sampler is explained in section 3.3 along with other Markov chain Monte Carlo algorithms. Resampling approach is for the case where a sample from h(a) is needed, but what is available is a sample o f random variates from a density g(a) absolutely continuous with respect to h(a) along with / ( a ) which is normalizable such that h(a) = f (a) / J / (a ) d a . Approximate resampling from h(a) can be done as follows: 1. Given or,., i = 1,...,», a sample from g , calculate = / ( a , ) / g ( a ,). m 2. Then calculate qt = ^ . y-i 3. Draw a* from the discrete distribution over { a r, , p l a c in g mass qt on ar Then a * is approximately distributed according to h with the approximation improving as n increases. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 17 3. MARKOV CHAIN MONTE CARLO METHODS 3.1 Monte Carlo Integration Monte Carlo integration evaluates expectation of a certain quantity by drawing samples form the distribution and then estimating the population mean by a sample mean ( Gilks et al. 1995 ): £ [ / W ] = - ! ■ £ /( * ,) . (3.1) n t-l If the samples are independent, laws of large numbers guarantee that the approximation can be made as accurate as possible by increasing n . However, drawing samples independently from distribution may not be feasible. In that case, it suffices if we can draw samples in correct proportions so that the ensuing distribution of the samples well approximates the original distribution. If we can generate a random variable which has such a sequence as the value of the current random variable depends only on its most recent value, then as time goes by, that sequence will eventually forget its initial value and the process will reach an invariant distribution. This is the stationary distribution of Markov chain and this kind of dependence is called Markov property, whose mathematical definition will be given in the next section. Thus, after a burn-in period of m samples, the samples will be dependent samples from the distribution that we are interested in. After discarding those samples during burn-in, we have ergodic average of Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 18 f — — Z / w - n in ± m m + | (3 .2 ) Convergence to the required expectation is ensured by the ergodic theorem. 3.2 Markov Chain and Stationary Distribution In the following sections, formal and mathematical definitions o f Markov chain and its properties are summarized. 3.2.1 Basic Concept and Definitions A family of variables which evolve as time passes are called stochastic processes. The stochastic process X is called Markov chain if it satisfies the Markov property for all n > 1 and s e S , where S is some countable set called the state space ( Grimmett and Stirzaker 1982; Medhi 1994). There is a simple dependence that the outcome of the n th trial depends directly on the (n - 1) st trial and only on it. If Xn = s, then the chain is said to be in the s th state at the n th step. The evolution of chain is described by its transition probabilities: which is the probability of the chain to move from state i to state j . In a more general case, m -step transition probability is defined as (3 .3 ) (3 .4 ) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 19 ( 3.5 ) It is the probability that, having started from the state i at the n th trial, the state j is reached in m steps. The transition matrix P = (ps) is the |5 |x|^| matrix of transition probabilities ( 3.4) and has the following properties: for all n , i , and j . For a homogeneous chain, the transition probabilities do not depend on specific time points. Henceforth, we will confine our consideration only to homogeneous chains. The period d(i) of a state i is defined to be the greatest common divisor of n for P„{n) > 0. It is the number of steps after which the state returns to itself. If d{i) = 1, the state i is said to be aperiodic. State i is called recurrent if (3 .6 ) (3 .7 ) j The chain X is called homogeneous if P(X,+ l = j\ X m =i) = P(Xl = j\XQ = /) (3 .8 ) P(X„ = i for some n > 1 | X0 - i) = 1, (3 .9 ) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 20 that is, having started from state i , the probability of returning to the same state is 1. If this probability is less than 1, then the state is called transient Let = P(Xi * j , X 2 * j \ X n = j\XQ = /) ( 3.10) be the probability that starting from i , the first visit to state j takes place at the n th step. Then, the mean recurrence time //f of a state is defined as 0 0 Mi = ' Z n f « W if i is recurrent (3.11) = oo if / is transient. ( 3.1 2 ) The recurrent state i is called positive-recurrent if //f < oo, or null-recurrent if fi{ - < x > . Lastly, a set C of states is called irreducible if i j for all i,j eC, where < -» is used to indicate that, having started from state /, there is a positive probability of the chain to visit state j , and vice versa. 3.2.2 Stationary Distribution In order to discuss the existence of stationary distribution of Markov chain, some of those terms defined in the previous section with respect to state and set of states need to be extended to transition matrix, and, equivalently, to chain. The most important terms are aperiodic, irreducible, and positive-recurrent. For the sake of generality, the definitions for non-negative matrices ( Seneta 1981) have been restated. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 21 Definition 3.1 An n x n non-negative matrix P is irreducible if for every pair i,j of its index set, there exists a positive integer m = m(i3 j) such that p~M ) > 0. An irreducible matrix is said to be cyclic ( periodic ) with period d , if the period of any one ( and so of each one) of its indices satisfies d > 1, and is said to be acyclic ( aperiodic ) if d = 1. Positive-recurrence of a chain ensures the existence of a stationary distribution. In other words, it implies that if a value is sampled from stationary distribution, then all subsequent iterates will also be distributed according to the stationary distribution. If a chain is irreducible, aperiodic, and positive-recurrent, then it has stationary distribution jc such that jt. > 0 for all j , and = 1 , ( 3.13 ) j 7lP = K . (3.14) If (3 .14) is iterated, xP2 = (nP)P = nP = it, so that nP" for all n > 0 . (3.15) When a Markov chain has a stationary distribution, it is also the limiting distribution of X m as n -> oo. For example, let’s consider the transition probabilities of p u = 05000, pl2 = 05000, p2l = 0.7000, and p2 2 = 03000. The transition matrix is P = 05000 05000 0.7000 03000 (3.16) It can be verified that this transition matrix satisfies the condition for existence o f stationary distribution. One way to find this stationary distribution is to keep on Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 22 multiplying the transition matrix by itself until we have constant limiting values for each element. That is, 0.6000 05600 05800 05800 05840 05824 05832 05835 05834 05833 05833 05833 Therefore, the stationary distribution of this Markov chain is n = [05833 0.4167]. In other words, if we keep on updating states according to the transition probabilities given as in (3.16 ), after a sufficient number of transitions we’ll be at state 1 with probability of 0.5833 and at state 2 with probability of 0.4167. Equation (3.15 ) can be easily verified by multiplying this stationary distribution with any of the transition matrices from ( 3.16) to (3.22 ) to obtain the resultant 2 x 2 matrix with k duplicated in each row. Analytically, stationary distribution is defined as * , = ! / / / „ (3.23) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0.4000 0.4400J 0.4200' 0.4120 J 0.4160' 0.4176 J 0.4168] 0.4165 J 0.4166] 0.4167 J ’ 0.4167] 0.4167J” (3-17) (3.18) (3.19) (3.20) (3.21) (3.22) 23 where ftt is given as in (3 .11 ). The existence of stationary distribution is closely related with the spectral radius of transition matrix. According to Perron-Frobenius theorem, transition matrix o f an irreducible and aperiodic chain has an eigenvalue of one, and the other eigenvalues are less than 1 in magnitude ( Seneta 1981 ). Therefore, if we write the multiple powers of transition matrix with eigenvectors and engenvalues using spectral theorem, those terms which should be multiplied by eigenvalues whose magnitudes are less than 1 vanish as the number of power becomes large. The only row that survives is the stationary distribution of the transition matrix. If a transition matrix is periodic with period d , it has d complex roots of unity. The example transition matrix of ( 3.16 ) has eigenvalues of I and - 0 .2 . 3.2.3 Construction of a Markov Chain In practical situations where we have to deal with infinitely many states, neither we can define probabilities to every possible move, nor can we verify that a certain Markov chain has stationary distribution by observing the limiting behavior of transition matrix or by an analytical analysis. Instead, in order to construct a Markov chain with stationary distribution, we make the transition probability satisfy the following reversibility condition ( Hastings 1970 ): 7T^. = 7tjPjj for all i and j . ( 3.24 ) This property guarantees that nP = P, for Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 24 Z * i P y = Z * J P fi = * I Z ^ = * > ’ ( 32 5 > / / / Hence, this k is the stationary distribution for P. In order to construct a Maikov chain with this method, we first define Ps =9sfif '*)< (3.26) P, = l - 2 > , . (3-27) where Q = {q^} is the transition matrix of an arbitrary Markov chain with states 0,1,...,5, and fi9 is given by (3 28 ) where s~ is a symmetric function of i and j , chosen so that 0 < f3~ < 1 for all i and j . It can be verified that (3.26) and (3.28 ) satisfy the reversibility condition ( 3. 24 ). One choice for s( > . in ( 3.28 ) is S ^ l + H illL * > l qH > 1 ), ( 3.2 9) k , / q.. n . I q.. = 1 + —-—— ( J- < 1). (3.30) Then, K . / q .. fi, = l ( - ^ - > 1 ) , (3.31) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. In order to simulate this process, we 1. Start with X(n) - i , and select a state j using Q , 2. Take X(n +1) = j with probability and X{n +1) = / with probability Very often, the stationary distribution to be realized is called target distribution, and the q(.\XK ) which the candidate samples come from is called proposal distribution. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 26 3J MCMC Algorithms The Metropolis algorithm ( Metropolis et al. 1953 ) is a classic example of applying the idea of sampling from the stationary distribution o f Markov chain. It considers only symmetrical proposals, having the proposal distribution of the form q(j\i) = q(i\j) for all i and j . Their problem was to calculate the properties of a substance which may be considered as composed of interacting individual molecules. Classical statistics was assumed, only two-body forces were considered, and the potential field of a molecule was taken to be spherically symmetric. They had to calculate the equilibrium value of F according to the equation - \ Fexp(— E / kT)d2N pd2Nq jexp(— E / kT)d2N pd2N q ’ K ' } where E is the potential energy of the system, k is Boltzmann’s constant, T is the absolute temperature o f the surroundings, and d2N p d 2Nq is a volume element in the 4N -dimensional phase space. In using Monte Carlo integration method, instead of choosing configurations randomly and then weighting them with e x p (-£ / kT) , they chose configurations with a probability exp(-E / kT) and weighted them evenly. Selection of configurations with a probability exp(— E / kT) was done as follows: 1. Move a particle to a new position. 2. Calculate the change in energy o f the system, AE . 3. If AE < 0 , allow the move and put the particle in new position. If AE > 0 , allow the move with probability exp(-AE / kT) . Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 27 For the purpose of taking averages, they always considered the particle was in a new configuration after each iteration. Then, after many moves the ensemble tends to a distribution which is proportional to exp(-E / kT). A special case of the Metropolis algorithm is random-walk Metropolis, for which q(J\i) = q(\i - y'|). The independence sampler is a Metropolis-Hastings algorithm whose proposal does not depend on / , i.e., q(J[i) = q(J) • Often it is more convenient and computationally efficient to divide the variables into components o f possibly differing dimensions, and then update these components one by one. This method is referred to as single-component Metropolis-Hastings. A special case of single-component Metropolis-Hastings is the Gibbs sampler ( Smith 1991 ). Gibbs sampling deals with a collection of random variables X x , X 2,..., X t , for which it is known that the joint density, [Xlt X 2>..., X t ] is uniquely determined by the full conditional densities [X,\Xr>r*s\, s = 1, 2 ,...,k . It is required that samples of X, can be generated straightforwardly and efficiently from [Xs\Xr,r * s] , given specified values of the conditioning variables, Xr, r * s . The Gibbs sampling algorithm proceeds as follows: 1. Given an arbitrary starting set of values X,(C \ ..., X { °], we draw A ,, < 1 ) from 2. Then draw X*" from [ * 2|* I (I),* j 0), . . . , ^ 0)] . 3. And so on up to drawing X from [Xt \X\x ) Ar J [ '_ > l] to complete one iteration. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 28 After t such iterations, we would arrive at (X f°y...,X];,}) . It is shown that under mild conditions ( Geman and Geman 1984 ), X™ —* - ► X, ~ [X,] as r co. (3.34) Thus, for t large enough, we can regard X ( p as a simulated observation from [Xs]. Replicating this process m times produces m independent and identically distributed k - tuples (Xf'},...,X u ) , j = l,...,m . For any s , the collection Xf(},...,X"} can be viewed as a simulated sample from [Xs] . The marginal density could then be estimated by the finite mixture density, m = = (3.35) m j.l Gibbs sampler candidates are always accepted and Gibbs sampling consists purely in sampling from full conditional distributions. MCMC sampler can be used in combination with other method into hybrid algorithm. One example is the regenerative independence Metropolis-Hastings algorithm where independence Metropolis-Hastings algorithm is combined with rejection sampling. Among the several MCMC algorithms mentioned here, independence Metropolis- Hastings algorithm and regenerative independence Metropolis-Hastings algorithm have been used in this study. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 29 3.4 Independence Metropolis-Hastings ( IMH ) Sampler The IMH sampler applies to cases where there exists independence between the current state and the candidate. 3.4.1 Independence and Acceptance Probability If independence between i and j exists, q~ = q. and qjt. = q .. This can be accomplished by drawing samples independently of the current state. In this case, we have from (3.31 ) and ( 3.32), If we correspond the posterior density p{a\z) to target density n and the prior density p(a) to proposal density q , we have according to Bayes’ Theorem ( 2.2 ), (3 .3 6 ) Jtjq., Ki lqi (3 .3 7 ) k p(a\z) _ Ka\z) 1 . 1 P(a ) | l(a\z)p(a)da' (3 .3 8 ) so that A,=l (3 .3 9 ) Kj\z) (/(>!£)<!) im im (3 .4 0 ) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 30 Thus, we don’t have to deal with the normalizing constant term in the denominator of ( 2.2 ), which is the main advantage of applying MCMC methods in Bayesian approach. This likelihood ratio is also the acceptance probability for the candidate samples generated independently from the prior distribution. 3.4.2 IMH Algorithm Below is the algorithm that accommodates the above idea and has been used and tested in application section o f this thesis: 1. Select a 0, initialize chain at a 0. 2. For n = do 3.-5. 3. Generate j from p(a) when currently at a n = i . 4. Generate u from C/(0,1), where C/(0,1) is uniform distribution between 0 and 1. 5. If u < - ^ 4 , then a„+ l = j , else a „ x = i . Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 31 3.5 Regenerative Independence Metropolis-Hastings ( RIMH ) Sampler Another algorithm, RIMH sampler, produces a regenerative Markov chain, which provides a way for error analysis for the Bayesian estimation. A process is called regenerative if there is a sequence of random times at which the process starts over independently and identically. The tours of the process between these times are independent and identically distributed ( iid ). One easy way to introduce regeneration into an MCMC sampler is to use a hybrid algorithm in which independence Metropolis- Hastings sampler o f previous section is combined in a mixture with an independence kernel where the proposal density is sampled by a rejection sampling algorithm. 3.5.1 Rejection Sampling and Acceptance Probability Suppose we wish to sample from a probability density function ( pdf) / but have a way to sample from the pdf g . Rejection sampling method retains the sampled values Y from g with a probability depending on K ( Ripley 1987). The procedure is as follows: 1. Generate Y from pdf g . 2. Take X = Y with probability h(Y), else iterate. Then ^fhas pdf proportional to gh. For P(Y< x and Tis accepted) = \lcag{y)h{y)dy, (3.41) and P(Y is accepted) = g(y)h(y)dy, ( 3.42) so Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 32 Ifrt-, g(y)h(y)dy P(Y<x\Y is accepted) = , (3.43 ) \Z,siy)hiy)dy which shows the accepted values have pdf ghfjgh. Since we want samples from f we take h — f / g M , provided that there exists a constant M such that f I g < M < oo. Then P(Y is accepted) = \g h = \ ( J I M) = (\l M )\ f = 1/ M, (3.44) so (3.43 ) can be rewritten as P(YSx\Y is accepted) ---- - 4 — = — 7— = / , (3.45) Igh gM] gh M]gh which verifies that X has pdf / . Employing the rejection sampling in Bayesian method requires relating pdf g to the prior distribution p(a) and h to the likelihood l(a\z). Then we can sample from / which corresponds to the posterior distribution p(a\z), for f = l L = ^ l r J f piaU), (3 .46) )gh I l(a\z)p(a)da according to ( 2 . 2 ). One problem with rejection sampling is that large value of M leads to an inefficient algorithm with many rejections. The art of using the rejection method is to find a suitable pdf g that matches f well and from which it is easy to sample. In order to introduce regeneration, we need to fix a particular state and consider the times at which the chain returns to this state ( Tierney 1996). In order to do this, we Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 33 * * define a set C such that l(a\z) <c , where o<c<l(a |z) and a = arg{max /(a|z)}, the maximum likelihood estimate of a . Then, with / as the current state and j as a candidate, we form acceptance probability of for i < = C, (3 .4 7 ) for i « C, j e C, ( 3.48 ) for i £ C, j e C. (3.49 ) Thus, whenever the current state of the chain is in C, any candidate is accepted, and the next state is generated as a draw from the distribution / , no matter what the particular state within C is. The visits to C therefore represent regeneration times. The constant c represents relative contribution from IMH sampler from rejection sampler. If c = 0, the algorithm reduces to IMH sampler, and if c = /(or |z) , it is purely rejection sampling scheme. The visits to the set C form a renewal process that can be used for a regenerative analysis. 3.5.2 Regenerative Process A regenerative process X has a sequence of regeneration points N x, that are random stopping times at which the future of the process is independent of its past ( Shedler 1993 ). The tours between regeneration points are then independent, and the sequence of regeneration points is a renewal process. Set 1 /(/|z) ’ K M l{i\z) ’ Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 34 = ( 3.5 0) which is the length o f each tour, and W ) = £ / ( * ,) ■ (3.51 ) the total observation on observed tour. Then the definition of regenerative process ensures that the sequence {(vk,Yk(f)):k > 1 } consists of iid random vectors. Define the quantity m j ) - = £ [/(* )]■ (3.52) and * » ) = - ! > , . (3.53) n *» 1 Y(r>) = - Y . Y l (f). (3.54) " t - l To obtain an asymptotic confidence interval, set Zk(f) = Yk(f) -9 ( f) v k. (3.55) We see that the k th cycle of the process {X{n):n > 0} determines the quantity Zk( f ) . Furthermore, ( 3.55 ) implies that £ [Z * (/)] = 0- (3-56) Then according to the central limit theorem for iid random variable with zero mean and finite variance, we have Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 35 ' r £ Z . ( / ) ~ JV(0,1) as * - » C O , ( 3.57 ) c r (/)V /i *-i where cr2 ( / ) = v a r[Z ,(/)]. With a little mathematical manipulation, (3.57) can be equivalently rewritten as ~ # (0,1) as n oo, ( 3.58 ) A where 0 (/i) is strongly consistent estimator for 0(f) , i.e., 0 (ji)« y (ii)/K ii)f (3 .5 9 ) and 0 (n) -»• # ( /) as n -> oo almost surely. Now define a weakly consistent estimator for var[Z, ( /) ] = a 2(f) as L r S ( n - « ( '> K ) i . (3 60 ) n _ * h that is, sf(n) -> <r2( / ) as n -> oo in probability. Then it follows that the interval (3.61) v{n)4n v(n)Jn is an asymptotic 100(1-2 y)% confidence interval for 0 ( f ) . The half-length of the confidence interval is 1 / 4n times a multiple zx _r of the estimate of the standard deviation constant s{f) / £ [v ,], so that as n increases, the length of the confidence interval converges to 0, and the midpoint converges to the true value. For 90% ( y = 0.5 ) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 36 confidence intervals, z,_r = 1.645, and for 95% (y = 0.025) confidence intervals, z,_r = 1.96. As a consequence o f ( 3.58 ), the number of tours required to obtain an asymptotic 100( 1- 2/)% confidence interval with a half-length that is 100<5% of 0(f) is approximately (z \ i-r 2 f <7(/) l I 8 In order to use this method, we need to first do a short pilot run to obtain estimates for <j(f) , 0 ( f ) > and £[ v, ], and determine n from (3.62 ). Then simulate the regenerative process again to get point and interval estimates using ( 3.59) and (3.61 ). 3.5.3 RIMH Algorithm Below is the algorithm that accommodates the above idea and has been implemented and used in application section of this thesis: * 1. Define a = arg{max. l(a\z)}. * 2. Select a constant c such that o<c< l(a \z) . 3. Let C represent the set of a such that l(a\z) < c . 4. Select a 0 in C as the initial value. 5. For n = , repeat 6.-10. 6 . Generate j from p(a) and u from U (0,1) until when currently at c Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 37 a, —i. 7. Generate v from C/(0,1). 8 . For i e C, then = J . 9. For / g C and j e C , i f V < l ( J \ z ) ' 1 1 1 6 1 1 = 7 ’ 6 l S e a ’ + l = ' ' 10. For / g C and j g C , /(y|z) . . ^ a .+ » = 7 . else a„+ I = i . Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 38 4. APPLICATION WITH BUSULFAN 4.1 General Description of the Drug Busulfan is usually administered as a high dose with cyclophosphamide in preparation for bone marrow transplantation ( BMT). However, the pharmacokinetics of the drug is highly variable and this variability can contribute to uncertainty in drug administration regimen. Low systemic exposures result in insufficient myeloablation in * the host to permit engraftment. High systemic exposure can cause severe toxicity such as veno-occlusive disease ( VOD ) of the liver. To achieve maximal targeted systemic exposure for each patient, individualized dose adjustment based on measured drug concentration is necessary. Therefore, it is important to accurately define the pharmacokinetics of the drug for safe and effective clinical administration. Commercially available busulfan formulation must be administered orally as tablets, and the variability in absorption is substantial. As an alternative, the parenteral formulation offers the potential to reduce both inter- and intra-patient variability in plasma concentrations. Moreover, a parenteral formulation avoids the imprecision of administering fractional tablets and vomiting commonly associated with high doses of oral administration. The variability in the busulfan absorption complicates pharmacokinetic monitoring because numerous blood samples are required in order to accurately characterize the Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 39 AUC. Bayesian estimation method which can incorporate previous patient data as the prior information is especially useful in this context. 4.2 Experiment and Simulation Protocol Busulfan is given to patients every 6 hours up to 16 times. The initial dose is parenteral, followed by a bolus dose. The third dose and on up to the last one are parenteral. Five 1-2 ml blood samples are taken at 1.5,2, 3,4, and 6 hours from the start of initial dose. These data are used to adjust the fourth dose. The amount of dose given is as follows: 1. Initial parenteral dose - 37.5 mg/m2 as one hour infusion. 2. Oral second dose - 37.5 mg/m2 rounded to nearest mg. The targeted exposure at the end of the third dose is concentrations between 0.2 and 0.3 (ig/ml. The targeted AUC is at the 60th percentile for the patient group, for it has been confirmed from previous study that there would be no serious incidents of VOD and other major regimen related toxicity with this level of AUC. In simulating the clinical situation, population simulation of 1000 patients was done first to obtain the necessary concentration data and to pick up specific subjects to investigate. The one-compartment model used is shown in Figure 1.1 and the model differential equations and the output equations are listed from ( 1.9 ) to ( 1.13 ). Busulfan dose was given as described above with the body surface area assumed to be 1 m2 , i.e., the initial and the third parenteral dose of 37.5 mg and the second oral dose o f 38 mg. The prior statistics for pharmacokinetic parameters came from previous clinical data Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 40 collected from 20 patients at St. Jude Children’s Research Hospital, Memphis, Tennessee. They are summarized in Tables 4.1 and 4.2. The drug concentration plot and AUC plot from the population simulation o f 1000 individuals are shown in Figures 4.1 and 4.2. The resultant statistics are listed in Tables 4.3 and 4.4. From the population simulation data, 3 subjects were chosen based on AUC, that is, 10th, 50th, and 90th percentile individuals were selected. Five simulated data from these subjects were used to apply Bayesian method using IMH and RIMH algorithms to estimate each patient’s pharmacokinetic parameter and to predict the dose concentration and AUC at 18 hours when the dose amount should be adjusted, along with. The measurement error model and Ke V Ka F Prior mean from patient data 0.299 29.19 1.52 0.7 Table 4.1 Prior mean for param eters from patient data. Ke V Ka F Ke V 0.0852* -0.881 18.282 Ka F -0.005 -2.559 0.7222 0.0253 0.142 Table 4.2 Prior covariance m atrix for parameters from patient data. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 41 I g Figure 4.1 Concentration plot from population simulation of 1000. Figure 4.2 AUC plot from population simulation of 1000. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 42 the variance model assumed are listed in ( 2.3 ) and ( 2.4 ). The 60th percentile AUC was found out to be 4.1579 p.g-hr/ml. Ke V Ka F Prior mean from population simulation 0.3019 28.58 1.521 0.6976 Table 4.3 Prior mean for parameters from population simulation. Ke V Ka F Ke 0.0831* V -0.871 17.092 Ka -0.0023 -2 .6 6 0.7182 F -0 .0 0 0 2 -0.004 0.031 0.1422 Table 4.4 Prior covariance matrix for parameters from population simulation. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Param eter SubjlO Subj50 Subj90 Ke 0.3119 0.3191 0.2951 V 44.5294 27.9277 18.3317 Ka 1.6460 1.5703 1.5009 F 0.7093 0.8691 0.9350 Table 4.5 Individual parameter data obtained from population simulation. Mean values. Time(hour) SubjlO Subj50 Subj90 1.5 0.4390 0.8517 1.4291 2 0.3757 0.7261 1.2331 3 0.2750 0.5277 0.9180 4 0.2013 0.3836 0.6835 6 0.1079 0.2026 0.3788 Table 4.6 Individual drug concentration data obtained from population simulation to be used in Bayesian estimation using MCMC methods. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 44 4 J Results from IMH Sampler Below are the results from Bayesian estimation using IMH sampler. 43.1 General Information A total o f 5000 samples were drawn through IMH sampler for each subject. Figure 4.3 shows the sequence o f samples for the case of Ke. The first 500 samples are shown together with the initial 102 discarded samples in consideration o f bum-in period. The flat portion o f the plot corresponds to the time when the present sample was retained without accepting candidate sample. It turned out that more transitions were associated for the subject with lower value of AUC. Table 4.7 summarizes the posterior sample mean and standard deviation of each parameter calculated from the 5000 samples. Posterior and prior population distribution of each parameter values along with concentration and AUC at 18 hours after the initial dose are plotted together in Figures 4.4.a to 4.4.c. There is not much difference in the means of Ke and Ka and those from the population simulation in Table 4.3. The similarity in the distribution for Ka can be explained by the fact that only the data from parenteral dose which does not require absorption process were used in estimation. Figure 4.4.b shows a decreasing trend in V and an increasing trend in F from SubjlO to Subj90. This is in general true, for less concentration would be associated with the individual with larger drug distribution volume and smaller F. The information in Figure 4.4.c has a significant implication in clinical situation. It gives a distribution of drug concentration and AUC at 18 hours after the initial dose when adjusted dose should be given with consideration of inter-patient Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. K e , S u b j9 0 K e,Sub)SO K e, SubjlO 45 0.6 0.4 0.2 -100 0 100 200 300 400 500 0.6 0.4 0.2 -100 0 100 200 300 400 500 0.6 0.4 0.2 -100 0 100 200 300 400 500 Figure 4.3 M arkov chains formed by IM H algorithm for Ke. Shown only up to 500 samples and the initial 102 discarded samples in consideration o f the bum-in period. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 46 SubjlO Mean Std SubjSO Mean Std Subj90 Mean Std Ke 0.2646 0.06935 0.3025 0.06618 0.3097 0.05983 V 45.20 14.21 24.92 7.078 15.82 3.821 Ka 1.237 0.5657 1.498 0.6856 1.896 0.8420 F 0.6393 0.1149 0.6930 0.1261 0.7745 0.1330 Table 4.7 Posterior mean and standard deviation calculated from 5000 samples with IM H algorithm. Com pare with Table 4.9. pharmacokinetic variability to obtain maximal targeted exposure of the drug. Since the targeted drug concentration is between 0 .2 and 0.3 pg/ml, it is clear that more amount of drug is necessary for Subj 10 to avoid ineffectiveness and less amount should be given to Subj90 to prevent possible toxic effect. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. p(Ka), p(Ka|z) p(Ke), p(Ka|z) 47 0.2 0.15 0.1 0.05 0.4 0.5 0.3 0.6 0.7 0.8 0.2 0.1 0.3 0.25 0.2 0.15 0.1 0.05 Ka(/hr) Figure 4.4.a Posterior distribution of Ke an d K a for the three individuals with population distribution. Result from IM H sampler. Compare with Figure 4.12.a. Solid: prior; dotted: SubjlO; dashdotted: S u b jS O ; dashed: Subj90. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 48 0.5 0.4 100 60 120 140 0.25 0.2 0.05 0.5 F Figure 4.4.b Posterior distribution of V and F for the three individuals with population distribution. Result from IMH sampler. Compare with Figure 4.12.b. Solid: prior; dotted: SubjlO; dashdotted: S u b jS O ; dashed: Subj90. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. p(AUC(t-18)), p(AUC(U18)|z) 49 S= 0.3 t 0.25 0.2 0.15 C O 0.1 0.05 0.8 (ugrtrt) 0.2 0.4 0.6 concentration (ug/tal) 1.2 1.4 0.5 0.4 0.2 0.1 AUC (ug-hr/ml) Figure 4.4.c Posterior distribution of drug concentration and AUC at t=18 hours for the three individuals with population distribution. Result from IM H sampler. Compare with Figure 4.12.c. Solid: prior; dotted: SubjlO; dashdotted: Subj50; dashed: Subj90. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 50 4.3.2 Convergence, M ixing, and Correlation Now that we have a meaningful result from IMH sampler, we need to verify the integrity o f the method. The first one is the issue of convergence. It is desired that the sequence should converge fast to the stationary distribution. Figure 4.5 shows how the chains converged to an identical distribution even with different initial values. Since it is reasonable and easy to choose the prior means for initial values, half, double, and quadruple o f the prior means were also tried. It was verified that it took only few samples for the chains to show identical sequences, meaning that they converged to an identical distribution. In this simulation, all of the four parameters were estimated. Therefore, all the components of variables were updated at each step o f the sequence. Hence, the rate of convergence was quite fast. The results for the cases when one or two of the parameters were fixed are not shown here, but the sampler took longer time to converge under those conditions. A chain may have converged to a stationary distribution, but still be mixing too slowly to be useful for Monte Carlo method. Mixing is the degree o f dependence between samples. Rapidly mixing means that the dependence decays rapidly as a function o f the distance between samples ( Geyer 1992 ). The rate of mixing can be measured by correlation. Figure 4.6 is the autocorrelation plots for a sample sequence. All of the three subjects are shown up to 50 sample lags. All shows a decreasing trend, but it turned out that the subject with extreme values mixes more slowly. For Subj90, the slowest case here, the degree o f autocorrelation became small enough after 30 lags, which can also be verified with the scatter plot in Figure 4.7 that the linear relationship has almost vanished. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 51 b 0.6 ■ g . i> 0.4 o 5* 0.2 -100 -90 -60 -70 -60 -50 -40 -30 -20 0 -10 10 0.6 0.4 * 0.2 -70 -60 -100 -90 -80 -50 -30 -20 -40 -10 0 10 0.6 -100 -90 -80 -70 -60 -50 -40 -20 -30 -10 0 1 0 0.6 0.4 * 0.2 -100 -90 -70 -80 -60 -50 -40 -30 -20 -10 0 10 Figure 4.5 Sequences of Ke started from different initial values for Subj90. Result from IM H sampler. They converged to an identical distribution during the buro-in period. Note the different starting points at y axis. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 52 0.5 -0.5 1 T T T T T T 0.5 •0.5 25 0.5 m -0.5 1 T T T T T T T T 0.5 u. -0.5 30 lag (samples) Figure 4.6 Autocorrelation plots to assess the rate of mixing. Result from IH M sam pler. Shown up to 50 lags. Dotted: SubjlO; dashdotted: Subj50; dashed: Subj90. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 53 In discarding a portion of the initial samples in consideration of the bum-in period, it is not necessary to throw away more samples than what is needed for the autocorrelation level to decay to a negligible level. Thus, discarding the first 102 samples in this simulation is legitimate. The prior joint densities for a few pairs o f parameters are shown in Figure 4.8. Those for Subj90 are shown in Figure 4.9. The degree of covariance between these selected parameters is in a good agreement. These pairs were selected based on the patient data, i.e., the pairs which were found to be correlated with each other. In addition to that, it turned out that there was a high degree of correlation between F and V in the posterior distribution which had not existed in the prior. This is due to the fact that only the concentration data from parenteral dose were used. Therefore, Kax 2( t) in ( 1.9 ) can be ignored, leaving only Kt as parameter, and F and V appear together in close relationship to each other in ( 1.12 ), the output equation. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 54 0.7 0.6 ‘0.5 *0.4 0 3 0.2 30 25 • • • • * _ * * 5 J • % +* • • ••*••• •• “T* >.* - _ • * o 20 > 1 5 10 • *•••*••• * * ; •• • i . • t * • • • T . . • \ • • • • * • * V % + 9 • 2 . • • • ^ * • • v * • « .,¥ • , * • • * • » | • • », • • • • • • • * % * • * . 5 0.2 0.3 0.4 0.5 Ke(n) 0.6 0.7 10 15 20 25 V(n) 30 8 6 2 0 2 8 4 6 0 Ka(n) 1.4 1.2 u. 0.8 0.6 0.4 0.6 0.4 0.8 1.2 F(n) Figure 4.7 Scatter plots for subj90 at the lag of 30. Result from IM H sampler. The linear dependence has become negligible at this lag. From the first 1000 samples. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 55 0.8 0.6 5 0.4 0.2 0 0.8 • •• • • • 0.6 5 0.4 • .• • • 'StE& li.-.' .^ ■ r d * . t » » V * # • * • * #.* * • • 0.2 0 - •"*.V 50 100 V 150 4 Ka 150 100 . . ?v •/. • 150 100 50 1.5 Figure 4.8 Selected prior joint densities ( p(V,Ke), p(Ka,Ke), p(Ka,V), and p(F,V)) from population simulation. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 56 Figure 4.9 Selected posterior joint densities ( p(V,Ke|z), p(Ka,Ke|z), p(Ka,V|z), and p(F,V|z)) calculated from 5000 samples for Subj90. Result from IM H sampler. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 57 4.3.3 Independent Samples To determine how long the chain should be run to obtain a certain accuracy o f the estimation, the method suggested by Raftery ( 1992 ) has been applied to the sequences for Subj90. The method consists o f first forming a binary 0-1 process depending on whether or not the samples from a pilot run fall within the bound of the quantile we want to estimate. Then constructing a new process consisting of every Ath iteration from the original 0-1 sequence gives us an approximate Markov chain for sufficiently large k. The number of bum-in period can be calculated using the result from the two-state Markov chain by making the intermediate transition matrix fall within a certain tolerance of its stationary distribution. To determine the total length of chain that should be run to obtain desired accuracy, we use the asymptotic normality of the estimator of cumulative distribution function by the binary sequence to make the target quantile be within a certain range we specified. The result from the pilot run of 1000 samples implied that the length of the sequence should be up to about 98,000 samples to estimate 0.5 quantile which is the median within 0.0125 quantile with probability of 95%. Therefore, 100,000 samples were generated. The suggestion on obtaining independent samples was to get every 15 to 17th samples to form a new sequence. But every 30th samples were selected to be more conservative and to be in accordance with the analysis on autocorrelation in section 4.3.2. Thus, the new sequence consisted of 3333 independent samples. Table 4.8 lists mean and standard deviation calculated from the three sequences. The slight difference between 5000 samples and the other two seems to result from the fact that median was estimated in the Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 58 cases of 100,000 samples and 3333 samples. The results from the two are not much different from each other, which can also be verified by the posterior distribution plots from Figures 4.10.a to 4.10.C. The degree of correlation between parameters in the 3333 independent samples in Figure 4.11 is similar to that from 5000 samples in Figure 4.9 and also that from 100,000 samples which is not shown here. Therefore, it is not necessary to keep all the 100,000 samples which require a lot of storage space. The inference made from the 3333 independent samples is equivalent in this case. 5000 Mean Std 1 0 0 ,0 0 0 Mean Std 3333 Mean Std Ke 0.3097 0.0598 0.3147 0.0586 0.3144 0.0587 V 15.82 3.821 15.52 3.948 15.60 3.990 Ka 1.896 0.8420 1.851 0.7922 1.844 0.7792 F 0.7745 0.1330 0.7627 0.1397 0.7652 0.1413 Cone, at t=18 0.4280 0.1184 0.4148 0.1134 0.4147 0.1115 Table 4.8 Posterior mean and standard deviation for Subj90. Calculated from 5000 samples, 100,000 samples, and independent 3333 samples originally from 100,000 samples with IM H sampler. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. P(K a). P(Ka|z) p(Ke), p(Ke|z) 59 0.2 0.15 0.1 0.05 0.6 0.4 0.5 0.8 0.2 0.7 0.1 0.3 0.25 0.2 0.15 0.1 0.05 Figure 4.10.a Posterior distribution of Ke and Ka for Subj90 with population distribution to consider independent samples. Result from IM H sampler. Solid: prior; dotted; 5000 samples; dashdotted: 100,000 samples; dashed: 3333 independent samples originally from 100,000 samples. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. P (F ),P (F iZ ) „ P(V).P(V|z) 60 0.5 0.4 0.1 100 120 140 V0) 0.2 0.1 0.05 - 0.5 F Figure 4. lO .b Posterior distribution of V and F for Subj90 with population distribution to consider independent samples. Result from IM H sampler. Solid: prior; dotted; 5000 samples; dashdotted: 100,000 samples; dashed: 3333 independent samples originally from 100,000 samples. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. PMlWW 6 1 Figure 4. lO.c Posterior distribution o f drug concentration at t = 18 hours for Subj90 with population distribution to consider independent samples. Result from IM H sampler. Solid: prior; dotted; 5000 samples; dashdotted: 100,000 samples; dashed: 3333 independent samples originally from 100,000 samples. 0.6 0.5 0.4 0.3 0.2 0.1 • • • 0.6 0.5 * * 0.4 * • * a S C 0 3 0.2 ' : • 0.1 10 20 V 30 40 4 Ka 40 30 > 2 0 10 0 4 Ka Figure 4.11 Selected posterior joint densities ( p(V,Ke|z), p(Ka,Ke|z), p(Ka, V|z), and p(F,V|z)) calculated from the 3333 independent samples for Subj90. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 62 4.4 Results from RIMH Sampler Below are the results of Bayesian estimation using RJMH sampler. 4.4.1 General Information The constant c to define the set C was chosen to be 05 x l{ a \z ) for all the three subjects. Likelihood value information was available from the previous IMH run of 5000 samples. Again, 5000 samples were drawn with RIMH sampler. Sample mean and standard deviation are summarized in Table 4.9. The results from the two samplers are similar, which verifies the consistency of these methods. Generally, sample variances are smaller in the results from RIMH sampler. The posterior and prior population distribution of each parameter and concentration and AUC at t = 18 hours are plotted together in Figures 4.12.a to 4.12.C. They also showed similarity with the results from IMH sampler in Figures 4.4.a to 4.4.c. SubjlO Mean Std SubjSO Mean Std Subj90 Mean Std Ke 0.2707 0.06310 0.3048 0.0564 0.3079 0.04673 V 45.78 12.43 24.22 5.794 15.48 3.455 Ka 1.173 0.5142 1.472 0.6257 1.851 0.7991 F 0.6357 0.1168 0.6935 0.1271 0.7664 0.1366 Table 4.9 Posterior mean and standard deviation calculated from 5000 samples with R IM H algorithm. Compare with Table 4.7. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. p(Ka), p(Ka|z) p(Ke), p(Ke|z) 63 0.2 0.15 0.1 0.05 0.2 0.1 0.4 0.7 0.8 0.5 0.6 Ka (/hr) 0.3 0.25 0.2 0.15 0.1 0.05 Figure 4.12.a Posterior distribution of Ke and Ka for the three individuals with population distribution. Result from R IM H sampler. Com pare with Figure 4.4.a. Solid: prior; dotted: SubjlO; dashdotted: S u b jS O ; dashed: Subj90. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 64 0.5 0.4 100 120 140 V O ) 0.25 0.2 0.05 0.5 F Figure 4.12.b Posterior distribution of V and F for the three individuals with population distribution. Result from RIM H sampler. Compare with Figure 4.4.b. Solid: prior; dotted: SubjlO; dashdotted: S u b jS O ; dashed: Subj90. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. p(AUC(Ul 8)), p(AUC(t-18)|z) p(oonc#-18)). p(oonc(t-18)|z) 65 0.1 0.6 concentration 0.2 0.4 0.8 concentration (ugftnl) 1.4 0.5 0.4 0.2 0.1 14 Figure 4.12.C Posterior distribution of drug concentration and AUC at t=18 hours for the three individuals with population distribution. Result from RIM H sampler. Compare with Figure 4.4.c. Solid: prior; dotted: SubjlO; dashdotted: Subj50; dashed: Subj90. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 66 4.4.2 Regeneration and Asymptotic Confidence Interval The tours o f the chain between regeneration times are independent and identically distributed. This eliminates the chore of throwing away a portion of initial samples in consideration o f the bum-in period as with IMH sampler. Standard errors can be computed using the methods based on iid observations. One way to diagnose a regenerative sampler is to examine the pattern of regeneration times graphically using the so-called scaled regeneration quantile ( SRQ ) plot ( Mykland et al. 1995 ). It is a plot of v, / vn against / / n . If the run length is long enough, then this plot should be close to a straight line because of the asymptotic uniform distribution o f regeneration. Deviations from this straight line imply that the observation period is not long enough for the sampler to have reached equilibrium. Figures 4.13.a and 4.13.b are the SRQ plots for Subj90 from 1000 samples and 5000 samples, respectively. There were 40 subchains in 1000 samples and 161 subchains in 5000 samples. Average lengths of the subchains were 24.82 and 30.75, respectively. There is not a serious deviation from straightness. Table 4.10 shows 95% asymptotic confidence intervals from different run lengths for V of Subj90. As it can be verified with the formula ( 3.61 ), the intervals becomes smaller as the number of tours increase. To decide how long the chain should be constructed to obtain a certain degree of accuracy, we can use ( 3.62 ) to determine the required tour length. Suppose we did a pilot run of 1000 samples. As a result, we have estimates 0<J) = 15.46, A <r(/) = 2232, Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 67 { S Figure 4.13.a SRQ plot for Subj90,40 subchains from 1000 samples. Result from RIMH sampler. i S Figure 4.13.b SRQ plot for Subj90,161 subchains from 5000 samples. Result from R IM H sampler. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 68 # Samples Lower bound Mean Upper bound # Subchains 1000 15.182 15.46 15.739 40 2500 15.268 15.47 15.676 87 5000 15.336 15.48 15.614 161 10000 15.363 15.46 15.565 323 Table 4.10 Asymptotic 95% confidence intervals for V of Subj90 for various sample lengths. E[vy] = 24.82. Therefore, if we want a confidence interval with a half-length o f 0.1 which corresponds to 6.53% of 6(f) = 15.46, the required number o f tours can be calculated according to the equation (3.62 ): f-g -Y f 2 2 3 2 Y — 304, i, 0.00653) U 5.46x24.82j which agrees with the case of 1 0 0 0 0 samples. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 69 4.4.3 Different c So far, the constant c to define the set C has been fixed to 03 x l(a\z ) . The value of l(a '\z) was approximated by the maximum of the likelihood values among the 5000 samples from the previous IMH run. In addition to 03 x l( a ’\z) , two more values, cd ef x l{a\z) with cd ef of 0.3 and 0.7, were tried for Subj90. Generally, smaller values of cd ef brought about longer subchains with reduced total number o f subchains, for they make it more difficult for the chain to return to the set C . The number o f subchains and the average length of them are listed in Table 4.11. Among the four parameters, the posterior distribution of Ke was the only sensitive one to different c . Drug concentration distribution was also affected from this variability. There was a very little difference in the shape of peaks in the cases of Ka, V, and F. For Ke and drug concentration at t = 18 hours, larger values of cdef resulted in smaller variance as shown in Table 4.12 and Figure 4.14. The plots for cdef = 0.3 were more similar to the results from IMH sampler, which conformed to the construction idea o f this hybrid sampling scheme. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 70 c d e f = 0 J cd ef - 0.5 cd ef = 0.7 # sabchalns Average length 130 37.58 161 30.75 187 26.68 Table 4.11 Number of subchains and average length of them calculated from 5000 samples with RIMH algorithm for Subj90 with different c. cdef — 0.3 Mean Std c d ef = 0.5 Mean Std cdef= 0.7 Mean Std Ke 0.3095 0.05266 0.3079 0.04673 0.3041 0.03858 V 15.50 3.685 15.48 3.455 15.58 3.375 Ka 1.858 0.7866 1.851 0.7991 1.848 0.7818 F 0.7635 0.1370 0.7664 0.1366 0.7690 0.1413 Cone, at t = 18 hours 0.4246 0 .1 0 2 0 0.4264 0.09075 0.4306 0.07531 Table 4.12 Posterior mean and standard deviation calculated from 5000 samples with R IM H algorithm for Subj90 with different c. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. p(conc(U18)), p(conc(t»18)|z) 71 0.15 Q- 0.1 0.05 0.7 0.5 0.6 0.1 0.2 0.4 0.8 Ke (ftv) 0.4 0.3 0.2 \\ 0.1 0.2 0.6 0.8 concentration (ugrtni) 0.4 1.4 Figure 4.14 Posterior distribution of Ke and drug concentration at t=18 hours for Subj90 with population distribution to consider different c. Result from R IM H sampler. Solid: prior; dotted: cdef=Q.3; dashdotted: cdef=0.S\ dashed: cdef=Q.l. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 72 5. DISCUSSION AND CONCLUSION With the help o f modem electronic computing power, Bayesian approaches can be successfully implemented and applied to many complex statistical problems. In this study, two methods incorporating Markov chain Monte Carlo methods, independence Metropolis-Hastings ( IMH) sampler and regenerative independence Metropolis- Hastings ( RIMH ) sampler, have been tested and applied in a simulated clinical situation. The information available through them is clinically meaningful and both of them produced consistent results. However, few techniques are available for quantitative characterization and analysis of their performance. Most of the methods that have been devised and used so far are oriented towards qualitative analysis based on graphical observation. The methods adopted in this work had been borrowed from the methods generally aimed for use with Gibbs sampler. But even with the few methods that had been tried here, both of the algorithms turned out to be consistent and useful methods in obtaining inferences on posterior distribution. MCMC methods offer a collection of samples from a distribution. Whether the samples can be taken to be a possible realization, especially in clinical environment, is not certain, even with the use o f resampling techniques. Sample based approaches enabled the investigator to overcome the difficulties o f parametric methods, but at the Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 73 same time, they leave a difficult task of how to analyze the sampling process and the final results. The IMH sampler is more straightforward to implement and evaluate than is the RIMH sampler. But the determination ofbum-in period and length of the chain depends on trial and error procedure. The RIMH algorithm produces iid sequences, which makes it easier to analyze the results in a sense. Yet, the criterion on how to choose the constant c, and what kind o f benefit can be expected in a certain situation need more investigation. MCMC methods have a great potential in that they can be combined with other sampling schemes to form a new algorithm. With a clever insight and understanding of the problem, they can provide a way for solving difficult problems in pharmacokinetic modeling and in other fields, too. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 74 REFERENCES Carlin, B. P., and Louis, T. (1996), Bayes and Empirical Bayes Methods fo r Data Analysis, London: Chapman & Hall. D’Argenio, D. Z., and Schumitzky, A. (1997), A D A P TII User’ s Guide: Pharmacokinetic/Pharmacodynamic Systems Analysis Software, Los Angeles: Biomedical Simulations Resource, University o f Southern California. Geman, S., and Geman, D. (1984), “Stochastic Relaxation, Gibbs Distributions and the Bayesian Restoration of Images,” IEEE Transactions on Pattern Analysis and Machine Intelligence, 6, 721-741. Geyer, C. J. (1992), “Practical Markov Chain Monte Carlo,” Statistical Science, 7,473- 511. Gilks, W. R., Richardson, S., and Spiegelhalter, D. J. (1995), “Introducing Markov Chain Monte Carlo,” in Markov Chain Monte Carlo in Practice, eds. W. R. Gilks, S. Richardson, and D. J. Spiegelhalter, London: Chapman & Hall, 1-19. Gladtke, E. (1988), “History o f Pharmacokinetics,” in Pharmacokinetics, eds. A. Pecile, and A. Rescigno, New York: Plenum Press, 1-9. Grimmett, G., and Stirzaker, D. (1982), Probability and Random Processes, Oxford: Clarendon Press; New York: Oxford University Press. Hastings, W. K. (1970), “Monte Carlo Sampling Methods using Markov Chains and Their Applications,” Biometrika, 57,97-109. Medhi, J. (1994), Stochastic Processes, 2nd ed., New York: John Wiley & Sons Inc. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E. (1953), “Equation of State Calculation by Fast Computing Machines,” The Journal o f Chemical Physics, 21, 1087-1092. Mykland, P., Tierney, L., and Yu, B. (1995), “Regeneration in Markov Chain Samplers,” Journal o f the American Statistical Association, 90, 233-241. Peck, C. C. (1993), “Rationale for the Effective Use of Pharmacokinetics and Pharmacodynamics in Early Drug Development,” in Integration o f Pharmacokinetics, Pharmacodynamics, and Toxicokinetics in Rational Drug Development, eds. A. Yacobi, J. P. Skelly, V. P. Shah, and L. Z. Benet, New York: Plenum Press, 1-5. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 75 Raftery, A. E., and Lewis, S. M. (1992), “How many intensions in the Gibbs Sampler?” in Bayesian Statistics 4, eds. J. M. Bernardo, J. O. Berger, A. P. Dawid, and A. F. M. Smith, Oxford: Oxford University Press, 763-773. Riegelman, S., Benet, L. Z., and Rowland, M. (1973), “Pharmacokinetics and Biopharmaceutics: A Definition o f Terms,” Journal o f Pharmacokinetics and Biopharmaceutics, I, 3-4. Ripley, B. D. (1987), Stochastic Simulation, New York: John Wiley & Sons. Seneta, E. (1981), Non-negative M atrices and M arkov Chains, 2nd ed., New York: Springer-Verlag. Shedler, G. S. (1993), Regenerative Stochastic Simulation, San Diego: Academic Press. Smith, A. F. M. (1991), “Bayesian Computational M ethods,” Philosophical Transactions o f the Royal Society o f London. A, 337,369-386. Tierney, L. (1996), “Introduction to General State-Space Markov Chain,” in Markov Chain Monte Carlo in Practice, eds. W. R. Gilks, S. Richardson, and D. J. Spiegelhalter, London: Chapman & Hall, 59-74. Welling, P. G. (1986), Pharmacokinetics: Processes and Mathematics, Washington, D.C.: American Chemical Society. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. IMAGE EVALUATION TEST TARGET (Q A -3 ) 150mm IIVW GE. In c 1653 East Main Street Rochester, NY 14609 USA Phone: 716/482-0300 F ax 716/288-5989 0 1993. Applied Image. Inc.. All Rights Reserved Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Bayesian inference using Markov chain Monte Carlo methods in pharmacokinetic /pharmacodynamic systems analysis

PDF

A user interface for the ADAPT II pharmacokinetic/pharmacodynamic systems analysis software under Windows 2000

PDF

Comparisons of deconvolution algorithms in pharmacokinetic analysis

PDF

Dynamics of the newly formed neuromuscular synapse

PDF

Head injury biomechanics: Quantification of head injury measures in rear-end motor vehicle collisions

PDF

A preliminary investigation to determine the effects of a crosslinking reagent on the fatigue resistance of the posterior annulus of the intervertebral disc

PDF

Comparing signal processing methods for spectral bio-imaging

PDF

Effects of prenatal cocaine exposure in quantitative sleep measures in infants

PDF

Biological materials investigation by atomic force microscope (AFM)

PDF

Propofol Effects On Eeg And Levels Of Sedation

PDF

Characteristic acoustics of transmyocardial laser revascularization

PDF

Cellular kinetic models of the antiviral agent (R)-9-(2-phosphonylmethoxypropyl)adenine (PMPA)

PDF

Construction and validation of multivariable population pharmacokinetic models: Utility of allometric forms to predict the pharmacokinetic disposition of gentamicin in pediatric patients with app...

PDF

Assessment of minimal model applicability to longitudinal studies

PDF

Destructive and non-destructive approaches for quantifying the effects of a collagen cross-linking reagent on the fatigue resistance of human intervertebral disc

PDF

A physiologic model of granulopoiesis

PDF

A model of upper airway dynamics in obstructive sleep apnea syndrome

PDF

Computation and validation of circulating blood volume with the indocyanine green dilution technique

PDF

A fluorescence microscopy study of quantum dots as fluorescent probes for brain tumor diagnosis

PDF

English phoneme and word recognition by nonnative English speakers as a function of spectral resolution and English experience

Asset Metadata

Creator Kang, Dongwoo (author)

Core Title Bayesian estimation using Markov chain Monte Carlo methods in pharmacokinetic system analysis

School Graduate School

Degree Master of Science

Degree Program Biomedical Engineering

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag engineering, biomedical,Health Sciences, Pharmacology,OAI-PMH Harvest,statistics

Language English

Contributor Digitized by ProQuest (provenance)

Advisor D'Argenio, David (committee chair), [illegible] (committee member), Kalaba, Robert E. (committee member)

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c16-25002

Unique identifier UC11336911

Identifier 1393173.pdf (filename),usctheses-c16-25002 (legacy record id)

Legacy Identifier 1393173.pdf

Dmrecord 25002

Document Type Thesis

Rights Kang, Dongwoo

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the au...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus, Los Angeles, California 90089, USA