equivalence: Checking optimality of a design with respect to locally D-optimal, minimax D-optimal and standardized maximin D-optimal criteria

Description

An approximate design $\xi$ is a probability measure defined on a user-selected design space $\chi$. Let $\Xi$ be the space of all such designs on $\chi$ and let $\xi$ be an approximate design with $k$ support points at $\boldsymbol{x}_1, \boldsymbol{x}_2, ..., \boldsymbol{x}_k$ from $\chi$ with corresponding weights $w_1, . . . ,w_k$, $\sum_{i=1}^{k} w_i = 1$. A design $\xi^*$ is minimax D-optimal among all designs on $\chi$ if and only if there exists a probability measure $\mu^*$ on $$A(\xi^*) = \left\{\nu \in \Theta \mid -log|M(\xi^*, \nu)| = \max_{\theta \in \Theta} -log|M(\xi^*, \theta)| \right\},$$ such that the following inequality holds for all $\boldsymbol{x} \in \chi$ $$c(\boldsymbol{x}, \mu^*, \xi^*) = \int_{A(\xi^*)} tr M^{-1}(\xi^*, \nu)I(\boldsymbol{x}, \nu)\mu^* d(\nu)-p \leq 0,$$ with equality at all support points of $\xi^*$. Here, $p$ is the number of model parameters. $c(\boldsymbol{x}, \mu^*, \xi^*)$ is called sensitivity function. The set $A(\xi^*)$ is sometimes called the answering set of $\xi^*$ and the measure $\mu^*$ is a subgradient of the non-differentiable criterion evaluated at $M(\xi^*,\nu)$. For standardized maximin D-optimal design the answering set $N(\xi^*)$ is $$N(\xi^*) = \left\{\boldsymbol{\nu} \in \Theta \mid \mbox{eff}_D(\xi^*, \boldsymbol{\nu}) = \min_{\boldsymbol{\theta} \in \Theta} \mbox{eff}_D(\xi^*, \boldsymbol{\theta}) \right\}. $$ where $\mbox{eff}_D(\xi, \boldsymbol{\theta}) = (\frac{|M(\xi, \boldsymbol{\theta})|}{|M(\xi_{\boldsymbol{\theta}}, \boldsymbol{\theta})|})^\frac{1}{p}$ and $\xi_\theta$ is the locally D-optimal design with respect to $\theta$. For locally optimal designs the answering set has only one element that is $\nu = \boldsymbol{\theta_0} $ and $\mu =1$.

Usage

equivalence(fimfunc, x, w, lx, ux, lp, up, type, n.seg = 6,
  maxeval_equivalence = 6000, locally = NULL, control_gosolnp = NULL,
  plot_3d = "lattice", plot_sensitivity = TRUE, ...)

Arguments

fimfunc

Fisher information matrix. Can be the name of the Fisher information matrix from FIM family functions available in this package as a character string or a function that returns the information matrix. See "Details" of mica.

a vector of design points. If the model has $n$ explanatory variables, let $x_{ij}$ be the $j$th component of the $$i$th design point. The argument x is $x = (x_{11}, x_{21},..., x_{k1},..., x_{1n}, x_{2n},... x_{kn})$. See "Examples" on how to set this argument when the design space does not have one dimension, e.g. is of two-dimension.

a vector of design weights

lower bound of the design space $\chi$

upper bound of the design space $\chi$

lower bound of the region of uncertainty $\Theta$. As same order as the argument param of fimfunc.

upper bound of the region of uncertainty $\Theta$. For the locally D-optimal design,lp and up must be fixed to the same values, i.e. lp = up.

type

a character string that shows the type of optimal design. "minimax" for a minimax optimal design, "standardized" for a standardized maximin D-optimal design and "locally" for a locally D-optimal design. When type = "locally", then lp must be set equal to up.

n.seg

number of starting points to find all local optima of the inner problem for constructing the answering set $A(\xi)$ or $N(\xi)$ is equal to (n.seg + 1)^p. Applicable only when type = "standardized" or type = "minimax". Defaults to 4. See "Details".

maxeval_equivalence

maximum number of evaluations (maxeval) that will be passed to optimization function directL to find the maximum of the sensitivity function required for calculating ELB. See "Details".

locally

a function that returns the value of determinant of FIM for the locally D-optimal design, i.e. $|M(\xi_{\bold{\theta}}, \bold{\theta})|$. Only required when type = "standardized". See "Details".

