mica: Imperialist Competitive Algorithm to find locally, minimax and standardized maximin D-optimal designs for nonlinear models

Description

Let $\Xi$ be the space of all approximate designs with $k$ 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 $\theta$ be the vector of unknown parameters. A minimax D-optimal design $\xi^*$ is defined by $$\xi^* = \arg \min_{\xi \in \Xi}\max_{\theta \in \Theta} -\log|M(\xi, \theta)|.$$ A standardized maximin D-optimal designs $\xi^*$ is defined by $$ \xi^* = \arg \max_{\xi \in \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 locally D-optimal designs $\xi^*$ is defined by $$\xi^* = \arg \min_{\xi \in \Xi} -\log|M(\xi, \theta_0)|.$$ where the minization is over $\Xi$ and $\theta_0$ is the user-given initial values for the unknown parameters.

Usage

mica(fimfunc, lx, ux, lp, up, iter, k, control = list(),
  control_gosolnp = list(), type, initial = NULL, locally = NULL, ...)

Arguments

fimfunc

the name of the FIM from available FIM functions in the package as a character string or the user-written function that returns the FIM as a matrix. See "Details".

lower bound of the design space $\chi$

upper bound of the design space $\chi$

lower bound of the region of uncertainty $\Theta$. As same order as the argument param of fimfunc.

upper bound of the region of uncertainty $\Theta$. For the locally D-optimal design,lp and up must be fixed to the same values, i.e. lp = up.

iter

maximum number of iterations.

number of design (support) points. Must be larger than the number of model parameters $p$ to avoid singularity of the FIM.

control

a list of contains the tuning parameters of ICA and the design problem. See "Details".

control_gosolnp

tuning parameters of the function gosolnp. Only when 1) type = "standardized" 2) fimfunc is a character string 3) locally D-optimal design is not available in a closed-form. See "Details" of equivalence for its components.

type

a character string that shows the type of optimal design. "minimax" for a minimax optimal design, "standardized" for a standardized maximin D-optimal design and "locally" for a locally D-optimal design. When type = "locally", then lp must be set equal to up.

initial

a matrix of the initial designs that will be used as initial countries in ICA. Every row is a design and concatenation of x and w. Default is NULL. See "Details" on how to define it when the design space is not of one dimension.

locally

a function that returns the value of determinant of FIM for the locally D-optimal design, i.e. $|M(\xi_{\bold{\theta}}, \bold{\theta})|$. Only required when type = "standardized". See "Details" of equivalence.

...

further arguments to be passed to the FIM function given by fimfunc.

Value

an object of class "ICA" that is a list containing three sub-lists:

arg

a list of design and algorithm parameters

evol

a list that stores the best design (the most powerful imperialist) in each iteration. The output design is stored in the last component of the list. evol[[iter]] contains:

`iter`		iteration.
`x`		design points
`w`		design weights
`min_cost`		cost (criterion value) of the best imperialist
`mean_cost`		mean of costs of all imperialists
`all_optima`		all optima of the inner problem. It is `NA` for a locally optimal design.
`all_optima_cost`		cost of all optima of the inner problem. It is `NA` for a locally optimal design.
`answering`		answering set. It is `NA` for a locally optimal design.
`answering_cost`		cost of each element of the answering set. It is `NA` for a locally optimal design.
`mu`		the probability measure on the answering set. It is `NA` for a locally optimal design.
`max_deriv`		maximum of the sensitivity function.
`ELB`		D-efficiency lower bound
`inner_param`		inner parameter

empires

a list of empires of the last iteration.

alg

a list with the following components:

`nfeval`		number of function evaluations (the function evaluations (calls) for checking the equivalence theorem is not counted)
`nlocal`		number of successful local search
`nrevol`		number of successful revolutions
`nimprove`		number of successful movements toward the imperialists in assimilation step
`convergence`		Why the algorithm has been stopped?

Given an iteration number, please note that all_optima, all_optima_cost, answering, answering_cost, mu, max_deriv and ELB are NA when the equivalence theorem was not requested for that iteration. For example, if equivalence_every = 100, then the equivalence theorem is only checked in every 100 iteration for the best designs.

