senslocally: Verifying Optimality of The Locally D-optimal Designs

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.

Vector of the design (support) points. See 'Details' of sensminimax for models with more than one predictors.

Vector of the corresponding design weights for x.

Vector of lower bounds for the predictors. Should be in the same order as predvars.

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.

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

Details

Let $\theta_0$ denotes the vector of initial estimates for the unknown parameters. A design $\xi^*$ is locally D-optimal among all designs on $\chi$ if and only if the following inequality holds for all $\boldsymbol{x} \in \chi$ $$c(\boldsymbol{x}, \xi^*, \theta_0) = tr M^{-1}(\xi^*, \theta_0)I(\boldsymbol{x}, \theta_0)-p \leq 0,$$ with equality at all support points of $\xi^*$. Here, $p$ is the number of model parameters. $c(\boldsymbol{x},\xi^*, \theta_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 $\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. See, for more details, Masoudi et al. (2017).

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

Run this code

# NOT RUN {
############################
# 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")

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

Run the code above in your browser using DataLab