ICAOD (version 0.9.8)

locally: Locally D-Optimal Designs

Description

Finds locally D-optimal designs for nonlinear models. It should be used when a vector of initial estimates is available for the unknown model parameters. Locally optimal designs may not be efficient when the initial estimates are far away from the true values of the parameters.

Usage

locally(formula, predvars, parvars, family = gaussian(), lx, ux, iter, k,
  inipars, fimfunc = NULL, ICA.control = list(),
  sens.minimax.control = list(), initial = NULL,
  npar = length(inipars), plot_3d = c("lattice", "rgl"), x = NULL)

Arguments

formula

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

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.

inipars

A vector of initial estimates for the unknown parameters. It must match parvars or the argument param of the function fimfunc, when provided.

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

Control parameters to verify the general equivalence theorem. For details, see the function sens.minimax.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 given, the sensitivity plot may be shifted below the y-axis. When NULL, it is set to length(inipars).

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.

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'?

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.

Details

Let \(M(\xi, \theta_0)\) be the Fisher information matrix (FIM) of a \(k-\)point design \(\xi\) and \(\theta_0\) be the vector of the initial estimates for the unknown parameters. A locally D-optimal design \(\xi^*\) minimizes over \(\Xi\) $$-\log|M(\xi, \theta_0)|.$$

One can adjust the tuning parameters in ICA.control to set a stopping rule based on the general equivalence theorem. See "Examples" below.

References

Masoudi E, Holling H, Wong W.K. (2017). Application of Imperialist Competitive Algorithm to Find Minimax and Standardized Maximin Optimal Designs. Computational Statistics and Data Analysis, 113, 330-345.

See Also

senslocally

Examples

# NOT RUN {
#################################
# Exponential growth model
################################
# See how we set stopping rule by adjusting 'stop_rule', 'checkfreq' and 'stoptol'
# It calls the 'senslocally' function every checkfreq = 50 iterations to
# calculate the ELB. if ELB is greater than stoptol = .95, then the algoithm stops.

# initializing by one iteration
res1 <- locally(formula = ~a + exp(-b*x), predvars = "x", parvars = c("a", "b"),
                lx = 0, ux = 1, inipars = c(1, 10),
                iter = 1, k = 2,
                ICA.control= ICA.control(rseed = 100,
                                         stop_rule = "equivalence",
                                         checkfreq = 20, stoptol = .95))
# }
# NOT RUN {
  # update the algorithm
  res1 <- iterate(res1, 150)
  #stops at iteration 21 because ELB is greater than .95
# }
# NOT RUN {
### fixed x, lx and ux are only required for equivalence theorem
# }
# NOT RUN {
res1.1 <- locally(formula = ~a + exp(-b*x), predvars = "x", parvars = c("a", "b"),
                  lx = 0, ux = 1, inipars = c(1, 10),
                  iter = 100, k = 3,
                  x = c(.25, .5, .75),
                  ICA.control= ICA.control(rseed = 100))
plot(res1.1)
# we can not have an optimal design using this x
# }
# NOT RUN {
################################
## two parameter logistic model
################################
res2 <- locally(formula = ~1/(1 + exp(-b *(x - a))),
                predvars = "x", parvars = c("a", "b"),
                family = binomial(), lx = -3, ux = 3,
                inipars = c(1, 3), iter = 1, k = 2,
                ICA.control= list(rseed = 100, stop_rule = "equivalence",
                                  checkfreq = 50, stoptol = .95))
# }
# NOT RUN {
  res2 <- iterate(res2, 100)
  # stops at iteration 51
# }
# NOT RUN {



################################
# A model with two predictors
################################
# mixed inhibition model
# }
# NOT RUN {
  res3 <- locally(formula =  ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),
                  predvars = c("S", "I"),
                  parvars = c("V", "Km", "Kic", "Kiu"),
                  family = gaussian(),
                  lx = c(0, 0), ux = c(30, 60),
                  k = 4,
                  iter = 300,
                  inipars = c(1.5, 5.2, 3.4, 5.6),
                  ICA.control= list(rseed = 100, stop_rule = "equivalence",
                                    checkfreq = 50, stoptol = .95))
  # stops at iteration 100
# }
# NOT RUN {

# }
# NOT RUN {
  # fixed x
  res3.1 <- locally(formula =  ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),
                  predvars = c("S", "I"),
                  parvars = c("V", "Km", "Kic", "Kiu"),
                  family = gaussian(),
                  lx = c(0, 0), ux = c(30, 60),
                  k = 5,
                  iter = 100,
                  x = c(20, 4, 20, 4, 10,  0, 0, 30, 3, 2),
                  inipars = c(1.5, 5.2, 3.4, 5.6),
                  ICA.control= list(rseed = 100))
# }