Details

fimfunc as a function must have three arguments: 1) design points x, 2) weights w and 3) model parameters param. The output should be of type matrix. Further parameters can be set, but should be passed by ... in mica (like parameter $s$ in power logistic model). See "Examples". fimfunc as a character string must be the name of the FIM functions available in this package. The implemented FIM functions are as follows:

`fiumfunc = "FIM_logistic"`	equivalent to	`fiumfunc = FIM_logistic`
`fiumfunc = "FIM_logistic_4par"`	equivalent to	`fiumfunc = FIM_logistic_4par`
`fiumfunc = "FIM_power_logistic"`	equivalent to	`fiumfunc = FIM_power_logistic`
`fiumfunc = "FIM_michaelis"`	equivalent to	`fiumfunc = FIM_michaelis`
`fiumfunc = "FIM_emax_3par"`	equivalent to	`fiumfunc = FIM_emax_3par`
`fiumfunc = "FIM_loglin"`	equivalent to	`fiumfunc = FIM_loglin`
`fiumfunc = "FIM_exp_2par"`	equivalent to	`fiumfunc = FIM_exp_2par`
`fiumfunc = "FIM_exp_3par"`	equivalent to	`fiumfunc = FIM_exp_3par`
`fiumfunc = "FIM_comp_inhibition"`	equivalent to	`fiumfunc = FIM_comp_inhibition`
`fiumfunc = "FIM_noncomp_inhibition"`	equivalent to	`fiumfunc = FIM_noncomp_inhibition`
`fiumfunc = "FIM_uncomp_inhibition"`	equivalent to	`fiumfunc = FIM_uncomp_inhibition`
`fiumfunc = "FIM_mixed_inhibition"`	equivalent to	`fiumfunc = FIM_mixed_inhibition`

See each information matrix for information about the parameters and the analytical locally optimal design (if available).

The control argument is a list that can supply any of the following components:

ncount: number of countries. Defaults to 40.
nimp: number of imperialists. Defaults to 10 percent of ncount.
assim_coeff: assimilation coefficient. Defaults to 4.
revol_rate: revolution rate. Defaults to 0.3
damp: damp ratio for revolution rate. revol_rate is decreased by damp in every iteration. Defaults to 0.99
zeta: a coefficient to find the 'total cost' of empires. Defaults to 0.1.
uniting_threshold: if the distance between two imperialists is less than the product of the uniting_threshold and the largest distance in the search space, then ICA unites these two empires. defaults to 0.02.
assim_strategy: a character strings that denotes the assimilation strategy; PICA for perturbed ICA and ICA for the original version. Defaults to PICA.
lsearch: logical. Perform a local search on imperialists in every iteration? Defaults to TRUE.
l: positive integer. the number of local search on each imperialist. Defaults to 2.
only_improve: logical, if TRUE, the new position in the assimilation step replaces the current position if the value of the criterion in the new position is better. . Defaults to TRUE.
stop_rule: a character string denotes stopping rule. When "maxiter", then ICA only stops when it reaches the maximum number of iterations. When "one_empire", then ICA stops if either all empires collapsed and one empire remains or it reaches the maximum number of iterations. When "equivalence", then ICA stops if either the Efficiency lower bound (ELB) of the current design is greater than stoptol or it reaches the maximum number of iterations.
stoptol: numeric between $0$ and $1$. The minimum ELB for the best current imperialist (best design) to stop the algorithm by equivalence theorem when stop_rule = "equivalence". Defaults to 0.99.
equivalence_every: a positive integer. Check and compute ELB in every equivalence_every iteration. Checking equivalence theorem in small intervals slows down the algorithm. Defaults to 200.
equal_weight: logical; whether the points should have equal weights. Defaults to FALSE.
sym: logical. Whether the design is symmetric around a point. If TRUE then sym_point must be given. Defaults to FALSE
sym_point: a vector of the same length as lx. The point that the design is symmetric around. Must have the same length as lx. See "Examples".
inner_space: a character string on how to deal with the parameter space (inner space): "continuous", "vertices" or "discrete". Defaults to "continuous". See "Details".
param_set: a matrix that denotes the fixed values for the parameters and is required when inner_space = "discrete". Each row of the matrix is the values of the components of the parameters, i.e. $\theta_{01}$ in the first row, $\theta_{02}$ in the second row and so on. The number of columns should be equal to length(lp). See "Examples".
inner_maxeval: maximum number of function evaluations for the continuous inner problem. It comes from the tuning parameters of directL. Its value should be large enough to not miss any global optima. It is only applicable for standardized maxinim and minimax optimal designs. Defaults to 600.
plot_deriv: logical. Should derivative be plotted whenever the equivalence theorem is checked? Defaults to TRUE.
plot_cost: logical. Should the evolution of ICA, i.e. mean cost of all imperialists and cost of the best imperialist, be plotted in every iteration? Defaults to TRUE; when type = "locally" is FALSE.
trace: logical. Should the best generated design (best imperialist) and corresponding algorithm parameters be printed in every iteration? Defaults to TRUE.
n.seg: a positive integer required when checking the equivalence theorem to construct the answering set. The number of initial starting points for the local optimizer (optim) to find all minima of the criterion on parameters space is equal to (n.seg + 1)^length(lp). Defaults to 4. See "Details" in equivalence.

