minimax {ICAOD}R Documentation

Minimax and Standardized Maximin D-Optimal Designs

Description

Finds minimax and standardized maximin D-optimal designs for linear and nonlinear models. It should be used when the user assumes the unknown parameters belong to a parameter region \Theta, which is called “region of uncertainty”, and the purpose is to protect the experiment from the worst case scenario over \Theta.

Usage

minimax(
  formula,
  predvars,
  parvars,
  family = gaussian(),
  lx,
  ux,
  lp,
  up,
  iter,
  k,
  n.grid = 0,
  fimfunc = NULL,
  ICA.control = list(),
  sens.control = list(),
  sens.minimax.control = list(),
  crt.minimax.control = list(),
  standardized = FALSE,
  initial = NULL,
  localdes = NULL,
  npar = length(lp),
  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.

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.

lp

Vector of lower bounds for the model parameters. Should be in the same order as parvars or param in the argument fimfunc.

up

Vector of upper bounds for the model parameters. Should be in the same order as parvars or param in the argument fimfunc. When a parameter is known (has a fixed value), its associated lower and upper bound values in lp and up must be set equal.

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.

n.grid

Only required when the parameter space is going to be discretized. The total number of grid points from the parameter space is n.grid^p. When n.grid > 0, optimal design protects the experimenter against the worst case scenario only over the grid points, and not over the continuous parameter space. The resulting designs may not be globally optimal. In some literature, this type of designs has been used as a compromise to the minimax type designs to avoid continuous optimization problem over the parameter space and simplify the minimax design problems. Especially when the design criterion is convex with respect to the given parameter space at every given design from the design space, the obtained design may also be globally optimal (because the maximum of a convex function is attained on the bounds, and here, are included in the grid points). See 'Details' of minimax.

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.

sens.minimax.control

Control parameters to construct the answering set required for verify the general equivalence theorem and calculating the ELB. For details, see the function sens.minimax.control.

crt.minimax.control

Control parameters to optimize the minimax or standardized maximin criterion at a given design over a continuous parameter space (when n.grid = 0). For details, see the function crt.minimax.control.

standardized

Maximin standardized design? When standardized = TRUE, the argument localdes must be given. Defaults to FALSE. See 'Details' of minimax.

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.

localdes

A function that takes the parameter values as inputs and returns the design points and weights of the locally optimal design. Required when standardized = "TRUE". 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 truly, the sensitivity (derivative) plot may be shifted below the y-axis. When NULL (default), it is set to length(lp).

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 minimax.

sensfunc

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

Details

Let \Xi be the space of all approximate designs with k design points (support points) at x_1, x_2, ..., x_k from the 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 \theta be the vector of unknown parameters. A minimax D-optimal design \xi^* minimizes over \Xi

\max_{\theta \in \Theta} -\log|M(\xi, \theta)|.

A standardized maximin D-optimal design \xi^* maximizes over \Xi

\inf_{\theta \in \Theta} \left[\left(\frac{|M(\xi, \theta)|}{|M(\xi_{\theta}, \theta)|}\right)^\frac{1}{p}\right],

where p is the number of model parameters and \xi_\theta is the locally D-optimal design with respect to \theta.

A minimax criterion (cost function or objective function) is evaluated at each design (decision variables) by maximizing the criterion over the parameter space. We call the optimization problem over the parameter space as inner optimization problem. Two different strategies may be applied to solve the inner problem at a given design (design points and weights):

  1. Continuous inner problem: we optimize the criterion over a continuous parameter space using the function nloptr. In this case, the tuning parameters can be regulated via the argument crt.minimax.control, when the most influential one is maxeval.

  2. Discrete inner problem: we map the parameter space to the grid points and optimize the criterion over a discrete parameter space. In this case, the number of grid points can be regulated via n.grid. This strategy is quite efficient (ans fast) when the maxima most likely attain the vertices of the continuous parameter space at any given design. The output design here protects the experiment from the worst scenario over the grid points.

The formula is used to automatically create the Fisher information matrix (FIM) for a linear or nonlinear model provided that the distribution of the response variable belongs to the natural exponential family. Function minimax also provides an option to assign a user-defined FIM directly via the argument fimfunc. In this case, the argument fimfunc takes a function that has three arguments as follows:

  1. x a vector of design points. For design points with more than one dimension, it is a concatenation of the design points, but dimension-wise. For example, let the model has three predictors (I, S, Z). Then, a two-point design is of the format \{\mbox{point1} = (I_1, S_1, Z_1), \mbox{point2} = (I_2, S_2, Z_2)\}. and the argument x is equivalent to x = c(I1, I2, S1, S2, Z1, Z2).

  2. w a vector that includes the design weights associated with x.

  3. param a vector of parameter values associated with lp and up.

The output must be the Fisher information matrix with number of rows equal to length(param). See 'Examples'.

Minimax optimal designs can have very different criterion values depending on the nominal set of parameter values. Accordingly, it is desirable to standardize the criterion and control for the potentially widely varying magnitude of the criterion (Dette, 1997). Evaluating a standardized maximin criterion requires knowing locally optimal designs. We strongly advise setting standardized = 'TRUE' only when analytical solutions for the locally D-optimal designs is available. In this case, for any initial estimate of the unknown parameters, an analytical solution for the locally optimal design, i.e, the support points x and the corresponding weights w, must be provided via the argument localdes. Therefore, depending on how the model is specified, localdes is a function with the following arguments (input).

This function must return a list with the components x and w (they match the same arguments in the function fimfunc). See 'Examples'.

The standardized D-criterion is equal to the D-efficiency and it must be between 0 and 1. However, in practice, when running the algorithm, it may be the case that the criterion takes a value larger than one. This may happen because the user-function that is given via localdes does not return the true (accurate) locally optimal designs for some requested initial estimates of the parameters from \Theta. In this case, the function minimax throw an error where the error message helps the user to debug her/his function.

Each row of initial is one design, i.e. a concatenation of values for design (support) points and the associated design weights. Let x0 and w0 be the vector of initial values with exactly the same length and order as x and w (the arguments of fimfunc). As an example, the first row of the matrix initial is equal to initial[1, ] = c(x0, w0). For models with more than one predictors, x0 is a concatenation of the initial values for design points, but dimension-wise. See the details of the argument fimfunc, above.

To verify the optimality of the output design by the general equivalence theorem, the user can either plot the results or set checkfreq in ICA.control to Inf. In either way, the function sensminimax is called for verification. Note that the function sensminimax always verifies the optimality of a design assuming a continues parameter space. See 'Examples'.

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

The crtfunc function must return the criterion value. 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 least criterion value) of each iteration. 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.
param Vector of parameters.
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 sensminimax for more details. It is given every ICA.control$checkfreq iterations and also the last iteration if ICA.control$checkfreq >= 0. Otherwise, NULL.

