Article pubs.acs.org/OPRD
Batch-to-Batch Variation: A Key Component for Modeling Chemical Manufacturing Processes Linas Mockus,*,† John J. Peterson,‡ Jose Miguel Lainez,§ and Gintaras V. Reklaitis† †
Purdue University, West Lafayette, Indiana 47907, United States GlaxoSmithKline Pharmaceuticals, Collegeville, Pennsylvania 19426, United States § University of Buffalo, Amherst, New York 14260, United States ‡
ABSTRACT: For chemical manufacturing processes, the chemical kinetics literature contains virtually no mention of quantitative models that involve batch-to-batch variation. Models for chemical process manufacturing quality improvement are being more carefully considered, particularly by the pharmaceutical industry and its regulators. This is evidenced in part by the recent ICH Q11 regulatory guidance on drug substance manufacture and quality. Quality improvement has been defined as a reduction in variation about a target. Hence the modeling of process variation plays an important role in quantifying quality improvement. Given that batch-to-batch variation is often a dominant source of process variation (usually exceeding measurement error variation), it is important for process models to incorporate such variance components. In this paper, we show how chemical kinetic models can incorporate batch-to-batch variation as well as measurement error. In addition, we show that these models can be used to quantify the reliability of meeting process specifications using Bayesian statistical methods. We also compare some different Bayesian computational approaches and recommend some software packages to aid with the computations.
■
INTRODUCTION For chemical processes, the laws of chemistry often allow engineers to construct mechanistic predictive models. While these models are typically meant to describe how a reaction profile changes over time, they tend to model the chemistry but not the entire process. The resulting chemical equation (or equations) are the backbone of the process model, but more is needed to describe the batch-to-batch variation that is observed each time a process is run. For example, slight perturbations in input material (and/or the ratio of input materials) that charge the reactor can induce some of the batch-to-batch variation that is observed. In some cases, uncontrolled variations in environmental conditions or manufacturing settings can also contribute to variations that are seen from one run to another. Fortunately, with the proper experimental design, data can be collected that will allow one to build an enhanced chemical kinetic prediction model that takes into account the effects of batch-to-batch variation as well as measurement error. For many chemical industries, particularly the pharmaceutical industry, quality improvement is an issue of growing importance among both producers and regulators. This is evidenced in part by the recent ICH Q11 regulatory guidance1 on drug substance manufacture and quality. While the notion of quality in ICH Q11 is rather broadly defined, it is important to note that quality improvement has been defined by some noted authors as a reduction in variation about a target.2 Hence the modeling of process variation plays an important role in quantifying quality improvement. Given that batch-to-batch variation is often a dominant source of process variation (usually exceeding measurement error variation), it is important for process models to incorporate such variance components. In fact, ICH Q11 lists “identification of potential sources of process variability” as a key aspect of quality-by-design (QbD). © XXXX American Chemical Society
Batch-to-batch variation is certainly one of those potential sources. In this paper, we show how chemical kinetic models can incorporate batch-to-batch variation as well as measurement error. In addition, we show that these models can be used to quantify the reliability of meeting process specifications using Bayesian statistical methods. We also compare some different Bayesian computational approaches and recommend some software packages to aid with the computations.
■
ROLE OF BATCH-TO-BATCH VARIATION IN QbD
Quality improvement can be thought of as the reduction in variation about a target.2 In addition, the larger the process variation, the more likely it is for a process to fail to meet manufacturing specifications. In fact, reduced process variation can even lower cycle time in manufacturing.3 As such, the ability to quantify and model variation is a key element in QbD. Unfortunately, many engineers tend to build models for their processes that do not properly include terms for batch-to-batch variation. In addition, experimental designs sometimes do not employ true replicates of the process under the same experimental conditions (e.g., same temperature, pressure, etc.). In such cases, information on batch-to-batch variation may be confounded with other experimental factors and not able to be modeled or directly observed. A recent review of experimental designs for chemical kinetic models does not appear to address the issue of batch-to-batch information.4 Special Issue: Application of ICH Q11 Principles to Process Development Received: July 28, 2014
A
dx.doi.org/10.1021/op500244f | Org. Process Res. Dev. XXXX, XXX, XXX−XXX
Organic Process Research & Development
Article
Ω is equal to zero.) Both the bj and eij random effects play an important role in modeling predictive values for future responses. More detailed and complex versions are possible,8 but the above model is quite flexible for many applications. The ability to fit the model form in eq 1 using a general procedure such as maximum likelihood estimation may require experiments involving replicated runs under the same experimental conditions so that information on the effects of batch-to-batch variation on the model can be obtained. In some cases, investigators will fit separate models over time (for different runs but under identical conditions) and compare the regression model parameters from each run to see how they are affected across batches. Not all model parameters will be equally affected by batch-to-batch variation, and in some cases investigators will choose to keep certain parameters fixed if there is little or no statistical evidence that batch-to-batch variation affects that model parameter. Model selection criteria such as Akaike’s Information Criterion (AIC) or Schwarz’ Bayesian Information Criterion (BIC)9 can be used to choose among models to avoid overfitting or underfitting the data. Once a model form, such as in eq 1, has been determined, it can be used to construct a predictive distribution to quantify the probability of a future response meeting specification limits. In some cases, the model may only exist in closed form as a differential equation. In such cases, numerical integration techniques will be needed to produce a prediction from the model. Using Bayesian statistical methods, a predictive distribution can be obtained from a posterior distribution of the model parameters.10−12 In accordance with Bayes’ theorem (eq 2), the posterior belief about the model parameters requires consideration of new data obtained from experiments (D) and the prior belief represented by the prior distribution p(ϕ) (e.g., a distribution of chemical kinetic parameters for the population). The probability that the experimental outcomes (D) are explained by the model, and a specific set of values for the unknown parameters (ϕ) is captured by the likelihood function (L(D| ϕ)).
Including true replicate experiments (even if at only the target set point) in experimental designs can help to locate and quantify the effects of batch-to-batch variation on the process. Proper model building to include terms that address batchto-batch variation can provide insight into how such variation affects the predictive model and hence the process itself. In the next section, we discuss how models can be built so that the random effects of batch-to-batch variation can be determined. Variations from run to run in a chemical manufacturing process can affect the chemical reaction in different ways over time. Even at the same manufacturing set point, batch-to-batch variation can induce slight perturbations in the time profile. However, such effects can be quantified using nonlinear mixedeffects models.5,6 Here, the phrase “mixed-effects models” refers to models that have both random and nonrandom parameters (in addition to the typical measurement random error term). It is the random parameters that induce the batch-to-batch variations that are observed in the time profiles. In addition, the issue of ICH Q8 design space7 has become an important topic for QbD, and is also discussed in ICH Q11 in the context of drug substance. ICH Q8 design space is the multidimensional combination and interaction of input variables and process parameters that have been “demonstrated to provide assurance of quality”. One way to demonstrate assurance of quality is to build (and validate) a model with an associated predictive distribution to compute (for each manufacturing condition) the probability of meeting quality specifications. Such a model may be greatly improved by including terms related to batch-to-batch variation. An ICH Q8 design space can be constructed using a predictive distribution generated from a nonlinear mixed-effects (mechanistic) model that includes terms related to batch-to-batch variation.
■
MODELING BATCH-TO-BATCH VARIATION Batch-to-batch variation will often appear to have more than just an additive effect. It may also affect each kinetic model parameter individually. As such, models of the form yij = f (ti , x j ; ϕ) + eij
p(ϕ|D) =
with the single random error term, eij, are often inadequate to model a chemical kinetic process from a stochastic perspective. Here, yij is the response at the ith time point for the jth experimental condition (batch), xj (e.g., temperature, pressure configuration), ti is the ith time point, ϕ is a vector of (unknown) model parameters, and eij is residual random error, typically modeled as having a normal distribution with mean zero and variance, σ2. To accommodate batch-to-batch variation, a more general class of models is needed. In most cases, the standard nonlinear mixed-effects model is able to model the various random influences incurred by variations that occur from one run to the next. The nonlinear, mixed effects model often has the form
yij = f (ti , x j ; ϕj ) + eij
L(D|ϕ)p(ϕ) p(D)
(2)
As kinetic data D becomes available, this information is processed to form the posterior distribution p(ϕ|D) for the unknown parameters using Bayes’ theorem. 1.1. Markov Chain Monte Carlo. Markov chain Monte Carlo (MCMC) methods have played an important role in making many complex Bayesian analyses computable.13 Markov chains are random processes where the current state value depends only upon the previous state value. Under proper conditions, a Markov chain process will converge to equilibrium, and the limiting values will have a particular distribution. Under a broad set of conditions, Markov chains can be used to sample from a posterior distribution of model parameters. Large samples from a posterior distribution can then be used to compute Monte Carlo estimates of probabilities of interest, such as the probability of future responses meeting specifications or the probability that a certain function of the model parameters will have a value in a desired range. In order for a Markov chain to reach the limiting distribution, it must be repeatedly run through many of its process states. This is often referred to as “burn-in”. Once burn-in has been
(1)
where ϕj = ϕ + bj. Here, bj is typically modeled as a random vector having a multivariate normal distribution with mean vector, 0, and variance−covariance matrix, Ω. The random effect bj is interpreted as the random batch effect for the jth batch or run. Note that the above model allows for batch-tobatch variation to affect some or all of the regression model parameters. (If the kth regression parameter is not affected by batch-to-batch variation, then the kth variance associated with B
dx.doi.org/10.1021/op500244f | Org. Process Res. Dev. XXXX, XXX, XXX−XXX
Organic Process Research & Development
Article
Figure 1. Schematic of the decomposition proposed for the variational Bayes’ approach.
achieved, the Markov chain has reached its equilibrium distribution. For Bayesian MCMC applications, one can sample from a posterior after burn-in has been achieved. However, samples from the posterior (which come from the Markov chain) are not necessarily independently distributed. In many cases, the samples are highly autocorrelated, meaning values “close” to each other in the Markov chain sequence are highly correlated, with the correlation dying off (slowly or quickly) as values grow further apart. Typically, software for Bayesian MCMC allows users to “thin” the Markov chain by discarding most values and keeping only every 10th, 50th, 100th value, and so forth. Such thinning can help to reduce the autocorrelation among the remaining values from the Markov chain. WinBUGS14,15 and Openbugs15 are popular software packages for implementing a type of MCMC called Gibbs sampling.16,17 However, hierarchical models with random effects, particularly nonlinear ones, are often associated with MCMC procedures whose Markov chains suffer from severe autocorrelation, even after much thinning. The unfortunate consequence of this severe autocorrelation is that there is noticeably less Monte Carlo information in such a sample from the posterior than there would be in a sample of the same size but where the elements are stochastically independent samples from the posterior. Some MCMC practitioners recommend no thinning, but rather the use of very large samples from the posterior if the sampled parameter values are highly autocorrelated. Unfortunately, this can cause practical problems associated with memory storage and computation time as millions of autocorrelated samples may be needed to obtain the information equivalent to a few thousand independent samples from the posterior. This is not to say that MCMC algorithms cannot be used to sample from the posterior for nonlinear, hierarchical models. In some cases, long burn-ins and much thinning (e.g., taking every 100th or 250th sample) will provide an ultimate subsample from the posterior that contain only mildly autocorrelated values that have an adequate amount of information for making Monte Carlo estimates. In fact, given the ease of use of WinBUGS or Openbugs, many investigators will start with an attempt to use MCMC by way of one of these software packages. However, problems with MCMC algorithm convergence or too much autocorrelation that cannot be reduced may necessitate the use of a non-MCMC approach to Bayesian computation. In the next two subsections, we discuss two such non-MCMC approaches, variational Bayes (VB)18 and generalized direct sampling.19 1.2. Variational Bayes. VB methods translate the Bayesian inference problem into an optimization problem by approximating the posterior distribution using a known distribution form (e.g., a family of Gaussian distributions). The methods are based on the separability of the logarithm of p(D) into two elements as shown in eqs 3−5.
ln p(D) = L(q) + KL(q||p)
L(q) = −
⎡
q(ϕ|θ)
(3)
⎤
∫ϕ q(ϕ|θ)ln⎢⎣ L(D|ϕ)p(ϕ) ⎥⎦dϕ
KL(q||p) =
(4)
⎡ q(ϕ|θ) ⎤
∫ϕ q(ϕ|θ) ln⎢⎣ p(ϕ|D) ⎥⎦dϕ
(5)
Here, q(ϕ|θ) is the parametric probability distribution that approximates the posterior distribution p(ϕ|D), and q is the set of parameters characterizing q (e.g., mean and covariance for a Gaussian distribution). If q(ϕ|θ) is free to be any probability function, then the maximum lower bound (L) is obtained when the Kullback−Leibler (KL) divergence is zero. This occurs when q(ϕ|θ) exactly matches the posterior distribution.20 Then, minimizing the KL divergence is equivalent to maximizing the lower bound (L). Therefore, the set of parameters q is determined as that which maximizes L. Laı ́nez et al.21 develop a decomposition strategy comprising three steps (Figure 1) to deal with the variational inference of models. 1.3. Generalized Direct Sampling (GDS). Generalized direct sampling provides a way to randomly sample from the posterior distribution but does not utilize MCMC. Instead, it makes clever use of rejection and slice sampling.22 It was first proposed by Braun and Damien.19 Braun has developed an R package, bayesGDS,23 for implementation of his approach. While computer run times for GDS may or may not be faster than MCMC (depending upon the software used, the manner of computational execution, and the nature of the problem being analyzed), GDS does provide stochastically independent samples from the posterior. This means that a sample of N values from a GDS procedure will have at least as much information than an MCMC sample of size N. In addition, it is easier to execute GDS in a parallel programming environment than for MCMC.19 Since GDS is a stochastic sampling algorithm (like MCMC), it is straightforward to compute Monte Carlo estimates of probabilities concerning functions of model parameters or posterior predictions. One particular advantage of GDS for chemical kinetic models with batch effects is that GDS tends to be computationally more efficient with nonlinear, mixed effect models than do many MCMC algorithms. It can be difficult to get convergence with MCMC algorithms for nonlinear regression.10 Another advantage for GDS is that it is straightforward to compute the marginal likelihood of the data (conditional only on the model form, not the model parameters). This allows one to compute posterior probabilities for several competing model forms. Such probabilities can be computed with MCMC but with much greater difficulty. C
dx.doi.org/10.1021/op500244f | Org. Process Res. Dev. XXXX, XXX, XXX−XXX
Organic Process Research & Development
■
Article
CHEMICAL KINETIC CASE STUDY This example comes from an experiment discussed in Crump et al.,24 who describe the influence of key reaction parameters on the development of a rate model which can be used to forecast starting material conversion independent of scale. Since the experimental design included three separate runs (center points) under identical conditions of temperature, pressure, and so forth, some information existed to allow one to model run-to-run variation in addition to measurement error over time. Figure 2 shows these three (center point) experimental
possible by the end of the reaction run. The vertical axis represents the proportion of remaining substrate. The horizontal axis represents the run time (in minutes). As one can easily see in Figure 2, the run-to-run variation is much greater than the within-run variation. The entire experiment included seven runs, the three shown in Figure 2 plus four additional runs at differing levels of experimental conditions. One of the goals of this experiment was to be able to predict the probability of the substrate being less than 1%, as a function of run time. In accordance with ICH Q11, this substrate was designated a critical quality attribute (CQA) as it is inversely proportional to the purity of the final drug substance. A predictive model that includes batch-to-batch variation terms as well as measurement error was used. This model is a nonlinear, mixed-effects kinetic model. The underlying kinetic model used has the following form, ln[Sub]i = μ0 − (k exp(− EaR /T )[Cat]pH /KH)(ti − t0) + ei 2
where [Sub]i is the measured substrate concentration as a function of time, ti (>t0) and ei is the (within run) measurement random error at time ti. There are four experimental factors: T (temperature, °C), [Cat] (catalyst concentration), pH2 (hydrogen partial pressure), and KH (Henry’s law constant). The unknown model parameters are μ0, the mean substrate at time t0, a kinetic rate constant, k, and EaR, where EaR = Ea/R and Ea is the activation energy and R is the gas constant. This model was then modified to include random batch effects perturbing the μ0and k model parameters. (The EaR parameter did not have any statistically significant variations from run to run, and hence no random effect was associated with that model parameter.) The chemical kinetic model,
Figure 2. Overlay plot of three experimental runs under identical conditions.
runs. In this reaction, the substrate in question is consumed over time. It is desired to have as little of this substrate as
Figure 3. Overlay plots of model parameter posterior distributions. D
dx.doi.org/10.1021/op500244f | Org. Process Res. Dev. XXXX, XXX, XXX−XXX
Organic Process Research & Development
Article
Figure 4. Overlay plots of r(t), reliability functions for three different experimental conditions.
r(t), or perform some kind of Bayesian predictive analysis needs access to a collection of different Bayesian computing tools, as the MCMC approaches may run into computational difficulties. In addition, for important QbD inference problems, it is prudent to execute more than one type of Bayesian algorithm so that we can compare the differences. This is because Bayesian calculations involving complex models (such as in eq 6) are virtually all approximate. In fact, for MCMC algorithms, the currently available convergence diagnostics11 provide no guarantee of convergence in general.19 As such, this paper considers three computationally different Bayesian modeling approaches: MCMC (Gibbs sampling), GDS, and VB. Using the example of this section, we illustrate the differences between these three Bayesian computational approaches. Using the above three different Bayesian computational approaches, we show overlay plots of the posterior distributions of the model parameters and then show an overlay plot of the computed Bayesian reliability function, r(t). As we can see from the posterior parameter density plots in Figure 3, the VB posterior provides an approximation which is tighter than the Monte Carlo algorithms: GDS and Gibbs sampling (via WinBUGS). The color coding is as follows: GDS is in red, Gibbs sampling is in blue, and VB is in green. However, the three corresponding Bayesian reliability function plots in Figure 4 (computed from the respective posterior predictive distributions) come surprisingly close to each other. However, this is just one example. We believe that the prudent modeler should consider more than one type of Bayesian algorithm for important inferences when dealing with complex models, particularly nonlinear mixed-effects models. The computations for these three examples were done on a laptop and took about 15 min each to complete.
modified to include the random batch effects, has the enriched hierarchical form, ln[Sub]ij = (μ0j + b1j) − {exp(ln k + b2j) exp(−EaR /Tj) × [Cat j]pH j /KHj}(tij − t0j) + eij 2
(6)
whereb1j and b2j are the random effects of the jth batch on the parameters μ0j and ln k, respectively. For the model in eq 6, it is assumed that the mean substrate value at time t0j is μ0j. Here, b1j and b2j are assumed to have independent normal distributions with means equal to zero and variances σ2b1 and σ2b2, respectively. Also, eij is the random measurement error for the ith time point within the jth batch. We assume that eij all have independent normal distributions with mean zero and variance equal to σ2. The j subscripts for the four experimental factors, T, [Cat], pH2, and KH, refer to the different values they take on for the jth experimental run. Crump et al.24 developed this model using the R package, nlme, which can compute maximum likelihood estimates9,25 and the associated AIC and BIC model selection statistics, which can be used to identify which kinetic model parameters should be enhanced to have random batch effects. They then used the R BRugs (Bayesian) package26 to compute the Bayesian reliability function, r(t) = Pr(Sub < 0.1|t, data), to assess how long of a run time, t, would be needed to have a high probability of the substrate being converted to less than 1%. However, the statistical details were omitted. Such modeling concepts and details are needed to use the R packages nlme and BRugs (or other statistical software) in performing such computations. In addition, we believe that the modeler who wishes to compute process reliabilities, such as E
dx.doi.org/10.1021/op500244f | Org. Process Res. Dev. XXXX, XXX, XXX−XXX
Organic Process Research & Development
■
Article
DISCUSSION Since quality improvement is intimately connected to reduction in variation about a target, process modeling for improving quality needs to take into account all key sources of process variation. As one can see from Figure 2, between-run variation can be substantially larger than within-run variation for chemical kinetic processes. Therefore, the modeling of chemical kinetic processes should consider models (such as in eq 6) which have random effects associated with batch-to-batch variation. This is particularly important for making inferences about process reliability based upon the posterior predictive distribution of future process responses. The distribution of such responses will, of course, depend strongly upon betweenrun as well as with-in run variation. Ignoring batch-to-batch variation can cause one to underestimate process variation, which in turn can bias QbD process reliability inferences. Fortunately, the Bayesian approach provides (in theory) a straightforward way to construct predictive distributions for simple or complex chemical kinetic models. However, the computational aspects are not fool proof. It is important to keep in mind that there are methods to assess MCMC convergence. Since these methods cannot guarantee convergence, it is prudent for important inference problems to try more than one algorithm to check for any differences among them. For modelers who do not have the time or expertise to program existing algorithms, it is good to know that there is a growing number of R packages and applications dedicated to Bayesian inference. In addition, some of these can provide Bayesian inference for chemical kinetic models with random batch effects. WinBUGS27 is a particularly flexible and powerful software application that can be run as a stand-alone application or called from the R using the R package R2WinBugs.28 Not only can WinBUGS address nonlinear mixed-effects models, but the WinBUGS add-on, WBDiff,29 can also model ordinary differential equations. The Variational Bayes algorithm is not readily available since it uses GAMS CONOPT or IPOPT solvers to solve the large Kullback−Leibler (KL) divergence optimization problem. There are plans of using a free version of ipopt30 and for it to be implemented as an R package. It is developed specifically to solve pharmacokinetics problems which may be formulated as a system of differential-algebraic equations (DAEs). We found it to work very well with smaller (up to 20 parameters) yet complicated multiple compartment pharmacokinetic models. This limit is imposed by sparse grids31 which are used for integration. In this paper we investigated the application of variational Bayes (VB) to hierarchical chemical kinetics models. A key advantage of the VB algorithm (over MCMC algorithms) is that it is based on optimization and is not affected by autocorrelation which is pronounced in hierarchical models. The GDS algorithm is available through the R package bayesGDS. While this package requires the user to construct an R function for the log posterior (up to an additive constant), in our experience it works well for nonlinear mixed-effects models. A key advantage of the GDS algorithm (over MCMC algorithms) is that it produces stochastically independent samples from the posterior. For some problems this can reduce overall computing time as compared to WinBUGS, with regard to obtaining a specified effective sample size.13
While we have discussed three different Bayesian algorithms and associated software for their implementation, still other algorithms and software are available for application to nonlinear mixed-effects models. The SAS/STAT package, v9.3,32 has a Bayesian procedure called PROC MCMC, which can sample from the posterior for nonlinear mixedeffects models. PROC MCMC uses an MCMC procedure called the random walk Metropolis algorithm. The STAN application33,34 can be called from R using the RStan package. STAN uses an MCMC variant of Hamiltonian Monte Carlo for computations that are generally more efficient than Gibbs sampling as found in WinBUGS. Bayesian researchers are continually producing new statistical algorithms and associated packaged software for their implementation, and there are further algorithmic approaches recently published, for which we do not have the space to discuss. So chemical process modelers can expect increasing access to tools for improved predictive modeling using key sources of variation for improved QbD and reliability-based inferences. We believe that Bayesian analysis is a powerful statistical tool for producing predictive distributions for modeling chemical kinetic processes. This is because it can be adapted to models that take into account batch-to-batch variation. Since it is still somewhat of a delicate computational tool, prudent analysis may require care and possibly utilization of two or more different algorithms. However, such algorithms are becoming more easily accessible by way of packaged software. We feel that Bayesian predictive methods should be an important part of a process modeler’s repertoire and will have effective application for chemical kinetic processes, particularly for QbD and/or risk assessment inference. While the Bayesian statistical and computational methodologies may seem like difficult hurdles to the uninitiated, we feel that the acquisition and use of such methods are worth it if one really wants to quantify risk associated with a process. Collaboration with engineers or statisticians skilled in Bayesian methods is an excellent way to climb the learning curve. Seeking advice on experimental design, with a view to including replications, is also important for acquiring data to support predictive models that quantify batch-to-batch variation and can be used for quantitative risk assessment.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS We thank Brian Crump and Rick Lewis for consultation and advice related to the chemical kinetics example in section on Chemical kinetic case study. In addition, we thank them for their innovative work in designing and executing a study that uses replication to improve the state-of-the-art for modeling of chemical kinetic processes.
■
REFERENCES
(1) ICH. ICH Harmonized Tripartite Guideline: Development and Manufacture of Drug Substances; International Conference on Harmonisation of Technical Requirements for Registration of Pharmaceuticals for Human Use Q11, Geneva, 2012. (2) Montgomery, D. C. Introduction to Statistical Quality Control; John Wiley & Sons, Inc.: New York, 2009. F
dx.doi.org/10.1021/op500244f | Org. Process Res. Dev. XXXX, XXX, XXX−XXX
Organic Process Research & Development
Article
(3) Nunnally, B. K.; McConnell, J. S. Six Sigma in the Pharmaceutical Industry − Understanding, Reducing, and Controlling Variation in Pharmaceuticals and Biologics; CRC Press: Boca Raton, FL, 2007. (4) Kitsos, C. P.; Kolovos, K. G. A Compilation of the D-Optimal Designs in Chemical Kinetics. Chem. Eng. Commun. 2013, 200 (2), 185−204. (5) Davidian, M.; Giltinan, D. M. Nonlinear Models for Repeated Measurement Data; Chapman and Hall: New York, 1995. (6) Vonesh, E. F.; Chinchilli, V. M. Linear and Nonlinear Models for the Analysis of Repeated Measurements; Marcel Dekker: New York, 1997. (7) ICH. ICH Harmonized Tripartite Guideline: Pharmaceutical Development; International Conference on Harmonisation of Technical Requirements for Registration of Pharmaceuticals for Human Use Q8, Geneva, 2005. (8) Davidian, M.; Giltinan, D. M. Nonlinear Models for Repeated Measurement Data: An Overview and Update. J. Agric., Biol., Environ. Stat. 2003, 8, 387−419. (9) Pinheiro, J.; Bates, D. Mixed-effect Models in S and S-PLUS; Springer: New York, 2000. (10) Christensen, R.; Johnson, W.; Branscum, A.; Hanson, T. E. Bayesian Ideas and Data Analysis − An Introduction for Scientists and Statisticians; Chapman & Hall/CRC Press: Boca Raton, FL, 2011. (11) Gelman, A.. Carlin, J. B.; Stern, H. S.; Rubin, D. B. Bayesian Data Analysis, 3rd ed.; Chapman and Hall/CRC: Boca Raton, FL, 2004. (12) Peterson, J. J. A Bayesian Approach to the ICH Q8 Definition of Design Space. J. Biopharm. Stat. 2008, 18, 959−975. (13) Gamerman, D.; Lopes, H. F. Markov Chain Monte Carlo, 2nd ed.; Chapman & Hall/CRC: Boca Raton, FL, 2006. (14) Ntzoufras, I. Bayesian Modeling Using WinBUGS; John Wiley & Sons, Inc.: Hoboken, NJ, 2009. (15) Lunn, D.; Jackson, C.; Best, N.; Thomas, A.; Spiegelhalter, D. The BUGS Book − A Practical Introduction to Bayesian Analysis; Chapman & Hall/CRC: Boca Raton, FL, 2013. (16) Geman, S.; Geman, D. Stochastic relaxation. Gibbs distributions and the Bayesian restoration of images. IEEE Trans. Pattern Anal. Machine Intell. 1984, 6, 721−741. (17) Gelfand, A. E.; Smith, A. F. M. Sampling-based approaches to calculating marginal densities. J. Am. Stat. Assoc. 1990, 85, 398−409. (18) Fox, C. W.; Roberts, S. J. A tutorial on variational Bayesian inference. Artif. Intell. Rev. 2012, 38, 85−95. (19) Braun, M.; Damien, P. Generalized Direct Sampling for Hierarchical Bayesian models. arXiV 2012, 1108.2245v3; http:// arxiv.org/abs/1108.2245v3. (20) Bishop, C. M. Pattern recognition and machine learning; Springer Science: New York, 2006. (21) Laínez, J. M.; Mockus, L.; Orcun, S.; Blau, G.; Reklaitis, G. V. A decomposition strategy for the variational inference of complex systems. Technometrics 2012, accepted. (22) Robert, C. P.; Casella, G. Monte Carlo Statistical Methods; Springer-Verlag: New York, NY, 1999. (23) Braun, M. bayesGDS: Functions to implement Generalized Direct Sampling, R package version 0.5.0, 2012. http://CRAN.Rproject.org/package=bayesGDS. (24) Crump, B. R.; Goss, C.; Lovelac, T.; Lewis, R.; Peterson, J. Influence of Reaction Parameters on the First Principles Reaction Rate Modeling of a Platinum and Vanadium Catalyzed Nitro Reduction. Org. Process Res. Dev. 2013, 17, 1277−1286. (25) Pinheiro, J.; Bates, D.; DebRoy, S.; Sarkar, D. nlme: Linear and Nonlinear Mixed Effects Models, R package version 3.1−118, 2014; http://cran.r-project.org/web/packages/nlme/nlme.pdf. (26) Thomas, A.; O’Hara, B.; Ligges, U.; Sturtz, S. Making BUGS Open. R News 2006, 6 (1), 12−17. (27) Spiegelhalter, D.; Thomas, A.; Best, N.; Lunn, D. WinBUGS User Manual v1.4, 2003; http://www.mrc-bsu.cam.ac.uk/wp-content/ uploads/manual14.pdf. (28) Sturtz, S.; Ligges, U.; Gelman, A. R2WinBUGS: A Package for Running WinBUGS from R. J. Stat. Software 2005, 12 (3), 1−16.
(29) Lunn, D. WinBUGS Differential Interface (WBDiff), 2004. http://www.winbugs-development.org.uk/wbdiff.html. (30) Wächter, A.; Biegler, L. T. On the Implementation of a PrimalDual Interior Point Filter Line Search Algorithm for Large-Scale Nonlinear Programming. Math. Programming 2006, 106 (1), 25−57. (31) Heiss, F.; Winschel, V. Likelihood approximation by numerical integration on sparse grids. J. Econometrics 2008, 144 (1), 62−80. (32) SAS. SAS/STAT 13.1 User’s Guide; SAS Institute Inc.: Cary, NC, 2013. (33) Stan. Stan Modeling Language - User’s Guide and Reference Manual, Stan Version 2.1.0, 2013; http://mc-stan.org/manual.html. (34) Stan. Stan: A C++ Library for Probability and Sampling, Version 2.1, 2013; http://mc-stan.org.
G
dx.doi.org/10.1021/op500244f | Org. Process Res. Dev. XXXX, XXX, XXX−XXX
