Description Usage Arguments Details Value Author(s) References See Also Examples
Generates an approximate utility function for generalised linear models and nonlinear regression models.
1 2 3 4 5 6  utilityglm(formula, family, prior,
criterion = c("D", "A", "E", "SIG", "NSEL", "SIGNorm", "NSELNorm"),
method = c("quadrature", "MC"), nrq)
utilitynlm(formula, prior, desvars, criterion = c("D", "A", "E", "SIG", "NSEL"),
method = c("quadrature", "MC"), nrq)

formula 
An argument providing a symbolic description of the model. For For 
family 
For 
prior 
An argument specifying the prior distribution. For For 
desvars 
For 
criterion 
An optional character argument specifying the utility function. There are currently seven utility functions implemented as follows:
The default value is 
method 
An optional character argument specifying the method of approximating the expected utility function. Current choices are 
nrq 
For 
Two utility functions are implemented.
Shannon information gain (SIG)
The utility function is
U^SIG(d) = f(θy,d)  f(θ),
where f(θy,d) and f(θ) denote the posterior and prior densities of the parameters θ, respectively.
Negative squared error loss (NSEL)
The utility function is
u^NSEL(d) =  (θ  E(θ y,d))^T(θ  E(θ y,d)),
where E(θ  y,d) denotes the posterior mean of θ.
In both cases the utility function is not available in closed form due to the analytical intractability of either the posterior distribution (for SIG) or the posterior mean (for NSEL). The acebayes
package implements two approximations to both utility functions. If criterion = "SIG"
or criterion = "NSEL"
then samplingbased Monte Carlo or importance sampling approximations will be employed. This was the original approach used by Overstall & Woods (2017). If criterion = "SIGNorm"
or criterion = "NSELNorm"
then approximations based on approximate normality of the posterior (Overstall et al., 2017) will be used.
The normal approximation to the posterior can be taken further leading to the approximation by some scalar function of the Fisher information matrix, I (θ;d), which only depends on θ (Chaloner & Verdinelli, 1995). In the case of SIG, the approximate utility is given by
u^D(d) = log  I(θ;d),
and the resulting design is typically called pseudoBayesian Doptimal. For NSEL, the approximate utility is given by
u^A(d) =  tr (I(θ;d)^(1))
with the resulting design termed pseudoBayesian Aoptimal. These designs are often used under the frequentist approach to optimal experimental design and so to complete the usual set, the following utility for finding a pseudoBayesian Eoptimal design is also implemented:
U^E(d) = min(e(I(θ;d))),
where e() denotes the function that calculates the eigenvalues of its argument.
The expected utilities can be approximated using Monte Carlo methods (method = "MC"
for all criteria) or using a deterministic quadrature method (method = "quadrature"
, implemented for the D, A and E criteria). The former approach approximates the expected utility via sampling from the prior. The latter approach uses a radialspherical integration rule (Monahan and Genz, 1997) and B[1]
specifies the number, n_r, of radial abscissas and B[2]
specifies the number, n_q, of random rotations. Larger values of n_r will produce more accurate, but also more computationally expensive, approximations. See Gotwalt et al. (2009) for further details.
For utilityglm
, note that the utility functions for SIG and NSEL are currently only implemented for logistic regression, i.e. family = binomial
, or Poisson regression, i.e. family = poisson(link = "log")
, whereas the utility functions for pseudoBayesian designs are implemented for generic GLM families.
For more details on the ACE algorithm, see Overstall & Woods (2017).
The function will return a list with the following components:
utility 
The utility function resulting from the choice of arguments. 
Antony M. Overstall A.M.Overstall@soton.ac.uk, David C. Woods, Maria Adamou & Damianos Michaelides
Chaloner, K. & Verdinelli, I. (1995). Bayesian experimental design: a review. Statistical Science, 10, 273304.
Gotwalt, C. M., Jones, B. A. & Steinberg, D. M. (2009). Fast computation of designs robust to parameter uncertainty for nonlinear settings. Technometrics, 51, 8895.
Monahan, J. and Genz, A. (1997). Sphericalradial integration rules for Bayesian computation,” Journal of the American Statistical Association, 92, 664674.
Overstall, A.M. & Woods, D.C. (2017). Bayesian design of experiments using approximate coordinate exchange. Technometrics, 59, 458470.
aceglm
, acenlm
, paceglm
, pacenlm
.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52  ## 1. This example uses utilityglm to generate the pseudoBayesian Doptimality
## approximate expected utility function using a Monte Carlo approximation.
low<c(3, 4, 5, 6, 2.5)
upp<c(3, 10, 11, 0, 3.5)
## Lower and upper limits of the uniform prior distributions.
prior<function(B){
t(t(6*matrix(runif(n=5*B),ncol=5))+low)}
## Create a function which specifies the prior. This function will return a
## B by 5 matrix where each row gives a value generated from the prior
## distribution for the model parameters.
ex < utilityglm(formula = ~x1+x2+x3+x4, family = binomial, prior = prior, method = "MC")
set.seed(1)
## Set seed for reproducibility.
n<6
## Specify the sample size (number of runs).
start.d<matrix( 2 * randomLHS(n = n,k = 4)  1,nrow = n,ncol = 4,
dimnames = list(as.character(1:n),c("x1", "x2", "x3", "x4")))
## Generate an initial design of appropriate dimension. The initial design is a
## Latin hypercube sample.
ex$utility(d = start.d, B = 10)
## Evaluate resulting approximate utility. Should get:
#[1] 13.98143 17.07772 19.88988 22.40720 15.27411 15.02717 16.17253 18.66600 13.75118
#[10] 21.83820
## 2. This example uses utilitynlm to generate the psuedoBayesian Aoptimality expected utility
## function using a quadrature approximation
low<c(0.01884, 0.298, 21.8)
upp<c(0.09884, 8.298, 21.8)
## Lower and upper limits of the uniform prior distributions. Note that the prior
## for the third element is a point mass.
prior2 < list(support = cbind(rbind(low, upp)))
colnames(prior2$support) < c("a", "b", "c")
## Specify a uniform prior with ranges given by low and upp
ex2 < utilitynlm(formula = ~ c * (exp(  a * t)  exp(  b *t)), prior = prior2,
desvars = "t")
n < 6
start.d < matrix(24 * randomLHS(n = n, k = 1), nrow = n)
colnames(start.d) < "t"
ex2$utility(d = start.d)
## 13.17817

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.