param is a vector of parameters that is the global minimum of the minimax criterion or the global maximum of the standardized maximin criterion over the parameter space, given the current x, w.

Note

For larger parameter space or model with more number of unknown parameters, it is always important to increase the value of ncount in ICA.control and optslist$maxeval in crt.minimax.control to produce very accurate designs.

Although standardized criteria have been preferred theoretically, in practice, they should be applied only when an analytical solution for the locally D-optimal designs is available for the model of interest. Otherwise, we encounter a three-level nested-optimization algorithm, which is very slow.

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.
Dette, H. (1997). Designing experiments with respect to 'standardized' optimality criteria. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 59(1), 97-110.
Masoudi E, Holling H, Wong WK (2017). Application of Imperialist Competitive Algorithm to Find Minimax and Standardized Maximin Optimal Designs. Computational Statistics and Data Analysis, 113, 330-345. <doi:10.1016/j.csda.2016.06.014>

See Also

sensminimax

Examples

########################################
# Two-parameter exponential growth model
########################################
res1 <- minimax (formula = ~a + exp(-b*x), predvars = "x", parvars = c("a", "b"),
                 lx = 0, ux = 1, lp = c(1, 1), up = c(1, 10),
                 iter = 1, k = 4,
                 ICA.control= ICA.control(rseed = 100),
                 crt.minimax.control = list(optslist = list(maxeval = 100)))
