Nothing
R/mosallocSTRS.R
In MOSAlloc: Constraint Multiobjective Sample Allocation
Defines functions
mosallocSTRS
Documented in
mosallocSTRS
# Filename: mosallocSTRS.R
#
# Date: 18.07.2025; modified: 10.09.2025, 16.01.2026
# Author: Felix Willems
# Contact: mail.willemsf+MOSAlloc@gmail.com
# (mail[DOT]willemsf+MOSAlloc[AT]gmail[DOT]com)
# Licensing: GPL-3.0-or-later
#
# Please report any bugs or unexpected behavior to
# mail.willemsf+MOSAlloc@gmail.com
# (mail[DOT]willemsf+MOSAlloc[AT]gmail[DOT]com)
#
#---------------------------------------------------------------------------
#
#' @title (Single-stage) stratified random sampling interface for functions
#' \code{mosalloc()} and \code{mosallocStepwiseFirst()}
#'
#' @description Interface for the functions \code{mosalloc()} and
#' \code{mosallocStepwiseFirst()} of the same package which allows for a
#' user-friendly data handling in the case of single-stage stratified random
#' sampling.
#'
#' @param X_var (type: \code{matrix})
#' A matrix containing stratum- (row) and variable- (column) specific
#' precision units (e.g. variances).
#' @param X_tot (type: \code{matrix})
#' A matrix containing stratum- (row) and variable- (column) specific totals.
#' @param N (type: \code{vector})
#' A vector of stratum sizes.
#' @param listD (type: \code{list})
#' A list of lists taking subpopulation- (domain/area) specific
#' arguments. \cr If \code{opts$sense == "max_precision"}, elements are lists
#' containing the following components
#' which correspond to one specific precision target each:
#' \cr \code{..$stratum_id} (type: \code{numeric}) A vector containing the
#' indices of the strata considered for the current restriction. The indices
#' must coincide with the row numbers of \code{X_var} and \code{X_tot}.
#' \cr \code{..$variate} (type: \code{character} or \code{numeric})
#' The column name or column index of \code{X_var} to be addressed.
#' \cr \code{..$measure} (type: \code{character} or \code{numeric})
#' A character ("relVAR" (relative variance) or "VAR" (variance)) or a numerical
#' value between 0 and 1 indicates a gradation coefficient that balances between
#' "relVar" (..$measure = 0) and "VAR" (..$measure = 1).
#' \cr \code{..$name} (type: \code{character}) The name of the subpopulation
#' (domain/area).
#' \cr If \code{opts$sense == "min_cost"}, elements are lists
#' containing the following components
#' which correspond to one specific cost target each:
#' \cr \code{..$stratum_id} (type: \code{numeric})
#' A vector containing the indices of the strata considered for the current
#' objective. The indices must coincide with the row numbers of \code{X_cost}.
#' \cr \code{..$c_type} (type: \code{character} or \code{numeric})
#' The column name or column index of \code{X_cost} to be addressed.
#' \cr \code{..$name} (type: \code{character})
#' The name of the subpopulation (domain/area).
#' @param listA (type: \code{list})
#' A list of lists taking subpopulation- (domain/area) specific
#' arguments. Elements are lists containing the following components
#' which in turn correspond to one specific precision restriction:
#' \cr \code{..$stratum_id} (type: \code{numeric}) A vector containing the
#' indices of the strata considered for the current restriction. The indices
#' must coincide with the corresponding row numbers of \code{X_var} and
#' \code{X_tot}.
#' \cr \code{..$variate} (type: \code{character} or \code{numeric})
#' The column name or column index of \code{X_var} to be addressed.
#' \cr \code{..$measure} (type: \code{character}) "RSE"
#' (relative standard error) or "VAR" (variance).
#' \cr \code{..$bound} (type: \code{numeric})
#' An upper bound to "RSE" (or "VAR").
#' \cr \code{..$name} (type: \code{character}) The name of the subpopulation
#' (domain/area).
#' @param listC (type: \code{list})
#' A list of lists taking subpopulation- (domain/area) specific
#' arguments. Elements are lists containing the following components which
#' correspond to one specific cost restriction:
#' \cr \code{..$stratum_id} (type: \code{numeric})
#' A vector containing the indices of the strata considered for the current
#' restriction. The indices must coincide with the row numbers of \code{X_var}
#' and \code{X_tot}.
#' \cr \code{..$c_coef} (type: \code{numeric}) A vector of length
#' \code{length(stratum_id)} containing the stratum-specific cost components
#' for the set of strata that is going to be bounded by above and/or below.
#' \cr \code{..$c_upper} (type: \code{numeric}) The cost upper bound value.
#' \code{NULL} if not present.
#' \cr \code{..$c_lower} (type: \code{numeric}) The cost lower bound value.
#' \code{NULL} if not present.
#' \cr \code{..$name} (type: \code{character}) The name of the subpopulation
#' (domain/area).
#' @param fpc (type: \code{logical})
#' A \code{TRUE} or \code{FALSE} statement indicating whether the finite
#' population correction should be considered.
#' @param l (type: \code{vector})
#' A vector of lower box constraints.
#' @param u (type: \code{vector})
#' A vector of upper box constraints.
#' @param opts (type: \code{list})
#' The options used by the algorithms:
#' \cr \code{$sense} (type: character) Sense of optimization
#' (default = \code{"max_precision"}, alternative \code{"min_cost"}).
#' \cr \code{$df} (type: function) The gradient of f (default = \code{NULL}).
#' \cr \code{$Hf} (type: function)
#' The Hesse matrix of f (default = \code{NULL}).
#' \cr \code{$method} (type: \code{character})
#' A character indicating scalarization method (default = \code{"WCM"},
#' alternative \code{"WSS"}), \code{$f} must be \code{NULL}.
#' \cr \code{$f} (type: function) Decision functional over the objective
#' vector (default = \code{NULL}). The weighted Chebyshev minimization
#' (\code{WCM}) minimizes for weighted maximum objective component. The
#' weighted sum scalarization (\code{WSS}) minimizes for the weighted sum.
#' For either case the weights are given by \code{$init_w}.
#' \cr \code{$init_w} (type: \code{numeric} or \code{matrix})
#' Preference weightings (default = 1; The weight for first objective component
#' must be 1).
#' \cr \code{$mc_cores} (type: \code{integer})
#' The number of cores for parallelizing
#' multiple input weightings stacked rowwise (default = \code{1L}).
#' \cr \code{$pm_tol} (type: \code{numeric}) The tolerance for the projection
#' method (default = \code{1e-5}).
#' \cr \code{max_iters} (type: \code{integer}) The maximum number of iterations
#' (default = \code{100L}).
#' \cr \code{$print_pm} (type: \code{logical}) A \code{TRUE} or \code{FALSE}
#' statement whether iterations of the projection method should be printed
#' (default = \code{FALSE}).
#' @param ForceOptimality (type: \code{logical})
#' A \code{TRUE} or \code{FALSE} statement indicating whether Pareto optimality
#' should be ensured. Default is \code{FALSE} (for most practical
#' problem formulations Pareto optimality is achieved by construction,
#' e.g. in the case of one overall cost constraint).
#' If \code{TRUE}, additional computation time is required. In this case,
#' \code{mosallocStepwiseFirst()} is used internally.
#' @param X_cost (type: \code{matrix})
#' A matrix containing stratum- (row) and type- (column) specific cost
#' coefficients associated with fixed cost. Types of cost might be, e.g.
#' '$ US', 'minutes', 'sample size', etc. Default is \code{NULL}. The argument
#' is required in the case of cost minimization (opt$sense == "min_cost").
#' @param X_fixed (type: \code{matrix})
#' A matrix containing stratum- (rows) and type- (columns) specific
#' cost coefficients associated with fixed cost. Default is \code{NULL}.
#' Used solely in the case of cost minimization (\code{opt$sense == "min_cost"}).
#'
#' @return A \code{mosaSTRS} object or a list of \code{mosaSTRS}
#' objects. A \code{mosaSTRS} object is a list containing the
#' following components:
#' @returns \code{$sense}
#' Sense of optimization; \code{max precision} or \code{min_cost}.
#' @returns \code{$method}
#' The method used, either weighted sum scalarization (WSS) or weighted
#' Chebyshev minimization (WCM).
#' @returns \code{$init_w}
#' The initial preference weightings (opts$init_w).
#' @returns \code{$opt_w}
#' The optimal weightings w.r.t. \code{opts$init_w} as \code{opts$init_w} might
#' not lead to Pareto optimality. \code{NULL} if \code{ForceOptimality = FALSE}.
#' @returns \code{$n_opt}
#' The vector of optimal sample sizes.
#' @returns \code{$objective}
#' The objective values, including the sensitivity and the RSE.
#' @returns \code{$precision}
#' A data frame corresponding to the precision constraints.
#' @returns \code{$cost}
#' A data frame corresponding to the cost constraints.
#' @returns \code{$problem_components}
#' A list containing the input data to the optimization problem.
#' @returns \code{$output_mosalloc}
#' A list of function returns of \code{mosalloc()}.
#' @returns If \code{opts$init_w} has multiple rows, the function returns a list
#' of \code{mosaSTRS} objects whose length equals the number of rows.
#'
#' @examples
#' # Artificial population of 50 568 business establishments and 5 business
#' # sectors (data from Valliant, R., Dever, J. A., & Kreuter, F. (2013).
#' # Practical tools for designing and weighting survey samples. Springer.
#' # https://doi.org/10.1007/978-1-4614-6449-5, Example 5.2 pages 133-9)
#'
#' # See also https://umd.app.box.com/s/9yvvibu4nz4q6rlw98ac/file/297813512360
#' # file: Code 5.3 constrOptim.example.R
#'
#' Nh <- c(6221, 11738, 4333, 22809, 5467) # stratum sizes
#' ch <- c(120, 80, 80, 90, 150) # stratum-specific cost of surveying
#'
#' # Revenues
#' mh.rev <- c(85, 11, 23, 17, 126) # mean revenue
#' Sh.rev <- c(170.0, 8.8, 23.0, 25.5, 315.0) # standard deviation revenue
#'
#' # Employees
#' mh.emp <- c(511, 21, 70, 32, 157) # mean number of employees
#' Sh.emp <- c(255.50, 5.25, 35.00, 32.00, 471.00) # std. dev. employees
#'
#' # Proportion of estabs claiming research credit
#' ph.rsch <- c(0.8, 0.2, 0.5, 0.3, 0.9)
#'
#' # Proportion of estabs with offshore affiliates
#' ph.offsh <- c(0.06, 0.03, 0.03, 0.21, 0.77)
#'
#' budget <- 300000 # overall available budget
#' n.min <- 100 # minimum stratum-specific sample size
#'
#' # Matrix containing stratum-specific variance components
#' X_var <- cbind(Sh.rev**2,
#' Sh.emp**2,
#' ph.rsch * (1 - ph.rsch) * Nh/(Nh - 1),
#' ph.offsh * (1 - ph.offsh) * Nh/(Nh - 1))
#' colnames(X_var) <- c("rev", "emp", "rsch", "offsh")
#'
#' # Matrix containing stratum-specific totals
#' X_tot <- cbind(mh.rev, mh.emp, ph.rsch, ph.offsh) * Nh
#' colnames(X_tot) <- c("rev", "emp", "rsch", "offsh")
#'
#' # Examples
#' #----------------------------------------------------------------------------
#' # Example 1: Univariate minimization of the variation of estimates for
#' # revenue subject to cost restrictions and precision restrictions to the
#' # relative standard error of estimates for the proportion of businesses with
#' # offshore affiliates. Additionally, there is one overall cost constraint and
#' # at least half of the provided budget must be spend to strata 1 to 3.
#'
#' # Specify objectives via listD
#' listD <- list(list(stratum_id = 1:5, variate = "rev", measure = "relVAR",
#' name = "pop"))
#'
#' # Specify precision constraints via listA
#' listA <- list(list(stratum_id = 1:5, variate = "offsh", measure = "RSE",
#' bound = 0.05, name = "pop"))
#'
#' # Specify cost constraints via listC
#' listC <- list(list(stratum_id = 1:5, c_coef = ch, c_lower = NULL,
#' c_upper = budget, name = "Overall"),
#' list(stratum_id = 1:3, c_coef = ch[1:3],
#' c_lower = 0.5 * budget, c_upper = NULL, name = "1to3"))
#'
#' # Specify stratum-specific box constraints
#' l <- rep(n.min, 5) # minimum sample size per stratum
#' u <- Nh # maximum sample size per stratum
#'
#' # Specify parameter for mosalloc (method = "WSS")
#' opts <- list(sense = "max_precision",
#' f = NULL, df = NULL, Hf = NULL,
#' method = "WSS", init_w = 1,
#' mc_cores = 1L, pm_tol = 1e-05,
#' max_iters = 100L, print_pm = FALSE)
#'
#' # Run mosallocSTRS with weighted sum scalarization (WSS)
#' resWSS <- mosallocSTRS(X_var, X_tot, Nh, listD, listA, listC,
#' fpc = TRUE, l, u, opts)
#'
#' summary(resWSS)
#'
#' # Specify parameter for mosalloc (method = "WCM")
#' opts = list(sense = "max_precision",
#' f = NULL, df = NULL, Hf = NULL,
#' method = "WCM", init_w = 1,
#' mc_cores = 1L, pm_tol = 1e-05,
#' max_iters = 100L, print_pm = FALSE)
#'
#' # Run mosallocSTRS with weighted Chebyshec minimization (WCM)
#' resWCM <- mosallocSTRS(X_var, X_tot, Nh, listD, listA, listC,
#' fpc = TRUE, l, u, opts)
#'
#' summary(resWCM)
#'
#' # The optimal sample sizes vector can also be obtained by
#' summary(resWCM)$n_opt
#'
#' # Hint: For univariate allocation problems 'WSS' and 'WCM' are equivalent!
#'
#' #----------------------------------------------------------------------------
#' # Example 2: Minimization of the maximum relative variation of estimates for
#' # the total revenue, the number of employee, the number of businesses claimed
#' # research credit and the number of businesses with offshore affiliates
#' # subject to one overall cost constraint and at least half of the provided
#' # budget must be spend to strata 1 to 3.
#'
#' # Specify objectives via listD
#' listD <- list(list(stratum_id = 1:5, variate = "rev", measure = "relVAR",
#' name = "pop"),
#' list(stratum_id = 1:5, variate = "emp", measure = "relVAR",
#' name = "pop"),
#' list(stratum_id = 1:5, variate = "rsch", measure = "relVAR",
#' name = "pop"),
#' list(stratum_id = 1:5, variate = "offsh", measure = "relVAR",
#' name = "pop"))
#'
#' # Specify cost constraints via listC
#' listC <- list(list(stratum_id = 1:5, c_coef = ch, c_lower = NULL,
#' c_upper = budget, name = "Overall"),
#' list(stratum_id = 1:3, c_coef = ch[1:3],
#' c_lower = 0.5 * budget, c_upper = NULL, name = "1to3"))
#'
#' # Specify stratum-specific box constraints
#' l <- rep(n.min, 5) # minimum sample size per stratum
#' u <- Nh # maximum sample size per stratum
#'
#' # Specify parameter for mosalloc (method = "WSS")
#' opts = list(sense = "max_precision",
#' f = NULL, df = NULL, Hf = NULL,
#' method = "WSS", init_w = 1,
#' mc_cores = 1L, pm_tol = 1e-05,
#' max_iters = 100L, print_pm = FALSE)
#'
#' # Run mosallocSTRS with weighted sum scalarization (WSS)
#' resWSS <- mosallocSTRS(X_var, X_tot, Nh, listD, NULL, listC,
#' fpc = TRUE, l, u, opts)
#'
#' summary(resWSS)
#'
#' # Specify parameter for mosalloc (method = "WCM")
#' opts = list(sense = "max_precision",
#' f = NULL, df = NULL, Hf = NULL,
#' method = "WCM", init_w = 1,
#' mc_cores = 1L, pm_tol = 1e-05,
#' max_iters = 100L, print_pm = FALSE)
#'
#' # Run mosallocSTRS with weighted Chebyshec minimization (WCM)
#' resWCM <- mosallocSTRS(X_var, X_tot, Nh, listD, NULL, listC,
#' fpc = TRUE, l, u, opts)
#'
#' summary(resWCM)
#'
#' # Since the WCM does not necessarily lead to Pareto optimal allocations,
#' # we might force this via a internal stepwise procedure by setting
#' # ForceOptimality = TRUE.
#'
#' resWCM_FO <- mosallocSTRS(X_var, X_tot, Nh, listD, NULL, listC,
#' fpc = TRUE, l, u, opts, ForceOptimality = TRUE)
#'
#' summary(resWCM_FO)
#'
#' #----------------------------------------------------------------------------
#' # Example 3: Example 2 with multiple sets of preference weightings.
#'
#' # Define a set of preference weightings, e.g.
#' w_1 <- c(1, 1, 1, 1)
#' w_2 <- c(1, 2, 2, 1)
#' w_3 <- c(1, 5, 5, 1)
#'
#' # Combine the weightings to a matrix stacked rowwise
#' w <- rbind(w_1, w_2, w_3)
#'
#' # Specify parameter for mosalloc() (method = "WCM"; not yet possible with WSS)
#' opts = list(sense = "max_precision",
#' f = NULL, df = NULL, Hf = NULL,
#' method = "WCM", init_w = w,
#' mc_cores = 1L, pm_tol = 1e-05,
#' max_iters = 100L, print_pm = FALSE)
#'
#' # Run mosallocSTRS with weighted Chebyshec minimization (WCM)
#' res <- mosallocSTRS(X_var, X_tot, Nh, listD, NULL, listC, fpc = TRUE,
#' l, u, opts)
#'
#' summary(res)
#'
#' #----------------------------------------------------------------------------
#' # Example 4: Minimization of survey cost subject to quality restrictions on
#' # subpopulation level.
#'
#' X_cost <- matrix(ch, 5, 1, dimnames = list(1:5,"cost"))
#'
#' # Specify cost objectives via listD
#' listD <- list(list(stratum_id = 1:5, c_type = "cost", name = "pop"))
#'
#' # Specify quailty restrictions via listD. Here: 5 % relative standard error
#' listA <- list(list(stratum_id = 1:2, variate = "rev", measure = "RSE",
#' bound = 0.05, name = "S1-2"),
#' list(stratum_id = 3:5, variate = "rev", measure = "RSE",
#' bound = 0.05, name = "S3-5"),
#' list(stratum_id = 1:2, variate = "emp", measure = "RSE",
#' bound = 0.05, name = "S1-2"),
#' list(stratum_id = 3:5, variate = "emp", measure = "RSE",
#' bound = 0.05, name = "S3-5"),
#' list(stratum_id = 1:2, variate = "rsch", measure = "RSE",
#' bound = 0.05, name = "S1-2"),
#' list(stratum_id = 3:5, variate = "rsch", measure = "RSE",
#' bound = 0.05, name = "S3-5"),
#' list(stratum_id = 1:2, variate = "offsh", measure = "RSE",
#' bound = 0.05, name = "S1-2"),
#' list(stratum_id = 3:5, variate = "offsh", measure = "RSE",
#' bound = 0.05, name = "S3-5"))
#'
#' # Specify cost constraints
#' listC <- NULL
#'
#' # Specify stratum-specific box constraints
#' l <- rep(n.min, 5) # minimum sample size per stratum
#' u <- Nh # maximum sample size per stratum
#'
#' # Specify parameters for mosalloc()
#' opts = list(sense = "min_cost",
#' f = NULL, df = NULL, Hf = NULL,
#' method = "WCM", init_w = 1,
#' mc_cores = 1L, pm_tol = 1e-05,
#' max_iters = 100L, print_pm = FALSE)
#'
#' # Run mosallocSTRS()
#' res <- mosallocSTRS(X_var, X_tot, Nh, listD, listA, NULL, fpc = TRUE,
#' l, u, opts, X_cost = X_cost)
#'
#' summary(res)
#'
#' # Optimal sample sizes
#' summary(res)$n_opt
#'
#' @importFrom utils capture.output
#'
#' @export
mosallocSTRS <- function(X_var, X_tot, N,
listD, listA = NULL, listC = NULL,
fpc = TRUE, l = 2, u = NULL, opts,
ForceOptimality = FALSE,
X_cost = NULL, X_fixed = NULL) {
message("\n----------------------------------------------------------\n")
# get number of strata in single-stage stratified random sampling (STRS)
H <- length(N)
# expand l and u if necessary
if (length(l) == 1) {
l <- rep(l, H)
}
if (length(u) == 1) {
l <- rep(u, H)
}
# check input
if (opts$sense == "min_cost") {
if (is.null(X_cost)) {
stop("X_cost not specified!")
} else {
if (is.null(X_fixed)) {
X_fixed <- X_cost - X_cost
}
Dc <- constructDobjCostSTRS(X_cost, X_fixed, listD)
}
} else {
if (length(unique(sapply(listD, function(L)L$measure))) != 1) {
message("CAUTION: ..$measures in listD varry among objectives!\n\n")
}
# Construct vector valued precision objective
Dc <- constructDobjPrecisionSTRS(X_var, X_tot, N, listD, fpc)
listD_RSE <- lapply(listD, function(L){L$measure = 0; L})
Dc_RSE <- constructDobjPrecisionSTRS(X_var, X_tot, N, listD_RSE)
}
if (is.null(opts$f)) {
method <- opts$method
} else {
stop("The possiblity to use decision functionals (opts$f) is not yet
implemented! Please set opts$f = NULL.")
}
# Construct precision constraints
if (is.null(listA)) {
Ac <- list(A = NULL, a = NULL)
} else {
Ac <- constructArestrSTRS(X_var, X_tot, N, listA, fpc)
}
# Construct cost constraints
if (is.null(listC)) {
Cc <- list(C = NULL, c = NULL)
} else {
Cc <- constructCrestrSTRS(H, listC)
}
if (is.null(dim(opts$init_w)) || dim(opts$init_w)[1] == 1) {
message(" mosalloc running...\n")
if (method == "WSS") {
# Weight objetives via opts$init_w and solve univariate problem
time_t0 <- proc.time()
if (is.null(dim(opts$init_w))) {
opts$init_w <- matrix(opts$init_w, 1, length(Dc$d))
}
WSS_weight <- opts$init_w
opts$init_w <- 1
solWSS <- mosalloc(WSS_weight %*% Dc$D,
as.vector(WSS_weight %*% Dc$d),
Ac$A, Ac$a, Cc$C, Cc$c, l, u, opts)
# Compute objective vector
if (opts$sense == "min_cost") {
obj <- Dc$D %*% sol$n + Dc$d
} else {
if (fpc == TRUE) {
obj <- Dc$D %*% (1 / solWSS$n - 1 / N)
} else {
obj <- Dc$D %*% (1 / solWSS$n)
}
}
# Obtain optimal weights for weighted Chebychev minimization
opts$init_w <- matrix(c(obj[1] / obj), nrow = 1)
if (ForceOptimality == TRUE) {
message("ForceOptimality is set to FALSE as Pareto optimality\n",
"will be achieved from using WSS!\n")
ForceOptimality <- FALSE
}
time_t1 <- proc.time()
}
if (ForceOptimality == FALSE) {
sol <- mosalloc(Dc$D, Dc$d, Ac$A, Ac$a, Cc$C, Cc$c, l, u, opts)
if (method == "WCM") {
Timing <- sol$Timing
rownames(Timing) <- ""
message(paste0(capture.output(print(Timing)), collapse = "\n"))
}
if (method == "WSS") {
Timing <- sol$Timing
Timing[1:2] <- Timing[1:2] + time_t1[3] - time_t0[3]
rownames(Timing) <- ""
message(paste0(capture.output(print(Timing)), collapse = "\n"))
}
message(" -> ECOSolveR statement: ", sol$Ecosolver$Ecoinfostring, "!\n")
if (any(sol$Sensitivity$D < 1e-7)) {
message("\nCaution! Preference weighting not optimal.\n",
"Verify init_w = ", paste0(capture.output(print(sol$J[1] / sol$J))),
"\n")
}
if (any(sol$Sensitivity$D < 1e-7) && sum(abs(sol$Qbounds - 1)) > 1e-8) {
message("Pareto optimality unclear! Please verify solver output.\n")
}
} else {
sol <- mosallocStepwiseFirst(Dc$D, Dc$d, Ac$A, Ac$a, Cc$C, Cc$c,
l, u, opts)
message(" -> ECOSolveR statement: ", sol$Ecosolver$Ecoinfostring, "!\n")
}
message("----------------------------------------------------------\n")
if (opts$sense == "min_cost") {
obj <- Dc$D %*% sol$n + Dc$d
outD <- data.frame(value = obj, sensitivity = sol$Sensitivity$D)
} else {
obj <- Dc$D %*% (1 / sol$n) - Dc$d
if (fpc == TRUE) {
out_RSE <- sqrt(Dc_RSE$D %*% (1 / sol$n - 1 / N))
} else {
out_RSE <- sqrt(Dc_RSE$D %*% (1 / sol$n))
}
outD <- data.frame(value = obj,
measure = sapply(listD, function(L)L$measure),
sensitivity = sol$Sensitivity$D,
RSE = out_RSE)
}
if (!is.null(Ac$A)) {
constA <- Ac$A %*% (1 / sol$n) - Ac$a + sapply(listA,
function(L)L$bound)**2
outA <- data.frame(value = constA,
bound = sapply(listA, function(L)L$bound),
measure = sapply(listA, function(L)L$measure),
sensitivity = sol$Sensitivity$A)
for (i in 1:nrow(outA)) {
if (outA[i, "measure"] == "RSE") {
outA[i, "value"] <- sqrt(outA[i, "value"])
}
}
} else {
outA <- NULL
}
if (!is.null(Cc$C)) {
outC <- data.frame(opt_cost = abs(Cc$C %*% sol$n),
bound = sapply(listC, function(L)c(L$c_lower, L$c_upper)),
sensitivity = sol$Sensitivity$C)
} else {
outC <- NULL
}
if (ForceOptimality == TRUE) {
init_w <- matrix(sol$w, 1)
rownames(init_w) <- c("WCM")
colnames(init_w) <- rownames(Dc$D)
out <- list(sense = opts$sense, method = method, init_w = init_w,
opt_w = as.vector(obj[1] / obj), n_opt = sol$n,
objectives = outD, precision = outA, cost = outC,
problem_components = list(objective = Dc,
precision_constr = Ac,
cost_constr = Cc,
box_constr = t(matrix(c(u, l),
ncol = 2,
dimnames = list(1:length(N), c("upper", "lower"))))),
output_mosalloc = sol)
} else {
if (method == "WSS") {
init_w <- rbind(WSS_weight, sol$w)
rownames(init_w) <- c("WSS", "WCM")
colnames(init_w) <- rownames(Dc$D)
} else {
init_w <- matrix(sol$w, 1)
rownames(init_w) <- c("WCM")
colnames(init_w) <- rownames(Dc$D)
}
out <- list(sense = opts$sense, method = method, init_w = init_w,
opt_w = NULL, n_opt = sol$n,
objectives = outD, precision = outA, cost = outC,
problem_components = list(objective = Dc,
precision_constr = Ac,
cost_constr = Cc,
box_constr = t(matrix(c(u, l),
ncol = 2,
dimnames = list(1:length(N), c("upper", "lower"))))),
output_mosalloc = sol)
}
class(out) <- "mosaSTRS"
return(out)
} else {
if (method == "WSS") {
stop("Parallelization of multiple preference weightings has not yet
been implemented for method == 'WSS'!")
}
message(" mosalloc running...\n Initial weight matrix:\n",
paste0(capture.output(print(opts$init_w)), collapse = "\n"),
"\n")
if (ForceOptimality == FALSE) {
message(" -> ECOSolveR statements:\n")
sol <- mosalloc(Dc$D, Dc$d, Ac$A, Ac$a, Cc$C, Cc$c, l, u, opts)
lapply(sol, function(L) {
message("Preference weighting:", paste0(capture.output(print(L$w))),
": ", L$Ecosolver$Ecoinfostring, "!\n")
if (any(L$Sensitivity$D < 1e-7)) {
message("\nCaution! Preference weighting not optimal.\n",
"Verify init_w = ", paste0(capture.output(print(L$J[1] / L$J))),
"\n")
}
if (any(L$Sensitivity$D < 1e-7) && sum(abs(L$Qbounds - 1)) > 1e-8) {
message("Pareto optimality unclear! Please verify solver output.\n")
}
message(" --- \n")
})
} else {
sol <- mosallocStepwiseFirst(Dc$D, Dc$d, Ac$A, Ac$a, Cc$C, Cc$c,
l, u, opts)
message(" -> ECOSolveR statements:\n")
lapply(sol, function(L) {
message("Preference weighting:", L$w, ": ", L$Ecosolver$Ecoinfostring,
"!\n", " --- \n")
})
}
message("----------------------------------------------------------\n\n")
Lout <- lapply(sol, function(Lsol) {
if (opts$sense == "min_cost") {
obj <- Dc$D %*% Lsol$n + Dc$d
outD <- data.frame(value = obj, sensitivity = Lsol$Sensitivity$D)
} else {
if (fpc == TRUE) {
obj <- Dc$D %*% (1 / Lsol$n - 1 / N)
} else {
obj <- Dc$D %*% (1 / Lsol$n)
}
out_RSE <- sqrt(Dc_RSE$D %*% (1 / Lsol$n - 1 / N))
outD <- data.frame(value = obj,
measure = sapply(listD, function(L)L$measure),
sensitivity = Lsol$Sensitivity$D,
RSE = out_RSE)
}
if (!is.null(Ac$A)) {
constA <- Ac$A %*% (1/Lsol$n) - Ac$a + sapply(listA, function(L)L$bound)
outA <- data.frame(value = constA,
bound = sapply(listA, function(L)L$bound),
measure = sapply(listA, function(L)L$measure),
sensitivity = Lsol$Sensitivity$A)
for (i in 1:nrow(outA)) {
if (outA[i, "measure"] == "RSE") {
outA[i, "value"] <- sqrt(outA[i, "value"])
}
}
} else {
outA <- NULL
}
if (!is.null(Cc)) {
outC <- data.frame(opt_cost = abs(Cc$C %*% Lsol$n),
bound = sapply(listC, function(L)c(L$c_lower, L$c_upper)),
sensitivity = Lsol$Sensitivity$C)
} else {
outC <- NULL
}
if (ForceOptimality == TRUE) {
init_w <- matrix(Lsol$w, 1)
rownames(init_w) <- c("WCM")
colnames(init_w) <- rownames(Dc$D)
out <- list(sense = opts$sense, method = method, init_w = init_w,
opt_w = as.vector(obj[1] / obj), n_opt = Lsol$n,
objectives = outD, precision = outA, cost = outC,
problem_components = list(objective = Dc,
precision_constr = Ac,
cost_constr = Cc,
box_constr = t(matrix(c(u, l),
ncol = 2,
dimnames = list(1:length(N), c("upper", "lower"))))),
output_mosalloc = Lsol)
} else {
init_w <- matrix(Lsol$w, 1)
rownames(init_w) <- c("WCM")
colnames(init_w) <- rownames(Dc$D)
out <- list(sense = opts$sense, method = method, init_w = init_w,
opt_w = NULL, n_opt = Lsol$n,
objectives = outD, precision = outA, cost = outC,
problem_components = list(objective = Dc,
precision_constr = Ac,
cost_constr = Cc,
box_constr = t(matrix(c(u, l),
ncol = 2,
dimnames = list(1:length(N), c("upper", "lower"))))),
output_mosalloc = Lsol)
}
class(out) <- "mosaSTRS"
return(out)
})
class(Lout) <- "mosaSTRS"
return(Lout)
}
}
Try the MOSAlloc package in your browser
Any scripts or data that you put into this service are public.
MOSAlloc documentation built on Feb. 14, 2026, 5:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions mosallocSTRS
Documented in mosallocSTRS
# Filename: mosallocSTRS.R
#
# Date: 18.07.2025; modified: 10.09.2025, 16.01.2026
# Author: Felix Willems
# Contact: mail.willemsf+MOSAlloc@gmail.com
# (mail[DOT]willemsf+MOSAlloc[AT]gmail[DOT]com)
# Licensing: GPL-3.0-or-later
#
# Please report any bugs or unexpected behavior to
# mail.willemsf+MOSAlloc@gmail.com
# (mail[DOT]willemsf+MOSAlloc[AT]gmail[DOT]com)
#
#---------------------------------------------------------------------------
#
#' @title (Single-stage) stratified random sampling interface for functions
#' \code{mosalloc()} and \code{mosallocStepwiseFirst()}
#'
#' @description Interface for the functions \code{mosalloc()} and
#' \code{mosallocStepwiseFirst()} of the same package which allows for a
#' user-friendly data handling in the case of single-stage stratified random
#' sampling.
#'
#' @param X_var (type: \code{matrix})
#' A matrix containing stratum- (row) and variable- (column) specific
#' precision units (e.g. variances).
#' @param X_tot (type: \code{matrix})
#' A matrix containing stratum- (row) and variable- (column) specific totals.
#' @param N (type: \code{vector})
#' A vector of stratum sizes.
#' @param listD (type: \code{list})
#' A list of lists taking subpopulation- (domain/area) specific
#' arguments. \cr If \code{opts$sense == "max_precision"}, elements are lists
#' containing the following components
#' which correspond to one specific precision target each:
#' \cr \code{..$stratum_id} (type: \code{numeric}) A vector containing the
#' indices of the strata considered for the current restriction. The indices
#' must coincide with the row numbers of \code{X_var} and \code{X_tot}.
#' \cr \code{..$variate} (type: \code{character} or \code{numeric})
#' The column name or column index of \code{X_var} to be addressed.
#' \cr \code{..$measure} (type: \code{character} or \code{numeric})
#' A character ("relVAR" (relative variance) or "VAR" (variance)) or a numerical
#' value between 0 and 1 indicates a gradation coefficient that balances between
#' "relVar" (..$measure = 0) and "VAR" (..$measure = 1).
#' \cr \code{..$name} (type: \code{character}) The name of the subpopulation
#' (domain/area).
#' \cr If \code{opts$sense == "min_cost"}, elements are lists
#' containing the following components
#' which correspond to one specific cost target each:
#' \cr \code{..$stratum_id} (type: \code{numeric})
#' A vector containing the indices of the strata considered for the current
#' objective. The indices must coincide with the row numbers of \code{X_cost}.
#' \cr \code{..$c_type} (type: \code{character} or \code{numeric})
#' The column name or column index of \code{X_cost} to be addressed.
#' \cr \code{..$name} (type: \code{character})
#' The name of the subpopulation (domain/area).
#' @param listA (type: \code{list})
#' A list of lists taking subpopulation- (domain/area) specific
#' arguments. Elements are lists containing the following components
#' which in turn correspond to one specific precision restriction:
#' \cr \code{..$stratum_id} (type: \code{numeric}) A vector containing the
#' indices of the strata considered for the current restriction. The indices
#' must coincide with the corresponding row numbers of \code{X_var} and
#' \code{X_tot}.
#' \cr \code{..$variate} (type: \code{character} or \code{numeric})
#' The column name or column index of \code{X_var} to be addressed.
#' \cr \code{..$measure} (type: \code{character}) "RSE"
#' (relative standard error) or "VAR" (variance).
#' \cr \code{..$bound} (type: \code{numeric})
#' An upper bound to "RSE" (or "VAR").
#' \cr \code{..$name} (type: \code{character}) The name of the subpopulation
#' (domain/area).
#' @param listC (type: \code{list})
#' A list of lists taking subpopulation- (domain/area) specific
#' arguments. Elements are lists containing the following components which
#' correspond to one specific cost restriction:
#' \cr \code{..$stratum_id} (type: \code{numeric})
#' A vector containing the indices of the strata considered for the current
#' restriction. The indices must coincide with the row numbers of \code{X_var}
#' and \code{X_tot}.
#' \cr \code{..$c_coef} (type: \code{numeric}) A vector of length
#' \code{length(stratum_id)} containing the stratum-specific cost components
#' for the set of strata that is going to be bounded by above and/or below.
#' \cr \code{..$c_upper} (type: \code{numeric}) The cost upper bound value.
#' \code{NULL} if not present.
#' \cr \code{..$c_lower} (type: \code{numeric}) The cost lower bound value.
#' \code{NULL} if not present.
#' \cr \code{..$name} (type: \code{character}) The name of the subpopulation
#' (domain/area).
#' @param fpc (type: \code{logical})
#' A \code{TRUE} or \code{FALSE} statement indicating whether the finite
#' population correction should be considered.
#' @param l (type: \code{vector})
#' A vector of lower box constraints.
#' @param u (type: \code{vector})
#' A vector of upper box constraints.
#' @param opts (type: \code{list})
#' The options used by the algorithms:
#' \cr \code{$sense} (type: character) Sense of optimization
#' (default = \code{"max_precision"}, alternative \code{"min_cost"}).
#' \cr \code{$df} (type: function) The gradient of f (default = \code{NULL}).
#' \cr \code{$Hf} (type: function)
#' The Hesse matrix of f (default = \code{NULL}).
#' \cr \code{$method} (type: \code{character})
#' A character indicating scalarization method (default = \code{"WCM"},
#' alternative \code{"WSS"}), \code{$f} must be \code{NULL}.
#' \cr \code{$f} (type: function) Decision functional over the objective
#' vector (default = \code{NULL}). The weighted Chebyshev minimization
#' (\code{WCM}) minimizes for weighted maximum objective component. The
#' weighted sum scalarization (\code{WSS}) minimizes for the weighted sum.
#' For either case the weights are given by \code{$init_w}.
#' \cr \code{$init_w} (type: \code{numeric} or \code{matrix})
#' Preference weightings (default = 1; The weight for first objective component
#' must be 1).
#' \cr \code{$mc_cores} (type: \code{integer})
#' The number of cores for parallelizing
#' multiple input weightings stacked rowwise (default = \code{1L}).
#' \cr \code{$pm_tol} (type: \code{numeric}) The tolerance for the projection
#' method (default = \code{1e-5}).
#' \cr \code{max_iters} (type: \code{integer}) The maximum number of iterations
#' (default = \code{100L}).
#' \cr \code{$print_pm} (type: \code{logical}) A \code{TRUE} or \code{FALSE}
#' statement whether iterations of the projection method should be printed
#' (default = \code{FALSE}).
#' @param ForceOptimality (type: \code{logical})
#' A \code{TRUE} or \code{FALSE} statement indicating whether Pareto optimality
#' should be ensured. Default is \code{FALSE} (for most practical
#' problem formulations Pareto optimality is achieved by construction,
#' e.g. in the case of one overall cost constraint).
#' If \code{TRUE}, additional computation time is required. In this case,
#' \code{mosallocStepwiseFirst()} is used internally.
#' @param X_cost (type: \code{matrix})
#' A matrix containing stratum- (row) and type- (column) specific cost
#' coefficients associated with fixed cost. Types of cost might be, e.g.
#' '$ US', 'minutes', 'sample size', etc. Default is \code{NULL}. The argument
#' is required in the case of cost minimization (opt$sense == "min_cost").
#' @param X_fixed (type: \code{matrix})
#' A matrix containing stratum- (rows) and type- (columns) specific
#' cost coefficients associated with fixed cost. Default is \code{NULL}.
#' Used solely in the case of cost minimization (\code{opt$sense == "min_cost"}).
#'
#' @return A \code{mosaSTRS} object or a list of \code{mosaSTRS}
#' objects. A \code{mosaSTRS} object is a list containing the
#' following components:
#' @returns \code{$sense}
#' Sense of optimization; \code{max precision} or \code{min_cost}.
#' @returns \code{$method}
#' The method used, either weighted sum scalarization (WSS) or weighted
#' Chebyshev minimization (WCM).
#' @returns \code{$init_w}
#' The initial preference weightings (opts$init_w).
#' @returns \code{$opt_w}
#' The optimal weightings w.r.t. \code{opts$init_w} as \code{opts$init_w} might
#' not lead to Pareto optimality. \code{NULL} if \code{ForceOptimality = FALSE}.
#' @returns \code{$n_opt}
#' The vector of optimal sample sizes.
#' @returns \code{$objective}
#' The objective values, including the sensitivity and the RSE.
#' @returns \code{$precision}
#' A data frame corresponding to the precision constraints.
#' @returns \code{$cost}
#' A data frame corresponding to the cost constraints.
#' @returns \code{$problem_components}
#' A list containing the input data to the optimization problem.
#' @returns \code{$output_mosalloc}
#' A list of function returns of \code{mosalloc()}.
#' @returns If \code{opts$init_w} has multiple rows, the function returns a list
#' of \code{mosaSTRS} objects whose length equals the number of rows.
#'
#' @examples
#' # Artificial population of 50 568 business establishments and 5 business
#' # sectors (data from Valliant, R., Dever, J. A., & Kreuter, F. (2013).
#' # Practical tools for designing and weighting survey samples. Springer.
#' # https://doi.org/10.1007/978-1-4614-6449-5, Example 5.2 pages 133-9)
#'
#' # See also https://umd.app.box.com/s/9yvvibu4nz4q6rlw98ac/file/297813512360
#' # file: Code 5.3 constrOptim.example.R
#'
#' Nh <- c(6221, 11738, 4333, 22809, 5467) # stratum sizes
#' ch <- c(120, 80, 80, 90, 150) # stratum-specific cost of surveying
#'
#' # Revenues
#' mh.rev <- c(85, 11, 23, 17, 126) # mean revenue
#' Sh.rev <- c(170.0, 8.8, 23.0, 25.5, 315.0) # standard deviation revenue
#'
#' # Employees
#' mh.emp <- c(511, 21, 70, 32, 157) # mean number of employees
#' Sh.emp <- c(255.50, 5.25, 35.00, 32.00, 471.00) # std. dev. employees
#'
#' # Proportion of estabs claiming research credit
#' ph.rsch <- c(0.8, 0.2, 0.5, 0.3, 0.9)
#'
#' # Proportion of estabs with offshore affiliates
#' ph.offsh <- c(0.06, 0.03, 0.03, 0.21, 0.77)
#'
#' budget <- 300000 # overall available budget
#' n.min <- 100 # minimum stratum-specific sample size
#'
#' # Matrix containing stratum-specific variance components
#' X_var <- cbind(Sh.rev**2,
#' Sh.emp**2,
#' ph.rsch * (1 - ph.rsch) * Nh/(Nh - 1),
#' ph.offsh * (1 - ph.offsh) * Nh/(Nh - 1))
#' colnames(X_var) <- c("rev", "emp", "rsch", "offsh")
#'
#' # Matrix containing stratum-specific totals
#' X_tot <- cbind(mh.rev, mh.emp, ph.rsch, ph.offsh) * Nh
#' colnames(X_tot) <- c("rev", "emp", "rsch", "offsh")
#'
#' # Examples
#' #----------------------------------------------------------------------------
#' # Example 1: Univariate minimization of the variation of estimates for
#' # revenue subject to cost restrictions and precision restrictions to the
#' # relative standard error of estimates for the proportion of businesses with
#' # offshore affiliates. Additionally, there is one overall cost constraint and
#' # at least half of the provided budget must be spend to strata 1 to 3.
#'
#' # Specify objectives via listD
#' listD <- list(list(stratum_id = 1:5, variate = "rev", measure = "relVAR",
#' name = "pop"))
#'
#' # Specify precision constraints via listA
#' listA <- list(list(stratum_id = 1:5, variate = "offsh", measure = "RSE",
#' bound = 0.05, name = "pop"))
#'
#' # Specify cost constraints via listC
#' listC <- list(list(stratum_id = 1:5, c_coef = ch, c_lower = NULL,
#' c_upper = budget, name = "Overall"),
#' list(stratum_id = 1:3, c_coef = ch[1:3],
#' c_lower = 0.5 * budget, c_upper = NULL, name = "1to3"))
#'
#' # Specify stratum-specific box constraints
#' l <- rep(n.min, 5) # minimum sample size per stratum
#' u <- Nh # maximum sample size per stratum
#'
#' # Specify parameter for mosalloc (method = "WSS")
#' opts <- list(sense = "max_precision",
#' f = NULL, df = NULL, Hf = NULL,
#' method = "WSS", init_w = 1,
#' mc_cores = 1L, pm_tol = 1e-05,
#' max_iters = 100L, print_pm = FALSE)
#'
#' # Run mosallocSTRS with weighted sum scalarization (WSS)
#' resWSS <- mosallocSTRS(X_var, X_tot, Nh, listD, listA, listC,
#' fpc = TRUE, l, u, opts)
#'
#' summary(resWSS)
#'
#' # Specify parameter for mosalloc (method = "WCM")
#' opts = list(sense = "max_precision",
#' f = NULL, df = NULL, Hf = NULL,
#' method = "WCM", init_w = 1,
#' mc_cores = 1L, pm_tol = 1e-05,
#' max_iters = 100L, print_pm = FALSE)
#'
#' # Run mosallocSTRS with weighted Chebyshec minimization (WCM)
#' resWCM <- mosallocSTRS(X_var, X_tot, Nh, listD, listA, listC,
#' fpc = TRUE, l, u, opts)
#'
#' summary(resWCM)
#'
#' # The optimal sample sizes vector can also be obtained by
#' summary(resWCM)$n_opt
#'
#' # Hint: For univariate allocation problems 'WSS' and 'WCM' are equivalent!
#'
#' #----------------------------------------------------------------------------
#' # Example 2: Minimization of the maximum relative variation of estimates for
#' # the total revenue, the number of employee, the number of businesses claimed
#' # research credit and the number of businesses with offshore affiliates
#' # subject to one overall cost constraint and at least half of the provided
#' # budget must be spend to strata 1 to 3.
#'
#' # Specify objectives via listD
#' listD <- list(list(stratum_id = 1:5, variate = "rev", measure = "relVAR",
#' name = "pop"),
#' list(stratum_id = 1:5, variate = "emp", measure = "relVAR",
#' name = "pop"),
#' list(stratum_id = 1:5, variate = "rsch", measure = "relVAR",
#' name = "pop"),
#' list(stratum_id = 1:5, variate = "offsh", measure = "relVAR",
#' name = "pop"))
#'
#' # Specify cost constraints via listC
#' listC <- list(list(stratum_id = 1:5, c_coef = ch, c_lower = NULL,
#' c_upper = budget, name = "Overall"),
#' list(stratum_id = 1:3, c_coef = ch[1:3],
#' c_lower = 0.5 * budget, c_upper = NULL, name = "1to3"))
#'
#' # Specify stratum-specific box constraints
#' l <- rep(n.min, 5) # minimum sample size per stratum
#' u <- Nh # maximum sample size per stratum
#'
#' # Specify parameter for mosalloc (method = "WSS")
#' opts = list(sense = "max_precision",
#' f = NULL, df = NULL, Hf = NULL,
#' method = "WSS", init_w = 1,
#' mc_cores = 1L, pm_tol = 1e-05,
#' max_iters = 100L, print_pm = FALSE)
#'
#' # Run mosallocSTRS with weighted sum scalarization (WSS)
#' resWSS <- mosallocSTRS(X_var, X_tot, Nh, listD, NULL, listC,
#' fpc = TRUE, l, u, opts)
#'
#' summary(resWSS)
#'
#' # Specify parameter for mosalloc (method = "WCM")
#' opts = list(sense = "max_precision",
#' f = NULL, df = NULL, Hf = NULL,
#' method = "WCM", init_w = 1,
#' mc_cores = 1L, pm_tol = 1e-05,
#' max_iters = 100L, print_pm = FALSE)
#'
#' # Run mosallocSTRS with weighted Chebyshec minimization (WCM)
#' resWCM <- mosallocSTRS(X_var, X_tot, Nh, listD, NULL, listC,
#' fpc = TRUE, l, u, opts)
#'
#' summary(resWCM)
#'
#' # Since the WCM does not necessarily lead to Pareto optimal allocations,
#' # we might force this via a internal stepwise procedure by setting
#' # ForceOptimality = TRUE.
#'
#' resWCM_FO <- mosallocSTRS(X_var, X_tot, Nh, listD, NULL, listC,
#' fpc = TRUE, l, u, opts, ForceOptimality = TRUE)
#'
#' summary(resWCM_FO)
#'
#' #----------------------------------------------------------------------------
#' # Example 3: Example 2 with multiple sets of preference weightings.
#'
#' # Define a set of preference weightings, e.g.
#' w_1 <- c(1, 1, 1, 1)
#' w_2 <- c(1, 2, 2, 1)
#' w_3 <- c(1, 5, 5, 1)
#'
#' # Combine the weightings to a matrix stacked rowwise
#' w <- rbind(w_1, w_2, w_3)
#'
#' # Specify parameter for mosalloc() (method = "WCM"; not yet possible with WSS)
#' opts = list(sense = "max_precision",
#' f = NULL, df = NULL, Hf = NULL,
#' method = "WCM", init_w = w,
#' mc_cores = 1L, pm_tol = 1e-05,
#' max_iters = 100L, print_pm = FALSE)
#'
#' # Run mosallocSTRS with weighted Chebyshec minimization (WCM)
#' res <- mosallocSTRS(X_var, X_tot, Nh, listD, NULL, listC, fpc = TRUE,
#' l, u, opts)
#'
#' summary(res)
#'
#' #----------------------------------------------------------------------------
#' # Example 4: Minimization of survey cost subject to quality restrictions on
#' # subpopulation level.
#'
#' X_cost <- matrix(ch, 5, 1, dimnames = list(1:5,"cost"))
#'
#' # Specify cost objectives via listD
#' listD <- list(list(stratum_id = 1:5, c_type = "cost", name = "pop"))
#'
#' # Specify quailty restrictions via listD. Here: 5 % relative standard error
#' listA <- list(list(stratum_id = 1:2, variate = "rev", measure = "RSE",
#' bound = 0.05, name = "S1-2"),
#' list(stratum_id = 3:5, variate = "rev", measure = "RSE",
#' bound = 0.05, name = "S3-5"),
#' list(stratum_id = 1:2, variate = "emp", measure = "RSE",
#' bound = 0.05, name = "S1-2"),
#' list(stratum_id = 3:5, variate = "emp", measure = "RSE",
#' bound = 0.05, name = "S3-5"),
#' list(stratum_id = 1:2, variate = "rsch", measure = "RSE",
#' bound = 0.05, name = "S1-2"),
#' list(stratum_id = 3:5, variate = "rsch", measure = "RSE",
#' bound = 0.05, name = "S3-5"),
#' list(stratum_id = 1:2, variate = "offsh", measure = "RSE",
#' bound = 0.05, name = "S1-2"),
#' list(stratum_id = 3:5, variate = "offsh", measure = "RSE",
#' bound = 0.05, name = "S3-5"))
#'
#' # Specify cost constraints
#' listC <- NULL
#'
#' # Specify stratum-specific box constraints
#' l <- rep(n.min, 5) # minimum sample size per stratum
#' u <- Nh # maximum sample size per stratum
#'
#' # Specify parameters for mosalloc()
#' opts = list(sense = "min_cost",
#' f = NULL, df = NULL, Hf = NULL,
#' method = "WCM", init_w = 1,
#' mc_cores = 1L, pm_tol = 1e-05,
#' max_iters = 100L, print_pm = FALSE)
#'
#' # Run mosallocSTRS()
#' res <- mosallocSTRS(X_var, X_tot, Nh, listD, listA, NULL, fpc = TRUE,
#' l, u, opts, X_cost = X_cost)
#'
#' summary(res)
#'
#' # Optimal sample sizes
#' summary(res)$n_opt
#'
#' @importFrom utils capture.output
#'
#' @export
mosallocSTRS <- function(X_var, X_tot, N,
listD, listA = NULL, listC = NULL,
fpc = TRUE, l = 2, u = NULL, opts,
ForceOptimality = FALSE,
X_cost = NULL, X_fixed = NULL) {
message("\n----------------------------------------------------------\n")
# get number of strata in single-stage stratified random sampling (STRS)
H <- length(N)
# expand l and u if necessary
if (length(l) == 1) {
l <- rep(l, H)
}
if (length(u) == 1) {
l <- rep(u, H)
}
# check input
if (opts$sense == "min_cost") {
if (is.null(X_cost)) {
stop("X_cost not specified!")
} else {
if (is.null(X_fixed)) {
X_fixed <- X_cost - X_cost
}
Dc <- constructDobjCostSTRS(X_cost, X_fixed, listD)
}
} else {
if (length(unique(sapply(listD, function(L)L$measure))) != 1) {
message("CAUTION: ..$measures in listD varry among objectives!\n\n")
}
# Construct vector valued precision objective
Dc <- constructDobjPrecisionSTRS(X_var, X_tot, N, listD, fpc)
listD_RSE <- lapply(listD, function(L){L$measure = 0; L})
Dc_RSE <- constructDobjPrecisionSTRS(X_var, X_tot, N, listD_RSE)
}
if (is.null(opts$f)) {
method <- opts$method
} else {
stop("The possiblity to use decision functionals (opts$f) is not yet
implemented! Please set opts$f = NULL.")
}
# Construct precision constraints
if (is.null(listA)) {
Ac <- list(A = NULL, a = NULL)
} else {
Ac <- constructArestrSTRS(X_var, X_tot, N, listA, fpc)
}
# Construct cost constraints
if (is.null(listC)) {
Cc <- list(C = NULL, c = NULL)
} else {
Cc <- constructCrestrSTRS(H, listC)
}
if (is.null(dim(opts$init_w)) || dim(opts$init_w)[1] == 1) {
message(" mosalloc running...\n")
if (method == "WSS") {
# Weight objetives via opts$init_w and solve univariate problem
time_t0 <- proc.time()
if (is.null(dim(opts$init_w))) {
opts$init_w <- matrix(opts$init_w, 1, length(Dc$d))
}
WSS_weight <- opts$init_w
opts$init_w <- 1
solWSS <- mosalloc(WSS_weight %*% Dc$D,
as.vector(WSS_weight %*% Dc$d),
Ac$A, Ac$a, Cc$C, Cc$c, l, u, opts)
# Compute objective vector
if (opts$sense == "min_cost") {
obj <- Dc$D %*% sol$n + Dc$d
} else {
if (fpc == TRUE) {
obj <- Dc$D %*% (1 / solWSS$n - 1 / N)
} else {
obj <- Dc$D %*% (1 / solWSS$n)
}
}
# Obtain optimal weights for weighted Chebychev minimization
opts$init_w <- matrix(c(obj[1] / obj), nrow = 1)
if (ForceOptimality == TRUE) {
message("ForceOptimality is set to FALSE as Pareto optimality\n",
"will be achieved from using WSS!\n")
ForceOptimality <- FALSE
}
time_t1 <- proc.time()
}
if (ForceOptimality == FALSE) {
sol <- mosalloc(Dc$D, Dc$d, Ac$A, Ac$a, Cc$C, Cc$c, l, u, opts)
if (method == "WCM") {
Timing <- sol$Timing
rownames(Timing) <- ""
message(paste0(capture.output(print(Timing)), collapse = "\n"))
}
if (method == "WSS") {
Timing <- sol$Timing
Timing[1:2] <- Timing[1:2] + time_t1[3] - time_t0[3]
rownames(Timing) <- ""
message(paste0(capture.output(print(Timing)), collapse = "\n"))
}
message(" -> ECOSolveR statement: ", sol$Ecosolver$Ecoinfostring, "!\n")
if (any(sol$Sensitivity$D < 1e-7)) {
message("\nCaution! Preference weighting not optimal.\n",
"Verify init_w = ", paste0(capture.output(print(sol$J[1] / sol$J))),
"\n")
}
if (any(sol$Sensitivity$D < 1e-7) && sum(abs(sol$Qbounds - 1)) > 1e-8) {
message("Pareto optimality unclear! Please verify solver output.\n")
}
} else {
sol <- mosallocStepwiseFirst(Dc$D, Dc$d, Ac$A, Ac$a, Cc$C, Cc$c,
l, u, opts)
message(" -> ECOSolveR statement: ", sol$Ecosolver$Ecoinfostring, "!\n")
}
message("----------------------------------------------------------\n")
if (opts$sense == "min_cost") {
obj <- Dc$D %*% sol$n + Dc$d
outD <- data.frame(value = obj, sensitivity = sol$Sensitivity$D)
} else {
obj <- Dc$D %*% (1 / sol$n) - Dc$d
if (fpc == TRUE) {
out_RSE <- sqrt(Dc_RSE$D %*% (1 / sol$n - 1 / N))
} else {
out_RSE <- sqrt(Dc_RSE$D %*% (1 / sol$n))
}
outD <- data.frame(value = obj,
measure = sapply(listD, function(L)L$measure),
sensitivity = sol$Sensitivity$D,
RSE = out_RSE)
}
if (!is.null(Ac$A)) {
constA <- Ac$A %*% (1 / sol$n) - Ac$a + sapply(listA,
function(L)L$bound)**2
outA <- data.frame(value = constA,
bound = sapply(listA, function(L)L$bound),
measure = sapply(listA, function(L)L$measure),
sensitivity = sol$Sensitivity$A)
for (i in 1:nrow(outA)) {
if (outA[i, "measure"] == "RSE") {
outA[i, "value"] <- sqrt(outA[i, "value"])
}
}
} else {
outA <- NULL
}
if (!is.null(Cc$C)) {
outC <- data.frame(opt_cost = abs(Cc$C %*% sol$n),
bound = sapply(listC, function(L)c(L$c_lower, L$c_upper)),
sensitivity = sol$Sensitivity$C)
} else {
outC <- NULL
}
if (ForceOptimality == TRUE) {
init_w <- matrix(sol$w, 1)
rownames(init_w) <- c("WCM")
colnames(init_w) <- rownames(Dc$D)
out <- list(sense = opts$sense, method = method, init_w = init_w,
opt_w = as.vector(obj[1] / obj), n_opt = sol$n,
objectives = outD, precision = outA, cost = outC,
problem_components = list(objective = Dc,
precision_constr = Ac,
cost_constr = Cc,
box_constr = t(matrix(c(u, l),
ncol = 2,
dimnames = list(1:length(N), c("upper", "lower"))))),
output_mosalloc = sol)
} else {
if (method == "WSS") {
init_w <- rbind(WSS_weight, sol$w)
rownames(init_w) <- c("WSS", "WCM")
colnames(init_w) <- rownames(Dc$D)
} else {
init_w <- matrix(sol$w, 1)
rownames(init_w) <- c("WCM")
colnames(init_w) <- rownames(Dc$D)
}
out <- list(sense = opts$sense, method = method, init_w = init_w,
opt_w = NULL, n_opt = sol$n,
objectives = outD, precision = outA, cost = outC,
problem_components = list(objective = Dc,
precision_constr = Ac,
cost_constr = Cc,
box_constr = t(matrix(c(u, l),
ncol = 2,
dimnames = list(1:length(N), c("upper", "lower"))))),
output_mosalloc = sol)
}
class(out) <- "mosaSTRS"
return(out)
} else {
if (method == "WSS") {
stop("Parallelization of multiple preference weightings has not yet
been implemented for method == 'WSS'!")
}
message(" mosalloc running...\n Initial weight matrix:\n",
paste0(capture.output(print(opts$init_w)), collapse = "\n"),
"\n")
if (ForceOptimality == FALSE) {
message(" -> ECOSolveR statements:\n")
sol <- mosalloc(Dc$D, Dc$d, Ac$A, Ac$a, Cc$C, Cc$c, l, u, opts)
lapply(sol, function(L) {
message("Preference weighting:", paste0(capture.output(print(L$w))),
": ", L$Ecosolver$Ecoinfostring, "!\n")
if (any(L$Sensitivity$D < 1e-7)) {
message("\nCaution! Preference weighting not optimal.\n",
"Verify init_w = ", paste0(capture.output(print(L$J[1] / L$J))),
"\n")
}
if (any(L$Sensitivity$D < 1e-7) && sum(abs(L$Qbounds - 1)) > 1e-8) {
message("Pareto optimality unclear! Please verify solver output.\n")
}
message(" --- \n")
})
} else {
sol <- mosallocStepwiseFirst(Dc$D, Dc$d, Ac$A, Ac$a, Cc$C, Cc$c,
l, u, opts)
message(" -> ECOSolveR statements:\n")
lapply(sol, function(L) {
message("Preference weighting:", L$w, ": ", L$Ecosolver$Ecoinfostring,
"!\n", " --- \n")
})
}
message("----------------------------------------------------------\n\n")
Lout <- lapply(sol, function(Lsol) {
if (opts$sense == "min_cost") {
obj <- Dc$D %*% Lsol$n + Dc$d
outD <- data.frame(value = obj, sensitivity = Lsol$Sensitivity$D)
} else {
if (fpc == TRUE) {
obj <- Dc$D %*% (1 / Lsol$n - 1 / N)
} else {
obj <- Dc$D %*% (1 / Lsol$n)
}
out_RSE <- sqrt(Dc_RSE$D %*% (1 / Lsol$n - 1 / N))
outD <- data.frame(value = obj,
measure = sapply(listD, function(L)L$measure),
sensitivity = Lsol$Sensitivity$D,
RSE = out_RSE)
}
if (!is.null(Ac$A)) {
constA <- Ac$A %*% (1/Lsol$n) - Ac$a + sapply(listA, function(L)L$bound)
outA <- data.frame(value = constA,
bound = sapply(listA, function(L)L$bound),
measure = sapply(listA, function(L)L$measure),
sensitivity = Lsol$Sensitivity$A)
for (i in 1:nrow(outA)) {
if (outA[i, "measure"] == "RSE") {
outA[i, "value"] <- sqrt(outA[i, "value"])
}
}
} else {
outA <- NULL
}
if (!is.null(Cc)) {
outC <- data.frame(opt_cost = abs(Cc$C %*% Lsol$n),
bound = sapply(listC, function(L)c(L$c_lower, L$c_upper)),
sensitivity = Lsol$Sensitivity$C)
} else {
outC <- NULL
}
if (ForceOptimality == TRUE) {
init_w <- matrix(Lsol$w, 1)
rownames(init_w) <- c("WCM")
colnames(init_w) <- rownames(Dc$D)
out <- list(sense = opts$sense, method = method, init_w = init_w,
opt_w = as.vector(obj[1] / obj), n_opt = Lsol$n,
objectives = outD, precision = outA, cost = outC,
problem_components = list(objective = Dc,
precision_constr = Ac,
cost_constr = Cc,
box_constr = t(matrix(c(u, l),
ncol = 2,
dimnames = list(1:length(N), c("upper", "lower"))))),
output_mosalloc = Lsol)
} else {
init_w <- matrix(Lsol$w, 1)
rownames(init_w) <- c("WCM")
colnames(init_w) <- rownames(Dc$D)
out <- list(sense = opts$sense, method = method, init_w = init_w,
opt_w = NULL, n_opt = Lsol$n,
objectives = outD, precision = outA, cost = outC,
problem_components = list(objective = Dc,
precision_constr = Ac,
cost_constr = Cc,
box_constr = t(matrix(c(u, l),
ncol = 2,
dimnames = list(1:length(N), c("upper", "lower"))))),
output_mosalloc = Lsol)
}
class(out) <- "mosaSTRS"
return(out)
})
class(Lout) <- "mosaSTRS"
return(Lout)
}
}
Try the MOSAlloc package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.