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
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
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
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:
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