R: Verifying Optimality of The Robust Designs

sensrobust {ICAOD}

R Documentation

Verifying Optimality of The Robust Designs

Description

It plots the sensitivity (derivative) function of the robust criterion at a given approximate (continuous) design and also calculates its efficiency lower bound (ELB) with respect to the optimality criterion. For an approximate (continuous) design, when the design space is one or two-dimensional, the user can visually verify the optimality of the design by observing the sensitivity plot. Furthermore, the proximity of the design to the optimal design can be measured by the ELB without knowing the latter. See, for more details, Masoudi et al. (2017).

Usage

sensrobust(
  formula,
  predvars,
  parvars,
  family = gaussian(),
  x,
  w,
  lx,
  ux,
  prob,
  parset,
  fimfunc = NULL,
  sens.control = list(),
  calculate_criterion = TRUE,
  plot_3d = c("lattice", "rgl"),
  plot_sens = TRUE,
  npar = dim(parset)[2],
  silent = FALSE,
  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.
`x`	Vector of the design (support) points. See 'Details' of `sensminimax` for models with more than one predictors.
`w`	Vector of the corresponding design weights for `x`.
`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`.
`prob`	A vector of the probability measure `\pi` associated with each row of `parset`.
`parset`	A matrix that provides the vector of initial estimates for the model parameters, i.e. support of `\pi`. 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`.
`sens.control`	Control Parameters for Calculating the ELB. For details, see `sens.control`.
`calculate_criterion`	Calculate the optimality criterion? See 'Details' of `sensminimax`.
`plot_3d`	Which package should be used to plot the sensitivity (derivative) function for models with two predictors. Either `"rgl"` or `"lattice"` (default).
`plot_sens`	Plot the sensitivity (derivative) function? Defaults to `TRUE`.
`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]`.
`silent`	Do not print anything? Defaults to `FALSE`.
`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 \Theta be the set initial estimates for the model parameters and \pi be a probability measure having support in \Theta. A design \xi^* is robust with respect to \pi if the following inequality holds for all \boldsymbol{x} \in \chi:

c(\boldsymbol{x}, \pi, \xi^*) = \int_{\pi} tr M^{-1}(\xi^*, \theta)I(\boldsymbol{x}, \theta)\pi(\theta) d(\theta)-p \leq 0,

with equality at all support points of \xi^*. Here, p is the number of model parameters.

ELB is a measure of proximity of a design to the optimal design without knowing the latter. Given a design, let \epsilon be the global maximum of the sensitivity (derivative) function over x \in \chi. ELB is given by

ELB = p/(p + \epsilon),

where p is the number of model parameters. Obviously, calculating ELB requires finding \epsilon and another optimization problem to be solved. The tuning parameters of this optimization can be regulated via the argument sens.minimax.control.

Value

an object of class sensminimax that is a list with the following elements:

type: Argument type that is required for print methods.
optima: A matrix that stores all the local optima over the parameter space. The cost (criterion) values are stored in a column named Criterion_Value. The last column (Answering_Set) shows if the optimum belongs to the answering set (1) or not (0). See 'Details' of sens.minimax.control. Only applicable for minimax or standardized maximin designs.
mu: Probability measure on the answering set. Corresponds to the rows of optima for which the associated row in column Answering_Set is equal to 1. Only applicable for minimax or standardized maximin designs.
max_deriv: Global maximum of the sensitivity (derivative) function (\epsilon in 'Details').
ELB: D-efficiency lower bound. Can not be larger than 1. If negative, see 'Note' in sensminimax or sens.minimax.control.
merge_tol: Merging tolerance to create the answering set from the set of all local optima. See 'Details' in sens.minimax.control. Only applicable for minimax or standardized maximin designs.
crtval: Criterion value. Compare it with the column Crtiterion_Value in optima for minimax and standardized maximin designs.
time: Used CPU time (rough approximation).

Note

Theoretically, ELB can not be larger than 1. But if so, it may have one of the following reasons:

max_deriv is not a GLOBAL maximum. Please increase the value of the parameter maxeval in sens.minimax.control to find the global maximum.
The sensitivity function is shifted below the y-axis because the number of model parameters has not been specified correctly (less value given). Please specify the correct number of model parameters via the argument npar.

Examples

# Verifying a robust design for the two-parameter logistic model
sensrobust(formula = ~1/(1 + exp(-b *(x - a))),
           predvars = c("x"),
           parvars = c("a", "b"),
           family = binomial(),
           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(0.260, 1, 1.739), w = c(0.275, 0.449, 0.275),
           lx = -5, ux = 5)


###################################
# user-defined optimality criterion
##################################
# When the model is defined by the formula interface
# Checking the A-optimality  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))
}

sensrobust(formula = ~1/(1 + exp(-b * (x-a))), predvars = "x",
           parvars = c("a", "b"), family = "binomial",
           crtfunc = Aopt,
           sensfunc = Aopt_sens,
           lx = -3, ux = 3,
           prob = c(.25, .5, .25),
           parset = matrix(c(-2, 0, 2, 1.25, 1.25, 1.25), 3, 2),
           x = c(-2.469, 0, 2.469), w = c(.317, .365, .317))
# not optimal. the optimal design has four points. see the last example in ?robust

[Package ICAOD version 1.0.1 Index]