bayes {ICAOD}

Bayesian D-Optimal Designs

Description

Finds (pseudo) Bayesian D-optimal designs for linear and nonlinear models. It should be used when the user assumes a (truncated) prior distribution for the unknown model parameters. If you have a discrete prior, please use the function robust.

Usage

bayes(
  formula,
  predvars,
  parvars,
  family = gaussian(),
  prior,
  lx,
  ux,
  iter,
  k,
  fimfunc = NULL,
  ICA.control = list(),
  sens.control = list(),
  crt.bayes.control = list(),
  sens.bayes.control = list(),
  initial = NULL,
  npar = NULL,
  plot_3d = c("lattice", "rgl"),
  x = NULL,
  crtfunc = NULL,
  sensfunc = NULL
)

Arguments

formula

A linear or nonlinear model formula. A symbolic description of the model consists of predictors and the unknown model parameters. Will be coerced to a formula if necessary.

predvars

A vector of characters. Denotes the predictors in the formula.

parvars

A vector of characters. Denotes the unknown parameters in the formula.

family

A description of the response distribution and the link function to be used in the model. This can be a family function, a call to a family function or a character string naming the family. Every family function has a link argument allowing to specify the link function to be applied on the response variable. If not specified, default links are used. For details see family. By default, a linear gaussian model gaussian() is applied.

prior

An object of class cprior. User can also use one of the functions uniform, normal, skewnormal or student to create the prior. See 'Details' of bayes.

lx

Vector of lower bounds for the predictors. Should be in the same order as predvars.

ux

Vector of upper bounds for the predictors. Should be in the same order as predvars.

iter

Maximum number of iterations.

k

Number of design points. Must be at least equal to the number of model parameters to avoid singularity of the FIM.

fimfunc

A function. Returns the FIM as a matrix. Required when formula is missing. See 'Details' of minimax.

ICA.control

ICA control parameters. For details, see ICA.control.

sens.control

Control Parameters for Calculating the ELB. For details, see sens.control.

crt.bayes.control

A list. Control parameters to approximate the integral in the Bayesian criterion at a given design over the parameter space. For details, see crt.bayes.control.

sens.bayes.control

A list. Control parameters to verify the general equivalence theorem. For details, see sens.bayes.control.

initial

A matrix of the initial design points and weights that will be inserted into the initial solutions (countries) of the algorithm. Every row is a design, i.e. a concatenation of x and w. Will be coerced to a matrix if necessary. See 'Details' of minimax.

npar

Number of model parameters. Used when fimfunc is given instead of formula to specify the number of model parameters. If not specified correctly, the sensitivity (derivative) plot may be shifted below the y-axis. When NULL (default), it will be set to length(parvars) or prior$npar when missing(formula).

plot_3d

Which package should be used to plot the sensitivity (derivative) function for two-dimensional design space. Defaults to "lattice".

x

A vector of candidate design (support) points. When is not set to NULL (default), the algorithm only finds the optimal weights for the candidate points in x. Should be set when the user has a finite number of candidate design points and the purpose is to find the optimal weight for each of them (when zero, they will be excluded from the design). For design points with more than one dimension, see 'Details' of sensminimax.

crtfunc

(Optional) a function that specifies an arbitrary criterion. It must have especial arguments and output. See 'Details' of bayes.

sensfunc

(Optional) a function that specifies the sensitivity function for crtfunc. See 'Details' of bayes.

Details

Let \Xi be the space of all approximate designs with k design points (support points) at x_1, x_2, ..., x_k from design space \chi with corresponding weights w_1, . . . ,w_k. Let M(\xi, \theta) be the Fisher information matrix (FIM) of a k-point design \xi and \pi(\theta) is a user-given prior distribution for the vector of unknown parameters \theta. A Bayesian D-optimal design \xi^* minimizes over \Xi

\int_{\theta \in \Theta} -\log|M(\xi, \theta)| \pi(\theta) d\theta.

An object of class cprior is a list with the following components:

• fn: Prior distribution as an R function with argument param that is the vector of the unknown parameters. See below.

• npar: Number of unknown parameters and is equal to the length of param.

• lower: Argument lower. It has the same length as param.

• upper: Argument upper. It has the same length as param.

A cprior object will be passed to the argument prior of the function bayes. The argument param in fn has the same order as the argument parvars when the model is specified by a formula. Otherwise, it is the same as the argument param in the function fimfunc.
The user can use the implemented priors that are uniform, normal, skewnormal and student to create a cprior object.