control_gosolnp

tuning parameters of the function gosolnp. Only when 1) type = "standardized" 2) fimfunc is a character string 3) locally D-optimal design is not available in a closed-form. See "Details" for the components.

plot_3d

a character strings shows which plotting package should be used to plot the two-dimensional sensitivity function. plot_3d = "lattice" to use package lattice and plot_3d = "rgl" to use package rgl. The package should be installed before use.

plot_sensitivity

logical, if TRUE, the sensitivity function will be plotted.

...

further argument to be passed to fimfunc.

Value

an object of class 'equivalence' that is a list contains:

type: argument type that is required for print methods.
all_optima: a matrix; all the local maxima of the minimax D-criterion for minimax optimal design problems or all the local minima of the standardized maximin D-criterion standardized maximin optimal design problems. Each row is one element of a set of parameters. For locally optimal designs it will be set to NA.
all_optima_cost: cost of each element of all_optima. NA for locally optimal designs.
answering: a matrix; answering set $A(\xi)$ or $N(\xi)$ chosen from all_optima based on a pre-chosen tolerance, say $0.005$. Each row is one element of the set. NA for locally optimal designs.
answering_cost: cost of each element of answering set. NA for locally optimal designs.
mu: the probability measure on the answering set. It is NA for a locally optimal design.
max_deriv: maximum of the sensitivity function
ELB: D-efficiency lower bound. If negative, probably a local maximum has been found. In this case, the value of maxeval_equivalence should be increased to find the global maximum.
crtval: criterion value. In minimax optimal design it is equal to the minimum of answering_cost and for the standardized maximin optimal design it is equal to the maximum of answering_cost. For the locally D-optimal design it is the value of the logarithm of the determinant of the FIM. See "Details".

Details

There is no theoretical rule on how to choose the number of points in $A(\xi^*)$ as support for the measure $\mu^*$ and they would have to be found by trial and error. To this end, we first find all the local maxima for optimization over $\Theta$ by a local search (L-BFGS-B) with different (n.seg + 1)^p starting points and then pick the ones nearest to the global minimum subject to a tolerance of $0.005$.

We measure the closeness of a design $\xi$ to the minimax optimal design using its minimax D-efficiency defined by $$d_{\mbox{eff}} = \left(\frac{\max_{\boldsymbol{\theta} \in \Theta} -\log|M(\xi_D, \boldsymbol{\theta})|}{\max_{\boldsymbol{\theta} \in \Theta} -\log|M(\xi, \boldsymbol{\theta})|}\right)^{1/p},$$ where $\xi_D$ is the minimax D-optimal design. We may obtain a D-efficiency Lower Bound (ELB) for the minimax D-efficiency of a design $\xi$ without knowing $\xi^*$. ELB is calculated by $p/(p + max_{\boldsymbol{x} \in \chi}c(\boldsymbol{x}, \mu, \xi))$, where $\mu^*$ is the probability measure defined on $A(\xi)$ that maximizes $c(x_\xi,\mu,\xi)$ over all probability measures $\mu$ and $x_\xi$ is any arbitrary support point of $\xi$. Using similar arguments, we can find the D-efficiency lower bound for the standardized maximin D-optimal and locally D-optimal designs.

For standardized maximin D-optimal designs, if fimfunc is a character strings from the list of available FIM, then the function locally is created automatically as follows:

if a closed-form is available for the locally D-optimal designs for the model of interest: function locally combines some algebraic functions to return the value of the determinant. See "Details" of each defined FIM for the formula.
if no closed-form is available for the locally D-optimal designs for the model of interest: function gosolnp is used to find the locally D-optimal design (within the class of minimally-supported and equally-weighted designs). Argument control_gosolnp is a list that contains some tuning parameters of the function gosolnp that are n.sim (default 500), n.restarts (default 1) and trace (default FALSE).

We note that if fimfunc is defined by the user, defining locally is obligatory. In this case function locally must have two arguments: 1) model parameters param 2) auxiliary that is for internal use.

To calculate the ELB, it is necessary to find the global maximum of the sensitivity function. if a local maximum is found then a negative value for ELB may be reported. In this case, increasing the value of maxeval_equivalence can fix the computational issue.