# The optimal design has 3 points, but we set k = 4 for illustration purpose to
#   show how the algorithm modifies the design by adjusting the weights
# The value of maxeval is changed to reduce the CPU time
## Not run: 
  res1 <- update(res1, 150)
  # iterating the algorithm up to 150 more iterations

## End(Not run)

res1 # print method
plot(res1) # Veryfying the general equivalence theorem

## Not run: 
  ## fixed x
  res1.1 <- minimax (formula = ~a + exp(-b*x), predvars = "x", parvars = c("a", "b"),
                     lx = 0, ux = 1, lp = c(1, 1), up = c(1, 10),
                     x = c(0, .5, 1),
                     iter = 150, k = 3, ICA.control= ICA.control(rseed = 100))
  # not optimal

## End(Not run)

########################################
# Two-parameter logistic model.
########################################
# A little playing with the tuning parameters
# The value of maxeval is reduced to 200 to increase the speed
cont1 <- crt.minimax.control(optslist = list(maxeval = 200))
cont2 <- ICA.control(rseed = 100, checkfreq = Inf, ncount = 60)

## Not run: 
  res2 <- minimax (formula = ~1/(1 + exp(-b *(x - a))), predvars = "x",
                   parvars = c("a", "b"),
                   family = binomial(), lx = -3, ux = 3,
                   lp = c(0, 1), up = c(1, 2.5), iter = 200, k = 3,
                   ICA.control= cont2, crt.minimax.control = cont1)
  print(res2)
  plot(res2)

## End(Not run)

############################################
# An example of a model with two predictors
############################################
# Mixed inhibition model
lower <- c(1, 4, 2, 4)
upper <- c(1, 5, 3, 5)
cont <- crt.minimax.control(optslist = list(maxeval = 100)) # to be faster
## Not run: 
  res3 <- minimax(formula =  ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),
                  predvars = c("S", "I"),
                  parvars = c("V", "Km", "Kic", "Kiu"),
                  lx = c(0, 0), ux = c(30, 60), k = 4,
                  iter = 100, lp = lower, up = upper,
                  ICA.control= list(rseed = 100),
                  crt.minimax.control = cont)

  res3 <- update(res3, 100)
  print(res3)
  plot(res3) # sensitivity plot
  res3$arg$time

## End(Not run)

# Now consider grid points instead of assuming continuous parameter space
# set n.grid to 5
## Not run: 
  res4 <- minimax(formula =  ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),
                  predvars = c("S", "I"),
                  parvars = c("V", "Km", "Kic", "Kiu"),
                  lx = c(0, 0), ux = c(30, 60),
                  k = 4, iter = 130, n.grid = 5, lp = lower, up = upper,
                  ICA.control= list(rseed = 100, checkfreq = Inf),
                  crt.minimax.control = cont)
  print(res4)
  plot(res4) # sensitivity plot

## End(Not run)

############################################
# Standardized maximin D-optimal designs
############################################
# Assume the purpose is finding STANDARDIZED designs
# We know from literature that the locally D-optimal design (LDOD)
# for this model has an analytical solution.
# The follwoing function takes the parameter as input and returns
# the design points and weights of LDOD.
# x and w are exactly similar to the arguments of 'fimfunc'.
# x is a vector and returns the design points 'dimension-wise'.
# see explanation of the arguments of 'fimfunc' in 'Details'.

LDOD <- function(V, Km, Kic, Kiu){
  #first dimention is for S and the second one is for I.
  S_min <- 0
  S_max <- 30
  I_min <- 0
  I_max <- 60
  s2 <- max(S_min, S_max*Km*Kiu*(Kic+I_min)/
              (S_max*Kic*I_min+S_max*Kic*Kiu+2*Km*Kiu*I_min+2*Km*Kiu*Kic))
  i3 <- min((2*S_max*Kic*I_min + S_max*Kic*Kiu+2*Km*Kiu*I_min+Km*Kiu*Kic)/
              (Km*Kiu+S_max*Kic), I_max)
  i4 <- min(I_min + (sqrt((Kic+I_min)*(Km*Kic*Kiu+Km*Kiu*I_min+
                                         S_max*Kic*Kiu+S_max*Kic*I_min)/
                            (Km*Kiu+S_max*Kic))), I_max )
  s4 <- max(-Km*Kiu*(Kic+2*I_min-i4)/(Kic*(Kiu+2*I_min-i4)), S_min)
  x <- c(S_max, s2, S_max, s4, I_min, I_min, i3, i4)
  return(list(x = x, w =rep(1/4, 4)))

}
formalArgs(LDOD)
## Not run: 
  minimax(formula =  ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),
          predvars = c("S", "I"),
          parvars = c("V", "Km", "Kic", "Kiu"),
          lx = c(0, 0), ux = c(30, 60),
          k = 4, iter = 300,
          lp = lower, up = upper,
          ICA.control= list(rseed = 100, checkfreq = Inf),
          crt.minimax.control = cont,
          standardized = TRUE,
          localdes = LDOD)

