POPT - Installation and user guide.
To run POPT and evaluate
sampling windows
Input files for running
POPT and evaluate sampling windows
(from WinPOPT User Guide)
A brief description of optimal design is provided here. For those readers who are not familiar with
optimal design an introduction is provided in: Monteleone
JPR and Duffull SB. Choice
of best design. In: Kimko HC and Duffull SB. Simulation for Designing
Clinical Trials. A pharmacokinetic-pharmacodynamic modeling perspective. Marcel
The purpose of optimal experimental design is to obtain data that
provides the most information about a model and its parameter values from a new
experiment. This contrasts with
modelling which proposes to provide the best, unbiased parameter estimates,
given a model and a design from which the current data arose.
Optimal design is a phrase used to describe a general set of design
methods that enable an optimum design to be constructed. These methods are often described as
information theoretic since they do not rely on the availability of real or
simulated data. However, it should be
noted that these methods are widely used in other industries and are gaining
popularity in drug development with numerous prospective applications of
optimally designed studies.
The general approach involves computation of the Fisher information
matrix (see bibliography for selected relevant readings). For nonlinear models, the Fisher information
matrix can be calculated analytically given any model, set of parameter values
and a candidate design. The inverse of
the Fisher information matrix is an estimate of the covariance matrix. The square roots of the diagonal elements of
the covariance matrix are the expected standard errors of the parameter
estimates (which is approximately equivalent to the output of NONMEM when
asking for $COV). The optimal design is
therefore a design that optimises a summary measure of the Fisher information
matrix. Generally, in pharmacokinetics
(PK) and pharmacokinetics-pharmacodynamics (PKPD) it
is usual to use the determinant of the matrix as the summary measure. Statistically it can be shown that maximizing
the determinant of the Fisher information matrix is the same as minimizing the
determinant of the covariance matrix.
Therefore the optimum design that maximises the determinant of the
Fisher information matrix will provide the smallest standard errors.
Since PK and PKPD models in populations are hierarchical and nonlinear
then the standard errors are approximate and asymptotic. The design is also dependent on the values of
the unknown parameters. It is common
therefore in optimal design to consider that the parameter values and perhaps
also that the model have some level of uncertainty.
Choosing the optimal design
The concept optimal design engenders the notion that there is one single
design which is better than all the rest and that this
design is in some manner the most desirable design. While it is true that the optimal design is
the best single design it is not necessarily the most desirable design. Indeed, there is, as yet, no proof that it is
even possible to locate the optimal design for a nonlinear mixed effects
model.
Similarly searching for the optimal design is a non-trivial
process. The determinant of the Fisher
information matrix over the design space is often a very convoluted surface and
many standard (local) search strategies, e.g. gradient methods, can become easily
trapped in local minima. It is usual
therefore to choose search methods which are considered to be global, or near
global procedures. Global means the
algorithm given enough time would look almost everywhere in the design
space. Local means that the algorithm
follows a specific path with little if any deviation to explore surrounding
areas. Two main advantages of global
search algorithms are 1) that it doesn’t matter where you start (not dependent
on initial estimates of the sampling times) and 2) they generally do not have
convergence problems. The principle
disadvantage is their very slow convergence.
WinPOPT has two (near) global search
algorithms. These are simulated
annealing (an asymptotically global search algorithm) and an exchange algorithm. Simulated annealing can be very slow but
provides good convergence properties.
The exchange algorithm is much quicker but has less rigorous convergence
properties. The exchange algorithm can
be used for quick exploration of the design space, but the final design should
be based on a run from simulated annealing.
It should be noted that simulated annealing and the exchange algorithm
always converge. However, if their
search is terminated, then the sampling times provided at this juncture are not
subject to bias. This contrasts with
modelling using maximum likelihood estimation methods (that most users will be
familiar with) where a lack of convergence implies that the parameters are not
unbiased best least squares estimators and therefore probably contain some
element of bias. So to recap, sampling
times gained during a non-converged optimisation run using WinPOPT
are just that, they are sampling times that provide the current estimates of
the standard errors of the parameters.
Interpreting the design
output
When exploring the design space using the determinant of the Fisher
information matrix there are various metrics that are provided in the
output. Not all have practical
significance.
|
Metric |
What is it? |
Practical use? |
|
Determinant |
The determinant is calculated mathematically from
the Fisher information matrix and is a summary measure of the volume of a
matrix. |
The optimal design will be the one that provides
the highest value of the determinant.
This metric has no practical utility and cannot be interpreted by
itself. |
|
Criterion |
The determinant to the power of 1 over the number
of parameters in the model. |
The ratio of 2 criterion for any given model can
be expressed as a fraction (or %) and is called efficiency; e.g. design A is
30% more efficient than design B. A
criterion by itself has no meaning. |
|
Eigenvalues |
These are the eigenvalues
of the Fisher information matrix. |
The ratio of the highest to lowest eigenvalue can be expressed as a condition number. The condition number is a measure of the
amenability of the matrix to mathematical computation. (Smaller numbers are better). |
|
Standard errors |
These are the asymptotic estimates of the
estimation standard deviation for a parameter. They are the square roots of the diagonal
elements of the inverse of the Fisher information matrix. |
Values less than some predefined limit (e.g. <
20% for fixed effects parameters) can be used to judge the utility of a
design |
From a practical perspective, the goodness of any particular design
should be judged on its ability to solve the problem at hand; i.e. provide
precise estimates of the parameters of interest and is able to be performed
clinically. For comparison between
designs then criterion values can be used as their ratio to provide an
efficiency rating.
POPT is a set of programs that have been written for MATLAB to help in
the design of clinical studies in which PK or PKPD data is to be collected for
the intention of developing models. POPT
has mostly been written and developed in MATLAB ver
6.x, but is backwardly and upwardly compatible with MATLAB ver
5.x and ver 7.x.
POPT does not require any specialised MATLAB toolboxes to be
available. The provision of toolboxes
will not improve the performance of POPT code as it is currently written, but
it is likely that some amendments could be made to increase performance.
The programs will run on any operating system platform (e.g. Windows,
Linux ...) under which MATLAB has been successfully installed. We have tried both Windows and Linux without
problems.
POPT is available free of charge and may be downloaded at
http://www.uq.edu.au/pharmacy/sduffull/popt.htm
POPT is subject to copyright.
Please refer to section xx for details of the copyright license.
There are no special instructions for installation of POPT. POPT is made available in a zipped format
which contains a directory of files "POPT Distribution copy" and a
subdirectory "Example files".
These directories can be unzipped to any location on your computer.
It is recommended that you then copy the contents of the directory
"POPT Distribution copy" to a new directory in your current work
space for MATLAB in order to ensure that you do not inadvertently change any of
the download files. From your Current Directory in MATLAB you can highlight and
run the appropriate programs.
The
POPT suite of programs BACK TO TOP
There are 21 objects in the main directory of POPT.
Programs for optimal design - DO NOT EDIT THESE
|
disp_2_scrn_sa.m |
graphical interim results of simulated annealing algorithm |
|
display_to_screen.m |
graphical interim results of exchange algorithm |
|
exch_setp.m |
initialization of the exchange algorithm |
|
exchange.m |
exchange algorithm |
|
exchange_with_v6_patch.m |
this is the
same as the exchange algorithm but is designed to work with version 6 of
MATLAB (or any version that does not have access to the subroutine unidrnd()). If the usual version of the exchange
algorithm does not work then rename exchange.m to
exchange_v7.m and rename exchange_with_v6_patch.m to exchange.m |
|
p_error.m |
checks for some user errors in POPT_INI.m |
|
p_expect.m |
computes expectation of det(information matrix) over upper
and lower bound |
|
p_nars.m |
searches for the max(det(information matrix)) by random
sampling |
|
p_output0.m |
produces output to screen and file for single model
optimizations |
|
p_output1.m |
produces output to screen and file for multiple model
optimizations |
|
p_sa_a.m |
simulated annealing algorithm for approximate design |
|
p_sa_dn.m |
simulated annealing algorithm for optimizing dose number
for samples |
|
p_sa_m8.m |
general simulated annealing algorithm |
|
pop_run.m |
sets up matrices for computation of the information matrix |
|
pop_runs.m |
evaluates information matrix and computes its determinant |
|
popt.m |
RUNS POPT (RUN THIS FILE) |
Programs for setting up optimal design - EDIT THESE AS
REQUIRED
|
P_MODEL.m |
enter your structural model here |
|
POPT_INI.m |
enter your design specification and model details here |
Programs for evaluating sampling windows - DO NOT EDIT
THESE
|
p_robutst_1.m |
RUNS sampling window estimation (RUN THIS FILE) |
|
p_robust_popt.m |
sets up sampling window estimation |
|
repeat_part_3.m |
resamples over
sampling windows |
Programs for evaluating sampling windows - EDIT AS
REQUIRED
|
POPT_INI_robust.m |
enter your sampling window and design specification and
model details here |
Miscellaneous
|
variable_names.txt |
a list of variable names |
|
P2NM |
Creates a NONMEM data file from the current design in POPT_INI.m |
To run POPT change MATLAB to the current directory that the POPT
programs are located and either type POPT at the command window prompt or open
POPT and click on the run button.
POPT will then run and produce output to a file which is echoed to the
screen. Outputs to files are appended -
so no previous runs will be overwritten.
This means that all runs will be appended one after the other and the
output file can get very large. The output file is generated automatically and
its naming is based on the optimization method, e.g. method0.res is the file
name when no optimization method is selected (i.e. method=0).
To evaluate sampling windows run P_ROBUST_1.
Input
files for running POPT BACK TO TOP
Example of POPT_INI_ROBUST file
POPT does not have an advanced or particularly user friendly method for
entering your models and designs. The
structural model must be written in the file P_MODEL.m
and the design component must be written in the file POPT_INI.m. Because of this inefficiency it is
recommended that you create a new directory (copying all POPT files) for each
new design project. The complete design
is printed out with the output files and therefore you will always have a
record of the design that you tried.
Examples of P_MODEL are provided for
- 2 compartment
model with first order input
- 1 compartment
model with zero order input
- 1 compartment
model with first order input defined by ODES
- multiple
models where a 1 compartment model and 2 compartment model are written
- multiple
response (PK and PD)
Each P_MODEL.m file contains a detailed
description of how to construct the structural model for each of these
occasions and further description is not provided here.
Examples of POPT_INI are provided for
- multiple
models where a 1 compartment model and 2 compartment model are written
- multiple
response (PK and PD)
POPT_INI is a relatively detailed (and therefore complex) input file. A description of each item is provided.
MULTIPLE MODELS: This situation arises when you want to
optimize a design for more than 1 model at the same time. This may be because
- you have 2
competing models and are not sure which is true,
- you have 2
drugs which you want to optimize the design for at the same time,
- you have 1
drug but 2 different sets of parameter values (e.g. for a paediatric and
adult population)
MULTIPLE RESPONSES: This situation arises when you have more than
1 response from a model. Consider the
framework that A à B à C, then C is dependent on B which is dependent
on A. The information matrix for C will
therefore also have information about parameters that describe the formation of
B from A as well as those that describe B à C. Situations where multiple
response models occur can include:
- parent,
metabolite PK models
- plasma, urine
PK models
- PKPD models
Example POPT_INI.m file. BACK
TO TOP
There are 5 parts to this file.
Part 1 deals with
model features
The description of the POPT_INI file will be divided into these
sections.
MODEL FEATURES Back to
this example
weighting=[];
Use this function only if you have either multiple models or multiple
responses. In this case you need to
enter a weighting value for each candidate model. The weightings need not add to 1, but this is
convenient. e.g.
use:
weighting = [0.5 0.5];
mr = 0;
POPT cannot determine the difference between a multiple response model
and multiple competing models. Use this
flag variable to indicate to POPT which it is.
If mr = 0 then this is a multiple model and if mr = 1 then it is a multiple response.
BETAI=[];
These are your prior values of the model parameters. If you have more than one set of parameters
(multiple models or multiple response models) then enter a row for each
model. You cannot have more than 2 rows
for a multiple response model. You can
have any number of rows for a multiple model.
BETAI must be symmetric. If there
are an uneven number of parameters between model 1 and model 2 then you need to
insert place holders, e.g. notation:
BETAI = [1
4 20 1 1
1 2 20
4 5];
in this example the first row has two 1's as
placeholders. Do not use 0's as
placeholders.
POP_PAR(1,:)=();
Enter the parameter names. This
variable is free text. Example use
POP_PAR(1,:)=(' Ka
CL V');
If you have 2 models then you will need two rows of names. The total length of the string vector must be
the same for each, example use
POP_PAR(1,:)=(' Ka
CL V ');
POP_PAR(2,:)=('
Ka CL V
Emax EC50');
here the spaces in POP_PAR(1,:) are
placeholders to make the POP_PAR matrix square.
FXFI=[];
Use this vector or matrix to fix none, one or more fixed effects
parameters (i.e. to remove them from consideration for optimization). Enter the position number of the parameter in
BETAI that you wish to fix. For a
single model an example use is
FXFI=[2];
this will fix parameter 2. If BETAI has more than 1 row then FXFI should also have
more than 1 row. You can use 0's as placeholders, example use
FXFI=[2
0];
in this syntax position 2 of the first row of BETAI is fixed but no
positions are fixed in the second row.
OI=[];
OI are the diagonal elements of the variance-covariance
matrix of random effects. Only the
diagonal elements can be entered, example use
OI=[0.25 0.1 0.1];
There needs to be an element in OI for every element
in BETAI. The value 0.00001 should be
used as placeholders for values that are not intended to be estimated - and
also FXI should indicate a fixed position.
FXI=[];
Use this vector or matrix to fix none, one or more variance of the
between subject random effects parameters (i.e. to remove them from
consideration for optimization). Enter
the position number of the parameter in OI that you wish to
fix. For a single model an example use
is
FXI=[1];
This fixes position 1 in the OI matrix.
BSVMOD=[];
This vector (or matrix if it is a multiple model) indicates the
structural form of the between subject variance-covariance matrix. A 0 indicates an additive (normal) model and
a 1 an exponential (lognormal) model.
Example use
BSVMODI=[1 1 1];
The dimensions of BSVMOD must equal the dimensions of OI which must equal
the dimensions of BETAI. You can use
either a 0 or 1 for a placeholder for BSVMOD - they are simply ignored if
not needed (which is determined by FXI).
sd_prior=[];
sd_addi=[];
These two terms are used to provide the prior values for the standard
deviation of the residual variance.
Either or both additive and proportional error terms can be added. It is always recommended to include an
additive term, even without prior evidence, as this ensures that the
optimization routine does not equate smaller values of the response variable
with greater precision
in its estimation, which can lead to the erroneous
situation where an infinitely small concentration is without error. An example of use is
sd_propi=[0.15];
sd_addi=[0.1];
If you have multiple models then an example of use is
sd_propi=[0.15
0.2];
sd_addi=[0.1
0.1];
These variables need to be the same number of rows as BETAI.
FXRI=[];
Use this vector or matrix to fix none, one or more variance of the
residual random effects parameters (i.e. to remove them from consideration for
optimization). Enter the position number
of the parameter that you wish to fix.
For a single model an example use is
FXRI=[2];
This fixes sd_addi. Whereas:
FXRI=[1];
fixes sd_propi
For a multiple model the syntax is not as obvious but to fix the
additive component from model 1 and the proportional component from model 2 you
would write:
FXRI=[2
1];
DESIGN FEATURES Back to
this example
ns =
maxs =
These variables are used with (method=6) and allow for an approximate
design to be constructed. ns is the total number of subjects in the study and maxs is the total number of samples allowed.
POPT will then determine the best arrangement of dividing the samples up
amongst the patients. Example syntax:
ns = 100;
maxs = 500;
The arrangement will be based on the allocation of samples per group
provided by NT. This method optimizes the number of patients per
group and the sampling times.
NUM_SUB=[];
This variable is used with all other methods (i.e. method <> 6) and
replaces ns and maxs. NUM_SUB is used to define
the number of subjects in each group.
Example use:
NUM_SUB=[100 100];
Where there are 100 patients allocated to group 1 and 100 patients
allocated to group 2. Each group has a
different design.
NT=[]
NT defines the number of sampling times per group. Example use:
NT=[3 3];
Indicates that group 1 has 3 sampling times as does group 2.
D=[];
This signifies the dose to be used for each group. Example use:
D=[100 100];
If you have multiple models (or a multiple response model) then dose
needs to be entered for each model. This
is especially important for multiple models where the models may represent
different drugs. Example use:
D=[100 100
100 100];
dose_zero_fix=;
This is a special feature that allows the user to set the information
matrix to a matrix of zeros when no dose is to be given. This is of importance
when model responses are present in the absence of dosing (e.g. a baseline
value). If the value is 1 then a dose of
zero yields an information matrix of zeros.
If the value is 0 then a dose of zero yields an information matrix
according to whatever components of the model can be estimated for this dose
level. Example use
Dose_zero_fix=0;
DI=[];
DI is the dose interval.
There must be a dose interval for every element in the dose matrix (D). Example use
DI=[24 24];
Example use for multiple models:
DI=[24 24
24 24];
DPM=[];
DPM are the initial estimates of the design points. In this case design points are sampling
times. This should be written as a row
vector, which can be expressed as a matrix (although still read as a column
vector) by the use of the
“…”, see example, the syntax
DPM=[0.1 1 4
...
2 6 12];
Is equivalent to
DPM=[0.1 1 4 2
4 12];
DPM must contain the same number of elements as the sum
of NT. All elements of DPM must fall between
the lower (LBA) and upper bounds (UBA) provided below.
DPE=[];
DPE are discrete potential values of the design variables
that may be chosen for optimization. These are only used with method=9 and is
optional. Example use
DPE=[0.01 0.5 1 2 3 4 5 6 7 8 9 10 11 12 18 24 …
0.01 0.5 1 2 3 4 5 6
7 8 9 10 11 12];
In this setting the exchange algorithm would select possible best
sampling times from only those provided in the list above. The row vector may be any length. Do not use
0 (zero) as it is a reserved variable in the exchange algorithm. It is usual not to enter these values, but
rather let the exchange algorithm select discrete values from within the
sampling space itself. Note in this example that the second sampling time
cannot exceed 12.
DPE_index=[];
This row vector provides the break points for DPE so that the
exchange algorithm can determine which set of discrete sampling times should be
used for each group. Example use (based
on DPE example above):
DPE_index=[16 14];
This indicates that group 1 has a total of 16 times from which to choose
but group 2 has only 14. The number of elements in DPE_index must equal the
number of elements in NT.
FX_DPM=[];
FX_DPM provides the index for any sampling times in DPM that
are fixed and are not permitted to be optimized. The default is not to fix any sampling times,
see example
FX_DPM=[];
Leaving the row vector blank is read by POPT as entering a row vector of
0’s the same length as DPM. To fix a single sampling time (e.g. the 12 hour
sampling time from group 2 only) then a 1 is used in the same position in the
row vector for FX_DPM as the position of 12 is in DPM. The other non-fixed elements in DPM are given a
0. Example use
FX_DPM=[0 0 0 …
0 0 1];
DN=[];
The dose number for which the samples are to be taken is provided by DN. For instance a DN value of 5 would
indicate that the sample should be taken on the 5th dose, i.e. on
the dose following 4 previous doses given at the specified dose interval. DN is a row vector that should
be of equal length as DPM. The default
value of DN is a row vector of 1’s, indicating that all samples are taken from the
first dose. Example use
DN=[];
If samples are to be taken from anything other than the first dose then
the dose number for each sample must be provided. The length of DN must be the same
as the length of DPM. Example use
DN=[4 4 4
…
4 4 1];
In this example all samples are taken from the 4th dose,
except the 3rd sampling time from group 2 which is taken on the 1st
dose.
MX_DN=[];
The dose number that any sample is taken can be optimized using method=7. In this setting MX_DN provides the
upper bound for the dose number that doses can be taken from. If it is allowed that samples are taken from
any dose up to steady state then it is recommended to set the MX_DN to be the dose
where approximately 5 half-lives have elapsed, otherwise the optimization
algorithm will have difficulty distinguishing the difference between a dose
after 10 half-lives vs a dose after 20
half-lives. Example use
MX_DN= [5 5 5 …
1 1 1];
In this example, the first group of subjects can have their samples
taken on any dose up to and including the 5th dose, while group 2
subjects will have all of their samples taken on the first dose.
UBA=[];
This is the upper bound on the optimized sampling times. The default is the upper bound equals the
dose interval. Example use
UBA=[];
Any elements in UBA cannot be greater than the dose interval. If some further sampling constraints remain,
e.g. a patient must leave a clinic at 6 hours after their dose then individual
upper bounds can be used. Example use
UBA=[24 24 24
…
12 12 12];
Here group 1 patients have an upper bound equal to the dose interval (DI), but group 1
patients are constrained to have their samples fall within the first 12 hours
post-dose only. UBA cannot be smaller
than LBA.
LBA=[];
This is equivalent to UBA but for the lower bound. The default lower bound is immediately
post-dose, provided by a row vector of zeros.
OPTIMIZATION
FEATURES Back
to this example
output_on=;
this option toggles GUI output. If output_on=1 then
the interim GUI is produced however if output_on=0
then the GUI screen output is set to off.
Optimization runs slightly slower with the GUI interim output feature
set to on. Example syntax
output_on=1;
prnt=;
prnt is the number of
iterations between the refresh of the graphical and tabulated results. The more frequently interim results are
echoed to the user the slower the algorithm will run. However, if the model is extremely complex
then up to an hour may elapse before 500 iterations have passed. In general, prnt can be set to 500
for simulated annealing (method=6, 7, 8) and 50 (method=9)for the exchange algorithm.
re_scale=;
When optimizing the user has the opportunity to re_scale the design variables so that they all have a value of
1.0. This may be of value if one design
variable has an initial value of 0.01 and another has a value of 1000, which
may result in numerical difficulties when attempting to change the
step-size. Example use to request
re-scaling
re_scale=1;
If no re-scaling is required then re_scale=0 can be used. It is not needed for the exchange algorithm (method=9), and generally
simulated annealing seems robust to scale issues.
constr=;
To force the optimization algorithm to arrange the sampling times per
group in a monotonic increasing manner such that sampling time n+1 >
sampling time n then set constr=1. This reduces the
search space and avoids samples being exchanged. This is not needed for
simulated annealing. Setting constr=0 avoids this
requirement.
dpm_constrain=
Further constraints to sampling times are available. Setting
|
dpm_constrain=0 |
no further
constraints. |
|
dpm_constrain=1 |
all sampling times are constrained to be the same for all
groups |
|
dpm_constrain=2 |
sampling times
are the same within a group for each occasion – user must define how to
divide samples between occasions (using first_occ). There
must be the same number of samples per occasion. |
|
dpm_constrain=3 |
constrains all
sampling times within a time period to be the same for all groups [only with method=8; user
must specify dpm_3_time]. When using this constraint
the user MUST ensure that there is the same number of sampling times for each
group within the time period specified by dpm_3_time. |
|
|
|
|
first_occ=; |
used in
conjunction with dpm_constrain=2. Defines how many samples constitute an
occasion in each group. For example if
there are 6 samples per group and samples 4, 5 and 6 are required to be
identical to samples 1, 2 & 3 then set first_occ=3.
This does not constrain samples 1, 2 and 3 from group 1 to equal samples 1, 2
and 3 from group 2. |
|
dpm_3_time=; |
defines the
time period that all samples from all groups should be the same. WARNING: you must ensure that ALL patients
have the same number of samples that will occur within the dpm_3_time
period. Do this by fixing boundaries
with UBA. WARNING: all samples within
this period should be exchangeable (as they are sorted) therefore DN must be
the same for every time - this cannot be tested for! |
lockout=;
To stop sampling times from occurring at the same time, or within a
defined time frame, set lockout to the minimal difference allowable. For example, if a pharmacodynamic response
takes 20 minutes to observe then two samples within a 20 minute period would
not be of use. In this case set lockout
to be, for example, lockout=0.5 (assuming time is in hours this gives 10
minutes extra).
OPTIMIZATION METHOD Back to this
example
There are currently 9 optimization methods available.
|
method=0 |
none |
|
method=1 |
Simplex (not available
with current version) |
|
method=2 |
non-adaptive random search |
|
method=3 |
simulated annealing
superseded by [method=8] |
|
method=4 |
grid (partitions=100 2 dimensions
only) |
|
method=5 |
compute expectation over
sampling times |
|
method=6 |
simulated annealing with
approximate design over treatment allocation |
|
method=7 |
simulated annealing with DN
optimized too [method=7] |
|
method=8 |
simulated annealing with
fixed times - this method supersedes [method=3] |
|
method=9 |
exchange algorithm |
|
|
|
|
|
For use with
non-adaptive random search |
|
Rwt= |
number of candidate sets
of parameter values |
|
|
|
|
|
For use with
method=5 |
|
n_exp= |
Provide the number of
evaluations of the information matrix from which to work out the expectation
of the determinant and standard errors, usually 1000. |
|
|
|
|
|
For use with
simulated annealing algorithm |
|
method_sa=0 |
normal temperature decline
(slowest – but more robust) |
|
method_sa=1 |
faster temperature decline
(geometric cooling) |
|
method_sa=2 |
fastest temperature
decline |
|
minv= |
Convergence criteria. This corresponds to the smallest stepsize which provides the number of significant digits
required. Usually set minv=0.0001 |
|
|
Initial value for
temperature. This value needs to be
adjusted. Generally |
|
|
|
|
|
For use with exchange
algorithm |
|
num_rand_arrangments |
This specifies the number
of times that the discrete variables are randomly rearranged between each set
of exchanges. This is used to help the algorithm out of local minima. Generally a value of 5 is sufficient. |
|
dpe_num |
This variable is used with
DPE is not provided. If discrete values of potential
sampling times are not provided by the user then the exchange algorithm will
provide its own discrete values. It does so by dividing the design space (UBA
– LBA) into a geometric series of discrete candidate
sampling times (that includes the upper (UBA) and lower bound (LBA) values).
The number of discrete sampling times is provided by dpe_num. Generally dpe_num=25
to 50 is fine. |
PRIOR INFORMATION Back to this
example
If prior information is available then this can be incorporated as a
fixed component of the information matrix which is specified by
MFD_prior=[];
This information matrix must have the same number of rows and columns as
the information matrix estimated in the current design. MFD_prior is then added to
the current information and the determinant of the sum of these matrices is
computed. The current sampling times are
therefore optimized for value added information.
The user has the opportunity to not use prior information by setting the
variable use_prior. If no prior information is available then use_prior should be set to zero. If prior information exists then the user can
decide whether to use it or not.
Example
POPT_INI_ROBUST BACK
TO TOP
There are 4 parts to this file
3) Setup
features for sampling windows
MODEL FEATURES FOR
SAMPLING WINDOWS
These are identical to those set out in the POPT_INI file
<GO TO POPT_INI Model features>.
DESIGN FEATURES FOR
WINDOWS Back to this
example
These are the same as those set out in the POPT_INI file
<GO TO POPT_INI Design features>
except for specification of the lower and upper bounds.
In this setting DPM should be set to the optimized sampling times from a
previous run. UBA should be the upper
boundary of the sampling window and LBA the lower boundary.
Lower and upper bounds must be chosen in order to represent the design
space from which to evaluate sampling windows.
The following suggested rules apply:
1)
all design points must have a specified upper and
lower bound
2)
no two bounds for sampling times should overlap for
samples in the same group (unless specified in rule 3).
3)
if 2 optimal sampling times have the same value then
it is recommended to set the upper and lower bounds to be the same
Some examples
For a 2 group design with 3 samples per group. Recall that the first 3 samples are for group
1 and the second 3 are for group 2. It
is OK to have sampling boundaries that overlap for samples that arise from
different groups. It is NOT OK to have
sampling boundaries that overlap for samples that arise from the same group
(except under suggested rule 3).
The following example satisfy all of the requirements
UBA = [0.5 2 8 3 6 12]
DPM = [0.1 1 5 2 4 8]
LBA = [0 0.5 2 0 3 6]
The following example does not satisfy all requirements and should not
be used
UBA = [1 2 8 3 6 12]
DPM = [0.1 1 5 2 4 8]
LBA = [0 0.5 2 0 3 6]
The following example shows how to deal with a replicated optimal sample
UBA = [2 2 8 3 6 12]
DPM = [1 1 5 2 4 8]
LBA = [0 0 2 0 3 6]
Specification of the values of LBA and UBA.
There are no right or wrong values to choose for LBA and UBA (as long as
the suggested rules are followed). If UBA and LBA are set to extreme values
(i.e. a long way apart) then it will take a long time to find the sampling
windows. If they are too close together
then your boundaries may be smaller than the windows and they may force
unnecessary constraints. It is usually
necessary to try a few runs first to find out how you should set your boundary
conditions. A suggested recipe:
1)
Try a set of LBA and UBA values, set the range to be
smaller for smaller values of the sampling time (sampling time windows are
generally proportional to the value of the sampling time)
2)
Run p_robust_1
3)
The 4th column in the output displays the fraction of
samples accepted. You should aim for a
value between 2 and 5%. If the value is
too small then constrict some of your sampling windows. If the value is too high then you have set
your sampling windows to be too constrictive.
SETUP FEATURES FOR
SAMPLING WINDOWS Back to this example
There is no current analytical solution to defining sampling windows for
nonlinear mixed effects models. In POPT
sampling windows are generated based on a 3 stage process.
1)
Generate random samples from a uniform distribution
with limits of UBA & LBA and accept those samples that have efficiencies
greater than a predefined percent.
2)
Re-sample from those samples that are accepted in
stage 1 to generate many different combinations of sampling times. Accept those that have an efficiency greater
than a pre-defined level. Continue this
to inflate the number of samples accepted to some predefined level.
3)
Re-sample a pre-defined number of times from those
samples accepted in stage 2. Accept any
that have an efficiency greater than some predefined value. Discard those that don’t. Repeat stage 3 until the fraction of those
that are discarded is lower than 10%.
Finally a joint sampling window is constructed as the percentile range
(defined by the user) of the sampling distributions from stage 3.
n_samples=;
This variable indicates how many samples that are generated at random
have to be accepted before stage 1 is completed. Typically n_samples = 1000;
win_eff_k =;
The minimum efficiency of any set of samples randomly generated in stage
1. Typically, win_eff_k = 0.90;
max_it=;
The maximum number of iterations.
Typically, max_it = 1e+6;
inflation =;
The number of re-sampled sets of sampling times to be generated in stage
2 from the joint samples generated from stage 1. Typically, inflation = 5000;
win_eff_r=;
The minimum relative efficiency of resampled sets of sampling times
generated in stage 3 at which level are accepted.
win_eff_r = 0.95;
lb_win=;
ub_win=;
P_ROBUST_1 reports the lower bound (lb_win) and upper bound
(ub_win) of the sampling windows as the percentiles of
the distribution for each sampling window.
The percentiles are typically, lb_win=5; & ub_win=95; indicating that P_ROBUST_1 will report the 5th
and 95th percentile of the distribution of the accepted samples from
stage 3.
PRIOR INFORMATION
FOR SAMPLING WINDOWS
These are identical to those set out in the POPT_INI file
<GO TO POPT_INI Prior Information>.
Interpreting Output BACK
TO TOP
When POPT is running with either method=8 or method=9 it will produce a
screen window like the following example.
This screen can be turned off by setting output_on = 0;.
For method=8 (simulated annealing)
- The upper
left panel represents the log of the determinant of the information
matrix.
- The upper middle
panel is the current sampling times (the y-axis is meaningless) presented
across the x-axis
- The upper
right panel is the current value of the temperature
- The lower
left panel is the step size (the amount by which the average sampling time
is allowed to deviate during any iteration of the algorithm. When this reaches 10-4 the algorithm
will have converged.
- The lower
middle panel is the fraction of moves that have been accepted away from
the optimum. This should generally
lie between 0.1 and 0.4
- The lower
right panel is the average period of time per iteration.
Output will also be numerically echoed to the screen as text. The column headings are the current values
of:
Column 1: Iteration number
Column 2: Temperature
Column 3: Maximum stepsize (stepsize for worst
estimated sampling time)
Column 4: Determinant
Column 5: Number of false
acceptances
For method=9 (exchange algorithm)
- The upper
left panel represents the log of the determinant of the information
matrix.
- The upper
middle panel is the current sampling times (the y-axis is meaningless)
presented across the x-axis
- The upper
right panel is the current number of successful exchanges made
- The lower
left panel is not used.
- The lower
middle panel is the number of random rearrangements of the discrete sample
space (usually the algorithm terminates after 5 rearrangements)
- The lower
right panel is the average period of time per iteration.
Output will also be numerically echoed to the screen as text. The column headings are the current values
of:
Column 1: Iteration Number
Column 2: Number of
re-randomizations
Column 3: Number of
successful exchanges
Column 4: Number of echoes to
the screen
Column 5: Determinant
Output file for
POPT Back to this example
Output files are created for all successful POPT runs. POPT creates an output file labelled methodx.res, where x is the method number used. For multiple models or multiple response
models POPT creates the file Methodx_comp.res. A break down of the output for a single model
single response is provided. The output
for a multiple model/multiple response is the same as the output for a single
model/response design except that the information matrix and standard errors are
computed for each model/response.
|
Variable |
Interpretation |
|
dose |
Dose per group |
|
NT |
Number of sampling times
per group |
|
NS |
Number of subjects per
group |
|
DN |
Dose number that the
sampling was taken from |
|
method |
Optimization method |
|
initial_starting_times |
Your initial estimates (DPM) |
|
sampling_times |
Optimized sampling times –
if you have chosen method=0 then these will equal your initial_sampling_times |
|
information_matrix |
The information matrix –
you can use this to get your prior information for subsequent runs. |
|
fixed_effect_parameters |
Values of the fixed
effects parameters (BETAI) |
|
POP_PAR |
The names of the
parameters |
|
FXFI |
Which parameters were fixed
– the numbers represent the position of the parameter vector (BETAI) |
|
standard_errors_fe |
Standard errors of the
fixed effects parameters |
|
standard_errors_percent_fe |
Standard errors of the
fixed effects parameters expressed as a percent |
|
BSV_parameters |
Values of the between
subject variance (OI) |
|
FXI |
Which parameters were
fixed – the numbers represent the position of the parameter vector (OI) |
|
standard_errors_bsv |
Standard errors of the BSV
parameters |
|
standard_errors_percent_bsv |
Standard errors of the BSV
parameters expressed as a percent |
|
BSV_model |
The variance model used
for BSV – 1 = exponential, 0 = additive |
|
residual_standard_deviation |
Values of the residual
standard deviation |
|
FXR |
Which values of the
residual standard deviation were fixed |
|
standard_errors_res |
Standard errors of the
residual standard deviation |
|
standard_errors_percent_res |
Standard errors of the
residual standard deviation expressed as a percent |
|
determinant |
The determinant of the
information matrix |
|
criterion |
The determinant to the
power of the inverse of the number of parameters (fixed effects + BSV + RUV)
in the model |
|
iterations |
The number of iterations used
by the optimization routine |
|
eigen_values |
The eigen
values of the information matrix. The
ratio of the highest to lowest numbers is often quoted as a measure of the |
|
time_taken_seconds |
Time taken to run the
problem and produce the output. If you
are using method=0 then echoing the output to the screen can take as long as
evaluating the information matrix. In
this case the screen output will include two times the first being the time
to evaluate the information matrix and the second the overall time. |
|
az |
Page separator. |
SCREEN OUTPUT FOR
P_ROBUST_1 Back
to this example
When P_ROBUST_1 is running it will produce a text screen output. The screen output changes depending on which
stage of the sampling design is being evaluated. For a description of the stages see section
<SET UP FEATURES FOR SAMPLING WINDOWS>.
Stage 1. Samples are generated randomly from ~U(LBA, UBA), until the total number of
accepted samples is 1000 (column 3). A
reasonable acceptance rate (column 4) is 0.02 to 0.1.
|
Stage |
Iteration number |
Number of accepted sets of
candidate samples |
Fraction of sets of
samples that are accepted |
For Stage 2 Samples are generated by resampling from the
stage 1 samples, until the number of accepted sets of candidate samples = 5000
(column 3).
|
Stage |
Iteration number |
Number of accepted sets of
candidate samples |
For Stage 3 Samples are generated by resampling from the
stage 2 samples until the proportion of rejected samples is < 0.1.
|
Stage |
Iteration number |
|
Fraction of sets of
samples that are accepted |
The stage number increments as stage 3 is repeated until the fraction
reduces to less than 0.1.
OUTPUT FILE FOR
P_ROBUST_1 Back to this example
P2NM = POPT too NONMEM
This program creates a NONMEM data file template suitable for
simulation-estimation runs in NONMEM.
This program calls POPT_INI.m and uses the
design provided in this file to create a NONMEM data file that has the
appropriate number of patients, doses, sampling times etc. The file is saved to NM_dat.csv. If the file NM_dat.csv
already exists then the file is overwritten.
The csv file is written from row 2, where the
first row is left blank in order for the user to manually input a header row if
desired. P2NM uses “0” as a “.”
Placeholder. The columns created are ID,
TIME, AMT, DV, MDV.
The NM data file is automatically started at ID=1. If a different starting ID value is required
then edit line 25
init_id_value=;
The value of init_id_value is taken as the first value of ID in the data
file and additional patients are given consecutive values. Example syntax init_id_value=2000 will
start the data file numbering at ID=2000.
Optimization algorithms BACK TO TOP
Exchange
The exchange algorithm implemented
in POPT is a simple k=1 algorithm, where on each iteration a single sampling
time is exchanged with one from a list of possible sampling times that provides
the same or a better value of the determinant.
The list of possible sampling times is either provided by the user or if
left blank is generated automatically by POPT.
If generated automatically POPT divides the dose interval into 50 (or
some other value as determined by the user) geometrically spaced samples and searches
over this discrete space. The search for
exchanges occurs by starting from the lower bound of the interval and
progressing through the list until an exchange occurs. If an exchange occurs then the search
re-starts at the lower bound of the search space again. When the search progresses through the list
one complete time without an exchange the algorithm randomizes the search space
and tries again. This process is repeated
until 5 re-randomizations occur (or as defined by the user). On some occasions this can result in renewed
exchanges.
Simulated Annealing
In this implementation of simulated
annealing the convergence criterion refers to the accuracy of the design points
(this is set by the user). A value of 0.01 means that all designs points are estimated
to a value of less than 0.01 significant digits of time units. Note that other implementations may set the
convergence to be a function of the current temperature. We have found this latter method to be
unreliable for searching over the surface of the determinant of a population
Fisher information matrix. We have found
the starting temperature of 1000 to be suitable for most PK and PKPD designs to
date. The initial temperature selected
should be one that provides few false acceptances (last column of interim
output). The choice of fast, faster and
slow is dependent on how sure you want to be that you have found a global
minimum. For most purposes fast or
faster are fine. The slower the process
the better the statistical quality of the convergence – but it can become very
slow. Don’t forget it is the utility of
the design that is located rather than finding the design with the absolute
highest possible determinant. Unlike
estimation methods, lack of convergence (or if the algorithm is terminated
during its run) does NOT indicate that the estimates of the sampling times are
in anyway biased.
In principle simulated annealing
works by search across the design space for regions of high probability. It determines these regions by improvements
in the objective function (the determinant in this case). It explores the design region by moving away
from the region of highest probability and accepting “false” candidates (i.e.
ones that are not as good as the current best region) at their probability
level of occurring using a Metropolis step.
The acceptance probability of these false samples becomes more difficult
as the temperature declines.
Copyright
Statement BACK TO TOP
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2 of
the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more
details.
POPT (C) copyright Stephen Duffull ver 3.0;
Oct, 5th 2003
PFIM_M (C) copyright Stephen Duffull ver 1.2(i) May; 28th 2003
TO MAINTAIN DESIRED FUNCTIONALITY, DO NOT EDIT THIS FILE
General Design
Atkinson AV. Donev AN. Ooptimum
Experimental Designs.
Box GEP. Lucas HL. Design of experiments in
non-linear situations. Biometrika 1959;46:77-90/
Cramér H. Mathematical
Methods of Statistics.
Draper NR. Hunter WG. Design of
experiments for parameter estimation in multiple response situations. Biometrika 1966;53:525-533
Population Fisher Information Matrix
Mentre F. Mallet A. Baccar
D. Optimal design in
random effects regression models. Biometrika. 1997;84:429-442
Retout S, Duffull S, Mentre F. Development and
implementation of the population Fisher information matrix for the evaluation
of population pharmacokinetic designs.
Comput Meth Prog Biomed 2001;65:141-151.
S. Retout and F. Mentre.
Further developments
of the Fisher information matrix in nonlinear mixed effects models with
evaluation in population pharmacokinetics. J. Biopharm. Stat. 13:209-227 (2003).
POPT
Hennig, S., Waterhouse, T.H., Bell, S.C.,
Duffull S. Waterhouse T, Eccleston J.
Some considerations on the design of population
pharmacokinetic studies. J Pharmacokinet Pharmacodyn 2005;32:441-457
Waterhouse TH, Redmann S, Duffull SB,
Eccleston JA. Optimal
design for model discrimination and parameter estimation for itraconazole
population pharmacokinetics in cystic fibrosis patients. J Pharmacokinet Pharmacodyn 2005;32:521-545
Green B, Duffull SB. Prospective
evaluation of a D-optimal designed population pharmacokinetic study. J Pharmacokinet Pharmacodyn 2003; 30:145-161
Duffull
SB, Retout S, Mentré
F. The use of simulated annealing for finding optimal
population designs. Comp. Methods Programs Biomed. 2002;69:25-35
Duffull SB, Mentré F, Aarons L. Optimal design of a
population pharmacodynamic experiment for ivabradine. Pharm Res 2001;18:83-9
Simulated Annealing
Duffull
SB. Retout S. and Mentre F. The use of simulated annealing
for finding optimal population designs. Comput. Meth. Prog. Biomed., 2002; 69:25-35
Goffe WL. Ferrier GD. Rogers J. Global optimisation
of statistical functions with simulated annealing. J Econom
1994:60:65-99