To verify the equivalence theorem of the output design, use plot function or change the value of the checkfreq in the argument ICA.control.

To increase the speed of the algorithm, change the value of the tuning parameters tol and maxEval via the argument crt.bayes.control when crt.bayes.control$method = "cubature". Similarly, this applies when crt.bayes.control$method = "quadrature". In general, if the CPU time matters, the user should find an appropriate speed-accuracy trade-off for her/his own problem. See 'Examples' for more details.

If some of the parameters are known and fixed, they should be set to their values via the argument paravars when the model is given by formula. In this case, the user must provide the number of parameters via the argument npar for verifying the general equivalence theorem. See 'Examples', Section 'Weibull', 'Richards' and 'Exponential' model.

crtfunc is a function that is used to specify a new criterion. Its arguments are:

• design points x (as a vector).

• design weights w (as a vector).

• model parameters as follows.

  • If formula is specified: they should be the same parameter specified by parvars. Note that crtfunc must be vectorized with respect to the parameters. The parameters enter the body of crtfunc as a vector with dynamic length.

  • If FIM is specified via the argument fimfunc: param that is a matrix where its row is a vector of parameters values.

• fimfunc is a function that takes the other arguments of crtfunc and returns the computed Fisher information matrices for each parameter vector. The output is a list of matrices.

The crtfunc function must return a vector of criterion values associated with the vector of parameter values. The sensfunc is the optional sensitivity function for the criterion crtfunc. It has one more argument than crtfunc, which is xi_x. It denotes the design point of the degenerate design and must be a vector with the same length as the number of predictors. For more details, see 'Examples'.

Value

an object of class minimax that is a list including three sub-lists:

arg

A list of design and algorithm parameters.

evol

A list of length equal to the number of iterations that stores the information about the best design (design with the minimum criterion value) of each iteration as follows: evol[[iter]] contains:

iter		Iteration number.
x		Design points.
w		Design weights.
min_cost		Value of the criterion for the best imperialist (design).
mean_cost		Mean of the criterion values of all the imperialists.
sens		An object of class 'sensminimax'. See below.

empires

A list of all the empires of the last iteration.

alg

A list with following information:

nfeval		Number of function evaluations. It does not count the function evaluations from checking the general equivalence theorem.
nlocal		Number of successful local searches.
nrevol		Number of successful revolutions.
nimprove		Number of successful movements toward the imperialists in the assimilation step.
convergence		Stopped by 'maxiter' or 'equivalence'?

method

A type of optimal designs used.

design

Design points and weights at the final iteration.

out

A data frame of design points, weights, value of the criterion for the best imperialist (min_cost), and Mean of the criterion values of all the imperialistsat each iteration (mean_cost).

The list sens contains information about the design verification by the general equivalence theorem. See sensbayes for more Details. It is only given every ICA.control$checkfreq iterations and also the last iteration if ICA.control$checkfreq >= 0. Otherwise, NULL.

References

Atashpaz-Gargari, E, & Lucas, C (2007). Imperialist competitive algorithm: an algorithm for optimization inspired by imperialistic competition. In 2007 IEEE congress on evolutionary computation (pp. 4661-4667). IEEE.
Masoudi E, Holling H, Duarte BP, Wong Wk (2019). Metaheuristic Adaptive Cubature Based Algorithm to Find Bayesian Optimal Designs for Nonlinear Models. Journal of Computational and Graphical Statistics. <doi:10.1080/10618600.2019.1601097>

See Also

sensbayes

Examples

#############################################
# Two parameter logistic model: uniform prior
#############################################
# set the unfirom prior
uni <- uniform(lower =  c(-3, .1), upper = c(3, 2))
# set the logistic model with formula
res1 <- bayes(formula = ~1/(1 + exp(-b *(x - a))),
              predvars = "x", parvars = c("a", "b"),
              family = binomial(), lx = -3, ux = 3,
              k =  5, iter = 1, prior = uni,
              ICA.control = list(rseed = 1366))

## Not run: 
  res1 <- update(res1, 500)
  plot(res1)

## End(Not run)
# You can also use your  Fisher information matrix (FIM) if you think it is faster!
## Not run: 
  bayes(fimfunc = FIM_logistic, lx = -3, ux = 3, k =  5, iter = 500,
        prior = uni, ICA.control = list(rseed = 1366))

## End(Not run)

