equivalence: Checking minimax, standardized maximin and locally D-optimality of a givne design by equivalence theorem

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 $x1, x2, ..., xk$ from $\chi$ with corresponding weights $w1, ..., wk$, $w1 + w2 + ...+ wk = 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 $x belong to \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(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)$.

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. When design space is multi-dimensional then x should be filled dimension by dimension. See "Examples".

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 unceratinty $\Theta$. The order of the dimension is the same as the order of the parameters in the argument paramof fimfunc.

upper bound of the region of unceratinty $\Theta$. If lp = up, then $\Theta = \theta_0 = lp$ and the design is checked for locally D-optimality, See also argument type.

type

a character strings; 'minimax' for minimax optimal design, 'standardized' for standardized maximin D-optimal design and 'locally' for locally D-optimal design.

n.seg

the number of starting points for finding all local optima to construct the answering set is (n.seg + 1)^p. Only for minimax and standardized maximin optimal designs. See "Details".

maxeval_equivalence

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

locally

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

control_gosolnp

tuning parameters of function gosolnp for models that their locally optimal design do not have an analytical solution and is find by gosolnp. Only required when type is set to 'standardized'. See "Details".

plot_3d

a character strings to show which packages should be used for plotting the sensitivity function when design space is of two dimension; "lattice" to use package lattice and "rgl" to use package rgl. The package should be installed before use.

plot_sensitivity

logical; sensitivity should be plotted? see "Details".

...

further argument to be passed to fimfunc.

Value

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

Details

For locally optimal designs the answering set has only one element that is $\nu = \theta_0$ and $\mu =1$. Thus, n.seg has no further use for locally optimal designs.

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

If $\chi$ is one or two dimensional, one may plot sensitivity function $c(x, \xi*, \mu*)$ versus $x (belongs to \chi)$ and visually inspect whether the graph meets the conditions in the equivalence theorem. If it does, the design $\xi*$ is minimax optimal; otherwise it is not optimal.

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. Using argument similar to Atwood (1969), we obtain a D-efficiency Lower Bound (DLB) for the minimax D-efficiency of a design $\xi$ without knowing $\xi*$. DLB is caluclated by $p/(p + maximum over \chi c(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$.

For standardazed 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 $ eff_D(\xi, \theta) = (|M(\xi, \theta)|/|M(\xi_\theta, \theta)|)^(1/p)$ and $\xi_\theta$ is the locally D-optimal design with respect to $\theta$. We measure the closeness of a design $\xi$ to the standardized maximin optimal design $\xi_D$ using its standardized maximin D-efficiency defined by $$ d_{\mbox{eff}} = \frac{\min_{\boldsymbol{\theta} \in \Theta} \mbox{eff}_D(\xi, \boldsymbol{\theta})}{\min_{\boldsymbol{\theta} \in \Theta} \mbox{eff}_D(\xi_D, \boldsymbol{\theta})}.$$ Similar to the minimax design, we can also find standardized maximin D-efficiency lower bound for the generated standardized maximin design and locally D-optimal design.

For standardized maximin D-optimal designs the function locally is created automatically when locally = NULL:

when fimfunc is a character strings from the available FIM and locally D-optimal design has a closed-form for the chosen model: locally is an algebraic function that returns the value of determinant of the FIM for the locally D-optimal design, i.e. $|M(\xi_\theta, \theta)|$. See "Details" of each defined FIM for the formula.
when fimfunc is a character strings but locally D-optimal design has no closed-form: gosolnp is used to find the locally D-optimal design (within the class of minimally-supported and equally-weighted designs). control_gosolnp is a list of some of the tunning parameters of gosolnp that are: n.sim (default 500), n.restarts (default 1) and trace (default FALSE).
when fimfunc is a user given function: same as the previous case.

User can also provide its own locally. In this case, args(locally) must be function (param, auxiliary), where param is the vector of parameters and auxiliary is an obligatory arguments for internal use when locally = NULL. Please see "Examples" on how to use locally.

For output max_deriv, if the local maximum is found (you can detect it from sensitivity plot) or DLB is negative, the value of maxeval_equivalence should be increased to return the global maximum.

The criterion value for locally D-optimal design is $$-\log|M(\xi, \boldsymbol{\theta_0} )|;$$ for minimax optimal design is $$\max_{\theta \in \Theta} -\log|M(\xi, \theta)|;$$ for standardized maximin optimal design is $$\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

#############################################################
## 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)
# ## End(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")
# ## End(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