## End(Not run)


################################################################
# Not necessary!
# The rest of the examples here are only for professional uses.
################################################################
# Imagine you have written your own FIM, say in Rcpp that is faster than
# the FIM created by the formula interface above.

###########################################
# An example of a model with two predictors
###########################################
# For example, th cpp FIM function for the mixed inhibition model is named:
formalArgs(FIM_mixed_inhibition)

# We should reparamterize the arguments to match the standard of the
# argument 'fimfunc' (see 'Details').
myfim <- function(x, w, param){
  npoint <- length(x)/2
  S <- x[1:npoint]
  I <- x[(npoint+1):(npoint*2)]
  out <- FIM_mixed_inhibition(S = S, I = I, w = w, param = param)
  return(out)
}
formalArgs(myfim)

# Finds minimax optimal design, exactly as before, but NOT using the
# formula interface.
## Not run: 
  res5 <- minimax(fimfunc = myfim,
                  lx = c(0, 0), ux = c(30, 60), k = 4,
                  iter = 100, lp = lower, up = upper,
                  ICA.control= list(rseed = 100),
                  crt.minimax.control = cont)
  print(res5)
  plot(res5) # sensitivity plot

## End(Not run)
#########################################
# Standardized maximin D-optimal designs
#########################################
# To match the argument 'localdes' when no formula inteface is used,
# we should reparameterize LDOD.
# The input must be 'param' same as the argument of 'fimfunc'
LDOD2 <- function(param)
  LDOD(V = param[1], Km = param[2], Kic = param[3], Kiu = param[4])

# compare these two:
formalArgs(LDOD)
formalArgs(LDOD2)
## Not run: 
  res6 <- minimax(fimfunc = myfim,
                  lx = c(0, 0), ux = c(30, 60), k = 4,
                  iter = 300, lp = lower, up = upper,
                  ICA.control= list(rseed = 100, checkfreq = Inf),
                  crt.minimax.control = cont,
                  standardized = TRUE,
                  localdes = LDOD2)
  res6
  plot(res6)

## End(Not run)

###################################
# user-defined optimality criterion
##################################
# When the model is defined by the formula interface
# A-optimal design for the 2PL model.
# the criterion function must have argument x, w fimfunc and the parameters defined in 'parvars'.
# use 'fimfunc' as a function of the design points x,  design weights w and
#  the 'parvars' parameters whenever needed.
Aopt <-function(x, w, a, b, fimfunc){
  sum(diag(solve(fimfunc(x = x, w = w, a = a, b = b))))
}
## the sensitivtiy function
# xi_x is a design that put all its mass on x in the definition of the sensitivity function
# x is a vector of design points
Aopt_sens <- function(xi_x, x, w, a, b, fimfunc){
  fim <- fimfunc(x = x, w = w, a = a, b = b)
  M_inv <- solve(fim)
  M_x <- fimfunc(x = xi_x, w = 1, a  = a, b = b)
  sum(diag(M_inv %*% M_x %*%  M_inv)) - sum(diag(M_inv))
}
## Not run: 
res7 <- minimax(formula = ~1/(1 + exp(-b * (x-a))), predvars = "x",
                parvars = c("a", "b"), family = "binomial",
                lx = -2, ux = 2,
                lp = c(-2, 1), up = c(2, 1.5),
                iter = 400, k = 3,
                crtfunc = Aopt,
                sensfunc = Aopt_sens,
                crt.minimax.control = list(optslist = list(maxeval = 200)),
                ICA.control = list(rseed = 1))
  plot(res7)