Each row of initial is one design, i.e. concatenation of $\bold{x} = (x_1,...x_k)$ and $\bold{w} = (w_1,...,w_k)$. The number of columns of initial is equal to $k \times n + k$, where $n$ is the number of model explanatory variables. For multi-dimensional design space, $\bold{x}$ must be given dimension by dimension. See description of argument x in equivalence.

In control, if inner_space = "continuous", then the inner problem is an optimization over $\Theta =$ (lp, up). If inner_space = "discrete", then the inner problem is a discrete optimization problem over a set of initial values for the parameters, i.e. $\Theta = \{\theta_{01}, \theta_{02},...\}$. In this case, the set of initial parameters should be given through param_set. If inner_space = "vertices" then the set of initial parameters are the vertices of $\Theta$. This should be set when the user is certain that the D-criterion is convex with respect to the parameter space for every given design. Please note that regardless of what inner_space is chosen, checking the equivalence theorem will be done on the continuous parameter space $\Theta$. See "Examples" on how to use this option.

For large parameter space or complex models it is important to increase ncount, inner_maxeval and n.seg (for checking the equivalence theorem).

Please note that the speed of mica highly depends on the inner_maxeval parameter in control when inner_space = "continuous".

References

Masoudi, E., Holling, H., & Wong, W. K. (in press). Application of imperialist competitive algorithm to find minimax and standardized maximin optimal designs. Computational Statistics & Data Analysis.

Examples

Run this code

# NOT RUN {
#######################################################################
## some examples for exponential model
# }
# NOT RUN {
# finding standardized maximin D-optimal design
res <- mica(fimfunc = "FIM_exp_2par", lx = 0, ux = 1, lp = c(1, 1), up = c(1, 5),
    iter = 100, k = 3, type = "standardized", control = list(seed = 215))
res <- iterate(res, 10)
plot(res)
# finding minimax D-optimal design
mica(fimfunc = "FIM_exp_2par", lx = 0, ux = 1, lp = c(1, 1), up = c(1, 5),
    iter = 100, k = 3, type = "minimax", control = list(seed = 215))
# }
# NOT RUN {
# finding locally D-optimal design. Please note that  'lp' and 'up' are equal
mica(fimfunc = "FIM_exp_2par", lx = 0, ux = 1, lp = c(2, 3), up = c(2, 3),
    iter = 40, k = 2, type = "locally", control = list(seed = 215))
# locally D-optimal design is x1 = lx, x2 = 1/lp[2]

# requesting an equally-weighted design, i.e w_1 = w_2 = ... w_k
res_loc <-mica(fimfunc = "FIM_exp_2par", lx = 0, ux = 1, lp = c(2, 3), up = c(2, 3),
              iter = 40, k = 2, type = "locally",
               control = list(seed = 215, equal_weight = TRUE))