# with fixed x
## Not run: 
  res1.1 <- bayes(formula = ~1/(1 + exp(-b *(x - a))),
                  predvars = "x", parvars = c("a", "b"),
                  family = binomial(), lx = -3, ux = 3,
                  k =  5, iter = 100, prior = uni,
                  x = c( -3, -1.5, 0,  1.5, 3),
                  ICA.control = list(rseed = 1366))
  plot(res1.1)
  # not optimal

## End(Not run)

# with quadrature formula
## Not run: 
  res1.2 <- bayes(formula = ~1/(1 + exp(-b *(x - a))),
                  predvars = "x", parvars = c("a", "b"),
                  family = binomial(), lx = -3, ux = 3,
                  k =  5, iter = 1, prior = uni,
                  crt.bayes.control = list(method = "quadrature"),
                  ICA.control = list(rseed = 1366))
  res1.2 <- update(res1.2, 500)
  plot(res1.2) # not optimal
  # compare it with res1 that was found by automatic integration
  plot(res1)

  # we increase the number of quadrature nodes
  res1.3 <- bayes(formula = ~1/(1 + exp(-b *(x - a))),
                  predvars = "x", parvars = c("a", "b"),
                  family = binomial(), lx = -3, ux = 3,
                  k =  5, iter = 1, prior = uni,
                  crt.bayes.control = list(method = "quadrature",
                                           quadrature = list(level = 9)),
                  ICA.control = list(rseed = 1366))
  res1.3 <- update(res1.3, 500)
  plot(res1.3)
  # by automatic integration (method = "cubature"),
  #  we did not need to worry about the number of nodes.

## End(Not run)
###############################################
# Two parameter logistic model: normal prior #1
###############################################
# defining the normal prior #1
norm1 <- normal(mu =  c(0, 1),
                sigma = matrix(c(1, -0.17, -0.17, .5), nrow = 2),
                lower =  c(-3, .1), upper = c(3, 2))
## Not run: 
  # initializing
  res2 <- bayes(formula = ~1/(1 + exp(-b *(x - a))), predvars = "x", parvars = c("a", "b"),
                family = binomial(), lx = -3, ux = 3, k =  4, iter = 1, prior = norm1,
                ICA.control = list(rseed = 1366))
  res2 <- update(res2, 500)
  plot(res2)

## End(Not run)

###############################################
# Two parameter logistic model: normal prior #2
###############################################
# defining the normal prior #1
norm2 <- normal(mu =  c(0, 1),
                sigma = matrix(c(1, 0, 0, .5), nrow = 2),
                lower =  c(-3, .1), upper = c(3, 2))
## Not run: 
  # initializing
  res3 <- bayes(formula = ~1/(1 + exp(-b *(x - a))), predvars = "x", parvars = c("a", "b"),
                family = binomial(), lx = -3, ux = 3, k =  4, iter = 1, prior = norm2,
                ICA.control = list(rseed = 1366))

  res3 <- update(res3, 700)
  plot(res3,
       sens.bayes.control = list(cubature = list(maxEval = 3000, tol = 1e-4)),
       sens.control = list(optslist = list(maxeval = 3000)))

## End(Not run)


######################################################
# Two parameter logistic model: skewed normal prior #1
######################################################
skew1 <- skewnormal(xi = c(0, 1),
                    Omega = matrix(c(1, -0.17, -0.17, .5), nrow = 2),
                    alpha = c(1, 0), lower =  c(-3, .1), upper = c(3, 2))
## Not run: 
  res4 <- bayes(formula = ~1/(1 + exp(-b *(x - a))), predvars = "x", parvars = c("a", "b"),
                family = binomial(), lx = -3, ux = 3, k =  4, iter = 700, prior = skew1,
                ICA.control = list(rseed = 1366, ncount = 60))
  plot(res4,
       sens.bayes.control = list(cubature = list(maxEval = 3000, tol = 1e-4)),
       sens.control = list(optslist = list(maxeval = 3000)))

## End(Not run)


######################################################
# Two parameter logistic model: skewed normal prior #2
######################################################
skew2 <- skewnormal(xi = c(0, 1),
                    Omega = matrix(c(1, -0.17, -0.17, .5), nrow = 2),
                    alpha = c(-1, 0), lower =  c(-3, .1), upper = c(3, 2))
## Not run: 
  res5 <- bayes(formula = ~1/(1 + exp(-b *(x - a))), predvars = "x", parvars = c("a", "b"),
                family = binomial(), lx = -3, ux = 3, k =  4, iter = 700, prior = skew2,
                ICA.control = list(rseed = 1366, ncount = 60))
  plot(res5,
       sens.bayes.control = list(cubature = list(maxEval = 3000, tol = 1e-4)),
       sens.control = list(optslist = list(maxeval = 3000)))

