locally {ICAOD} | R Documentation |
Locally D-Optimal Designs
Description
Finds locally D-optimal designs for linear and 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.control = list(),
initial = NULL,
npar = length(inipars),
plot_3d = c("lattice", "rgl"),
x = NULL,
crtfunc = NULL,
sensfunc = NULL
)
Arguments
formula |
A linear or nonlinear model |
predvars |
A vector of characters. Denotes the predictors in the |
parvars |
A vector of characters. Denotes the unknown parameters in the |
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 |
lx |
Vector of lower bounds for the predictors. Should be in the same order as |
ux |
Vector of upper bounds for the predictors. Should be in the same order as |
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 |
fimfunc |
A function. Returns the FIM as a |
ICA.control |
ICA control parameters. For details, see |
sens.control |
Control Parameters for Calculating the ELB. For details, see |
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 |
npar |
Number of model parameters. Used when |
plot_3d |
Which package should be used to plot the sensitivity (derivative) function for two-dimensional design space. Defaults to |
x |
A vector of candidate design (support) points.
When is not set to |
crtfunc |
(Optional) a function that specifies an arbitrary criterion. It must have especial arguments and output. See 'Details' of |
sensfunc |
(Optional) a function that specifies the sensitivity function for |
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.
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
.
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
Examples
#################################
# 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 <- update(res1, 150)
#stops at iteration 21 because ELB is greater than .95
## End(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,
x = c(.25, .5, .75),
ICA.control= ICA.control(rseed = 100))
plot(res1.1)
# we can not have an optimal design using this x
## End(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 <- update(res2, 100)
# stops at iteration 51
## End(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
## End(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),
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))
## 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))
}
res4 <- locally(formula = ~1/(1 + exp(-b * (x-a))), predvars = "x",
parvars = c("a", "b"), family = "binomial",
lx = -3, ux = 3, inipars = c(1, 1.25),
iter = 1, k = 2,
crtfunc = Aopt,
sensfunc = Aopt_sens,
ICA.control = list(checkfreq = Inf))
## Not run:
res4 <- update(res4, 50)
## 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 param 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))
}
res4.1 <- locally(fimfunc = FIM_logistic,
lx = -3, ux = 3, inipars = c(1, 1.25),
iter = 1, k = 2,
crtfunc = Aopt2,
sensfunc = Aopt_sens2,
ICA.control = list(checkfreq = Inf))
## Not run:
res4.1 <- update(res4.1, 50)
plot(res4.1)
## End(Not run)
# locally c-optimal design
# example from Chaloner and Larntz (1989) Figure 3
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))
}
res4.2 <- locally(formula = ~1/(1 + exp(-b * (x-a))), predvars = "x",
parvars = c("a", "b"), family = "binomial",
lx = -1, ux = 1, inipars = c(0, 7),
iter = 1, k = 2,
crtfunc = c_opt, sensfunc = c_sens,
ICA.control = list(rseed = 1, checkfreq = Inf))
## Not run:
res4.2 <- update(res4.2, 100)
## End(Not run)