The criterion value for locally D-optimal design is $$-\log|M(\xi, \boldsymbol{\theta_0} )|;$$ for minimax optimal design is (global maximum) $$\max_{\theta \in \Theta} -\log|M(\xi, \theta)|;$$ for standardized maximin optimal design is (global minimum) $$\inf_{\theta \in \Theta} \left[\left(\frac{|M(\xi, \theta)|}{|M(\xi_{\theta}, \theta)|}\right)^\frac{1}{p}\right].$$

References

Atwood, C. L. (1969). Optimal and efficient designs of experiments. The Annals of Mathematical Statistics, 1570-1602.

Examples

Run this code

# NOT RUN {
#############################################################
## check locally optimality: lp = up and type = "locally"
inipar <- c(2, 3)
equivalence (fimfunc = "FIM_logistic",
            x = c(1.485526, 2.51446 ),
            w = c(.5, .5),
            lx = -5, ux = 5,
            lp = inipar, up = inipar,
            type = "locally")

##############################################################################
## standardized maximin D-optimal design does not depend on theta0 and theta1,
##  so we fix them locally D-optimal design has a closed-form which is defined
##  internally
equivalence (fimfunc = "FIM_loglin",
            x = c(0, 4.2494, 17.0324, 149.9090),
            w = c(0.3204, 0.1207, 0.2293, 0.3296),
            lx = 0, ux = 150,
            lp = c(2, 2, 1), up = c(2, 2, 15),
            type = "standardized")

##############################################################################
# }
# NOT RUN {
## there is no analytical solution for locally optimal design for this model
## gosolnp automatically will be used to find the locally
##  optimal design in the denominator of standardized criterion.
##  Becasue, it is two-level nested optimization
##  (first level on parameter space) and second level on design space)
##  it takes so long to find 'all_optima' and construct 'answerign' set.
equivalence (fimfunc = "FIM_power_logistic",
            x = c(-4.5515, 0.2130, 2.8075),
            w = c(0.4100, 0.3723, 0.2177),
            lx = -5, ux = 5,
            lp = c(0, 1), up = c(3, 1.5),
            type = "standardized",
            s = .2)
# }
# NOT RUN {

############################################################################
### when a design point is of two dimension
# }
# NOT RUN {
equivalence (fimfunc = "FIM_mixed_inhibition",
            x = c(3.4614, 4.2801, 30, 30, 0, 3.1426, 0, 4.0373 ),
            w = rep(1/4, 4),
            lx = c(0, 0), ux = c(30, 60),
            lp = c(7, 4, 2, 4), up = c(7, 5, 3, 5),
            type = "standardized")
## here the design points are x1 = c(3.4614, 0), x2 = c(4.2801, 3.1426),
## x3 = c(30, 0), x4 = c(30, 4.0373)
## using package rgl (rgl must be installed for plot)
equivalence (fimfunc = "FIM_mixed_inhibition",
            x = c(3.4614, 4.2801, 30, 30, 0, 3.1426, 0, 4.0373 ),
            w = rep(1/4, 4),
            lx = c(0, 0), ux = c(30, 60),
            lp = c(7, 4, 2, 4), up = c(7, 5, 3, 5),
            type = "standardized", plot_3d = "rgl")

equivalence (fimfunc = "FIM_comp_inhibition",
            x = c(3.4432, 30, 30, 0, 0, 18.8954),
            w = rep(1/3, 3),
            lx = c(0, 0), ux = c(30, 60),
            lp = c(7, 4, 2), up = c(7, 5, 3),
            type = "standardized")
# }
# NOT RUN {
##########################################################################
##########################################################################
## defining function 'locally'
locally_det<- function(param,  auxiliary){
 ## param is the vector of theta = (theta0, theta1, theta2)
 ux <- 0
 lx <- 150
 xstar <- (ux + param[3]) * (lx + param[3]) * (log(ux + param[3]) -
  log(lx + param[3]))/(ux - lx) - param[3]
 denominator <- det(FIM_loglin(x = c(lx, xstar, ux) , w = rep(1/3, 3), param = param))
 return(denominator)
}
equivalence (fimfunc = "FIM_loglin",
            x = c(0, 4.2494, 17.0324, 149.9090),
            w = c(0.3204, 0.1207, 0.2293, 0.3296),
            lx = 0, ux = 150,
            lp = c(2, 2, 1), up = c(2, 2, 15),
            locally = locally_det,
            type = "standardized")
# }

Run the code above in your browser using DataLab