mosalloc: Multiobjective sample allocation for constraint multivariate...
In MOSAlloc: Constraint Multiobjective Sample Allocation
mosalloc R Documentation
Multiobjective sample allocation for constraint multivariate and
multidomain optimal allocation in survey sampling
Description
Computes solutions to standard sample allocation problems
under various precision and cost restrictions. The input data is
transformed and parsed to the Embedded COnic Solver (ECOS) from the
'ECOSolveR' package. Multiple survey purposes can be optimized
simultaneously through a weighted Chebyshev minimization. Note that in
the case of multiple objectives, mosalloc() does not necessarily
lead to Pareto optimality. This highly depends on the problem structure.
A strong indicator for Pareto optimality is when the weighted objective
values given by Dbounds are constant over all objective
components or when all components of Qbounds equal 1.
In addition, mosalloc() can handle twice-differential convex
decision functionals (in which case Pareto optimality is ensured).
mosalloc() returns dual variables, enabling a detailed sensitivity
analysis.
Usage
mosalloc(
D,
d,
A = NULL,
a = NULL,
C = NULL,
c = NULL,
l = 2,
u = NULL,
opts = list(sense = "max_precision", f = NULL, df = NULL, Hf = NULL, init_w = 1,
mc_cores = 1L, pm_tol = 1e-05, max_iters = 100L, print_pm = FALSE)
)
Arguments
D
(type: matrix)
The objective matrix. A matrix of either precision or cost units.
d
(type: vector)
The objective vector. A vector of either fixed precision components
(e.g. finite population corrections) or fixed costs.
A
(type: matrix)
A matrix of precision units for precision constraints.
a
(type: vector)
The right-hand side vector of the precision constraints.
C
(type: matrix)
A matrix of cost coefficients for cost constraints
c
(type: vector)
The right-hand side vector of the cost constraints.
l
(type: vector)
A vector of lower box constraints.
u
(type: vector)
A vector of upper box constraints.
opts
(type: list)
The options used by the algorithms:
$sense (type: character) Sense of optimization
(default = "max_precision", alternative "min_cost").
$f (type: function) Decision functional over the objective
vector (default = NULL).
$df (type: function)
The gradient of f (default = NULL).
$Hf (type: function)
The Hesse matrix of f (default = NULL).
$init_w (type: numeric or matrix)
Preference weightings (default = 1; The weight for first objective
component must be 1).
$mc_cores (type: integer)
The number of cores for parallelizing multiple input weightings stacked
rowwise (default = 1L).
$pm_tol (type: numeric)
The tolerance for the projection method (default = 1e-5).
max_iters (type: integer)
The maximum number of iterations (default = 100L).
$print_pm (type: logical) A TRUE or FALSE
statement whether iterations of the projection method should be printed
(default = FALSE).
Value
mosalloc() returns a list containing the following
components:
$w The initial preference weighting opts$init_w.
$n The vector of optimal sample sizes.
$J The optimal objective vector.
$Objective The objective value with respect to decision
functional f. NULL if opts$f = NULL.
$Utopian The component-wise univariate optimal
objective vector. NULL if opts$f = NULL.
$Normal The vector normal to the Pareto frontier at
$J.
$dfJ The gradient of opts$f evaluated at $J.
$Sensitivity The dual variables of the objectives and
constraints.
$Qbounds The Quality bounds of the Lorentz cones.
$Dbounds The weighted objective constraints ($w * $J).
$Scalepar An internal scaling parameter.
$Ecosolver A list of ECOSolveR returns including:
...$Ecoinfostring The info string of
ECOSolveR::ECOS_csolve().
...$Ecoredcodes The redcodes of ECOSolveR::ECOS_csolve().
...$Ecosummary Problem summary of ECOSolveR::ECOS_csolve().
$Timing Run time info.
$Iteration A list of internal iterates. NULL if
opts$f = NULL.
Note
Precision optimization (opts$sense == "max_precision",
opts$f == NULL)
The mathematical problem solved is
\min_{n, z, t}\{t: Dz-d\leq w^{-1}t, Az\leq a, Cn\leq c,
1\leq n_iz_i\, \forall i, l\leq n \leq u\}
with \textbf{1}\leq l\leq u \leq N, where
-
n,z\in\mathbb{R}^m and t\in\mathbb{R} are the
optimization variables,
-
D\in\{M\in\mathbb{R}^{k_D\times m}: \sum_i M_{ij}>0\, \forall j,
\sum_j M_{ij}>0\, \forall i\} a matrix of nonnegative objective precision
components with k_D the number of precision objectives
(the number of variables of interest),
-
d\in\mathbb{R}^{k_D} the objective right hand-side (RHS), e.g.
the finite population correction (fpc),
-
A\in\mathbb{R}^{k_A\times m} a nonnegative precision matrix with
k_A the number of precision constraints,
a\in \mathbb{R}^{k_A} the corresponding RHS, e.g. fpc + (upper bound to
the coefficient of variation)^2,
-
C\in\mathbb{R}^{k_C\times m} a cost matrix with k_C the number
of cost constraint,
-
c\in\mathbb{R}^{k_C} the corresponding RHS,
-
l,u\in\mathbb{R}^m the bounds to the sample size vector n,
-
N\in\mathbb{N}^m the vector of population sizes, and
-
w a given strictly positive preference weighting.
Special cases of this formulation are
Neyman-Tschuprow allocation (Neyman, 1934 and Tschuprow, 1923):
\min_n \Big\{\sum_{h=1}^H
\Big(\frac{N_h^2S_h^2}{n_h}-N_hS_h^2\Big):\sum_{h=1}^H n_h \leq c\Big\}
\quad \Leftrightarrow \quad \min_{n, z, t}\{t: Dz-d\leq t, Cn\leq c,
1\leq n_iz_i\, \forall i\}
with D = (N_1^2S_1^2,\dots,N_H^2S_H^2),
d = \sum_{h=1}^H N_hS_h^2, C = (1,\dots,1) and c a
maximum sample size. Here, H is the number of strata, N_h
the size of stratum h and S_h^2 the variance of the variable
of interest in stratum h.
box-constrained optimal allocation:
\min_{n, z, t}\{t: Dz-d\leq t, Cn\leq c,
1\leq n_iz_i\, \forall i, l\leq n \leq u\}
with
D = (N_1^2S_1^2,\dots,N_H^2S_H^2),
d = \sum_{h=1}^H N_hS_h^2, C = (1,\dots,1) and c a
maximum sample size (cf. Srikantan, 1963 and Münnich et al., 2012).
Here, l and u are bounds to the optimal sample size vector.
cost and precision constrained univariate optimal allocation:
\min_{n, z, t}\{t: Dz-d\leq t, Az\leq a, Cn\leq c,
1\leq n_iz_i\, \forall i, l\leq n \leq u\}
with k_D = 1
(cf. Willems, 2025, Chapter 3).
multivariate optimal allocation with weighted sum scalarization:
\min_{n, z, t}\{t: w^\top Dz-w^\top d\leq t, Az\leq a, Cn\leq c,
1\leq n_iz_i\, \forall i, l\leq n \leq u\}
where w\in\mathbb{R}^{k_D}
is a strictly positive preference weighting (cf. Folks and Antle, 1965,
and Rupp, 2018). Note that for this case the problem reduces to cost and
precision constrained univariate optimal allocation. Solutions are ensured
to be optimal in the Pareto sense.
box-constraint two-stage cluster sampling:
\quad \quad \quad\min_{n_\textbf{I},n_\textbf{II}}
\Big\{\Big(\frac{N_\textbf{I}^2
S_\textbf{I}^2}{n_\textbf{I}}-N_\textbf{I}S_\textbf{I}^2\Big) +
\frac{N_\textbf{I}}{n_\textbf{I}}\sum_{j=1}^{N_\textbf{I}}
\Big(\frac{N_{\textbf{II}j}^2S_{\textbf{II}j}^2}{n_{\textbf{II}j}}-
N_{\textbf{II}j}S_{\textbf{II}j}^2\Big):
c_{\textbf{I}}n_\textbf{I} + \frac{n_\textbf{I}}{N_\textbf{I}}
\sum_{j=1}^{N_\textbf{I}}c_{\textbf{II}j}n_{\textbf{II}j} \leq
c_\textrm{max},\hspace{3.5cm}
\hspace{9.1cm} l_\textbf{I}\leq n_\textbf{I}\leq u_\textbf{I},
l_\textbf{II}\leq n_\textbf{II}\leq u_\textbf{II} \Big\}
\Leftrightarrow \quad \min_{n, z, t}\{t: Dz-d\leq t, Cn\leq c,
1\leq n_iz_i\, \forall i, l\leq n \leq u\}\hspace{4cm}
with
D = (N_\textbf{I}^2S_\textbf{I}^2 -
N_\textbf{I}\sum_{j=1}^{N_\textbf{I}} N_{\textbf{II}j}S_{\textbf{II}j}^2,
N_\textbf{I}N_{\textbf{II}1}^2 S_{\textbf{II}1}^2,\dots,
N_\textbf{I}N_{\textbf{II}N_\textbf{I}}^2S_{\textbf{II}N_\textbf{I}}^2),
d = N_\textbf{I}S_\textbf{I}^2,
C = [C_1, C_2], where
C_1=(c_\textbf{I},l_\textbf{II}^\top,-u_\textbf{II}^\top)^\top and
C_2 = [N_\textbf{I}^{-1}c_\textbf{II}^\top;-\textbf{I};\textbf{I}]
(\textbf{I} is the identity matrix),
c = (c_\textrm{max},0,\dots,0)^\top,
where l = (l_\textbf{I},l_\textbf{I}l_\textbf{II}^\top)^\top and
u = (u_\textbf{I},u_\textbf{I}u_\textbf{II}^\top)^\top
(cf. Willems, 2025, Chapter 3).
Here, N_\textbf{I} is the number of
clusters, N_{\textbf{II}j} the size of cluster j,
S_\textbf{I}^2 the between cluster variance and
S_{\textbf{II}j}^2 the within cluster variances
of the variable of interest. Furthermore, c_\textrm{max} is a
maximum expected cost, c_\textbf{I} a variable cost for sampling
one cluster, and c_{\textbf{II}j} a variable cost for sampling
one unit in cluster j. The optimal number of clusters to be drawn
and the optimal sample sizes are given through n = (n_\textbf{I},
n_\textbf{I}n_{\textbf{II}1},\dots, n_\textbf{I}n_{\textbf{II}
N_\textbf{I}})^\top.
For the special cases above, solutions are unique and, thus, Pareto optimal.
For the general multiobjective problem formulation this is not the case.
However, a strong indicator for uniqueness of solutions is
n_iz_i = 1\, \forall i (Qbounds) or
Dz-d = w^{-1}t (Dbounds). Uniqueness can be ensured via a
stepwise procedure implemented in mosallocStepwiseFirst().
Precision optimization (opts$sense == "max_precision",
opts$f ==f, opts$f == \nabla f,
opts$f == Hf)
The mathematical problem solved is
\min_{n, z}\{f(Dz-d): Az\leq a, Cn\leq c,
1\leq n_iz_i\, \forall i, l\leq n \leq u\}
with components as specified above and where
f:\mathbb{R}^{k_D}\rightarrow \mathbb{R}, x \mapsto f(x)
is a twice-differentiable convex decision functional.
E.g. a p-norm f(x) = \lVert x \rVert_p with
p\in\mathbb{N}.
Cost optimization (opts$sense == "min_cost")
The mathematical problem solved is
\min_{n, z, t}\{t: Dn-d\leq\textbf{1}t, Az\leq a, Cn\leq c,
1\leq n_iz_i\, \forall i, l\leq n \leq u\}
with 1\leq l\leq u \leq N. Hence, the only difference to
precision optimization is the type of objective constraint
Dn-d\leq\textbf{1}t.
Special cases of this formulation are
the cost optimal allocation (possibly multivariate, i.e. k_A\geq 2):
\min_{n, z, t}\{t: Dn-d\leq\textbf{1}t, Az\leq a, 1\leq n_iz_i\,
\forall i\}
where D^\top is a vector of stratum-specific sampling
cost and d some fixed cost.
References
See:
Folks, J.L., Antle, C.E. (1965). Optimum Allocation of Sampling Units to
Strata when there are R Responses of Interest. Journal of the American
Statistical Association, 60(309), 225-233.
\Sexpr[results=rd]{tools:::Rd_expr_doi("10.1080/01621459.1965.10480786")}.
Münnich, R., Sachs, E., Wagner, M. (2012). Numerical solution of optimal
allocation problems in stratified sampling under box constraints.
AStA Advances in Statistical Analysis, 96, 435-450.
\Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/s10182-011-0176-z")}.
Neyman, J. (1934). On the Two Different Aspects of the Representative
Method: The Method of Stratified Sampling and the Method of Purposive
Selection. Journal of the Royal Statistical Society, 97(4), 558–625.
Tschuprow, A.A. (1923). On the Mathematical Expectation of the Moments of
Frequency Distribution in the Case of Correlated Observations. Metron,
2(3,4), 461-493, 646-683.
Rupp, M. (2018). Optimization for Multivariate and Multi-domain Methods
in Survey Statistics (Doctoral dissertation). Trier University.
\Sexpr[results=rd]{tools:::Rd_expr_doi("10.25353/UBTR-8351-5432-14XX")}.
Srikantan, K.S. (1963). A Problem in Optimum Allocation.
Operations Research, 11(2), 265-274.
Willems, F. (2025). A Framework for Multiobjective and Uncertain Resource
Allocation Problems in Survey Sampling based on Conic Optimization
(Doctoral dissertation). Trier University.
\Sexpr[results=rd]{tools:::Rd_expr_doi("10.25353/ubtr-9200-484c-5c89")}.
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
# Examples
#----------------------------------------------------------------------------
# Example 1: Minimization of the variation of estimates for revenue subject
# to cost restrictions and precision restrictions to the coefficient of
# variation of estimates for the proportion of businesses with offshore
# affiliates.
l <- rep(n.min, 5) # minimum sample size per stratum
u <- Nh # maximum sample size per stratum
C <- rbind(ch,
ch * c(-1, -1, -1, 0, 0))
c <- c(budget, # Maximum overall survey budget
- 0.5 * budget) # Minimum overall budget for strata 1-3
# We require at maximum 5 % relative standard error for estimates of
# proportion of businesses with offshore affiliates
A <- matrix(ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2,
nrow = 1)
a <- sum(ph.offsh * (1 - ph.offsh) * Nh**2/(Nh - 1)
)/sum(Nh * ph.offsh)**2 + 0.05**2
D <- matrix(Sh.rev**2 * Nh**2, nrow = 1) # objective variance components
d <- sum(Sh.rev**2 * Nh) # finite population correction
opts <- list(sense = "max_precision",
f = NULL, df = NULL, Hf = NULL,
init_w = 1,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 100L, print_pm = FALSE)
sol <- mosalloc(D = D, d = d, A = A, a = a, C = C, c = c, l = l, u = u,
opts = opts)
# Check solution statement of the internal solver to verify feasibility
sol$Ecosolver$Ecoinfostring # [1] "Optimal solution found"
# Check constraints
c(C[1, ] %*% sol$n) # [1] 3e+05
c(C[2, ] %*% sol$n) # [1] -150000
c(sqrt(A %*% (1 / sol$n) - A %*% (1 / Nh))) # 5 % rel. std. err.
#----------------------------------------------------------------------------
# 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 cost restrictions
l <- rep(n.min, 5) # minimum sample size ber stratum
u <- Nh # maximum sample size per stratum
C <- rbind(ch, ch * c(-1, -1, -1, 0, 0))
c <- c(budget, - 0.5 * budget)
A <- NULL # no precision constraint
a <- NULL # no precision constraint
# Precision components (Variance / Totals^2) for multidimensional objective
D <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
d <- as.vector(D %*% (1 / Nh)) # finite population correction
opts <- list(sense = "max_precision",
f = NULL, df = NULL, Hf = NULL,
init_w = 1,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 100L, print_pm = FALSE)
sol <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
# Obtain optimal objective value
sol$J # [1] 0.0017058896 0.0004396972 0.0006428475 0.0017058896
# Obtain corresponding normal vector
sol$Normal # [1] 6.983113e-01 1.337310e-11 1.596167e-11 3.016887e-01
# => Revenue and offshore affiliates are dominating the solution with a
# ratio of approximately 2:1 (sol$Normal[1] / sol$Normal[4])
#----------------------------------------------------------------------------
# Example 3: Example 2 with preference weighting
w <- c(1, 3.85, 3.8, 1.3) # preference weighting
l <- rep(n.min, 5) # minimum sample size ber stratum
u <- Nh # maximum sample size per stratum
C <- rbind(ch, ch * c(-1, -1, -1, 0, 0))
c <- c(budget, - 0.5 * budget)
A <- NULL # no precision constraint
a <- NULL # no precision constraint
D <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
d <- as.vector(D %*% (1 / Nh))
opts <- list(sense = "max_precision",
f = NULL, df = NULL, Hf = NULL,
init_w = w,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 100L, print_pm = FALSE)
mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
#----------------------------------------------------------------------------
# Example 4: Example 2 with multiple preference weightings for simultaneous
# evaluation
w <- matrix(c(1.0, 1.0, 1.0, 1.0, # matrix of preference weightings
1.0, 3.9, 3.9, 1.3,
0.8, 4.2, 4.8, 1.5,
1.2, 3.5, 4.8, 2.0,
2.0, 1.0, 1.0, 2.0), 5, 4, byrow = TRUE)
w <- w / w[,1] # rescale w (ensure the first weighting to be one)
l <- rep(n.min, 5) # minimum sample size ber stratum
u <- Nh # maximum sample size per stratum
C <- rbind(ch, ch * c(-1, -1, -1, 0, 0))
c <- c(budget, - 0.5 * budget)
A <- NULL # no precision constraint
a <- NULL # no precision constraint
D <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
d <- as.vector(D %*% (1 / Nh))
opts <- list(sense = "max_precision",
f = NULL, df = NULL, Hf = NULL,
init_w = w,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 100L, print_pm = FALSE)
sols <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
lapply(sols, function(sol){sol$Qbounds})
#----------------------------------------------------------------------------
# Example 5: Example 2 where a weighted sum scalarization of the objective
# components is minimized
l <- rep(n.min, 5) # minimum sample size ber stratum
u <- Nh # maximum sample size per stratum
C <- matrix(ch, nrow = 1)
c <- budget
A <- NULL # no precision constraint
a <- NULL # no precision constraint
# Objective variance components
D <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
d <- as.vector(D %*% (1 / Nh)) # finite population correction
# Simple weighted sum as decision functional
wss <- c(1, 1, 0.5, 0.5) # preference weighting (weighted sum scalarization)
Dw <- wss %*% D
dw <- as.vector(Dw %*% (1 / Nh))
opts <- list(sense = "max_precision",
f = NULL, df = NULL, Hf = NULL,
init_w = 1,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 1000L, print_pm = FALSE)
# Solve weighted sum scalarization (WSS) via mosalloc
sol_wss <- mosalloc(D = Dw, d = dw, C = C, c = c, l = l, u = u, opts = opts)
# Obtain optimal objective values
J <- D %*% (1 / sol_wss$n) - d
# Reconstruct solution via a weighted Chebyshev minimization
wcm <- J[1] / J
opts = list(sense = "max_precision",
f = NULL, df = NULL, Hf = NULL,
init_w = matrix(wcm, 1),
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 1000L, print_pm = FALSE)
sol_wcm <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
# Compare solutions
rbind(t(J), sol_wcm$J)
# [,1] [,2] [,3] [,4]
# [1,] 0.00155645 0.0004037429 0.0005934474 0.001327165
# [2,] 0.00155645 0.0004037429 0.0005934474 0.001327165
rbind(sol_wss$n, sol_wcm$n)
# [,1] [,2] [,3] [,4] [,5]
# [1,] 582.8247 236.6479 116.7866 839.5988 841.4825
# [2,] 582.8226 236.6475 116.7871 839.5989 841.4841
rbind(wss, sol_wcm$Normal / sol_wcm$Normal[1])
# [,1] [,2] [,3] [,4]
#wss 1 1.0000000 0.5000000 0.5000000
# 1 0.9976722 0.4997552 0.4997462
#----------------------------------------------------------------------------
# Example 6: Example 1 with two subpopulations and a p-norm as decision
# functional
l <- rep(n.min, 5) # minimum sample size per stratum
u <- Nh # maximum sample size per stratum
C <- rbind(ch, ch * c(-1, -1, -1, 0, 0))
c <- c(budget, - 0.5 * budget)
# At maximum 5 % relative standard error for estimates of proportion of
# businesses with offshore affiliates
A <- matrix(ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2,
nrow = 1)
a <- sum(ph.offsh * (1 - ph.offsh) * Nh**2/(Nh - 1)
)/sum(Nh * ph.offsh)**2 + 0.05**2
D <- rbind((Sh.rev**2 * Nh**2)*c(0,0,1,1,0),
(Sh.rev**2 * Nh**2)*c(1,1,0,0,1))# objective variance components
d <- as.vector(D %*% (1 / Nh)) # finite population correction
# p-norm solution
p <- 5 # p-norm
opts <- list(sense = "max_precision",
f = function(x) sum(x**p),
df = function(x) p * x**(p - 1),
Hf = function(x) diag(p * (p - 1) * x**(p - 2)),
init_w = 1,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 1000L, print_pm = TRUE)
sol <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
c(sol$Normal/sol$dfJ)/mean(c(sol$Normal/sol$dfJ))
# [1] 0.9999972 1.0000028
#----------------------------------------------------------------------------
# Example 7: Example 2 with p-norm as decision functional and only one
# overall cost constraint
l <- rep(n.min, 5) # minimum sample size ber stratum
u <- Nh # maximum sample size per stratum
C <- matrix(ch, nrow = 1)
c <- budget
A <- NULL # no precision constraint
a <- NULL # no precision constraint
# Objective precision components
D <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
d <- as.vector(D %*% (1 / Nh)) # finite population correction
# p-norm solution
p <- 5 # p-norm
opts <- list(sense = "max_precision",
f = function(x) sum(x**p),
df = function(x) p * x**(p - 1),
Hf = function(x) diag(p * (p - 1) * x**(p - 2)),
init_w = 1,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 1000L, print_pm = TRUE)
sol <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
c(sol$Normal/sol$dfJ)/mean(c(sol$Normal/sol$dfJ))
# [1] 1.0014362 0.9780042 1.0197807 1.0007789
#----------------------------------------------------------------------------
# Example 8: Minimization of sample sizes subject to precision constraints
l <- rep(n.min, 5) # minimum sample size ber stratum
u <- Nh # maximum sample size per stratum
# We require at maximum 4.66 % relative standard error for the estimate of
# total revenuee, 5 % for the number of employees, 3 % for the proportion of
# businesses claiming research credit, and 3 % for the proportion of
# businesses with offshore affiliates
A <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
a <- as.vector(A%*%(1 / Nh) + c(0.0466, 0.05, 0.03, 0.03)**2)
# We do not consider any additional sample size or cost constraints
C <- NULL # no cost constraint
c <- NULL # no cost constraint
# Since we minimize the sample size, we define D and d as follows:
D <- matrix(1, nrow = 1, ncol = length(Nh)) # objective cost components
d <- as.vector(0) # vector of possible fixed cost
opts <- list(sense = "min_cost", # Sense of optimization is survey cost
f = NULL,
df = NULL,
Hf = NULL,
init_w = 1,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 100L, print_pm = TRUE)
sol <- mosalloc(D = D, d = d, A = A, a = a, l = l, u = u, opts = opts)
sum(sol$n) # [1] 2843.219
sol$J # [1] 2843.219
#----------------------------------------------------------------------------
#----------------------------------------------------------------------------
# Note: Sample size optimization for two-stage cluster sampling can be
# reduced to the structure of optimal stratified random samplin when
# considering expected costs. Therefore, mosalloc() can handle such
# designs. A benefit is that mosalloc() allows relatively complex
# sample size restrictions such as box constraints for subsampling.
# Optimal sample sizes at secondary stages have to be reconstructed
# from sol$n.
#
# Example 9: Optimal number of primary sampling units (PSU) and secondary
# sampling units (SSU) in 2-stage cluster sampling.
set.seed(1234)
pop <- data.frame(value = rnorm(100, 100, 35),
cluster = sample(1:4, 100, replace = TRUE))
CI <- 36 # Sampling cost per PSU
CII <- 10 # Average sampling cost per SSU
NI <- 4 # Number of PSUs
NII <- table(pop$cluster) # PSU/cluster sizes
S2I <- var(by(pop$value, pop$cluster, sum)) # between PSU variance
S2II <- by(pop$value, pop$cluster, var) # within PSU variances
D <- matrix(c(NI**2 * S2I - NI * sum(NII * S2II), NI * NII**2 * S2II), 1)
d <- as.vector(NI * S2I) # = D %*% (1 / c(NI, NI * NII))
C <- cbind(c(CI, rep(2, NI), -NII),
rbind(rep(CII / NI, 4), -diag(4), diag(4)))
c <- as.vector(c(500, rep(0, 8)))
l <- c(2, rep(4, 4))
u <- c(NI, NI * NII)
opts <- list(sense = "max_precision",
f = NULL,
df = NULL,
Hf = NULL,
init_w = 1,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 100L, print_pm = TRUE)
sol <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
# Optimum number of clusters to be drawn
sol$n[1] # [1] 2.991551
# Optimum number of elements to be drawn within clusters
sol$n[-1] / sol$n[1] # [1] 12.16454 11.60828 15.87949 12.80266
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.
| mosalloc | R Documentation |
Multiobjective sample allocation for constraint multivariate and multidomain optimal allocation in survey sampling
Description
Computes solutions to standard sample allocation problems
under various precision and cost restrictions. The input data is
transformed and parsed to the Embedded COnic Solver (ECOS) from the
'ECOSolveR' package. Multiple survey purposes can be optimized
simultaneously through a weighted Chebyshev minimization. Note that in
the case of multiple objectives, mosalloc() does not necessarily
lead to Pareto optimality. This highly depends on the problem structure.
A strong indicator for Pareto optimality is when the weighted objective
values given by Dbounds are constant over all objective
components or when all components of Qbounds equal 1.
In addition, mosalloc() can handle twice-differential convex
decision functionals (in which case Pareto optimality is ensured).
mosalloc() returns dual variables, enabling a detailed sensitivity
analysis.
Usage
mosalloc(
D,
d,
A = NULL,
a = NULL,
C = NULL,
c = NULL,
l = 2,
u = NULL,
opts = list(sense = "max_precision", f = NULL, df = NULL, Hf = NULL, init_w = 1,
mc_cores = 1L, pm_tol = 1e-05, max_iters = 100L, print_pm = FALSE)
)
Arguments
D |
(type: |
d |
(type: |
A |
(type: |
a |
(type: |
C |
(type: |
c |
(type: |
l |
(type: |
u |
(type: |
opts |
(type: |
Value
mosalloc() returns a list containing the following
components:
$w The initial preference weighting opts$init_w.
$n The vector of optimal sample sizes.
$J The optimal objective vector.
$Objective The objective value with respect to decision
functional f. NULL if opts$f = NULL.
$Utopian The component-wise univariate optimal
objective vector. NULL if opts$f = NULL.
$Normal The vector normal to the Pareto frontier at
$J.
$dfJ The gradient of opts$f evaluated at $J.
$Sensitivity The dual variables of the objectives and
constraints.
$Qbounds The Quality bounds of the Lorentz cones.
$Dbounds The weighted objective constraints ($w * $J).
$Scalepar An internal scaling parameter.
$Ecosolver A list of ECOSolveR returns including:
...$Ecoinfostring The info string of
ECOSolveR::ECOS_csolve().
...$Ecoredcodes The redcodes of ECOSolveR::ECOS_csolve().
...$Ecosummary Problem summary of ECOSolveR::ECOS_csolve().
$Timing Run time info.
$Iteration A list of internal iterates. NULL if
opts$f = NULL.
Note
Precision optimization (opts$sense == "max_precision",
opts$f == NULL)
The mathematical problem solved is
\min_{n, z, t}\{t: Dz-d\leq w^{-1}t, Az\leq a, Cn\leq c,
1\leq n_iz_i\, \forall i, l\leq n \leq u\}
with \textbf{1}\leq l\leq u \leq N, where
-
n,z\in\mathbb{R}^mandt\in\mathbb{R}are the optimization variables, -
D\in\{M\in\mathbb{R}^{k_D\times m}: \sum_i M_{ij}>0\, \forall j, \sum_j M_{ij}>0\, \forall i\}a matrix of nonnegative objective precision components withk_Dthe number of precision objectives (the number of variables of interest), -
d\in\mathbb{R}^{k_D}the objective right hand-side (RHS), e.g. the finite population correction (fpc), -
A\in\mathbb{R}^{k_A\times m}a nonnegative precision matrix withk_Athe number of precision constraints,a\in \mathbb{R}^{k_A}the corresponding RHS, e.g. fpc + (upper bound to the coefficient of variation)^2, -
C\in\mathbb{R}^{k_C\times m}a cost matrix withk_Cthe number of cost constraint, -
c\in\mathbb{R}^{k_C}the corresponding RHS, -
l,u\in\mathbb{R}^mthe bounds to the sample size vectorn, -
N\in\mathbb{N}^mthe vector of population sizes, and -
wa given strictly positive preference weighting.
Special cases of this formulation are
Neyman-Tschuprow allocation (Neyman, 1934 and Tschuprow, 1923):
\min_n \Big\{\sum_{h=1}^H \Big(\frac{N_h^2S_h^2}{n_h}-N_hS_h^2\Big):\sum_{h=1}^H n_h \leq c\Big\} \quad \Leftrightarrow \quad \min_{n, z, t}\{t: Dz-d\leq t, Cn\leq c, 1\leq n_iz_i\, \forall i\}with
D = (N_1^2S_1^2,\dots,N_H^2S_H^2),d = \sum_{h=1}^H N_hS_h^2,C = (1,\dots,1)andca maximum sample size. Here,His the number of strata,N_hthe size of stratumhandS_h^2the variance of the variable of interest in stratumh.box-constrained optimal allocation:
\min_{n, z, t}\{t: Dz-d\leq t, Cn\leq c, 1\leq n_iz_i\, \forall i, l\leq n \leq u\}with
D = (N_1^2S_1^2,\dots,N_H^2S_H^2),d = \sum_{h=1}^H N_hS_h^2,C = (1,\dots,1)andca maximum sample size (cf. Srikantan, 1963 and Münnich et al., 2012). Here,landuare bounds to the optimal sample size vector.cost and precision constrained univariate optimal allocation:
\min_{n, z, t}\{t: Dz-d\leq t, Az\leq a, Cn\leq c, 1\leq n_iz_i\, \forall i, l\leq n \leq u\}with
k_D = 1(cf. Willems, 2025, Chapter 3).multivariate optimal allocation with weighted sum scalarization:
\min_{n, z, t}\{t: w^\top Dz-w^\top d\leq t, Az\leq a, Cn\leq c, 1\leq n_iz_i\, \forall i, l\leq n \leq u\}where
w\in\mathbb{R}^{k_D}is a strictly positive preference weighting (cf. Folks and Antle, 1965, and Rupp, 2018). Note that for this case the problem reduces to cost and precision constrained univariate optimal allocation. Solutions are ensured to be optimal in the Pareto sense.box-constraint two-stage cluster sampling:
\quad \quad \quad\min_{n_\textbf{I},n_\textbf{II}} \Big\{\Big(\frac{N_\textbf{I}^2 S_\textbf{I}^2}{n_\textbf{I}}-N_\textbf{I}S_\textbf{I}^2\Big) + \frac{N_\textbf{I}}{n_\textbf{I}}\sum_{j=1}^{N_\textbf{I}} \Big(\frac{N_{\textbf{II}j}^2S_{\textbf{II}j}^2}{n_{\textbf{II}j}}- N_{\textbf{II}j}S_{\textbf{II}j}^2\Big): c_{\textbf{I}}n_\textbf{I} + \frac{n_\textbf{I}}{N_\textbf{I}} \sum_{j=1}^{N_\textbf{I}}c_{\textbf{II}j}n_{\textbf{II}j} \leq c_\textrm{max},\hspace{3.5cm}\hspace{9.1cm} l_\textbf{I}\leq n_\textbf{I}\leq u_\textbf{I}, l_\textbf{II}\leq n_\textbf{II}\leq u_\textbf{II} \Big\}\Leftrightarrow \quad \min_{n, z, t}\{t: Dz-d\leq t, Cn\leq c, 1\leq n_iz_i\, \forall i, l\leq n \leq u\}\hspace{4cm}with
D = (N_\textbf{I}^2S_\textbf{I}^2 - N_\textbf{I}\sum_{j=1}^{N_\textbf{I}} N_{\textbf{II}j}S_{\textbf{II}j}^2, N_\textbf{I}N_{\textbf{II}1}^2 S_{\textbf{II}1}^2,\dots, N_\textbf{I}N_{\textbf{II}N_\textbf{I}}^2S_{\textbf{II}N_\textbf{I}}^2),d = N_\textbf{I}S_\textbf{I}^2,C = [C_1, C_2], whereC_1=(c_\textbf{I},l_\textbf{II}^\top,-u_\textbf{II}^\top)^\topandC_2 = [N_\textbf{I}^{-1}c_\textbf{II}^\top;-\textbf{I};\textbf{I}](\textbf{I}is the identity matrix),c = (c_\textrm{max},0,\dots,0)^\top, wherel = (l_\textbf{I},l_\textbf{I}l_\textbf{II}^\top)^\topandu = (u_\textbf{I},u_\textbf{I}u_\textbf{II}^\top)^\top(cf. Willems, 2025, Chapter 3). Here,N_\textbf{I}is the number of clusters,N_{\textbf{II}j}the size of clusterj,S_\textbf{I}^2the between cluster variance andS_{\textbf{II}j}^2the within cluster variances of the variable of interest. Furthermore,c_\textrm{max}is a maximum expected cost,c_\textbf{I}a variable cost for sampling one cluster, andc_{\textbf{II}j}a variable cost for sampling one unit in clusterj. The optimal number of clusters to be drawn and the optimal sample sizes are given throughn = (n_\textbf{I}, n_\textbf{I}n_{\textbf{II}1},\dots, n_\textbf{I}n_{\textbf{II} N_\textbf{I}})^\top.
For the special cases above, solutions are unique and, thus, Pareto optimal.
For the general multiobjective problem formulation this is not the case.
However, a strong indicator for uniqueness of solutions is
n_iz_i = 1\, \forall i (Qbounds) or
Dz-d = w^{-1}t (Dbounds). Uniqueness can be ensured via a
stepwise procedure implemented in mosallocStepwiseFirst().
Precision optimization (opts$sense == "max_precision",
opts$f ==f, opts$f == \nabla f,
opts$f == Hf)
The mathematical problem solved is
\min_{n, z}\{f(Dz-d): Az\leq a, Cn\leq c,
1\leq n_iz_i\, \forall i, l\leq n \leq u\}
with components as specified above and where
f:\mathbb{R}^{k_D}\rightarrow \mathbb{R}, x \mapsto f(x)
is a twice-differentiable convex decision functional.
E.g. a p-norm f(x) = \lVert x \rVert_p with
p\in\mathbb{N}.
Cost optimization (opts$sense == "min_cost")
The mathematical problem solved is
\min_{n, z, t}\{t: Dn-d\leq\textbf{1}t, Az\leq a, Cn\leq c,
1\leq n_iz_i\, \forall i, l\leq n \leq u\}
with 1\leq l\leq u \leq N. Hence, the only difference to
precision optimization is the type of objective constraint
Dn-d\leq\textbf{1}t.
Special cases of this formulation are
the cost optimal allocation (possibly multivariate, i.e.
k_A\geq 2):\min_{n, z, t}\{t: Dn-d\leq\textbf{1}t, Az\leq a, 1\leq n_iz_i\, \forall i\}where
D^\topis a vector of stratum-specific sampling cost anddsome fixed cost.
References
See:
Folks, J.L., Antle, C.E. (1965). Optimum Allocation of Sampling Units to Strata when there are R Responses of Interest. Journal of the American Statistical Association, 60(309), 225-233. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1080/01621459.1965.10480786")}.
Münnich, R., Sachs, E., Wagner, M. (2012). Numerical solution of optimal allocation problems in stratified sampling under box constraints. AStA Advances in Statistical Analysis, 96, 435-450. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/s10182-011-0176-z")}.
Neyman, J. (1934). On the Two Different Aspects of the Representative Method: The Method of Stratified Sampling and the Method of Purposive Selection. Journal of the Royal Statistical Society, 97(4), 558–625.
Tschuprow, A.A. (1923). On the Mathematical Expectation of the Moments of Frequency Distribution in the Case of Correlated Observations. Metron, 2(3,4), 461-493, 646-683.
Rupp, M. (2018). Optimization for Multivariate and Multi-domain Methods in Survey Statistics (Doctoral dissertation). Trier University. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.25353/UBTR-8351-5432-14XX")}.
Srikantan, K.S. (1963). A Problem in Optimum Allocation. Operations Research, 11(2), 265-274.
Willems, F. (2025). A Framework for Multiobjective and Uncertain Resource Allocation Problems in Survey Sampling based on Conic Optimization (Doctoral dissertation). Trier University. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.25353/ubtr-9200-484c-5c89")}.
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
# Examples
#----------------------------------------------------------------------------
# Example 1: Minimization of the variation of estimates for revenue subject
# to cost restrictions and precision restrictions to the coefficient of
# variation of estimates for the proportion of businesses with offshore
# affiliates.
l <- rep(n.min, 5) # minimum sample size per stratum
u <- Nh # maximum sample size per stratum
C <- rbind(ch,
ch * c(-1, -1, -1, 0, 0))
c <- c(budget, # Maximum overall survey budget
- 0.5 * budget) # Minimum overall budget for strata 1-3
# We require at maximum 5 % relative standard error for estimates of
# proportion of businesses with offshore affiliates
A <- matrix(ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2,
nrow = 1)
a <- sum(ph.offsh * (1 - ph.offsh) * Nh**2/(Nh - 1)
)/sum(Nh * ph.offsh)**2 + 0.05**2
D <- matrix(Sh.rev**2 * Nh**2, nrow = 1) # objective variance components
d <- sum(Sh.rev**2 * Nh) # finite population correction
opts <- list(sense = "max_precision",
f = NULL, df = NULL, Hf = NULL,
init_w = 1,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 100L, print_pm = FALSE)
sol <- mosalloc(D = D, d = d, A = A, a = a, C = C, c = c, l = l, u = u,
opts = opts)
# Check solution statement of the internal solver to verify feasibility
sol$Ecosolver$Ecoinfostring # [1] "Optimal solution found"
# Check constraints
c(C[1, ] %*% sol$n) # [1] 3e+05
c(C[2, ] %*% sol$n) # [1] -150000
c(sqrt(A %*% (1 / sol$n) - A %*% (1 / Nh))) # 5 % rel. std. err.
#----------------------------------------------------------------------------
# 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 cost restrictions
l <- rep(n.min, 5) # minimum sample size ber stratum
u <- Nh # maximum sample size per stratum
C <- rbind(ch, ch * c(-1, -1, -1, 0, 0))
c <- c(budget, - 0.5 * budget)
A <- NULL # no precision constraint
a <- NULL # no precision constraint
# Precision components (Variance / Totals^2) for multidimensional objective
D <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
d <- as.vector(D %*% (1 / Nh)) # finite population correction
opts <- list(sense = "max_precision",
f = NULL, df = NULL, Hf = NULL,
init_w = 1,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 100L, print_pm = FALSE)
sol <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
# Obtain optimal objective value
sol$J # [1] 0.0017058896 0.0004396972 0.0006428475 0.0017058896
# Obtain corresponding normal vector
sol$Normal # [1] 6.983113e-01 1.337310e-11 1.596167e-11 3.016887e-01
# => Revenue and offshore affiliates are dominating the solution with a
# ratio of approximately 2:1 (sol$Normal[1] / sol$Normal[4])
#----------------------------------------------------------------------------
# Example 3: Example 2 with preference weighting
w <- c(1, 3.85, 3.8, 1.3) # preference weighting
l <- rep(n.min, 5) # minimum sample size ber stratum
u <- Nh # maximum sample size per stratum
C <- rbind(ch, ch * c(-1, -1, -1, 0, 0))
c <- c(budget, - 0.5 * budget)
A <- NULL # no precision constraint
a <- NULL # no precision constraint
D <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
d <- as.vector(D %*% (1 / Nh))
opts <- list(sense = "max_precision",
f = NULL, df = NULL, Hf = NULL,
init_w = w,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 100L, print_pm = FALSE)
mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
#----------------------------------------------------------------------------
# Example 4: Example 2 with multiple preference weightings for simultaneous
# evaluation
w <- matrix(c(1.0, 1.0, 1.0, 1.0, # matrix of preference weightings
1.0, 3.9, 3.9, 1.3,
0.8, 4.2, 4.8, 1.5,
1.2, 3.5, 4.8, 2.0,
2.0, 1.0, 1.0, 2.0), 5, 4, byrow = TRUE)
w <- w / w[,1] # rescale w (ensure the first weighting to be one)
l <- rep(n.min, 5) # minimum sample size ber stratum
u <- Nh # maximum sample size per stratum
C <- rbind(ch, ch * c(-1, -1, -1, 0, 0))
c <- c(budget, - 0.5 * budget)
A <- NULL # no precision constraint
a <- NULL # no precision constraint
D <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
d <- as.vector(D %*% (1 / Nh))
opts <- list(sense = "max_precision",
f = NULL, df = NULL, Hf = NULL,
init_w = w,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 100L, print_pm = FALSE)
sols <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
lapply(sols, function(sol){sol$Qbounds})
#----------------------------------------------------------------------------
# Example 5: Example 2 where a weighted sum scalarization of the objective
# components is minimized
l <- rep(n.min, 5) # minimum sample size ber stratum
u <- Nh # maximum sample size per stratum
C <- matrix(ch, nrow = 1)
c <- budget
A <- NULL # no precision constraint
a <- NULL # no precision constraint
# Objective variance components
D <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
d <- as.vector(D %*% (1 / Nh)) # finite population correction
# Simple weighted sum as decision functional
wss <- c(1, 1, 0.5, 0.5) # preference weighting (weighted sum scalarization)
Dw <- wss %*% D
dw <- as.vector(Dw %*% (1 / Nh))
opts <- list(sense = "max_precision",
f = NULL, df = NULL, Hf = NULL,
init_w = 1,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 1000L, print_pm = FALSE)
# Solve weighted sum scalarization (WSS) via mosalloc
sol_wss <- mosalloc(D = Dw, d = dw, C = C, c = c, l = l, u = u, opts = opts)
# Obtain optimal objective values
J <- D %*% (1 / sol_wss$n) - d
# Reconstruct solution via a weighted Chebyshev minimization
wcm <- J[1] / J
opts = list(sense = "max_precision",
f = NULL, df = NULL, Hf = NULL,
init_w = matrix(wcm, 1),
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 1000L, print_pm = FALSE)
sol_wcm <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
# Compare solutions
rbind(t(J), sol_wcm$J)
# [,1] [,2] [,3] [,4]
# [1,] 0.00155645 0.0004037429 0.0005934474 0.001327165
# [2,] 0.00155645 0.0004037429 0.0005934474 0.001327165
rbind(sol_wss$n, sol_wcm$n)
# [,1] [,2] [,3] [,4] [,5]
# [1,] 582.8247 236.6479 116.7866 839.5988 841.4825
# [2,] 582.8226 236.6475 116.7871 839.5989 841.4841
rbind(wss, sol_wcm$Normal / sol_wcm$Normal[1])
# [,1] [,2] [,3] [,4]
#wss 1 1.0000000 0.5000000 0.5000000
# 1 0.9976722 0.4997552 0.4997462
#----------------------------------------------------------------------------
# Example 6: Example 1 with two subpopulations and a p-norm as decision
# functional
l <- rep(n.min, 5) # minimum sample size per stratum
u <- Nh # maximum sample size per stratum
C <- rbind(ch, ch * c(-1, -1, -1, 0, 0))
c <- c(budget, - 0.5 * budget)
# At maximum 5 % relative standard error for estimates of proportion of
# businesses with offshore affiliates
A <- matrix(ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2,
nrow = 1)
a <- sum(ph.offsh * (1 - ph.offsh) * Nh**2/(Nh - 1)
)/sum(Nh * ph.offsh)**2 + 0.05**2
D <- rbind((Sh.rev**2 * Nh**2)*c(0,0,1,1,0),
(Sh.rev**2 * Nh**2)*c(1,1,0,0,1))# objective variance components
d <- as.vector(D %*% (1 / Nh)) # finite population correction
# p-norm solution
p <- 5 # p-norm
opts <- list(sense = "max_precision",
f = function(x) sum(x**p),
df = function(x) p * x**(p - 1),
Hf = function(x) diag(p * (p - 1) * x**(p - 2)),
init_w = 1,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 1000L, print_pm = TRUE)
sol <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
c(sol$Normal/sol$dfJ)/mean(c(sol$Normal/sol$dfJ))
# [1] 0.9999972 1.0000028
#----------------------------------------------------------------------------
# Example 7: Example 2 with p-norm as decision functional and only one
# overall cost constraint
l <- rep(n.min, 5) # minimum sample size ber stratum
u <- Nh # maximum sample size per stratum
C <- matrix(ch, nrow = 1)
c <- budget
A <- NULL # no precision constraint
a <- NULL # no precision constraint
# Objective precision components
D <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
d <- as.vector(D %*% (1 / Nh)) # finite population correction
# p-norm solution
p <- 5 # p-norm
opts <- list(sense = "max_precision",
f = function(x) sum(x**p),
df = function(x) p * x**(p - 1),
Hf = function(x) diag(p * (p - 1) * x**(p - 2)),
init_w = 1,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 1000L, print_pm = TRUE)
sol <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
c(sol$Normal/sol$dfJ)/mean(c(sol$Normal/sol$dfJ))
# [1] 1.0014362 0.9780042 1.0197807 1.0007789
#----------------------------------------------------------------------------
# Example 8: Minimization of sample sizes subject to precision constraints
l <- rep(n.min, 5) # minimum sample size ber stratum
u <- Nh # maximum sample size per stratum
# We require at maximum 4.66 % relative standard error for the estimate of
# total revenuee, 5 % for the number of employees, 3 % for the proportion of
# businesses claiming research credit, and 3 % for the proportion of
# businesses with offshore affiliates
A <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
a <- as.vector(A%*%(1 / Nh) + c(0.0466, 0.05, 0.03, 0.03)**2)
# We do not consider any additional sample size or cost constraints
C <- NULL # no cost constraint
c <- NULL # no cost constraint
# Since we minimize the sample size, we define D and d as follows:
D <- matrix(1, nrow = 1, ncol = length(Nh)) # objective cost components
d <- as.vector(0) # vector of possible fixed cost
opts <- list(sense = "min_cost", # Sense of optimization is survey cost
f = NULL,
df = NULL,
Hf = NULL,
init_w = 1,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 100L, print_pm = TRUE)
sol <- mosalloc(D = D, d = d, A = A, a = a, l = l, u = u, opts = opts)
sum(sol$n) # [1] 2843.219
sol$J # [1] 2843.219
#----------------------------------------------------------------------------
#----------------------------------------------------------------------------
# Note: Sample size optimization for two-stage cluster sampling can be
# reduced to the structure of optimal stratified random samplin when
# considering expected costs. Therefore, mosalloc() can handle such
# designs. A benefit is that mosalloc() allows relatively complex
# sample size restrictions such as box constraints for subsampling.
# Optimal sample sizes at secondary stages have to be reconstructed
# from sol$n.
#
# Example 9: Optimal number of primary sampling units (PSU) and secondary
# sampling units (SSU) in 2-stage cluster sampling.
set.seed(1234)
pop <- data.frame(value = rnorm(100, 100, 35),
cluster = sample(1:4, 100, replace = TRUE))
CI <- 36 # Sampling cost per PSU
CII <- 10 # Average sampling cost per SSU
NI <- 4 # Number of PSUs
NII <- table(pop$cluster) # PSU/cluster sizes
S2I <- var(by(pop$value, pop$cluster, sum)) # between PSU variance
S2II <- by(pop$value, pop$cluster, var) # within PSU variances
D <- matrix(c(NI**2 * S2I - NI * sum(NII * S2II), NI * NII**2 * S2II), 1)
d <- as.vector(NI * S2I) # = D %*% (1 / c(NI, NI * NII))
C <- cbind(c(CI, rep(2, NI), -NII),
rbind(rep(CII / NI, 4), -diag(4), diag(4)))
c <- as.vector(c(500, rep(0, 8)))
l <- c(2, rep(4, 4))
u <- c(NI, NI * NII)
opts <- list(sense = "max_precision",
f = NULL,
df = NULL,
Hf = NULL,
init_w = 1,
mc_cores = 1L, pm_tol = 1e-05,
max_iters = 100L, print_pm = TRUE)
sol <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
# Optimum number of clusters to be drawn
sol$n[1] # [1] 2.991551
# Optimum number of elements to be drawn within clusters
sol$n[-1] / sol$n[1] # [1] 12.16454 11.60828 15.87949 12.80266
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.