bayescomp: Bayesian Compound DP-Optimal Designs
In ICAOD: Optimal Designs for Nonlinear Models via ICA

bayescomp

R Documentation

Bayesian Compound DP-Optimal Designs

Description

Finds compound Bayesian DP-optimal designs that meet the dual goal of parameter estimation and increasing the probability of a particular outcome in a binary response model. A compound Bayesian DP-optimal design maximizes the product of the Bayesian efficiencies of a design \xi with respect to D- and average P-optimality, weighted by a pre-defined mixing constant 0 \leq \alpha \leq 1.

Usage

bayescomp(
  formula,
  predvars,
  parvars,
  family = binomial(),
  prior,
  alpha,
  prob,
  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")
)

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`.
`alpha`	A value between 0 and 1. Compound or combined DP-criterion is the product of the efficiencies of a design with respect to D- and average P- optimality, weighted by `alpha`.
`prob`	Either `formula` or a `function`. When function, its argument are `x` and `param`, and they are the same as the arguments in `fimfunc`. `prob` as a function takes the design points and vector of parameters and returns the probability of success at each design points. See 'Examples'.
`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"`.

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, \pi(\theta) is a user-given prior distribution for the vector of unknown parameters \theta and p(x_i, \theta) is the ith probability of success given by x_i in a binary response model. A compound Bayesian DP-optimal design maximizes over \Xi

\int_{\theta \in \Theta} \frac{\alpha}{q}\log|M(\xi, \theta)| + (1- \alpha) \log \left( \sum_{i=1}^k w_ip(x_i, \theta) \right) \pi(\theta) d\theta.

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 its method component is equal to "cubature". 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.

Value

An object of class 'sensminimax'. See sensminimax for structure details.

References

McGree, J. M., Eccleston, J. A., and Duffull, S. B. (2008). Compound optimal design criteria for nonlinear models. Journal of Biopharmaceutical Statistics, 18(4), 646-661.

Examples

##########################################################################
# DP-optimal design for a logitic model with two predictors: with formula
##########################################################################
p <- c(1, -2, 1, -1)
myprior <- uniform(p -1.5, p + 1.5)
myformula1 <- ~exp(b0+b1*x1+b2*x2+b3*x1*x2)/(1+exp(b0+b1*x1+b2*x2+b3*x1*x2))
res1 <- bayescomp(formula = myformula1,
                  predvars = c("x1", "x2"),
                  parvars = c("b0", "b1", "b2", "b3"),
                  family = binomial(),
                  lx = c(-1, -1), ux = c(1, 1),
                  prior = myprior, iter = 1, k = 7,
                  prob = ~1-1/(1+exp(b0 + b1 * x1 + b2 * x2 + b3 * x1 * x2)),
                  alpha = .5, ICA.control = list(rseed = 1366),
                  crt.bayes.control = list(cubature = list(tol = 1e-4, maxEval = 1000)))


 
  res1 <- update(res1, 1000)
  plot(res1, sens.bayes.control = list(cubature = list(tol = 1e-3, maxEval = 1000)))
  # or use quadrature method
  plot(res1, sens.bayes.control= list(method = "quadrature"))


##########################################################################
# DP-optimal design for a logitic model with two predictors: with fimfunc
##########################################################################
# The function of the Fisher information matrix for this model is 'FIM_logistic_2pred'
# We should reparameterize it to match the standard of the argument 'fimfunc'
 
myfim <- function(x, w, param){
  npoint <- length(x)/2
  x1 <- x[1:npoint]
  x2 <- x[(npoint+1):(npoint*2)]
  FIM_logistic_2pred(x1 = x1,x2 = x2, w = w, param = param)
}

## The following function is equivalent to the function created
# by the formula: ~1-1/(1+exp(b0 + b1 * x1 + b2 * x2 + b3 * x1 * x2))
# It returns probability of success given x and param
# x = c(x1, x2) and param = c()

myprob <- function(x, param){
  npoint <- length(x)/2
  x1 <- x[1:npoint]
  x2 <- x[(npoint+1):(npoint*2)]
  b0 <- param[1]
  b1 <- param[2]
  b2 <- param[3]
  b3 <- param[4]
  out <- 1-1/(1+exp(b0 + b1 * x1 + b2 * x2 + b3 * x1 * x2))
  return(out)
}

res2 <- bayescomp(fimfunc = myfim,
                  lx = c(-1, -1), ux = c(1, 1),
                  prior = myprior, iter = 1000, k = 7,
                  prob = myprob, alpha = .5,
                  ICA.control = list(rseed = 1366))
  plot(res2, sens.bayes.control = list(cubature = list(maxEval = 1000, tol = 1e-4)))
  # quadrature with 6 nodes (default)
  plot(res2, sens.bayes.control= list(method = "quadrature"))

ICAOD documentation built on Feb. 3, 2026, 5:07 p.m.