senslocally: Verifying Optimality of The Locally D-optimal Designs

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

Description

It plots the sensitivity (derivative) function of the locally D-optimal 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

senslocally(
  formula,
  predvars,
  parvars,
  family = gaussian(),
  x,
  w,
  lx,
  ux,
  inipars,
  fimfunc = NULL,
  sens.control = list(),
  calculate_criterion = TRUE,
  plot_3d = c("lattice", "rgl"),
  plot_sens = TRUE,
  npar = length(inipars),
  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.
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.
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 length(inipars).
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 θ_0 denotes the vector of initial estimates for the unknown parameters. A design ξ* is locally D-optimal among all designs on χ if and only if the following inequality holds for all x belong to χ

c(x, ξ*, θ_0) = tr M^-1(ξ*, θ0)I(x, θ_0)-p <= 0,

with equality at all support points of ξ*. Here, p is the number of model parameters. c(x,ξ*, θ_0) is called sensitivity or derivative function.

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

ELB = p/(p + ε),

where p is the number of model parameters. Obviously, calculating ELB requires finding ε and another optimization problem to be solved. The tuning parameters of this optimization can be regulated via the argument sens.minimax.control. See, for more details, Masoudi et al. (2017).

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 (ε 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.

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.

Examples

############################
# Exponential growth model
############################
# Verifying optimailty of a locally D-optimal design
senslocally(formula = ~a + exp(-b*x),
            predvars = "x", parvars = c("a", "b"),
            x = c(.1, 1), w = c(.5, .5),
            lx = 0, ux = 1, inipars = c(1, 10))


##############################
# A model with two predictors
##############################
x0 <- c(30, 3.861406, 30, 4.600633, 0, 0, 5.111376, 4.168798)
w0 <- rep(.25, 4)
senslocally(formula =  ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),
            predvars = c("S", "I"),
            parvars = c("V", "Km", "Kic", "Kiu"),
            x = x0, w = w0,
            lx = c(0, 0), ux = c(30, 60),
            inipars = c(1.5, 5.2, 3.4, 5.6))
## Not run: 
  # using package rgl for 3d plot:
  res<- senslocally(formula =  ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),
                    predvars = c("S", "I"),
                    parvars = c("V", "Km", "Kic", "Kiu"),
                    x = x0, w = w0,
                    lx = c(0, 0), ux = c(30, 60),
                    inipars = c(1.5, 5.2, 3.4, 5.6),
                    plot_3d = "rgl")


## End(Not run)

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

senslocally(formula = ~1/(1 + exp(-b * (x-a))), predvars = "x",
            parvars = c("a", "b"), family = "binomial",
            inipars = c(0, 1.5),
            crtfunc = Aopt,
            lx = -2, ux = 2,
            sensfunc = Aopt_sens,
            x = c(-1,  1), w = c(.5, .5))
# not optimal

ICAOD documentation built on Oct. 23, 2020, 6:40 p.m.