## End(Not run)

###############################################
# Two parameter logistic model: t student prior
###############################################
# set the prior
stud <- student(mean =  c(0, 1), S   = matrix(c(1, -0.17, -0.17, .5), nrow = 2),
                df = 3, lower =  c(-3, .1), upper = c(3, 2))
## Not run: 
  res6 <- bayes(formula = ~1/(1 + exp(-b *(x - a))), predvars = "x", parvars = c("a", "b"),
                family = binomial(), lx = -3, ux = 3, k =  5, iter = 500, prior = stud,
                ICA.control = list(ncount = 50, rseed = 1366))
  plot(res6)

## End(Not run)
# not bad, but to find a very accurate designs we increase
# the ncount to 200 and repeat the optimization
## Not run: 
  res6 <- bayes(formula = ~1/(1 + exp(-b *(x - a))),
                predvars = "x", parvars = c("a", "b"),
                family = binomial(), lx = -3, ux = 3, k =  5, iter = 1000, prior = stud,
                ICA.control = list(ncount = 200,  rseed = 1366))
  plot(res6)

## End(Not run)


##############################################
# 4-parameter sigmoid Emax model: unform prior
##############################################
lb <- c(4, 11, 100, 5)
ub <- c(8, 15, 130, 9)
## Not run: 
  res7 <- bayes(formula = ~ theta1 + (theta2 - theta1)*(x^theta4)/(x^theta4 + theta3^theta4),
                predvars = c("x"), parvars = c("theta1", "theta2", "theta3", "theta4"),
                lx = .001, ux = 500, k = 5, iter = 200, prior = uniform(lb, ub),
                ICA.control = list(rseed = 1366, ncount = 60))
  plot(res7,
       sens.bayes.control = list(cubature = list(maxEval = 500, tol = 1e-3)),
       sens.control = list(optslist = list(maxeval = 500)))

## End(Not run)

#######################################################################
# 2-parameter Cox Proportional-Hazards Model for type one cenosred data
#######################################################################
# The Fisher information matrix is available here with name FIM_2par_exp_censor1
# However, we should reparameterize the function to match the standard of the argument 'fimfunc'
myfim <- function(x, w, param)
  FIM_2par_exp_censor1(x = x, w = w, param = param, tcensor = 30)
## Not run: 
  res8 <- bayes(fimfunc = myfim, lx = 0, ux = 1, k = 4,
                iter = 1, prior = uniform(c(-11, -11), c(11, 11)),
                ICA.control = list(rseed = 1366))

  res8 <- update(res8, 200)
  plot(res8,
       sens.bayes.control = list(cubature = list(maxEval = 500, tol = 1e-3)),
       sens.control = list(optslist = list(maxeval = 500)))

## End(Not run)


#######################################################################
# 2-parameter Cox Proportional-Hazards Model for random cenosred data
#######################################################################
# The Fisher information matrix is available here with name FIM_2par_exp_censor2
# However, we should reparameterize the function to match the standard of the argument 'fimfunc'
myfim <- function(x, w, param)
  FIM_2par_exp_censor2(x = x, w = w, param = param, tcensor = 30)
## Not run: 
  res9 <- bayes(fimfunc = myfim, lx = 0, ux = 1, k = 2,
                iter = 200, prior = uniform(c(-11, -11), c(11, 11)),
                ICA.control = list(rseed = 1366))
  plot(res9,
       sens.bayes.control = list(cubature = list(maxEval = 100, tol = 1e-3)),
       sens.control = list(optslist = list(maxeval = 100)))

## End(Not run)

#################################
# Weibull model: Uniform prior
################################
# see Dette, H., & Pepelyshev, A. (2008).
# Efficient experimental designs for sigmoidal growth models.
# Journal of statistical planning and inference, 138(1), 2-17.

## See how we fixed a some parameters in Bayesian designs
## Not run: 
  res10 <- bayes(formula = ~a - b * exp(-lambda * t ^h),
                 predvars = c("t"),
                 parvars = c("a=1", "b=1", "lambda", "h=1"),
                 lx = .00001, ux = 20,
                 prior = uniform(.5, 2.5), k = 5, iter = 400,
                 ICA.control = list(rseed = 1366))
  plot(res10)

## End(Not run)

