ICAOD: robust – R documentation

Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!

robust

Robust D-Optimal Designs

Description

Finds Robust designs or optimal in-average designs for linear and nonlinear models. It is useful when a set of different vectors of initial estimates along with a discrete probability measure are available for the unknown model parameters. It is a discrete version of bayes.

Usage

robust(
  formula,
  predvars,
  parvars,
  family = gaussian(),
  lx,
  ux,
  iter,
  k,
  prob,
  parset,
  fimfunc = NULL,
  ICA.control = list(),
  sens.control = list(),
  initial = NULL,
  npar = dim(parset)[2],
  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`.
`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.
`prob`	A vector of the probability measure π associated with each row of `parset`.
`parset`	A matrix that provides the vector of initial estimates for the model parameters, i.e. support of π. Every row is one vector (`nrow(parset) == length(prob)`). See 'Details'.
`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`.
`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 `dim(parset)[2]`.
`plot_3d`	Which package should be used to plot the sensitivity (derivative) function for models with two predictors. Either `"rgl"` or `"lattice"` (default).
`x`	Vector of the design (support) points. See 'Details' of `sensminimax` for models with more than one predictors.
`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 Θ be a set of initial estimates for the unknown parameters. A robust criterion is evaluated at the elements of Θ weighted by a probability measure π as follows:

B(ξ, Π) = intergation over Θ Ψ(ξ, θ)π(θ) dθ.

A robust design ξ* maximizes B(ξ, π) over the space of all designs.

When the model is given via formula, columns of parset must match the parameters introduced in parvars. Otherwise, when the model is introduced via fimfunc, columns of parset must match the argument param in fimfunc.

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 sensrobust is called for verification. One can also 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.

Note

When a continuous prior distribution for the unknown model parameters is available, use bayes. When only one initial estimates of the unknown model parameters is available (Θ has only one element), use locally.

Examples

# Finding a robust design for the two-parameter logistic model
# See how we set a stopping rule.
# The ELB is computed every checkfreq = 30 iterations
# The optimization stops when the ELB is larger than stoptol = .95
res1 <- robust(formula = ~1/(1 + exp(-b *(x - a))),
               predvars = c("x"), parvars = c("a", "b"),
               family = binomial(),
               lx = -5, ux = 5, prob = rep(1/4, 4),
               parset = matrix(c(0.5, 1.5, 0.5, 1.5, 4.0, 4.0, 5.0, 5.0), 4, 2),
               iter = 1, k =3,
               ICA.control = list(stop_rule = "equivalence",
                                  stoptol = .95, checkfreq = 30))

## Not run: 
  res1 <- update(res1, 100)
  # stops at iteration 51

## End(Not run)


## Not run: 
  res1.1 <- robust(formula = ~1/(1 + exp(-b *(x - a))),
                   predvars = c("x"), parvars = c("a", "b"),
                   family = binomial(),
                   lx = -5, ux = 5, prob = rep(1/4, 4),
                   parset = matrix(c(0.5, 1.5, 0.5, 1.5, 4.0, 4.0, 5.0, 5.0), 4, 2),
                   x = c(-3, 0, 3),
                   iter = 150, k =3)
  plot(res1.1)
  # not optimal

## 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))
}

res2 <- robust(formula = ~1/(1 + exp(-b * (x-a))), predvars = "x",
               parvars = c("a", "b"), family = "binomial",
               lx = -3, ux = 3,
               iter = 1, k = 4,
               crtfunc = Aopt,
               sensfunc = Aopt_sens,
               prob = c(.25, .5, .25),
               parset = matrix(c(-2, 0, 2, 1.25, 1.25, 1.25), 3, 2),
               ICA.control = list(checkfreq = 50, stoptol = .999,
                                  stop_rule = "equivalence",
                                  rseed = 1))
## Not run: 
  res2 <- update(res2, 500)

## 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))
}


res3 <- robust(formula = ~1/(1 + exp(-b * (x-a))), predvars = "x",
               parvars = c("a", "b"), family = "binomial",
               lx = -1, ux = 1,
               parset = matrix(c(0, 7, .2, 6.5), 2, 2, byrow = TRUE),
               prob = c(.5, .5),
               iter = 1, k = 3,
               crtfunc = c_opt, sensfunc = c_sens,
               ICA.control = list(rseed = 1, checkfreq = Inf))

## Not run: 
  res3 <- update(res3, 300)

## End(Not run)

ICAOD

Optimal Designs for Nonlinear Statistical Models by Imperialist Competitive Algorithm (ICA)

v1.0.1

GPL (>= 2)

Authors

Ehsan Masoudi [aut, cre], Heinz Holling [aut], Weng Kee Wong [aut], Seongho Kim [ctb]

Initial release