overstallwoods {acebayes} | R Documentation |
Functions implementing the examples of Overstall & Woods (2017).
Description
This suite of functions implement the examples in Overstall & Woods (2017).
Usage
######## Compartmental model #################################
utilcomp18bad(d, B)
optdescomp18bad(type = "ACE")
utilcomp15bad(d, B)
optdescomp15bad()
utilcomp15sig(d, B)
optdescomp15sig()
utilcomp15sigDRS(d, B)
optdescomp15sigDRS()
######## Logistic regression model ###########################
utillrbad(d, B)
optdeslrbad(n, type = "ACE")
utillrsig(d, B)
inideslrsig(n, rep)
optdeslrsig(n)
utilhlrbad(d, B)
optdeshlrbad(n)
utilhlrsig(d, B)
inideshlrsig(n, rep)
optdeshlrsig(n)
utillrbaa(d, B)
optdeslrbaa(n)
utillrnsel(d, B)
inideslrnsel(n, rep)
optdeslrnsel(n)
optdeshlrbaa(n)
utilhlrbaa(d, B)
utilhlrnsel(d, B)
inideshlrnsel(n, rep)
optdeshlrnsel(n)
######## Beetle mortality experiment #########################
utilbeetle(d, B)
optdesbeetle(n)
######## Linear model ########################################
utillinmod(d, B)
optdeslinmod(n, type = "ACE")
##############################################################
Arguments
d |
An |
n |
The number of runs ins the experiment. |
rep |
A scalar integer in the set |
B |
A scalar integer specifying the Monte Carlo sample size. |
type |
An optional character argument specifying which design to return. For For |
Details
This section provides details on the examples considered and the functions implemented in acebayes
.
Compartmental model
Compartmental models are used in Pharmokinetics to study how materials
flow through an organism. A drug is administered to an individual or animal and then the
amount present at a certain body location is measured at a set of n
pre-determined sampling
times (in hours). There is one design variable: sampling time. Therefore the design matrix d
is
an n
by 1 matrix with elements controlling the n
sampling times, i.e. the number of factors is
k=1
.
In Overstall & Woods (2017), two different compartmental model examples are considered. The first (in the Supplementary Material)
comes from Atkinson et al (1993) and Gotwalt et al (2009) where there are n = 18
sampling times and interest
lies in finding a Bayesian D-optimal design. The functions whose name includes "comp18"
refer to this example.
The second example (in Section 3.2) comes from Ryan et al (2014), where there are n = 15
sampling times and
the ultimate interest lies in finding an optimal design under the Shannon information gain utility. Also considered is the
Bayesian D-optimal design. The functions whose name includes "comp15"
refer to this example. Note that Ryan et al (2014) used a dimension reduction scheme (DRS) to find optimal designs. The function whose name is suffixed by "DRS"
refer to this situation.
Logistic regression model
In Section 3.3 of Overstall & Woods (2017), a first-order logistic regression model in k = 4
factors and n
runs is considered. Woods et al (2006) and Gotwalt et al (2009) considered generating Bayesian D-optimal designs for n = 16
and n = 48
. Overstall & Woods (2017) extended this example by considering Bayesian A-optimal, Shannon information gain (SIG) and negative squared error loss (NSEL) utility functions, a range of number of runs from 6 to 48, and "random effects" to form a hierarchical logistic regression model.
Beetle mortality experiment
Overstall & Woods (2017, Section 3.4) considers generating an optimal design for a follow-up experiment. The original design and data (Bliss, 1935) involves administering different doses of poison to N = 8
groups of beetles. The number of beetles that die in each group are recorded. Six different models are considered formed from the combination of three link functions and two linear predictors (following the analysis of O'Hagan & Forster, 2004). Interest lies in the quantity known as lethal dose 50 (LD50) which is the dose required to kill 50% of the beetles and is a function of the model parameters for a given model. Consider finding an optimal design for estimating LD50 under the negative squared error loss (NSEL) function for n
new doses of poison (i.e. k = 1
factor). The prior distribution is equivalent to the posterior distribution arising from the original data and includes model uncertainty.
Linear model
In the supplementary material, Overstall & Woods (2017) considers finding D-optimal designs for a second-order (i.e. k = 2
) response surface in n=6,7,8,9
runs. Note that the D-optimal design is equivalent to the optimal design under Shannon information gain and a non-informative prior distribution.
The expected utility function in this case is available in closed form, i.e. it does not require approximation. Box & Draper (1971) found optimal designs analytically for the number of runs considered here. Overstall & Woods (2017) use this example to demonstrate the efficacy of the ACE algorithm.
Value
Compartmental model
For the example in the Supplementary Material;
-
The function
utilcomp18bad
will return a vector of lengthB
where each element is the value of the Bayesian D-optimal utility function (i.e. the log determinant of the Fisher information matrix) evaluated at a sample of sizeB
generated from the prior distribution of the model parameters. -
The function
optdescomp18bad
will return an 18 by 1 matrix giving the optimal design (specified by the argumenttype
). The elements will be scaled to be in the interval [-1, 1], i.e. a -1 corresponds to a sampling times of 0 hours, and 1 corresponds to a time of 24 hours.
For the example in Section 3.2;
-
The function
utilcomp15bad
will return a vector of lengthB
where each element is the value of the Bayesian D-optimal utility function (i.e. the log determinant of the Fisher information matrix) evaluated at a sample of sizeB
generated from the prior distribution of the model parameters. -
The function
optdescomp15bad
will return an 15 by 1 matrix giving the optimal design (found using ACE) under Bayesian D-optimality. The elements will be scaled to be in the interval [-1, 1], i.e. a -1 corresponds to a sampling times of 0 hours, and 1 corresponds to a time of 24 hours. -
The function
utilcomp15sig
will return a vector of lengthB
where each element is the value of the SIG utility function evaluated at a sample of sizeB
generated from the joint distribution of model parameters and unobserved responses. -
The function
optdescomp15sig
will return an 18 by 1 matrix giving the optimal design (found using ACE) under the SIG utility. The elements will be scaled to be in the interval [-1, 1], i.e. a -1 corresponds to a sampling times of 0 hours, and 1 corresponds to a time of 24 hours. -
The function
utilcomp15sigDRS
will return a vector of lengthB
where each element is the value of the SIG utility function (where a DRS has been used) evaluated at a sample of sizeB
generated from the joint distribution of model parameters and unobserved responses. Here the Beta DRS (see Overstall & Woods, 2017) has been used sod
should be a 2 by 1 matrix containing the positive beta parameters. -
The function
optdescomp15sigDRS
will return a 2 by 1 matrix giving the optimal design (found using ACE) under the SIG utility, where a DRS has been used. The elements correspond to the parameters of a beta distribution.
Logistic regression model
A function whose name includes "lr"
refers to standard logistic regression, whereas "hlr"
refers to hierarchical logistic regression. Under standard logistic regression the possible values for the argument n
can be any even integer between 6 and 48. For hierarchical logistic regression, n
can be any integer divisible by 6 between 12 and 48. The function name also indicates the utility function:
"bad"
Bayesian D-optimal"baa"
Bayesian A-optimal"sig"
Shannon information gain"nsel"
Negative squared error loss
The functions prefixed by "util"
will return a vector of length B
where each element is the utility function evaluated at a sample generated from the prior distribution of model parameters (for Bayesian D- and A-optimality) or the joint distribution of model parameters and unobserved responses (for SIG and NSEL).
The functions prefixed by "optdes"
will return an n
by k = 4
matrix giving the optimal design found by ACE. The designs given by this function are those reported on in Overstall & Woods (2017). The function optdeslrbad
will also return designs (for n = 16, 48
) found by Woods et al (2006) and Gotwalt et al (2009) by specifying the type
argument appropriately.
The functions prefixed by "inides"
will return an n
by k = 4
matrix giving an initial design for ACE to find the optimal designs under the SIG and NSEL utility functions. These are 20 designs found using ACE under approximations to the Bayesian A- and D-optimal utility functions, respectively. The argument rep
specifies which of these 20 designs to use.
Beetle mortality experiment
The function utilbeetle
will return a vector of length B
where each element is the value of the utility function for a sample generated from the joint distribution of the model parameters, model and unobserved responses.
The function optdesbeetle
will return an n
by 1 matrix giving the optimal design under the NSEL utility function (found using ACE) for estimating the LD50. The elements will be scaled to be in the interval [-1, 1], where -1 corresponds to a dose of 1.6907, 0 to a dose of 1.7873 and 1 to a dose of 1.8839. The designs given by this function are those reported in Overstall & Woods (2017) for n
= 1, ..., 10.
Linear model
The function utillinmod
will return a vector of length B
where each element is a realisation of a stochastic approximation to the expected utility.
The function optdeslinmod
will return an n
by 2 matrix giving the D-optimal design (specified by the argument type
). If type = "ACE"
, the designs returned by this function are those found using the ACE algorithm and are reported in the Supplementary Material of Overstall & Woods (2017), and if type = "BoxDraper"
, the designs returned are the exact D-optimal designs.
Author(s)
Antony M. Overstall A.M.Overstall@soton.ac.uk, David C. Woods, Maria Adamou & Damianos Michaelides
References
Atkinson, A., Chaloner, K., Herzberg, A., & Juritz, J. (1993). Experimental Designs for Properties of a Compartmental Model. Biometrics, 49, 325-337.
Bliss, C. (1935). The calculation of the dosage-mortality curve. Annals of Applied Biology, 22, 134-167.
Box, M. & Draper, N. (1971). Factorial designs, the |F^T F|
criterion and some related matters.
Techometrics, 13, 731-742.
Gotwalt, C., Jones, B. & Steinberg, D. (2009). Fast Computation of Designs Robust to Parameter Uncertainty for Nonlinear Settings. Technometrics, 51, 88-95.
O'Hagan, A, & Forster, J.J. (2004). Kendall's Advanced Theory of Statistics, Volume 2B: Bayesian Inference. 2nd edition. John Wiley & Sons.
Overstall, A.M. & Woods, D.C. (2017). Bayesian design of experiments using approximate coordinate exchange. Technometrics, 59, 458-470.
Ryan, E., Drovandi, C., Thompson, M., Pettitt, A. (2014). Towards Bayesian experimental design for nonlinear models that require a large number of sampling times. Computational Statistics and Data Analysis, 70, 45-60.
Woods, D.C., Lewis, S., Eccleston, J., Russell, K. (2006). Designs for Generalized Linear Models With Several Variables and Model Uncertainty. Technometrics, 48, 284-292.
See Also
Examples
######## Compartmental model #################################
set.seed(1)
## Set seed for reproducibility
d<-optimumLHS(n = 18, k = 1) * 2 - 1
## Generate an 18-run design.
u<-utilcomp18bad(d = d, B = 20000)
## Calculate the D-optimal utility function for a
## sample of size 20000.
u[1:5]
## Look at first 5 elements.
#[1] 14.283473 10.525390 4.126233 7.061498 12.793245
d0<-optdescomp18bad( )
u0<-utilcomp18bad(d = d0, B = 20000)
## Optimal design found by ACE and calculate the D-optimal
## utility function for a sample of size 20000.
u0[1:5]
## Look at first 5 elements.
#[1] 15.04721 15.37141 16.84287 14.06750 14.01523
mean(u)
mean(u0)
## Calculate expected Bayesian D-optimal utility.
d<-matrix(runif(2), ncol = 1)
## Generate two beta parameters.
u<-utilcomp15sigDRS(d = d, B = 5)
u
## Calculate the SIG utility function for a
## sample of size 5.
#[1] 17.652044 4.878998 19.919559 22.017760 5.600473
######## Logistic regression model ###########################
set.seed(1)
## Set seed for reproducibility
d<-optimumLHS(n = 16,k = 4) * 2 - 1
## Generate an 16-run design.
u<-utillrbad(d = d, B = 20000)
## Calculate the D-optimal utility function for a
## sample of size 20000.
u[1:5]
## Look at first 5 elements.
#[1] -11.630683 -5.748912 -9.554413 -10.150132 -7.940938
d0<-optdeslrbad(16)
u0<-utillrbad(d = d0, B = 20000)
## Optimal design found by ACE and calculate the D-optimal
## utility function for a sample of size 20000.
u0[1:5]
## Look at first 5 elements.
#[1] -4.644116 -2.411431 -4.999891 -2.906558 -2.282687
mean(u)
mean(u0)
## Calculate expected Bayesian D-optimal utility.
#[1] -9.38253
#[1] -2.992012
######## Beetle mortality experiment #########################
set.seed(1)
## Set seed for reproducibility
d<-optimumLHS(n = 10, k = 1)*2-1
## Generate a design of 10 doses with elements in [-1, 1]
utilbeetle(d = d, B = 5)
## Calculate the utility function for a sample of size 5.
#-4.720491e-06 -1.198955e-05 -1.681380e-05 -3.123498e-06 -1.412722e-05
d0<-optdesbeetle(10)
d0
## Print out optimal design from Overstall & Woods (2017) for 10 doses
0.5*( d0 + 1)*( 1.8839 - 1.6907 ) + 1.6907
## On original scale.
# [,1]
# [1,] 1.769957
# [2,] 1.769520
# [3,] 1.768641
# [4,] 1.777851
# [5,] 1.768641
# [6,] 1.769520
# [7,] 1.777851
# [8,] 1.765997
# [9,] 1.768641
#[10,] 1.768641
######## Linear model ########################################
set.seed(1)
## Set seed for reproducibility
d<-cbind(rep(c( -1, 0, 1), each = 3), rep(c( -1, 0, 1), 3))
## Create a 9-run design which is the true D-optimal design
utillinmod(d = d, B = 5)
## Calculate the approximation to the true expected D-optimal utility
## function for a sample of size 5.
#[1] 7.926878 8.736976 7.717704 10.148613 8.882840
d0<-optdeslinmod(9)
## Optimal D-optimal design found using ACE
X<-cbind(1, d, d^2, d[,1] * d[,2])
X0<-cbind(1, d0, d0^2, d0[,1] * d0[,2])
## Calculate model matrices under both designs
detX<-determinant( t(X) %*% X)$modulus[1]
detX0<-determinant( t(X0) %*% X0)$modulus[1]
## Calculate true expected D-optimal utility function for both designs
100 * exp( 0.2 * ( detX0 - detX ))
## Calculate D-efficiency of ACE design.
# 99.93107