#################################
# Weibull model: Normal prior
################################
norm3 <- normal(mu = 1, sigma = .1, lower = .5, upper = 2.5)
res11 <- bayes(formula = ~a - b * exp(-lambda * t ^h),
               predvars = c("t"),
               parvars = c("a=1", "b=1", "lambda", "h=1"),
               lx = .00001, ux = 20, prior = norm3, k = 4, iter = 1,
               ICA.control = list(rseed = 1366))

## Not run: 
  res11 <- update(res11, 400)
  plot(res11)

## End(Not run)

#################################
# Richards model: Normal prior
#################################
norm4 <- normal(mu = c(1, 1), sigma = matrix(c(.2, 0.1, 0.1, .4), 2, 2),
                lower = c(.4, .4), upper = c(1.6, 1.6))
## Not run: 
  res12 <- bayes(formula = ~a/(1 + b * exp(-lambda*t))^h,
                 predvars = c("t"),
                 parvars = c("a=1", "b", "lambda", "h=1"),
                 lx = .00001, ux = 10,
                 prior = norm4,
                 k = 5, iter = 400,
                 ICA.control = list(rseed = 1366))
  plot(res12,
       sens.bayes.control = list(cubature = list(maxEval = 1000, tol = 1e-3)),
       sens.control = list(optslist = list(maxeval = 1000)))
  ## or we can use the quadrature formula to plot the derivative function
  plot(res12,
       sens.bayes.control = list(method = "quadrature"),
       sens.control = list(optslist = list(maxeval = 1000)))


## End(Not run)

#################################
# Exponential model: Uniform prior
#################################
## Not run: 
res13 <- bayes(formula = ~a + exp(-b*x), predvars = "x",
               parvars = c("a = 1", "b"),
               lx = 0.0001, ux = 1,
               prior = uniform(lower = 1, upper = 20),
               iter = 300, k = 3,
               ICA.control= list(rseed = 100))
  plot(res13)

## End(Not run)

#################################
# Power logistic model
#################################
# See, Duarte, B. P., & Wong, W. K. (2014).
# A Semidefinite Programming based approach for finding
# Bayesian optimal designs for nonlinear models
uni1 <- uniform(lower = c(-.3, 6, .5), upper = c(.3, 8, 1))
## Not run: 
  res14 <- bayes(formula = ~1/(1 + exp(-b *(x - a)))^s, predvars = "x",
                 parvars = c("a", "b", "s"),
                 lx = -1, ux = 1, prior = uni1, k = 5, iter = 1)
  res14 <- update(res14, 300)
  plot(res14)

## End(Not run)

############################################################################
# A two-variable generalized linear model with a gamma distributed response
############################################################################
lb <- c(.5, 0, 0, 0, 0, 0)
ub <- c(2, 1, 1, 1, 1, 1)
myformula1 <- ~beta0+beta1*x1+beta2*x2+beta3*x1^2+beta4*x2^2+beta5*x1*x2
## Not run: 
  res15 <- bayes(formula = myformula1,
                 predvars = c("x1", "x2"), parvars = paste("beta", 0:5, sep = ""),
                 family = Gamma(),
                 lx = rep(0, 2), ux = rep(1, 2),
                 prior = uniform(lower = lb, upper = ub),
                 k = 7,iter = 1, ICA.control = list(rseed = 1366))
  res14 <- update(res14, 500)
  plot(res14,
       sens.bayes.control = list(cubature = list(maxEval = 5000, tol = 1e-4)),
       sens.control = list(optslist = list(maxeval = 3000)))

## End(Not run)

#################################
# Three parameter logistic model
#################################
## Not run: 
sigma1 <- matrix(-0.1, nrow = 3, ncol = 3)
diag(sigma1) <- c(.5, .4, .1)
norm5 <- normal(mu =  c(0, 1, .2), sigma = sigma1,
                lower =  c(-3, .1, 0), upper = c(3, 2, .7))
  res16 <- bayes(formula = ~ c + (1-c)/(1 + exp(-b *(x - a))), predvars = "x",
                 parvars = c("a", "b", "c"),
                 family = binomial(), lx = -3, ux = 3,
                 k =  4, iter = 500, prior = norm5,
                 ICA.control = list(rseed = 1366, ncount = 50),
                 crt.bayes.control = list(cubature = list(maxEval = 2500, tol = 1e-4)))
  plot(res16,
       sens.bayes.control = list(cubature = list(maxEval = 3000, tol = 1e-4)),
       sens.control = list(optslist = list(maxeval = 3000)))
  # took 925 second on my system

## End(Not run)