######################################################################################################*
######################################################################################################*
#'@title Calculates Relative Efficiency for Locally Optimal Designs
#' @description
#' Given a vector of initial estimates for the parameters, this function calculates the D-and PA- efficiency of a design \eqn{\xi_1} with respect to a design \eqn{\xi_2}.
#' Usually, \eqn{\xi_2} is an optimal design.
#'
#'
#' @details For a known \eqn{\theta_0}, relative D-efficiency is
#' \deqn{exp(\frac{log|M(\xi_1, \theta_0)| - log|M(\xi_2, \theta_0)|}{npar})}{
#' exp\{(log|M(\xi_1, \theta_0)| - log|M(\xi_2, \theta_0)|)/npar\}.}
#' The relative P-efficiency is
#' \deqn{\exp(\log(\sum_{i=1}^k w_{1i}p(x_{1i}, \theta_0) - \log(\sum_{i=1}^k w_{2i}p(x{2_i}, \theta_0))}{
#' exp(log (\sum w1_i p(x1_i, \theta_0) - log(\sum w2_i p(x2_i, \theta_0)),
#' }
#' where \eqn{x_2}{x1} and \eqn{w_2}{w1} are usually the support points and the corresponding weights of the optimal design, respectively.
#'
#' The argument \code{x1} is the vector of design points.
#' For design points with more than one dimension (the models with more than one predictors),
#' it is a concatenation of the design points, but \strong{dimension-wise}.
#' For example, let the model has three predictors \eqn{(I, S, Z)}.
#' Then, a two-point optimal design has the following points:
#' \eqn{\{\mbox{point1} = (I_1, S_1, Z_1), \mbox{point2} = (I_2, S_2, Z_2)\}}{{point1 = (I1, S1, Z1), point2 = (I2, S2, Z2)}}.
#' Then, the argument \code{x1} is equal to
#' \code{x = c(I1, I2, S1, S2, Z1, Z2)}.
#'
#' @export
#' @inheritParams locallycomp
#' @param x1 Vector of design (support) points of \eqn{\xi_1}. See 'Details' of \code{\link{leff}}.
#' @param w1 Vector of corresponding design weights for \code{x1}.
#' @param x2 Vector of design (support) points of the optimal design (\eqn{\xi_2}). Similar to \code{x1}.
#' @param w2 Vector of corresponding design weights for \code{x2}.
#' @param type A character. \code{"D"} denotes the D-efficiency and \code{"PA"} denotes the average P-efficiency.
#' @return A value between 0 and 1.
#' @references McGree, J. M., Eccleston, J. A., and Duffull, S. B. (2008). Compound optimal design criteria for nonlinear models. Journal of Biopharmaceutical Statistics, 18(4), 646-661.
#' @example inst/examples/leff_examples.R
leff <- function(formula,
predvars,
parvars,
family = gaussian(),
inipars,
type = c("D", "PA"),
fimfunc = NULL,
x2, w2, x1, w1,
npar = length(inipars),
prob = NULL){
## bayesian D-efficiency
### relative efficieny of x with respect to x2
if (!(type[1] %in% c("D", "PA")))
stop("'type' must be 'D' or 'PA'")
npred <- length(x1)/length(w1)
fimfunc_formula <- check_common_args(fimfunc = fimfunc, formula = formula,
predvars = predvars, parvars = parvars,
family = family, lx =rep(0, npred), ux = rep(1, npred), iter = 1, k = 1,
paramvectorized = FALSE,
prior = NULL, x = NULL,
user_crtfunc = NULL,
user_sensfunc = NULL)
if (!missing(formula)){
if (length(inipars) != length(parvars))
stop("length of 'inipars' is not equal to the length of 'parvars'")
}
if(missing(formula)){
# to handle ...
fimfunc2 <- function(x, w, param)
fimfunc(x = x, w = w, param = param)
} else{
fimfunc2 <- fimfunc_formula$fimfunc_formula ## can be vectorized with respect to parameters!
}
if (type[1] == "D")
releff <- (det2(fimfunc2(x = x1, w = w1, param = inipars))/
det2(fimfunc2(x = x2, w = w2, param = inipars)))^(1/npar)
if (type[1] == "PA"){
if (is.null(prob))
stop("'prob' must be given for P-optimality")
if (is.formula(prob)){
prob <- create_prob(prob = prob, predvars = predvars, parvars = parvars)
}else{
if (!is.function(prob))
stop("'prob' must be either a function or a formula")
if (!formalArgs(prob) %in% c("x", "param"))
stop("arguments of 'prob' must be 'x' and 'param'")
}
releff <- sum(w1 * prob(x = x1, param = inipars))/sum(w2 * prob(x = x2, param = inipars))
}
return(releff)
}
######################################################################################################*
######################################################################################################*
#' @title Calculates Relative Efficiency for Bayesian Optimal Designs
#' @description
#' Given a prior distribution for the parameters, this function calculates the Bayesian D-and PA- efficiency of a design \eqn{\xi_1} with respect to a design \eqn{\xi_2}.
#' Usually, \eqn{\xi_2} is an optimal design.
#' This function is especially useful for investigating the robustness of Bayesian optimal designs under different prior distributions (See 'Examples').
#'
#' @inheritParams leff
#' @inheritParams bayes
#' @details
#' See Masoudi et al. (2018) for formula details (the paper is under review and will be updated as soon as accepted).
#'
#' The argument \code{x1} is the vector of design points.
#' For design points with more than one dimension (the models with more than one predictors),
#' it is a concatenation of the design points, but \strong{dimension-wise}.
#' For example, let the model has three predictors \eqn{(I, S, Z)}.
#' Then, a two-point optimal design has the following points:
#' \eqn{\{\mbox{point1} = (I_1, S_1, Z_1), \mbox{point2} = (I_2, S_2, Z_2)\}}{{point1 = (I1, S1, Z1), point2 = (I2, S2, Z2)}}.
#' Then, the argument \code{x} is equal to
#' \code{x = c(I1, I2, S1, S2, Z1, Z2)}.
#'
#'
#' @export
#' @example inst/examples/beff_examples.R
beff <- function(formula,
predvars,
parvars,
family = gaussian(),
prior,
fimfunc = NULL,
x2, w2, x1, w1,
crt.bayes.control = list(),
npar = NULL,
type = c("D", "PA"),
prob = NULL){
## bayesian D-efficiency
### relative efficieny of x with respect to x2
if (!(type[1] %in% c("D", "PA")))
stop("'type' must be 'D' or 'PA'")
if (is.null(npar)){
if (!missing(formula))
npar <- length(parvars)
else
npar <- prior$npar
}
## only to pass the check_common_eargs
npred <- length(x1)/length(w1)
fimfunc_formula <- check_common_args(fimfunc = fimfunc, formula = formula,
predvars = predvars, parvars = parvars,
family = family, lx =rep(0, npred), ux = rep(1,npred),
iter = 1, k = length(w1),
paramvectorized = FALSE, prior = prior,
x = NULL, user_crtfunc = NULL,
user_sensfunc = NULL)
if(missing(formula)){
# to handle ...
fimfunc2 <- function(x, w, param){
fimfunc(x = x, w = w, param = param)
#fimfunc(x = x, w = w, param = param,...)
}
} else{
#if (length(prior$lower) != length(parvars))
if (length(prior$lower) != fimfunc_formula$num_unknown_param)
stop("length of 'prior$lower' is not equal to the number of unknown (not fixed) parameters")
# fim_localdes <- fimfunc_formula$fimfunc_formula
fimfunc2 <- fimfunc_formula$fimfunc_formula ## can be vectorized with respect to parameters!
}
# do not change its position
if (!is.null(crt.bayes.control$method))
if (crt.bayes.control$method == "quadrature")
warning("The 'quadrature' method is not available for 'beff'. 'method' will be coerced to 'cubature'")
crt.bayes.control <- do.call("crt.bayes.control", crt.bayes.control)
if (type[1] == "D")
cr_integrand <- function(param, x, w){
# bcrfunc1 <- apply(param, 2,
# FUN = function(col_par)-det2(fimfunc2(x = x, w = w, param = col_par), logarithm = TRUE)) * prior$fn(t(param))
bcrfunc1 <- apply(param, 2,
FUN = function(col_par)-det2(fimfunc2(x = x, w = w, param = col_par), logarithm = TRUE)) * prior$fn(t(param))
dim(bcrfunc1) <- c(1, length(bcrfunc1))
return(bcrfunc1)
}
if (type[1] == "PA"){
if (is.formula(prob)){
prob <- create_prob(prob = prob, predvars = predvars, parvars = parvars)
}else{
if (!is.function(prob))
stop("'prob' must be either a function or a formula when type is 'PA'")
if (!formalArgs(prob) %in% c("x", "param"))
stop("arguments of 'prob' must be 'x' and 'param'")
}
cr_integrand <- function(param, x, w){
bcrfunc1 <- (apply(param, 2, function(col_par) -log(sum(w * prob(x = x, param = col_par))))) * prior$fn(t(param))
dim(bcrfunc1) <- c(1, length(bcrfunc1))
return(bcrfunc1)
}
}
crfunc_bayesian <- function(x, w, maxEval, tol) {
out <- hcubature(f = cr_integrand, lowerLimit = prior$lower,
upperLimit = prior$upper,
vectorInterface = TRUE,
x = x, w = w, tol = tol, maxEval = maxEval)
val <- out$integral
return(list(val = val, fneval = out$functionEvaluations))
}
if (type[1] == "D")
releff<- exp((crfunc_bayesian(x = x2, w = w2, maxEval = crt.bayes.control$cubature$maxEval, tol = crt.bayes.control$cubature$tol)$val-
crfunc_bayesian(x = x1, w = w1, maxEval = crt.bayes.control$cubature$maxEval, tol = crt.bayes.control$cubature$tol)$val)/npar)
# if (type[1] == "PA")
# releff<- (crfunc_bayesian(x = x, w = w, maxEval = crt.bayes.control$cubature$maxEval, tol = crt.bayes.control$cubature$tol)$val/
# crfunc_bayesian(x = x2, w = w2, maxEval = crt.bayes.control$cubature$maxEval, tol = crt.bayes.control$cubature$tol)$val)
if (type[1] == "PA")
releff <- exp(crfunc_bayesian(x = x2, w = w2, maxEval = crt.bayes.control$cubature$maxEval, tol = crt.bayes.control$cubature$tol)$val - crfunc_bayesian(x = x1, w = w1, maxEval = crt.bayes.control$cubature$maxEval, tol = crt.bayes.control$cubature$tol)$val)
# releff<- crfunc_bayesian_D(x = x2, w = w2, maxEval = crt.bayes.control$cubature$maxEval, tol = crt.bayes.control$cubature$tol)$val/
# crfunc_bayesian_D(x = x, w = w, maxEval = crt.bayes.control$cubature$maxEval, tol = crt.bayes.control$cubature$tol)$val
return(releff)
}
######################################################################################################*
######################################################################################################*
######################################################################################################*
######################################################################################################*
#'@title Calculates Relative Efficiency for Minimax Optimal Designs
#' @description
#' Given a parameter space for the unknown parameters, this function calculates the D-efficiency of a design \eqn{\xi_1} with respect to a design \eqn{\xi_2}.
#' Usually, \eqn{\xi_2} is an optimal design.
#'
#'
#' @details
#' See Masoudi et al. (2017) for formula details.
#'
#' The argument \code{x1} is the vector of design points.
#' For design points with more than one dimension (the models with more than one predictors),
#' it is a concatenation of the design points, but \strong{dimension-wise}.
#' For example, let the model has three predictors \eqn{(I, S, Z)}.
#' Then, a two-point optimal design has the following points:
#' \eqn{\{\mbox{point1} = (I_1, S_1, Z_1), \mbox{point2} = (I_2, S_2, Z_2)\}}{{point1 = (I1, S1, Z1), point2 = (I2, S2, Z2)}}.
#' Then, the argument \code{x} is equal to
#' \code{x = c(I1, I2, S1, S2, Z1, Z2)}.
#'
#' @export
#' @inheritParams minimax
#' @param x1 Vector of design (support) points of \eqn{\xi_1}. See 'Details' of \code{\link{leff}}.
#' @param w1 Vector of corresponding design weights for \code{x}.
#' @param x2 Vector of design (support) points of the optimal design (\eqn{\xi_2}). Similar to \code{x}.
#' @param w2 Vector of corresponding design weights for \code{x2}.
#' @return A value between 0 and 1.
#' @example inst/examples/meff_examples.R
meff <- function(formula,
predvars,
parvars,
family = gaussian(),
lp, up,
fimfunc = NULL,
x2, w2, x1, w1,
standardized = FALSE,
localdes = NULL,
crt.minimax.control = list(),
npar = length(lp)){
## minimax D-efficiency
### relative efficieny of x with respect to x2
if (!is.logical(standardized))
stop("'standardized' must be logical")
npred <- length(x1)/length(w1)
if (standardized){
if (!missing(formula))
localdes_par <- create_localdes(parvars = parvars, localdes = localdes) else
localdes_par <- localdes
}else
localdes_par <- NULL
fimfunc_formula <- check_common_args(fimfunc = fimfunc, formula = formula,
predvars = predvars, parvars = parvars,
family = family, lx =rep(0, npred), ux = rep(1, npred),
iter = 1, k = 1,
paramvectorized = FALSE,
prior = NULL, x = NULL,
user_crtfunc = NULL,
user_sensfunc = NULL)
if (!missing(formula)){
if (length(lp) != length(parvars))
stop("length of 'lp' is not equal to the length of 'parvars'")
}
if(missing(formula)){
# to handle ...
fimfunc2 <- function(x, w, param)
fimfunc(x = x, w = w, param = param)
#fimfunc(x = x, w = w, param = param,...)
} else{
fimfunc2 <- fimfunc_formula$fimfunc_formula ## can be vectorized with respect to parameters!
}
crt.minimax.control <- do.call("crt.minimax.control", crt.minimax.control)
if (!standardized){
crfunc <- function(x, w, param){
minimax_crfunc <- -det2(fimfunc2(x = x, w = w, param = param), logarithm = TRUE) + 5000 * (sum(w) - 1) ^ 2
return(minimax_crfunc)}
}else{
crfunc <- function(x, w, param) {
localdes_res <- localdes_par(param = param)
denominator <- det2(fimfunc2(x = localdes_res$x, w = localdes_res$w, param = param), logarithm = FALSE)
numerator <- det2(fimfunc2(x = x, w = w, param = param), logarithm = FALSE)
eff <- (numerator/denominator)
if (is.nan(eff))
stop("The criterion (D-efficiency) value is 'NaN' for ", paste("iniparam = c(", paste(round(param, 5), collapse = ","), ").",sep = ""),
"\nCheck the Fisher information matrix, number of design points, lx, ux and .... for possible unmatched design parameters.")
if (round(eff, 7) > 1)
stop("Efficiency value ", eff,
" is greater than one. This results in a wrong conclusion that the non-optimal design is better than the true optimal design when ",
paste("iniparam = c(", paste(round(param, 5), collapse = ","), ").",sep = ""),
"\n Probably 'localdes' does not return the true locally optimal designfor iniparam.")
if (npar %% 2 != 0) {
eff <- (eff)^(1/npar)
}else{
eff <- ifelse(eff < 0, 0,(eff)^(1/npar))
}
return(-eff)
}
}
any_fixed <- sapply(1:length(lp), function(i) lp [i] == up[i]) # is a vector
if (any(any_fixed)){
is_fixed <- TRUE
fixedpar_id <- which(any_fixed)
fixedpar <- lp[fixedpar_id]
lp_nofixed <- lp[-fixedpar_id]
up_nofixed <- up[-fixedpar_id]
## warning: lp and up are channged here if 'any_fix == TRUE'
}else{
fixedpar <- NA
fixedpar_id <- NA
is_fixed <- FALSE
lp_nofixed <- lp
up_nofixed <- up
}
if (is_fixed){
crfunc2 <- function(param, x, w, fixedpar = NA, fixedpar_id = NA){
param1 <- insert_fixed(param = param, fixedpar = fixedpar, fixedpar_id = fixedpar_id)
out <- crfunc(param = param1, x = x, w = w)
return(out)
}
}else{
crfunc2 <- function(param, x, w, fixedpar = NA, fixedpar_id = NA){
out <- crfunc(param = param, x = x, w = w)
return(out)
}
}
vertices_inner <- make_vertices(lower = lp_nofixed, upper = up_nofixed)
optim_nloptr <- function(fn, lower, upper, w, x, fixedpar, fixedpar_id, standardized){
if (is.null(crt.minimax.control$x0))
x0 <- (lower + upper)/2 else
x0 <- crt.minimax.control$x0
out_nloptr <- nloptr::nloptr(x0= x0, eval_f = fn, lb = lower, ub = upper,
opts = crt.minimax.control$optslist,
x = x, w = w, fixedpar = fixedpar,
fixedpar_id = fixedpar_id)
out <- find_on_points(fn = fn,
points = vertices_inner,
x = x, w = w,
fixedpar = fixedpar,
fixedpar_id = fixedpar_id)
minima <- out$minima
minima_nloptr <- c(out_nloptr$solution, out_nloptr$objective)
#minima_nloptr <- c(out_nloptr$par, out_nloptr$value)
minima <- rbind(minima, minima_nloptr)
if (standardized)
crt_val <- -min(minima[, dim(minima)[2]]) else
crt_val <- max(minima[, dim(minima)[2]])
return(crt_val)
}
crt_x1 <- optim_nloptr(fn = crfunc2, lower = lp_nofixed,
upper = up_nofixed, x = x1, w = w1,
fixedpar_id = fixedpar_id, fixedpar = fixedpar,
standardized = standardized)
# optimal
crt_x2 <- optim_nloptr(fn = crfunc2, lower = lp_nofixed,
upper = up_nofixed, x = x2, w = w2,
fixedpar_id = fixedpar_id, fixedpar = fixedpar,
standardized = standardized)
# (crt_x1/crt_x2)^(1/npar) when log = FALSE
#0.011237
if (!standardized)
releff <- exp((crt_x2 - crt_x1)/npar) else
releff <- crt_x1/crt_x2
return(releff)
}
######################################################################################################*
######################################################################################################*
######################################################################################################*
######################################################################################################*
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.