## End(Not run)
# with grid points
res7.1 <- minimax(formula = ~1/(1 + exp(-b * (x-a))), predvars = "x",
                  parvars = c("a", "b"), family = "binomial",
                  lx = -2, ux = 2,
                  lp = c(-2, 1), up = c(2, 1.5),
                  iter = 1, k = 3,
                  crtfunc = Aopt,
                  sensfunc = Aopt_sens,
                  n.grid = 9,
                  ICA.control = list(rseed = 1))
## Not run: 
  res7.1 <- update(res7.1, 400)
  plot(res7.1)

## End(Not run)

# When the FIM of the model is defined directly via the argument 'fimfunc'
# the criterion function must have argument x, w fimfunc and param.
# use 'fimfunc' as a function of the design points x,  design weights w and
#  the 'parvars' parameters whenever needed.
Aopt2 <-function(x, w, param, fimfunc){
  sum(diag(solve(fimfunc(x = x, w = w, param = param))))
}
## the sensitivtiy function
# xi_x is a design that put all its mass on x in the definition of the sensitivity function
# x is a vector of design points
Aopt_sens2 <- function(xi_x, x, w, param, fimfunc){
  fim <- fimfunc(x = x, w = w, param = param)
  M_inv <- solve(fim)
  M_x <- fimfunc(x = xi_x, w = 1, param = param)
  sum(diag(M_inv %*% M_x %*%  M_inv)) - sum(diag(M_inv))
}
## Not run: 
res7.2 <- minimax(fimfunc = FIM_logistic,
                  lx = -2, ux = 2,
                  lp = c(-2, 1), up = c(2, 1.5),
                  iter = 1, k = 3,
                  crtfunc = Aopt2,
                  sensfunc = Aopt_sens2,
                  crt.minimax.control = list(optslist = list(maxeval = 200)),
                  ICA.control = list(rseed = 1))
  res7.2 <- update(res7.2, 200)
  plot(res7.2)

## End(Not run)
# with grid points
res7.3 <- minimax(fimfunc = FIM_logistic,
                  lx = -2, ux = 2,
                  lp = c(-2, 1), up = c(2, 1.5),
                  iter = 1, k = 3,
                  crtfunc = Aopt2,
                  sensfunc = Aopt_sens2,
                  n.grid = 9,
                  ICA.control = list(rseed = 1))
## Not run: 
  res7.3 <- update(res7.2, 200)
  plot(res7.3)

## End(Not run)


# robust c-optimal design
# example from Chaloner and Larntz (1989), Figure 3, but robust design
c_opt <-function(x, w, a, b, fimfunc){
  gam <- log(.95/(1-.95))
  M <- fimfunc(x = x, w = w, a = a, b = b)
  c <- matrix(c(1, -gam * b^(-2)), nrow = 1)
  B <- t(c) %*% c
  sum(diag(B %*% solve(M)))
}

c_sens <- function(xi_x, x, w, a, b, fimfunc){
  gam <- log(.95/(1-.95))
  M <- fimfunc(x = x, w = w, a = a, b = b)
  M_inv <- solve(M)
  M_x <- fimfunc(x = xi_x, w = 1, a = a, b = b)
  c <- matrix(c(1, -gam * b^(-2)), nrow = 1)
  B <- t(c) %*% c
  sum(diag(B %*% M_inv %*% M_x %*%  M_inv)) - sum(diag(B %*% M_inv))
}

## Not run: 
res8 <- minimax(formula = ~1/(1 + exp(-b * (x-a))), predvars = "x",
                parvars = c("a", "b"), family = "binomial",
                lx = -1, ux = 1,
                lp = c(-.3, 6), up = c(.3, 8),
                iter = 500, k = 3,
                crtfunc = c_opt, sensfunc = c_sens,
                ICA.control = list(rseed = 1, ncount = 100),
                n.grid = 12)
  plot(res8)

## End(Not run)
[Package ICAOD version 1.0.1 Index]