# }
# NOT RUN {
res_loc <- iterate(res_loc, 10) ## update the result
plot(res_loc)
# using symetric option for the logisitic model
mica(fimfunc = "FIM_logistic", lx = -5, ux = 5, lp = c(0, 1), up = c(3.5 , 1.25),
    iter = 100, k = 5, control = list(rseed = 215, sym = TRUE,
     sym_point = (0 + 3.5)/2),type = "minimax")
#######################################################################
# 2PL model
mica(fimfunc = "FIM_logistic", lx = -5, ux = 5, lp = c(-1, 1), up = c(1 , 2),
    iter = 100, k = 3, control = list(rseed = 215), type = "minimax")

# an example on how to supply 'fimfunc' with a function
logistic_fim <- function(x, w, param){
 a <- param[1]
 b <- param[2]
 constant <- 1/(1 + exp(-b * (x - a)))
 constant <- constant * (1 - constant)
 A <-  sum(w * b^2 * constant)
 B <- sum(-w * b * (x - a) * constant)
 C <- sum(w * ((x -a)^2) * constant)
 mat <- matrix(c(A, B, B, C), 2, 2)
 return(mat)
}

mica(fimfunc = logistic_fim, lx = -5, ux = 5, lp = c(-1, 1), up = c(1 , 2),
    iter = 100, k = 3,  control = list(rseed = 215), type = "minimax")
     ## is the same when 'fimfunc = "FIM_logistic'

mica(fimfunc = logistic_fim, lx = -5, ux = 5, lp = c(-1, 1), up = c(-1 , 1),
    iter = 100, k = 3, control = list(rseed = 215), type = "locally")
#######################################################################

#######################################################################
## how to use inner_space option in control list

#### Enzyme kinetic models. examples for3D plots
mica(fimfunc = "FIM_comp_inhibition", lx = c(0, 0), ux = c(30, 60),
    lp = c(7, 4, 2), up = c(7, 5, 3), k =3, type = "standardized",
    iter = 300, control = list(rseed = 215, inner_maxit = 300,
                               stop_rule = "equivalence"))

## setting the parameter space as only the points on the vertices
mica(fimfunc = "FIM_comp_inhibition", lx = c(0, 0), ux = c(30, 60),
    lp = c(7, 4, 2), up = c(7, 5, 3), k =3, type = "standardized",
    iter = 300, control = list(rseed = 215, inner_space = "vertices",
                               stop_rule = "equivalence"))

## every row is one of the vertices of Theta
param_set <- matrix(c(7, 4, 2, 7, 5, 2, 7, 4, 3, 7, 5, 3),
                   ncol = 3, nrow = 4, byrow = TRUE)
res <-mica(fimfunc = "FIM_comp_inhibition", lx = c(0, 0), ux = c(30, 60),
          lp = c(7, 4, 2), up = c(7, 5, 3), k =3, type = "standardized",
          iter = 300, control = list(rseed = 215,inner_space = "discrete",
                                     stop_rule = "equivalence",
                                      param_set = param_set))


#######################################################################

#######################################################################
## optimal designs for the 1Pl model
mica(fimfunc = "FIM_logistic_1par", lx = 0, ux = 5,
     lp = 2, up = 2, k = 3, iter = 100, type = "locally")

lx <- -.5
ux <- .5
ux - lx <= 2 * log(2 + sqrt(3))
mica(fimfunc = "FIM_logistic_1par", lx = lx, ux = ux,
    lp = -1, up = 1, k = 1, iter = 10,
    type = "standardized")



lx <- -2
ux <- 2
ux - lx <= 2 * log(2 + sqrt(3))
mica(fimfunc = "FIM_logistic_1par", lx = lx, ux = ux,
    lp = -1, up = 1, k = 1, iter = 10,
   type = "standardized")
#######################################################################

# }
# NOT RUN {
# }

Run the code above in your browser using DataLab