Bayesian Applications to QualitybyDesign and Assay Development John Peterson PhD Director Statistical Sciences Group GlaxoSmithKline Pharmaceuticals Collegeville Pennsylvania USA NonClinical Statistics Conference ID: 532579
Download Presentation The PPT/PDF document "Short Course on" is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.
Slide1
Short Course onBayesian Applications to Quality-by-Design and Assay Development
John Peterson, Ph.D.Director, Statistical Sciences GroupGlaxoSmithKline Pharmaceuticals Collegeville, Pennsylvania, USA
Non-Clinical Statistics Conference, Brugge, Belgium, 8 October, 2014
1
john.peterson@gsk.comSlide2
2OUTLINE
Some Bayesian preliminaries
What is QbD? Some Basic Concepts The Flaw of Averages vs. Predictive Distributions
Process Optimization using Bayesian Predictive Distributions
A Bayesian Approach to ICH Q8 Design Space & Scale-up
Bayesian Monte Carlo Studies for USP Test Assessment Assay Development & Dissolution
Some Computational Recommendations
Acknowledgements
References Slide3
3
3
3
3
The Bayesian Paradigm
P
osterior
distribution of
unknown model parameters
prior distribution of
unknown model parameters
Bayes’ rule
weighting function
x
=Slide4
4
4
4
4
4
The Bayesian Paradigm
P
osterior
distribution of
unknown model parameters
prior distribution of
unknown model parameters
Bayes’ rule
weighting function
x
=Slide5
5
5
5
5
5
5
The Bayesian Paradigm
P
osterior
distribution of
unknown model parameters
prior distribution of
unknown model parameters
Bayes’ rule
weighting function
x
=
A
simulation procedure
known as “Markov Chain Monte Carlo” (
MCMC
)
can be used to
sample from
the posterior distribution of unknown model parameters
given experimental data
and
prior
information.
Other computational procedures can sometimes be used such as : numerical integration,
generalized direct sampling, Laplace approximation approaches (e.g. INLA). Slide6
6Three Practical Benefits of Bayesian Statistics
Provides a scientific way to incorporate informative prior information.
- Probability theory and risk assessment can be used to build informative priors. Priors can be used to provide a certain degree of “regularization”. - Priors can be used to induce necessary boundaries (e.g. positive variance components).
- Priors can be used to stabilize parameter estimates (e.g. ridge regression)
- Priors on model forms can be used to obtain “model averaging”, which improves
predictions. (This has been successfully applied in data mining competition.)
Bayesian methodology provides a direct method for getting a predictive distribution.
- For complex models (often needed for pharmaceutical process modeling) the
frequentist approach may prove difficult here (e.g. nonlinear mixed-effect model). Slide7
Prior distributions can be used to incorporate prior information
about a process (in a scientific way ….via Bayes Rule) - This can be done by using a distribution that maps out the
relative possibilities of parameter values. - Or, a prior can simply designate parameter boundaries. (e.g. a gamma distribution prior for a parameter that must be between positive.)
7
Prior Distributions
q
s
0
shape=0.1
rate=0.1Slide8
8
Prior distributions can be weak or non-informative e.g. using a flat prior for the mean in a normal-theory model.
Sometimes priors are used to induce a certain degree of “regularization”. - For example we may choose a multivariate normal prior for a vector of regression coefficients (intercept omitted) that has a mean vector
of zero. - This will tend to “pull” the posterior a bit towards the zero vector and
thereby prevent regression coefficient estimates from becoming too large in situations of multi-collinearity. (Ridge regression can be
thought of a Bayesian procedure with respect to a zero-centered prior for the regression coefficients.)
8
Prior DistributionsSlide9
9
9
In some cases, care must be used in choosing priors:
(
i) For small sample sizes or weak information about a model parameter.
(ii) For informative priors (iii) Even for non-informative priors in certain situations.
For details see: Seaman, J. W. III, et al. (2012). “Hidden Dangers of Specifying
Noninformative Priors”,
The American Statistician, 66(2), 77-84.
9
Prior Distributions
Seaman et al. recommend that if one is interested a particular function of the
model parameters, then one should use the chosen prior to observe the
induced prior for that function (e.g. by Monte Carlo simulation). Slide10
10Posterior Predictive Distributions
Predictive distributions is why I became a Bayesian statistician.
Process optimization and QbD is why I desire predictivedistributions.
Why? Slide11
11
11
Posterior Predictive Distributions
Predictive distributions is why I became a Bayesian statistician.
Process optimization and QbD is why I desire predictive
distributions.
Why?
Because predictive distributions are a very good way to
quantify process capability! Slide12
12
12
12
Posterior Predictive Distributions
Predictive distributions
is why I became a Bayesian statistician.
Process optimization and
QbD
is why I desire predictive
distributions. Why?
Because predictive distributions are a very good way to quantify process capability!
A process can be optimized, but that does not imply that it is “capable”, i.e. likely to meet specifications. Slide13
13
13
13
13
Posterior Predictive Distributions
Posterior distribution
for model parameters
0
standard
normal
distribution
Posterior predictive
distributionSlide14
14
14
14
14
14
Posterior Predictive Distributions
Posterior distribution
for model parameters
0
standard
normal
distribution
Posterior predictive
distribution
Sampling from the posterior predictive distribution, we can then compute
as a measure of process capability. (Here,
S
is the specification interval.) Slide15
15What is
QbD? Quality-by-Design (QbD
): Quality by Design (QbD) is a concept put forth by Joseph
Juran. (See
Juran, J. M. (1992) Juran on Quality by Design
)
Juran believed that quality could be designed into products and processes.
He called this idea “Quality by Design”
QbD has been applied in many industries, most notably the automotive industry.
Recently, the US Food and Drug Administration has adopted QbD principles for drugmanufacture.
The FDA QbD imperative is outlined in its report “Pharmaceutical Quality for the 21st Century: A Risk-Based Approach.”
Quote from Juran on Quality by Design: “Organizations that have adopted such methodsof quantification (process capability) have significantly outperformed those which havenot.”Slide16
16
16
Risk is an important concept for the ICH Q8-Q10 Guidances.
ICH Q8 - The word ‘risk’ or ‘risk-based’ occurs 35 times. ICH Q9 - The word ‘risk’ or ‘risk-based’ occurs 290 times. ICH Q10 - The word ‘risk’ or ‘risk-based’ occurs 34 times.
The growing importance of process reliability and risk assessment
for pharmaceutical manufacturing
The FDA continues to push its initiative on
“Pharmaceutical Quality for the 21st Century”
The ICH Quality (Q) Guidances are supported by the FDA. - They outline (at a high level) what the pharmaceutical companies should do to improve the quality of the manufacturing of their products.
Such improved quality strategies are referred to as “Quality by Design” or QbD. Currently, QbD is “optional” for Rx sponsors, but there appears to be increasing pressure to incorporate QbD as a strategy in Rx submissions. Slide17
17
Some Basic Quality Concepts
Process capability – is a measure of the inherent variability of a process measure about its target.Quality improvement – reduction in variation of a process measure about its target.
Quality improvement can therefore be thought of as an improvement
in process capability.
Quality concepts for a single quality measure: Slide18
18
18
Some
Basic
Quality Concepts
Process capability – is a measure of the inherent variability of a process measure about its target.
Quality improvement – reduction in variation of a process measure about its target.
Quality improvement can therefore be thought of as an improvement
in process capability.
Quality concepts for a single quality measure:
But what about products and processes involving multiple quality criteria?
What do we do with multiple variance (and covariance) measures? - Add them up? Find some function of a variance-covariance matrix?Slide19
19
19
Some Basic Quality
Concepts
Process capability
– is the inherent (multivariate) distribution of the process measures about a (vector) target.
Quality improvement – shrinking of the (multivariate) distribution of the
process measures about a (vector) target
Quality concepts for a multiple
quality measures: Quality improvement can therefore be thought of as an improvement
in process capability. Slide20
20
20
20
Some Basic Quality
Concepts
Process capability
– is the inherent (multivariate)
distribution
of the
process measures about a (vector) target.
Quality improvement – shrinking of the (multivariate) distribution of the process measures about a (vector) target
Quality concepts for a multiple quality measures:
Quality improvement can therefore be thought of as an improvementin process capability.
Generalizing from variances to
distributions
helps us to think more clearly
about dealing with single or multiple quality endpoints. Slide21
21
A natural process capability index for both univariate and multivariatequality measures is the proportion of the distribution of qualitymeasures that are within the quality specification limits.
This process capability index can be used as a desirability function for process optimization. Mathematically, we can express this as:
where Y
is a vector of quality measures, S is a specification region, and
x is a vector of process conditions.
p
(x
)
is a function of all of the means and variance components associatedwith the process. The expression for
p(x
) involves the posterior predictive distribution for
Y. Some Basic Quality
ConceptsSlide22
22
22
Some confusion about process capability
Process capability indices have been criticized as being incapable of capturing all aspects of a process’s capability of meeting specifications because they are a single number.
It has been argued that “process capability” should be assessed using various aspects of the distribution of a process’s responses, instead of using a single metric (
Pignatiello and Ramberg, 1993)
I suspect that this confusion arises from not recognizing the difference between “exploratory analysis” and “confirmatory analysis”.
I believe that one should explore the nature of a process’s response distribution, but in the end one must choose some optimization criterion.
If we are to optimize a process, we should optimize it to be as capable aspossible. Hence, the use of a process capability index makes sense for
process optimization (Plante (2001), Jeang (2010) ). But care is required.Slide23
23
23
23
Some Points about Process Capability Assessment
Process Capability Assessment is important.
- Just because your process mean is on target does not necessarily imply that
your process will meet specifications a high proportion of the time.
Process Capability Assessment is statistically more difficult than other
quality improvement methods, particularly for multiple
responses. (This is why statisticians are important here!)
It is important to have the process in control as much as possible forprocess capability assessment. We need to assess process capability in both exploratory and confirmatory
ways.
Exploratory:
- Explore the process means relative to specifications.
- Identify and measure process variance components and correlations.
- Assess response distribution at various operating conditions.
Confirmatory:
- Compute process capability index (after process is in control) .
- Validate process capability dynamically during SPC phase. Slide24
24
Common Cause vs. Special Cause Variation Some Basic Quality Concepts
Common Cause Variation: - Common cause variation is the process variation due to “common causes”. - Common causes are associated with the usual historically known, quantifiable variation
of a process. - Typically, it is variation that can be modelled.
Special Cause Variation:
- Special cause variation is the process variation due to “special causes” (sometimes referred to as “assignable causes”).
- Special causes are typically come in the form of surprises or are outside the historical experience base.
- Some examples of special causes may be: operator error, machine malfunction, a surprisingly poor batch of raw material, or a change in a process condition that
was not previously recognized as influential. Slide25
25This process capability index as a desirability function
We assume that
Y
can be modeled as a (stochastic) function of process
factors,
x.
For example: Note that
We can optimize a process (under the above model) by maximizing
p
(x) with respect to the factors in
x. p(
x) forms a desirability function for process optimization that is easy tointerpret and avoids flaws associated with mean response surface optimization.Slide26
26
Why is a
multivariate reliability approach needed? (Accounting for correlation among the responses…a simple example)
Suppose we have a process with four key responses,
Y
1, Y2
, Y3
, Y
4
For simplicity, let’s assume that Let
Consider If S = I , then But if
then and if Slide27
27
27
The Flaw of Averages:
Why We Underestimate Risk in the Face of Uncertainty by Dr. Sam Savage
The Flaw of Averages vs. Predictive Distributions
Average depth of river is 3 feet.Slide28
28
28
28
The Flaw of Averages vs. Predictive Distributions
Inferences about population
means
are rarely sufficient to
quantify good quality.
Recall:
Quality improvement – reduction in variation of a process measure about its target.
Or more generally…. Quality improvement – shrinking of the (multivariate) distribution
of the process measures about a (vector) target Slide29
29
29
The Flaw of Averages vs. Predictive Distributions
All classical textbooks on response surface methodology state a
response surface model for a process in the basic form:
where:
x
is a vector of process factors (e.g. temperature, pressure, etc.)
and
e
is a “statistical error” or “random error” term, where E(
e
)=0.
In these textbooks it is stated that is a “mechanistic model”
with a vector of unknown model parameters
b
).
In cases where the form of is unknown, it is recommended
that one approximate by a (second-order) Taylor series. Slide30
30
30
What is the truth?: mean models vs. distribution models
Given that the basic response surface model is represented in
classical textbooks asit is not too surprising that practitioners think (in the back of their
minds) that a response surface model is: observed process response = truth + measured error.
or possibly
observed process response = approx. truth + measured error.
Since I have seen some authors refer to the“true mean model” with regard to process optimization. Slide31
31
31
Consider the following thought experiment…
Suppose you want to optimize a complex manufacturing process.
Suppose further that you are able to acquire instruments to measure the
process response that are extremely accurate (producing truly negligible measurement error).
If you were to keep the process factors constant, but stop and started the
process several times using different batches of raw materials each time, would you expect to measure the exact same process response values?
NO! Slide32
32
32
Consider the following real experiment…
where
y
t
is the measured amount of substrate remaining after time
t
.
The parameters in
red
are unknown model parameters and the factors in
blue
are
process conditions of temperature, amount of catalyst, etc.
We want to be able to run the
experiment long enough to remove
most all of the substrate (less than 0.1% say)Slide33
33
Process modeling: Scientific modeling is the backbone… but is there more?
33
Chemistry equation:Slide34
34
Process modeling: Scientific modeling is the backbone… but is there more?
34
Chemistry equation:
Three replicate batches run under the
same
experimental conditions:
Same chemistry but different
profiles! Slide35
35
35
35
Traditional mechanistic model equation for prediction
Stochastic-mechanistic model equation that models batch-to-batch variation
Measurement
random error
Batch-to-batch
random effects
Predictive
distribution
for
risk assessment
Possible drivers:
Reactor is charged with
slightly different amounts
of raw materials.
Proportion of catalyst is
slightly different.
Changes in environmental
conditions.
Subtle operator differences
… Slide36
36
36
36
36
36
Traditional mechanistic model equation for prediction
Stochastic-mechanistic model equation that models batch-to-batch variation
Measurement
random error
Batch-to-batch
random effects
Predictive
distribution
for
risk assessment
Why is variation important
to quality improvement ?
Because quality improvement
can be defined as
“reduction in
variation about a target”Slide37
37
37
Predictive Distributions and Multiple Response Process Optimizatio
nSlide38
38
38
Predictive Distributions and Multiple Response Process OptimizationSlide39
39
39
Predictive Distributions and Multiple Response Process Optimization
Y = f
(
x
,
z
,
e
|
q
)
Variation due to
estimating unknown
model parameters.
Variation due to
common
-
cause
error variability
Multivariate predictive
distribution of quality
responses
Variation due to noisy
control variables, input
raw materials, etc.
This quality response
distribution is from
a much better
x
-
point
in the process.
Multivariate
predictive
model
A
very
reliable process
y2=friability
85%
60%
80%
100%
0%
99%
mean
Probability of meeting both
specifications is well above
0.99
!
Process
control
values
y1=Percent dissolved
(at 30 min.)
h
3%Slide40
Multiple response surface optimization:Beware the “sweet spot” for ICH Q8 Design Space
The classical textbooks in response surface methodology (RSM)refer to “overlapping mean” response surfaces as a way to optimize processes with multiple response types.
Popular point & click packages such as Design Expert and
SAS/JMP will produce such a “sweet spot” graph with a few clicks
of your mouse. Box & Draper
Myers, Montgomery, Anderson-Cook
Box, Hunter & HunterSlide41
41
41
ICH Q8 Regulatory Guidance: The “Design Space” Concept
ICH (
International
Congress on Harmonization) is an international body of representatives from the US, Europe, and Japan for the purposes of harmonization
of technical requirements for registration of pharmaceuticals for human use. The ICH has a set of quality guidelines, all with the prefix Q. ICH Q8 addresses
“pharmaceutical development” (i.e. development of the actual pill, liquid,
injectible drug product)
ICH Q8 contains a key concept called “design space”. (This is not a statistical term, but refers to the set of all manufacturing recipes that should produce acceptable product quality results.)
ICH Q8 defines the “design space” as: The multidimensional combination and interaction of input variables (e.g., materialattributes) and process parameters that have been demonstrated to provide assurance of quality.Slide42
42
42
42
ICH Q8 Regulatory Guidance: The “Design Space” Concept
ICH Q8 defines the “design space” as:
The multidimensional combination and interaction of input variables (e.g., material
attributes) and process parameters that have been demonstrated to provide
assurance of quality.
A “design space” example from the ICH Q8 regulatory guidance document:
&
Contour plot of “dissolution”
Contour plot of “friability”
Design space (white region)Slide43
43
43
43
43
ICH Q8 Regulatory Guidance: The “Design Space” Concept
ICH Q8 defines the “design space” as:
The multidimensional combination and interaction of input variables (e.g., material
attributes) and process parameters that have been demonstrated to provide
assurance of quality.
A “design space” example from the ICH Q8 regulatory guidance document:
Specification Region
Dissolution Spec. Limits
Friability spec. limits
Process
distribution
Process
distribution
Process
distributionSlide44
44
Overlapping Mean Design
Space Situation
with Three Response Types
y1 = tablet disintegration time,
y2 = friability, and y3 = hardness
y1
y2
y3
Three-dimensional distribution of responses at a point on the boundary of an
overlapping means design space.Red points are “out of spec”.
From Peterson, J. and Lief, K. (2010) “The ICH Q8 Definition of Design Space: A Comparison of the Overlapping Means and the Bayesian Predictive Approaches”, Statistics in Biopharmaceutical Research, 2, 249-259.
85% Out of Spec
44
For this process, the overlapping means
Design Space from the
Design Expert
package
harbored process configurations with a low
probability of meeting all three quality specifications.
The worst probability was about 15%.Slide45
45
Example
Number of factorsfor the DS
Number of responses for the DS
Minimum posterior
probability of acceptance
Maximum posterior
probability of acceptance
1
4
3
0.15
0.84
2
4
4
0.23
0.72
3
3
3
0.34
1
4
3
*
2
0.26
0.74
5
4
3
0.11
0.33
6
3
4
0.11
0.84
*
Mixture experiment. Factors sum to 1.
From Peterson, J. and Lief, K. (2010) “The ICH Q8 Definition of Design Space: A Comparison of the
Overlapping Means and the Bayesian Predictive Approaches”,
Statistics in Biopharmaceutical Research
,
2
, 249-259.
Examples from Six Different (Real Data) Experiments
Results Corresponding to Overlapping Mean Design SpacesSlide46
46A
Distribution Approach to Process OptimizationBuild your process model taking into account all key sources of variation
(e.g. batch-to-batch, measurement error) Design an experiment to gather data that will enable you to compute aposterior distribution for the unknown model parameters.
Using the posterior distribution, compute a posterior predictive distribution for your process.
Compute the probability of meeting specifications (or some other utility measure) as a function of process factors.
For example: Optimize
Optimize the process using the utility measure in 4. above.
(Here,
Y
is a vector of response types and
x
is a vector of process factors)Slide47
47
Example: An intermediate stage of a multi-stage route of manufacture for anActive Pharmaceutical Ingredient (API).Measurements: Four controllable quality factors (x’s) were used in a designed experiment.
(x1=‘catalyst’, x2= ‘temperature’, x3=‘pressure’, x4=‘run time’.) A (face centered) Central Composite Design (CCD) was employed. (It was a Full Factorial + axial points + 6 center points [30 runs])
Four quality-related response variables, Y ’s, were measured.
(These were three side products and purity measure for the final API.) Y1= ‘Starting material Isomer’, Y2=‘Product Isomer’, Y3=‘Impurity #1 Level’,
Y4=‘Overall Purity measure’
Quality Specification limits: Y1<=0.15%, Y2<=2%, Y3<=3.5%, Y4>=95%.Multidimensional Acceptance region,
Overlapping Means vs. Bayesian Reliability Approach to Design Space:
An Example
– due to Greg Stockdale, GSK.Slide48
48
Overlapping Means vs. Bayesian Reliability Approach to Design Space:An Example – due to Greg Stockdale, GSK.
Prediction Models:
Model Terms
Response
x1
x2
x3
x4
x11
x22
x33
x44
x12
x13
x14
x23
x24
x34
SM Isomer
D
D
D
D
D
D
D
D
Prod Isomer
D
D
D
D
D
Impurity
D
D
D
Purity
D
D
D
Temperature = x1 Pressure = x2 Catalyst Amount = x3 Reaction time = x4Slide49
49
Design Space Table of Computed Reliabilities1for the API (sorted by joint
probability)
Temp
Pressure
Catalyst
Rxntime
Joint
Prob
SM
Isomer
Prob
Prod
Isomer
Prob
Impurity
Prob
Purity
Prob
35
60
6
3
0.752
1
0.9985
0.8435
0.79
32.5
60
7
3
0.743
1
0.9995
0.7875
0.8295
37.5
60
6
3
0.7375
0.9995
0.9995
0.7855
0.8255
32.5
60
6.5
3
0.737
1
0.9975
0.821
0.7845
30
60
7.5
3
0.7335
1
0.9995
0.7775
0.8175
37.5
60
6.5
3
0.725
1
1
0.7485
0.845
35
60
6.5
3
0.7225
1
1
0.77
0.812
32.5
60
6
3
0.7195
1
0.9955
0.864
0.7415
30
60
7
3
0.717
1
0.999
0.8075
0.759
32.5
60
7.5
3
0.716
1
1
0.734
0.859
37.5
60
5.5
3
0.7145
1
0.993
0.8065
0.7565
35
60
7
3
0.712
1
1
0.731
0.8555
[1] This is only a small portion of the Monte Carlo output
.
Optimal Reaction Conditions
Marginal Probabilities
Note that the largest probability of meeting specifications is only about 0.75Slide50
50
Posterior Predicted Reliability with
Temp=20 to 70, Catalyst=2 to 12, Pressure=60, Rxntime=3.0
Catalyst
Temp
20
30
40
50
60
70
2
4
6
8
10
12
Pressure
Rxntime
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Overlapping
Mean
Contours from analysis of each response individually.
= Design Space
Contour plot
of
p
(
x
) equal to
Prob
(
Y
is in
A
given
x
& data).
The region inside the
red ellipse is the
design space.
x
1
=
x
2
=
But this
x
-point (in the yellow
“sweet spot”)
has a probability of only
0.23
!
This
x
-point (in the yellow sweet spot)
has only a probability of
0.75 .Slide51
51How was this model fit?
Prediction Models:
Model Terms
Response
x1
x2
x3
x4
x11
x22
x33
x44
x12
x13
x14
x23
x24
x34
SM Isomer
D
D
D
D
D
D
D
D
Prod Isomer
D
D
D
D
D
Impurity
D
D
D
Purity
D
D
DSlide52
52
How was this model fit? … Using a SUR ModelSUR=Seemingly Unrelated Regressions
52
Prediction Models:
Prediction Models in Matrix Form:Slide53
53
53
Fitting the SUR Model
53
Prediction Models in Matrix Form:
The
frequentist
fit: Slide54
54
54
54
Fitting the SUR Model
54
Prediction Models in Matrix Form:
The
frequentist
fit:
These two forms
strongly suggest
Gibbs sampling
for a Bayesian analysis! Slide55
55A Bayesian Analysis for the SUR Model
Consider a non-informative prior approach:
Gibbs sampling:
Notes:
See John
Geweke’s
book:
Contemporary Bayesian Econometrics and Statistics
for details.
The
rsurGibbs
function in the
bayesm
R package uses Gibbs sampling for
the SUR model.
A similar analysis can also be done with BUGS with weak or informative priors. Slide56
56How do We Sample Posterior Predictive Values
from the SUR Model?Slide57
57A
Preposterior Analysis Suppose that is not large enough.
The spread of the posterior predictive distribution depends upon both the process variation and the model parameter uncertainty.
If the sample size for the experiment was small, it may be that the model parameter uncertainty is playing a significant factor in the size
of
By “simulating” additional data from the model, one can assess how muchdata is likely to be needed to substantially reduce the model parameteruncertainty and thereby increase
Experimental regionSlide58
58
58
A
Preposterior Analysis for the SUR Model
Gibbs sampling:
Increase the rows of
X
to
correspond to additional
experimental runs.
Increase
N
to
correspond to
the additional rows
of
X.
Note: In some cases, parametric bootstrap simulations will be needed to generate
the additional data. Simulating from the posterior predictive distribution may not
work! For details see:
See
Peterson, J. J.,
Miró-Quesada, G., and del Castillo, E.
(2009), “A Bayesian Reliability
Approach to Multiple Response Optimization with Seemingly Unrelated Regression Models”,
Quality Technology and Quality Management,
6
(4), 353-369.Slide59
Design Space Scale-upfor the Maeda et al. (2012) Study
Process: Lubricant Blending process for theophylline tablets. Optimization factors: x1=Froude number, x2=blending timeResponses: Y1=“compression rate of powder mixture”
Y2=“tablet hardness (N)” Y3=“dissolution percentage at 30 min.”Specification limits: Y1<=0.32, Y2>=30, Y3>=75
59Slide60
Design Space Scale-upfor the Maeda et al. (2012) Study
Small-scale design: 32 factorial. Large-scale design: three replications at x1=0.36, x2=21 plus experimental runs at:
x1=0.20, x2=2 x1=0.40, x2=5860Slide61
Design Space Scale-upfor the Maeda et al. (2012) Study
61Optimization factors: x1=Froude number, x2=blending timeResponses: Y1=“compression rate of powder mixture” Y2=“tablet hardness (N)”
Y3=“dissolution percentage at 30 min.”Small-scale Model: Seemingly Unrelated Regressions (SUR) model.Slide62
62Design Space Scale-up
for the Maeda et al. (2012) Study
62
Optimization factors: x1=Froude number, x2=blending time
Responses: Y1=“compression rate of powder mixture”
Y2=“tablet hardness (N)” Y3=“dissolution percentage at 30 min.”Small-scale Model: Seemingly Unrelated Regressions (SUR) model.Slide63
63Design Space Scale-up
for the Maeda et al. (2012) Study
63
Optimization factors: x1=Froude number, x2=blending time
Responses: Y1=“compression rate of powder mixture”
Y2=“tablet hardness (N)” Y3=“dissolution percentage at 30 min.”Large-scale Model: A SUR-like model
is a additive adjustment coefficient.
is a multiplicative adjustment coefficient.
j
= 1,2,3Slide64
64
64
Design Space Scale-up
for the Maeda et al. (2012) Study
64
Optimization factors: x1=Froude number, x2=blending time
Responses: Y1=“compression rate of powder mixture”
Y2=“tablet hardness (N)”
Y3=“dissolution percentage at 30 min.”
Large-scale Model: A Reduced SUR-like model
is a additive adjustment coefficient.
is a multiplicative adjustment coefficient.
j
= 1,2,3Slide65
65
65
65
Design Space Scale-up
for the Maeda et al. (2012) Study
65
S
is the specification region.
simulated from the posterior
N
= 1,000 simulations
for some reliability level
R
.Slide66
66
66
66
66
Design Space Scale-up
for the Maeda et al. (2012) Study
66
Data from the Small-scale & Large-scale Experiments:
X1 X2
Response
Types
Froude
no.
Blending
Time
Y1
Y2
Y3
0.10 2 0.340 71.7 90.2
0.10 30 0.320 51.6 85.0
0.10 58 0.289 37.8 81.1
0.25 2 0.345 67.3 89.1
0.25 30 0.286 39.4 81.1
0.25 58 0.242 24.2 75.6
0.40 2 0.339 60.3 86.3
0.40 30 0.255 31.0 77.4 0.40 58 0.240 21.9 67.8 0.36 21 0.286 25.1 75.3
0.36 21 0.277 29.7 75.8 0.36 21 0.270 27.5 77.2 0.10 2 0.340 53.8 85.0 0.40 58 0.236 18.7 67.8
Large-scale data
Small-scale data
Specification limits: Y1<=0.32 Y2>=30 Y3>=75Slide67
67
67
67
67
67
Design Space Scale-up
for the Maeda et al. (2012) Study
67
Data from the Small-scale & Large-scale Experiments:
X1 X2
Response
Types
Froude
no.
Blending
Time
Y1
Y2
Y3
0.10 2 0.340 71.7 90.2
0.10 30 0.320 51.6 85.0
0.10 58 0.289 37.8 81.1
0.25 2 0.345 67.3 89.1
0.25 30 0.286 39.4 81.1
0.25 58 0.242 24.2 75.6
0.40 2 0.339 60.3 86.3
0.40 30 0.255 31.0 77.4
0.40 58 0.240 21.9 67.8
0.36 21 0.286 25.1 75.3 0.36 21 0.277 29.7 75.8 0.36 21 0.270 27.5 77.2
0.10 2 0.340 53.8 85.0 0.40 58 0.236 18.7 67.8 Large-scale data
Small-scale data
Specification limits: Y1<=0.32 Y2>=30
Y3>=75
0.10 2 0.340 71.7 90.2 0.10 2 0.340 53.8 85.0
0.40 58 0.240 21.9 67.8 0.40 58 0.236 18.7 67.8
Small-scale data
Large-scale data
Small-scale data
Large-scale dataSlide68
68Design Space Scale-up
for the Maeda et al. (2012) StudyModel fits for small-scale experiments:
Multivariate
Wilks
-Shapiro normality test on SUR residuals:
p=0.014
Multivariate normal scores plot
Scatter plot matrix of SUR residuals
No outliers detected.Slide69
69
69
Design Space Scale-up
for the Maeda et al. (2012) Study
Model fits for small-scale to large-scale analysis:
No outliers detected.
SS=“small-scale”
Large-scale
prediction
Large-scale
Design
Points X1 X2 Y1 (Y1-pred.) Y2 (Y2-pred.) Y3 (Y3-pred.) 0.36 21 0.270 (0.285) 27.5 (29.8) 77.2 (76.8)
0.10 2 0.340 (0.335) 53.8 (51.4) 85.0 (84.2)
0.40 58 0.236 (0.217) 18.7 (13.5) 67.8 (66.6)
Specification limits:
Y1<=0.32
Y2>=30
Y3>=75Slide70
70
70
70
70
Design Space Scale-up
for the Maeda et al. (2012) Study
70
S
is the specification region.
simulated from the posterior
N
= 1,000 simulations
for some reliability level
R
.Slide71
71
= 3 points
= small-scale point
= large-scale point
Both small & large scale information used
0.92 probability
0.80 probabilitySlide72
72The Posterior Predictive Approach Easily Solves the
Multivariate Robust Parameter Design Problem
72
72
.
and use
p
(
x
1
,…,
x
k
) to optimize the process.
See
Miró
-Quesada, del Castillo, Peterson (2004
)
,
Journal of Applied Statistics
for details.
Sometimes we have process factors that we can control in a lab setting but are
noisy in a manufacturing plant.
If such noise factors interact with other completely controllable factors, we may
be able to dampen the effects of the noise factors to improve the probability of
meeting specifications.
Suppose where
x
1
,…,
xk are control factors and
Xk+1,…, Xn are (random) noise factors. Slide73
73
73
Advantages of the Multivariate Bayesian Predictive Approach
to Process Optimization and Design Space
It provides a quantifiable assessment of process capability for a single
operating target or a region of process operating conditions. It considers the uncertainty of all of the model parameters.
It models the correlations among the response types.
It is easy to add noise variables.
It allows for a “preposterior” analysis to see how much further data would
reduce the model parameter uncertainty. It can be easily adapted to special desirability or cost functions. It has Bayesian “flexibility”
- Informative prior information can be used if desired. - Bayesian regularization can be used (e.g. to induce model stability or parameter constraints) - Predictive distributions for complex (e.g. nonlinear mixed effects) models can
obtained. Slide74
74Some Cautions on the Multivariate Bayesian Predictive Approach to Process Optimization and Design Space
Process capability surface is statistically more difficult to
estimate than a mean response surface. - By its definition, process capability is sensitive to distributional assumptions. - Of course, this is not a reason to avoid inference for process capability.
In most all cases, process capability cannot be quantified to depend
additionally on special cause variation.
- For example, the “probability of meeting specification” in most cases cannotbe easily modeled to take into account issues related to machine failure,
unexpected contamination in raw materials, mistakes by workers, etc.
- As such, the
overall probability of being out of spec. will be somewhat larger
than that based upon the common cause variation of the process. - Nonetheless, actions can be taken to lessen the extent of special cause
variation (better machine maintenance, better raw material acquisition & screening, better staff training).Slide75
75How Can We Deal with these Cautions?
Use good experimental design.
Use sensitivity analyses for regression models, not just priors. - If you can, fit more than one model form and compare the results.
Make sure clients and managers understand the difference between
common cause and special cause variation….particularly if the
estimated process capability is large.
If the estimated process capability is too small, first consider a “preposterior
” analysis, particularly for experiments with small
sample sizes. -
It may be that by simulating additional data, one can show that the process capability can be increased to acceptable levels by additional experiments. - If additional data still does not increase process capability far enough,
further process variation reduction and/or mean improvement may be needed. Slide76
76Some Future Challenges for Process Optimization
using a Predictive Distribution Approach Processes with many response types and associatedspecification limits.
- In biopharmaceutical manufacturing it is common to have about a dozen or more different response types, each with their own specification limits that need to be met. - Often, too little experimental runs will be available with which to check distributional assumptions in a clear way. -
Current research underway involving use of a special desirability function to
reduce the dimensionality of the problem. Processes that have many latent variables.
- Modeling of raw material properties can help with building more realistic predictive models for pharmaceutical manufacturing.
- However, it is still “early days” for Bayesian predictive modeling using these types of models. Slide77
77
77
Some Future Challenges for Process Optimization
using a Predictive Distribution Approach
Functional responses
- Some manufacturing processes produces a functional trace or profile as a o measured response.
del Castillo, E.,
Colosimo
, B. M., and
Alshraided, H., (2012), “Bayesian Modeling and Optimization of Functional Responses Affected by Noise Factors”, Journal of Quality Technology 44, 117-135.Slide78
78
78
Summary
Classical response surface methodology (based on means) has had a profound
influence on process optimization…but more needs to be done!
The “more” involves thinking about processes as entities which producedistributions
of results over time.
In other words, process optimization needs to be a
distribution oriented
endeavor. “Quality improvement” has been defined in a nutshell as “the reduction in variation about a target”. This can be generalized to “the shrinking of a multivariate distribution about a vector of quality target values”.
The Bayesian posterior predictive distribution can be an effective tool for quality improvement, particularly for processes that involve complex modeling. Slide79
79Some Further Applications Slide80
80
Bayesian Monte Carlo Studies for USP Test Assessment
Many USP tests involve multi-stage algorithms. - Here it is not obvious what the probability of acceptance or rejection would be for processes with assumed means and variance components (within and between batch
variances). Slide81
81
81
Bayesian Monte Carlo Studies for USP Test Assessment
Many USP tests involve multi-stage algorithms.
-
Here it is not obvious what the probability of acceptance or rejection would be for processes with assumed means and variance components (within and between batch
variances). Suppose we have
b batches of data
- Let A
(X) = 1 if the batch is accepted,
0 if rejected., where A is a USP data algorithm. - Assumed model: - Consider the reliability measure:
- Based upon data, X, from several batches we can compute a posterior distribution
for to make process qualification decisions about a process. - For example we could compute the posterior probability
batch 1,……...…..…, batch
bSlide82
82
82
Bayesian Monte Carlo Studies for USP Test Assessment
For further details see:
LeBlond , D. and Mockus, L. (2014) “The Posterior Probability of Passing a Compendial
Standard, Part 1: Uniformity of Dosage Units”, Statistics in Biopharmaceutical Research,
DOI: 10.1080/19466315.2014.928231
Using Monte Carlo simulations to generate data from given numbers of batches, we can determine the number of batches (of a given size) needed to have a
specified degree of certainty for process qualification. - This determination would make use of for simulated values
of x.
A less conservative approach could involve the use of instead of Note: Of course, prior information about the product quality endpoint (content uniformity,
dissolution, etc.) can be used to possibly reduce the number of batches needed for process qualification. Slide83
83
83
Assay Development
Assessing Assay Parallelism over a Pre-specified Interval , [
xL, xU]
- The reference standard and test sample are similar if they contain the same effective constituent. - Under such condition, the test sample behaves as a dilution
(or concentration) of the reference standard.
Graphically, this line is a copy of shifted by a fixed amount on the log-concentration axis.
In other words, there exists a number such that for all x.
This calibration constant is commonly known as log-relative potency of the test sample.
x
L
x
USlide84
84
84
84
Assay Development
Assessing Assay Parallelism over a Pre-specified Interval , [
x
L
,
xU]
If two lines are ”close” to parallel
over a finite interval, then there
should exist a shift that will bring the lines close together over that interval.
As such, the criterion below
could be used to quantify assay
parallelism in a practical way, over
a finite interval. Slide85
85
85
85
85
Assay Development
Assessing Assay Parallelism over a Pre-specified Interval , [
x
L
,
x
U
]
Therefore, one can assess assay parallelism (for linear assay profiles) by computing
For
sigmoidal
assay profiles, one can assess parallelism by computing
For details see:
Novick, S. J., Yang, H., and Peterson, J. J., (2012) “A Bayesian Approach to Parallelism in Testing in Bioassay”.
Statistics in Biopharmaceutical Research
,
4
(4) 357-374. Slide86
86
86
86
Assay Development
Assessing Assay Parallelism over a Pre-specified Interval , [
x
L
,
xU
] - Note, for sigmoidal assay profiles, some people might prefer to assess differences after a logistic transformation has been done so that comparisons
will be between linear forms. Slide87
87
87
87
87
Assay Development
Assessing Assay Parallelism over a Pre-specified Interval , [
x
L
,
x
U
] - Note, for sigmoidal assay profiles, some people might prefer to assess
differences
after
a logistic transformation has been done so that comparisons
will be between
linear
forms.
But the transformed data depends upon the (unknown) parameters
A
and
B
! Slide88
88
88
88
88
88
Assay Development
Assessing Assay Parallelism over a Pre-specified Interval , [
x
L
,
x
U
]
- Note, for
sigmoidal
assay profiles, some people might prefer to assess
differences
after
a logistic transformation has been done so that comparisons
will be between
linear
forms.
But the transformed data depends upon the (unknown) parameters
A
and
B
!
This does not stop the Bayesian approach! Slide89
89
89
89
89
89
89
Assay Development
Assessing Assay Parallelism over a Pre-specified Interval , [
x
L
,
x
U
]
- Note, for
sigmoidal
assay profiles, some people might prefer to assess
differences
after
a logistic transformation has been done so that comparisons
will be between
linear
forms.
Compute instead: Slide90
90
90
Assay Development
Assay ruggedness
- Bayesian design space modeling ideas can be used here as well.
- Suppose one had several factors, say x
1
, x
2,…, x
8 , in a “screening design”over a small (robustness) factor region, X.
- Suppose also that one had one (or more) process response types:
- One could compute
For further details see:
Peterson, J. J. and Yahyah, M., (2009) "A Bayesian Design Space Approach to Robustness and
System Suitability for Pharmaceutical Assays and Other Processes",
Statistics in
Biopharmaceutical Research
1
(4), 441-449.
Specification
regionSlide91
91
91
Assay Validation
The “closeness” of a result
X
to the unknown true value of the sample, T, is simultaneously linked to both the size of the bias and precision of the method.
The “total error concept” in analytical method validation: A Bayesian perspective.
Suppose we have an assay response,
X
, and a true (gold standard) value, T.
Suppose further that
Note: Consider
This is equivalent to
For further details see:
Boulanger, B. et al. (2007), “Risk management for analytical methods based
on the total error concept: Conciliating the objectives of the pre-study and in-study validation phases”,
Chemometrics
and Intelligent Laboratory Systems
86, 198–207.Slide92
The f2 statistic is commonly used to make decisions about equivalence of referenceand test batches.
92
Dissolution
Consider the dissolution data for a test and reference
batch:
Reference
Reference
Test
TestSlide93
93
93
Dissolution
The f
2
statistic is commonly used to make decisions
about equivalence of reference and test batches.
F
2 has been proposed as a population-based analogue of f
2. Slide94
94
94
94
Dissolution
Guidance for Industry: SUPAC-MR:
Modified Release Solid Oral Dosage Forms. Scale-Up and
Postapproval
Changes: Chemistry, Manufacturing and Controls; In Vitro Dissolution Testing and In Vivo Bioequivalence Documentation.
“An f2 value between 50 and 100 suggests the two dissolution profiles are similar. Also, the
average difference at any dissolution sampling time point should not be greater than 15% between
the changed drug product and the
biobatch
or marketed batch (unchanged drug product)dissolution profiles.”
f
2
=52.3
16%
So we might want to consider: Slide95
95
95
95
95
Dissolution
f
2
=52.3
16%
However,
frequentist
inference related to
could be difficult to implement.
But, if we have the posterior for the model parameters, the computation of
is straightforward (at least from a Monte Carlo approach).Slide96
96
Some Computational Recommendations Generally, the most difficult part of a Bayesian analysisis computing the posterior distribution of model parameters.
Posterior predictive distributions and associated risk probabilities are more straightforward to compute. As a statistical consultant, it is advisable to have one (or
more) backup strategies for statistical computing. Sooner or later,
WinBUGS or PROC MCMC will fail to produce adequate results!
Slide97
97
97
Some Computational Recommendations
“Access” to other Bayesian applications or specialty R
packages can be useful.
In addition to the software, “access” can mean - Working knowledge of the software.
or - Access to someone who is knowledgeable and willing to help to help you get up to speed quickly.
But it can be useful to know, and often use, 2 or 3 different
Bayesian applications. - MCMC algorithms can be delicate and it may be prudent to
check any important results using a different algorithm. Slide98
98
98
Some Computational Recommendations
Some alternative Bayesian applications to
WinBUGS, jags, or
PROC MCMC in SAS:
STAN - Uses the ‘No-U-Turn Sampler’ algorithm.
- Generally faster convergence for analysis situations where
BUGS will have problems. - Can be invoked in R using the
RStan package. Generalized Direct Sampling (GDS)
- Does not use MCMC, but rather a sophisticated rejection algorithm. - GDS produces independent samples from the posterior! - Requires finding the posterior mode.
- There is a R package for this: bayesGDSSlide99
99
99
99
Some Computational Recommendations
Some alternative Bayesian applications to
WinBUGS
, jags, or PROC MCMC in SAS:
The R package,
MCMCglmm
-
glmm
stands for “generalized linear mixed-effect models” - It can also sample from the posterior for multivariate
normal mixed-effect models. - It can sample from the posterior for binary, multinomial, and Poisson models. Slide100
100Acknowledgements
Steven Novick Harry Yang
Stan Altan David LeBlond Yan Shen Bruno Boulanger Mohammad Yahyah
Kevin LiefSlide101
101Some Useful Bayesian Books
Christensen, R., Johnson, W., Branscum, A., and Hanson, T. E. (2011).
Bayesian Ideas and Data Analysis: An Introduction for Scientists and Statisticians, Chapman & Hall, CRC, Boca Raton, FL. del Castillo, E. (2007), Process Optimization - A Statistical Approach, Springer,
New York, NY. (Chapters 11 & 12 discuss Bayesian analysis for process optimization)Krusche
, J. K. (2010), Doing Bayesian Data Analysis: A Tutorial with R and BUGS, Academic Press, Oxford, UK.
Lunn, D., Jackson, C., Best, N., Thomas, A., Spiegelhalter, D., (2013
), The BUGS Book – A Practical Introduction to Bayesian Analysis, CRC Press, Boca Raton, FL.
Ntzoufras, I. (2009). Bayesian Modeling Using
WinBUGS, John Wiley and Sons, Inc. Hoboken, NJ.
Gellman, A. and Hill, J, (2006), Data Analysis Using Regression and Multilevel/Hierarchical Models, Cambridge University Press, NY, NY. Slide102
102Concerned about becoming a Bayesian?
You might consider reading the recent book: The Theory That Would Not Die: How Bayes' Rule Cracked the Enigma Code,
Hunted Down Russian Submarines, and Emerged Triumphant from Two Centuries of Controversy by Sharon Bertsch McGrayne. Slide103
103
Bernardo, J. M. and Irony, T. Z. (1996), "A General Multivariate Bayesian Process Capability Index", The Statistician, 45, 487-502.del Castillo, E. (2007), Process Optimization - A Statistical Approach, Springer, New York, NY.
(This is a good book for Bayesian process optimization) del Castillo, E., Colosimo, B. M., and Alshraided, H., (2012), “Bayesian Modeling and Optimization of
Functional Responses Affected by Noise Factors”, Journal of Quality Technology 44, 117-135.
Duncan, A. J. (1986). Quality Control and Industrial Statistics. Irwin, Homewood, IL.
Flaig, J. J. (1999), “Process Capability Sensitivity Analysis”, Quality Engineering, 11(4),
587 - 592 Fu, Z., Leighton, J., Cheng, A., Appelbaum, E., and Aon, J. C. (2012) “Optimization of a
Saccharomyces cerevisiae fermentation process for production of a therapeutic recombinant protein using a multivariate Bayesian approach”,
Biotechnology Progress (to appear).
Gelman, A. and Hill, J., (2007) Data Analysis Using Regression and Multilevel/Hierarchical Models, Cambridge University Press, New York, NY. (This is a good book for Bayesian mixed effect modeling.)
ReferencesSlide104
104
104
References
(continued)
Jeang, A. (2010), “Optimal Process Capability Analysis”, International Journal of
Production Research, 48(4), 957-989.LeBlond , D. and Mockus, L. (2014) “The Posterior Probability of Passing a Compendial
Standard, Part 1: Uniformity of Dosage Units”, Statistics in Biopharmaceutical Research,DOI: 10.1080/19466315.2014.928231
Maeda, J., Suzuki, T., and
Takayamab, K. (2012) “Novel Method for Constructing a Large-ScaleDesign Space in Lubrication Process by Using Bayesian Estimation Based on the Reliability of
a Scale-Up Rule”, Chem. Pharm. Bull. 60(9) 1155–1163. Miró
-Quesada, G., del Castillo, E., and Peterson, J.J., (2004) "A Bayesian Approach for Multiple Response Surface Optimization in the Presence of Noise variables", Journal of Applied Statistics, 31, 251-270Montgomery, D. C. (2009),
Introduction to Statistical Quality Control (6th ed), John Wiley & Sons,
Inc., Hoboken, NJ.Ng, S. H. (2010), “A Bayesian Model-Averaging Approach for Multiple-Response Optimization”, Journal of Quality Technology, 42, 52-68. Peterson, J. J., (2004) "A Posterior Predictive Approach to Multiple Response Surface
Optimization", Journal of Quality Technology, 36, 139-153.Slide105
References (continued)
105Peterson, J. J. (2006), "A Review of Bayesian Reliability Approaches to Multiple Response Surface Optimization", Chapter 12 of
Bayesian Statistics for Process Monitoring, Control, and Optimization , pp 269-290, (eds. Colosimo, B. M and del Castillo, E.) Chapman and Hall/CRCPress Inc.
Peterson, J. J. (2007), “A Review of Bayesian Reliability Approaches to Multiple Response
Surface Optimization”, in Bayesian Process
Monitorng , Control & Optimization,
Chapman & Hall/CRC.Peterson, J. J. (2008), “A Bayesian Approach to the ICH Q8 Definition of Design Space”,
Journal of Biopharmaceutical Statistics
, 18, 959-975.
Peterson, J. J., Miró-Quesada, G., and del Castillo, E. (2009), “A Bayesian Reliability
Approach to Multiple Response Optimization with Seemingly Unrelated Regression Models”, Quality Technology and Quality Management, 6(4), 353-369.
Peterson, J. J. (2009) “What Your ICH Q8 Design Space Needs: A Multivariate Predictive Distribution”, Pharmaceutical Manufacturing, 8(10) 23-28.
Peterson, J. J. and Lief, K. (2010) "The ICH Q8 Definition of Design Space: A Comparison of theOverlapping Means and the Bayesian Predictive Approaches", Statistics in Biopharmaceutical
Research,
2
,249-259.Slide106
106References
(continued)Plante, R. D. (2001), “Process Capability: A Criterion for Optimizing Multiple Response
Product and Process Design”, IIE Transactions, 33, 497-509.Polansky, A.M. (2001), “A Smooth Nonparametric Approach to Process Capability”, Technometrics, 43(2), 199-211.
Rajagopal
, R. and del Castillo,, E. (2005), “Model-Robust Process Optimization Using Bayesian Model Averaging”,
Technometrics, 47, 152-163
Rajagopal, R., del Castillo, E., Peterson, J. J. (2005), "Model and Distribution-Robust Process
Optimization with Noise Factors", Journal of Quality Technology
37, 210-222.
Rajagopal, R. and del Castillo, E. (2006), “A Bayesian approach for multiple criteria decision making with applications in Design for Six Sigma”,
Journal of the Operational Research Society, 1–12.Robinson, T.J., Anderson-Cook, C.M. and Hamada, M.S. (2009). “Bayesian Analysis of
Split-Plot Experiments with Non-Normal Responses for Evaluating Non- Standard Performance Criteria”. Technometrics, 51, pp. 56-65.
Savage, S. (2009) The Flaw of Averages - Why We Underestimate Risk in the Face of Uncertainty , John Wiley and Sons, Inc. , Hoboken, NJ.Slide107
107
Stockdale, G. and Cheng, Aili (2009), “Finding Design Space and a Reliable OperatingRegion using a Multivariate Bayesian Approach with Experimental Design”, Quality Technology and Quantitative Management
6(4),391-408
References (continued)