# Robust D-Optimal Designs Finds Robust designs or optimal in-average designs for linear and nonlinear models. It is useful when a set of different vectors of initial estimates along with a discrete probability measure are available for the unknown model parameters. It is a discrete version of `bayes`. ```r robust( formula, predvars, parvars, family = gaussian(), lx, ux, iter, k, prob, parset, fimfunc = NULL, ICA.control = list(), sens.control = list(), initial = NULL, npar = dim(parset)[2], plot_3d = c("lattice", "rgl"), x = NULL, 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. - `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`. - `iter`: Maximum number of iterations. - `k`: Number of design points. Must be at least equal to the number of model parameters to avoid singularity of the FIM. - `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`. - `ICA.control`: ICA control parameters. For details, see `ICA.control`. - `sens.control`: Control Parameters for Calculating the ELB. For details, see `sens.control`. - `initial`: A matrix of the initial design points and weights that will be inserted into the initial solutions (countries) of the algorithm. Every row is a design, i.e. a concatenation of `x` and `w`. Will be coerced to a `matrix` if necessary. See 'Details' of `minimax`. - `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]`. - `plot_3d`: Which package should be used to plot the sensitivity (derivative) function for models with two predictors. Either `"rgl"` or `"lattice"` (default). - `x`: Vector of the design (support) points. See 'Details' of `sensminimax` for models with more than one predictors. - `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`. ## Returns an object of class `minimax` that is a list including three sub-lists: - **`arg`**: A list of design and algorithm parameters. - **`evol`**: A list of length equal to the number of iterations that stores the information about the best design (design with least criterion value) of each iteration. `evol[[iter]]` contains: |||| |:--|:--|:--| |`iter`||Iteration number.| |`x`||Design points.| |`w`||Design weights.| |`min_cost`||Value of the criterion for the best imperialist (design).| |`mean_cost`||Mean of the criterion values of all the imperialists.| |`sens`||An object of class `'sensminimax'` . See below.| |`param`||Vector of parameters.| - **`empires`**: A list of all the empires of the last iteration. - **`alg`**: A list with following information: |||| |:--|:--|:--| |`nfeval`||Number of function evaluations. It does not count the function evaluations from checking the general equivalence theorem.| |`nlocal`||Number of successful local searches.| |`nrevol`||Number of successful revolutions.| |`nimprove`||Number of successful movements toward the imperialists in the assimilation step.| |`convergence`||Stopped by `'maxiter'` or `'equivalence'` ?| - **`method`**: A type of optimal designs used. - **`design`**: Design points and weights at the final iteration. - **`out`**: A data frame of design points, weights, value of the criterion for the best imperialist (min_cost), and Mean of the criterion values of all the imperialistsat each iteration (mean_cost). The list `sens` contains information about the design verification by the general equivalence theorem. See `sensminimax` for more details. It is given every `ICA.control$checkfreq` iterations and also the last iteration if `ICA.control$checkfreq >= 0`. Otherwise, `NULL`. `param` is a vector of parameters that is the global minimum of the minimax criterion or the global maximum of the standardized maximin criterion over the parameter space, given the current `x`, `w`. ## Details Let $\Theta$ be a set of initial estimates for the unknown parameters. A robust criterion is evaluated at the elements of $\Theta$ weighted by a probability measure $\pi$ as follows: $$ B(\xi, \pi) = \int_{\Theta}|M(\xi, \theta)|\pi(\theta) d\theta.B(\xi, \Pi) = intergation over \Theta \Psi(\xi, \theta)\pi(\theta) d\theta. $$ A robust design $\xi*$ maximizes $B(\xi, \pi)$ over the space of all designs. When the model is given via `formula`, columns of `parset` must match the parameters introduced in `parvars`. Otherwise, when the model is introduced via `fimfunc`, columns of `parset` must match the argument `param` in `fimfunc`. To verify the optimality of the output design by the general equivalence theorem, the user can either `plot` the results or set `checkfreq` in `ICA.control` to `Inf`. In either way, the function `sensrobust` is called for verification. One can also adjust the tuning parameters in `ICA.control` to set a stopping rule based on the general equivalence theorem. See 'Examples' below. ## Note When a continuous prior distribution for the unknown model parameters is available, use `bayes`. When only one initial estimates of the unknown model parameters is available ($\Theta$ has only one element), use `locally`. ## Examples ```r # Finding a robust design for the two-parameter logistic model # See how we set a stopping rule. # The ELB is computed every checkfreq = 30 iterations # The optimization stops when the ELB is larger than stoptol = .95 res1 <- robust(formula = ~1/(1 + exp(-b *(x - a))), predvars = c("x"), parvars = c("a", "b"), family = binomial(), lx = -5, ux = 5, 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), iter = 1, k =3, ICA.control = list(stop_rule = "equivalence", stoptol = .95, checkfreq = 30)) ## Not run: res1 <- update(res1, 100) # stops at iteration 51 ## End(Not run) ## Not run: res1.1 <- robust(formula = ~1/(1 + exp(-b *(x - a))), predvars = c("x"), parvars = c("a", "b"), family = binomial(), lx = -5, ux = 5, 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(-3, 0, 3), iter = 150, k =3) plot(res1.1) # not optimal ## End(Not run) ################################### # user-defined optimality criterion ################################## # When the model is defined by the formula interface # A-optimal design 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)) } res2 <- robust(formula = ~1/(1 + exp(-b * (x-a))), predvars = "x", parvars = c("a", "b"), family = "binomial", lx = -3, ux = 3, iter = 1, k = 4, crtfunc = Aopt, sensfunc = Aopt_sens, prob = c(.25, .5, .25), parset = matrix(c(-2, 0, 2, 1.25, 1.25, 1.25), 3, 2), ICA.control = list(checkfreq = 50, stoptol = .999, stop_rule = "equivalence", rseed = 1)) ## Not run: res2 <- update(res2, 500) ## End(Not run) # robust c-optimal design # example from Chaloner and Larntz (1989), Figure 3, but robust design c_opt <-function(x, w, a, b, fimfunc){ gam <- log(.95/(1-.95)) M <- fimfunc(x = x, w = w, a = a, b = b) c <- matrix(c(1, -gam * b^(-2)), nrow = 1) B <- t(c) %*% c sum(diag(B %*% solve(M))) } c_sens <- function(xi_x, x, w, a, b, fimfunc){ gam <- log(.95/(1-.95)) M <- fimfunc(x = x, w = w, a = a, b = b) M_inv <- solve(M) M_x <- fimfunc(x = xi_x, w = 1, a = a, b = b) c <- matrix(c(1, -gam * b^(-2)), nrow = 1) B <- t(c) %*% c sum(diag(B %*% M_inv %*% M_x %*% M_inv)) - sum(diag(B %*% M_inv)) } res3 <- robust(formula = ~1/(1 + exp(-b * (x-a))), predvars = "x", parvars = c("a", "b"), family = "binomial", lx = -1, ux = 1, parset = matrix(c(0, 7, .2, 6.5), 2, 2, byrow = TRUE), prob = c(.5, .5), iter = 1, k = 3, crtfunc = c_opt, sensfunc = c_sens, ICA.control = list(rseed = 1, checkfreq = Inf)) ## Not run: res3 <- update(res3, 300) ## End(Not run) ``` ## See Also `bayes` `sensrobust`

robust function