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 $x1, x2, ..., xk$ from design space $\chi$ with corresponding weights $w1, . . . ,wk$. 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)|.$$ Throughout, the terms ``inner problem'' and ``outer problem'' are used for optimization over $\Theta$ and $\Xi$, respectively. 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 paramters 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 $\theta0$ is the 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

FIM as a function or character. As a character, it must be the name of the FIM function that are available, e.g. as.character(substitute(FIM_logistic)). As a function it must return the information matrix. See "Details".

lower bound of the design space $\chi$.

upper bound of the design space $\chi$.

lower bound of the region of unceratinty $\Theta$. The order of the dimension is the same as the order of the parameters in the argument paramof fimfunc.

upper bound of the region of unceratinty $\Theta$. If lp = up, then $\Theta = \theta0$ and the genearted design is locally D-optimal design. See "Examples".

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 control parameters. See "Details".

control_gosolnp

tuning parameters of function gosolnp for models that their locally optimal design do not have an analytical solution and is find by gosolnp. Only required when type is set to 'standardized'. See "Details" of equivalence.

type

a character strings; "minimax" for minimax optimal design, "standardized" for standardized maximin D-optimal design and "locally" for locally D-optimal design. When "locally", then lp must be set equal to up.

initial

a matrix of user intial countries or a vector of a country that will be inserted into the initial countries of ICA. See "Details" .

locally

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

...

further arguments to be passed to the FIM function corresponding to fimfunc. For power logisitc model when fimfunc is equal to "FIM_power_logistic", the value of s must be given here.

Value

an object of class "ICA" that is a list contains another three lists:

The design in best is the same best imperialist generated in last iteration and stored in evol.

Details

fimfunc as a function must have three arguments: 1) desig 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 defined in this package. The implemented FIM functions are given 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`

Setting fimfunc as a character strings only results in using the internal defined locally for standardized maximin optimal designs. 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:

Each row of intial is one design, i.e. concatenation of $x = (x1,...xk)$ and $w = (w1,...,wk)$. 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 = {\theta01, \theta02, ....}$. In this case, the set of initial parameters should be given through param_set. If inner_space = "vertices" then the set of intial 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, checking the equivalence theorem is done on 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". equivalence_every and l in local search are the other factors that impact the CPU time.

From Section "Value", note that all_optima, all_optima_cost, answering, answering_cost, mu, max_deriv and DLB are NA when the equivalence theorem was not requested for that iteration by equivalence_every in control. For example, if equivalence_every = 100, then the equivalence theorem is only check for best designs in iteration $100$, $200$, $300$ and so on. From Section "Value", inner_param is equal to $arg max -log|M(\xi, \theta)|,$ for minimax and $ = arg inf {|M(\xi, \theta)| / |M(\xi_\theta, \theta)|}^p,$ for standardized maximin D-optimal designs.

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

#######################################################################
## 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))
# ## End(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",
#                                countries = 100, nimperialists = 10))
# 
# ## 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",
#                                countries = 100, nimperialists = 10))
# 
# ## 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", countries = 100,
#                                      nimperialists = 10, 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")
# #######################################################################
# 
# ## End(Not run)

Run the code above in your browser using DataLab