Nothing
Optimum Sample Allocation in Stratified Sampling with `stratallo`
In stratallo: Optimum Sample Allocation in Stratified Sampling
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
\def\R{{\mathbb R}}
\def\x{\mathbf x}
\def\n{\mathbf n}
\newcommand{\texteq}{\mathrm}
The goal of stratallo package is to provide implementations of the efficient
algorithms that solve a classical problem in survey methodology - an optimum
sample allocation in stratified sampling. In this context, the classical problem
of optimum sample allocation is the Tschuprow-Neyman's sense [@Tschuprow; @neyman].
It is formulated as determination of a vector of strata sample sizes that
minimizes the variance of the stratified $\pi$ estimator of the population
total of a given study variable, under constraint on total sample size. More
specifically, the algorithms provided in this package are valid given that the
variance of the stratified estimator is of the following generic form:
$$
V_{st}(\n) = \sum_{h=1}^{H} \frac{A_h^2}{n_h} - A_0,
$$
where $H$ denotes total number of strata, $\n = (n_h)_{h \in {1,\ldots,H}}$ is
the allocation vector with strata sample sizes, and population parameters
$A_0,\, A_h > 0,\, h = 1,\ldots,H$, do not depend on the $x_h,\, h = 1,\ldots,H$.
The allocation problem mentioned, can be further complemented by imposing lower
or upper bounds on sample sizes is strata.
Among stratified estimators and stratified sampling designs that jointly give
rise to a variance of the above form, is the so called stratified $\pi$
estimator of the population total with stratified simple random sampling
without replacement design, which is one of the most basic and commonly used
stratified sampling designs. This case yields
$A_0 = \sum_{h = 1}^H N_h S_h^2$, $A_h = N_h S_h,\, h = 1,\ldots,H$, where
$S_h$ denotes stratum standard deviation of study variable and $N_h$ is the
stratum size (see e.g. @sarndal, Result 3.7.2, p.103).
A minor modification of the classical optimum sample allocation problem leads
to the minimum cost allocation. This problem lies in the determination of a
vector of strata sample sizes that minimizes total cost of the survey, under
assumed fixed level of the stratified $\pi$ estimator's variance. As in the case
of the classical optimum allocation, the problem of minimum cost allocation can
be complemented by imposing upper bounds on sample sizes in strata.
Package stratallo provides two user functions:
opt()optcost()
that solve sample allocation problems briefly characterized above as well as the
following helpers functions:
var_st()var_st_tsi()asummary()ran_round()round_oric().
Functions var_st() and var_st_tsi() compute a value of the variance $V_{st}$.
The var_st_tsi() is a simple wrapper of var_st() that is dedicated for the
case when $A_0 = \sum_{h = 1}^H N_h S_h^2$ and $A_h = N_h S_h,\, h = 1,\ldots,H$.
asummary() creates a data.frame object with summary of the allocation.
Functions ran_round() and round_oric() are the rounding functions that can
be used to round non-integers allocations (see section Rounding, below).
The package comes with three predefined, artificial populations with 10, 507
and 969 strata. These are stored under pop10_mM, pop507 and pop969
objects, respectively.
Minimization of the variance with opt() function
The opt() function solves the following three problems of the optimum sample
allocation, formulated in the language of mathematical optimization. User of
opt() can choose whether the solution computed will be for
Problem 1, Problem 2 or Problem 3. This is achieved with the proper
use of m and M arguments of the function. Also, if required, the inequality
constraints can be removed from the optimization problem. For more details, see
the help page for opt() function.
Problem 1 (one-sided upper bounds)
Given numbers $n > 0,\, A_h > 0,\, M_h > 0$, such that
$M_h \leq N_h,\, h = 1,\ldots,H$, and $n \leq \sum_{h=1}^H M_h$,
\begin{align}
\underset{\x \in \R_+^H}{\texteq{minimize ~\,}} & \quad f(\x) = \sum_{h=1}^H \tfrac{A_h^2}{x_h} \
\texteq{subject ~ to} & \quad \sum_{h=1}^H x_h = n \
& \quad x_h \leq M_h, \quad{h = 1,\ldots,H,}
\end{align}
where $\x = (x_h)_{h \in {1,\ldots,H}}$.
There are four different algorithms available to use for Problem 1,
RNA (default), SGA, SGAPLUS, COMA. All these algorithms, except
SGAPLUS, are described in detail in @wesolowski2021. The SGAPLUS is defined
in @wojciak2019 as Sequential Allocation (version 1) algorithm.
Examples
library(stratallo)
Define example population.
N <- c(3000, 4000, 5000, 2000) # Strata sizes.
S <- c(48, 79, 76, 16) # Standard deviations of a study variable in strata.
A <- N * S
n <- 190 # Total sample size.
Tschuprow-Neyman allocation (no inequality constraints).
xopt <- opt(n = n, A = A)
xopt
sum(xopt) == n
# Variance of the st. estimator that corresponds to the optimum allocation.
var_st_tsi(xopt, N, S)
One-sided upper bounds.
M <- c(100, 90, 70, 80) # Upper bounds imposed on the sample sizes in strata.
all(M <= N)
n <= sum(M)
# Solution to Problem 1.
xopt <- opt(n = n, A = A, M = M)
xopt
sum(xopt) == n
all(xopt <= M) # Does not violate upper-bounds constraints.
# Variance of the st. estimator that corresponds to the optimum allocation.
var_st_tsi(xopt, N, S)
Problem 2 (one-sided lower bounds)
Given numbers $n,\, A_h > 0,\, m_h > 0$, such that
$m_h \leq N_h,\, h = 1,\ldots,H$, and $n \geq \sum_{h=1}^H m_h$,
\begin{align}
\underset{\x \in \R_+^H}{\texteq{minimize ~\,}} & \quad f(\x) = \sum_{h=1}^H \tfrac{A_h^2}{x_h} \
\texteq{subject ~ to} & \quad \sum_{h=1}^H x_h = n \
& \quad x_h \geq m_h, \quad{h = 1,\ldots,H,}
\end{align}
where $\x = (x_h)_{h \in {1,\ldots,H}}$.
The optimization Problem 2 is solved by the LRNA that in principle is
based on the RNA and it is introduced in @wojciak2023.
Examples
m <- c(50, 120, 1, 2) # Lower bounds imposed on the sample sizes in strata.
n >= sum(m)
# Solution to Problem 2.
xopt <- opt(n = n, A = A, m = m)
xopt
sum(xopt) == n
all(xopt >= m) # Does not violate lower-bounds constraints.
# Variance of the st. estimator that corresponds to the optimum allocation.
var_st_tsi(xopt, N, S)
Problem 3 (box constraints)
Given numbers $n,\, A_h > 0,\, m_h > 0,\, M_h > 0$, such that
$m_h < M_h \leq N_h,\, h = 1,\ldots,H$, and $\sum_{h=1}^H m_h \leq n \leq \sum_{h=1}^H M_h$,
\begin{align}
\underset{\x \in \R_+^H}{\texteq{minimize ~\,}} & \quad f(\x) = \sum_{h=1}^H \tfrac{A_h^2}{x_h} \
\texteq{subject ~ to} & \quad \sum_{h=1}^H x_h = n \
& \quad m_h \leq x_h \leq M_h, \quad{h = 1,\ldots,H,}
\end{align}
where $\x = (x_h)_{h \in {1,\ldots,H}}$.
The optimization Problem 3 is solved by the RNABOX which is a new algorithm
proposed by the authors of this package and it will be published soon.
Examples
m <- c(100, 90, 500, 50) # Lower bounds imposed on sample sizes in strata.
M <- c(300, 400, 800, 90) # Upper bounds imposed on sample sizes in strata.
n <- 1284
n >= sum(m) && n <= sum(M)
# Optimum allocation under box constraints.
xopt <- opt(n = n, A = A, m = m, M = M)
xopt
sum(xopt) == n
all(xopt >= m & xopt <= M) # Does not violate any lower or upper bounds constraints.
# Variance of the st. estimator that corresponds to the optimum allocation.
var_st_tsi(xopt, N, S)
Minimization of the total cost with optcost() function
The optcost() function solves the following minimum total cost allocation
problem, formulated in the language of mathematical optimization.
Problem 4
Given numbers $A_h > 0,\, c_h > 0,\, M_h > 0$, such that
$M_h \leq N_h,\, h = 1,\ldots,H$, and $V \geq \sum_{h=1}^H \tfrac{A_h^2}{M_h} - A_0$,
\begin{align}
\underset{\x \in \R_+^H}{\texteq{minimize ~\,}} & \quad c(\x) = \sum_{h=1}^H c_h x_h \
\texteq{subject ~ to} & \quad \sum_{h=1}^H \tfrac{A_h^2}{x_h} - A_0 = V \
& \quad x_h \leq M_h, \quad{h = 1,\ldots,H,}
\end{align}
where $\x = (x_h)_{h \in {1,\ldots,H}}$.
The algorithm that solves Problem 4 is based on the LRNA and it is
described in @wojciak2023.
Examples
A <- c(3000, 4000, 5000, 2000)
A0 <- 70000
unit_costs <- c(0.5, 0.6, 0.6, 0.3) # c_h, h = 1,...4.
M <- c(100, 90, 70, 80)
V <- 1e6 # Variance constraint.
V >= sum(A^2 / M) - A0
xopt <- optcost(V = V, A = A, A0 = A0, M = M, unit_costs = unit_costs)
xopt
sum(A^2 / xopt) - A0 == V
all(xopt <= M)
Rounding
stratallo comes with 2 functions: ran_round() and round_oric() that can
be used to round non-integer allocations.
Examples
m <- c(100, 90, 500, 50)
M <- c(300, 400, 800, 90)
n <- 1284
# Optimum, non-integer allocation under box constraints.
xopt <- opt(n = n, A = A, m = m, M = M)
xopt
xopt_int <- round_oric(xopt)
xopt_int
Installation
You can install the released version of stratallo package from
CRAN with:
install.packages("stratallo")
Note on finite precision arithmetic
Consider the following example
N <- c(3000, 4000, 5000, 2000)
S <- c(48, 79, 76, 17)
a <- N * S
n <- 190
xopt <- opt(n = n, A = A) # which after simplification is (n / sum(a)) * a
xopt
and note that
sum(xopt) == n
which results from the fact that
options(digits = 22)
sum(xopt)
sum((n / sum(A)) * A) == n # mathematically, it should be TRUE!
References
Try the stratallo package in your browser
Any scripts or data that you put into this service are public.
stratallo documentation built on Nov. 27, 2023, 1:07 a.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.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
\def\R{{\mathbb R}} \def\x{\mathbf x} \def\n{\mathbf n} \newcommand{\texteq}{\mathrm}
The goal of stratallo package is to provide implementations of the efficient algorithms that solve a classical problem in survey methodology - an optimum sample allocation in stratified sampling. In this context, the classical problem of optimum sample allocation is the Tschuprow-Neyman's sense [@Tschuprow; @neyman]. It is formulated as determination of a vector of strata sample sizes that minimizes the variance of the stratified $\pi$ estimator of the population total of a given study variable, under constraint on total sample size. More specifically, the algorithms provided in this package are valid given that the variance of the stratified estimator is of the following generic form:
$$ V_{st}(\n) = \sum_{h=1}^{H} \frac{A_h^2}{n_h} - A_0, $$ where $H$ denotes total number of strata, $\n = (n_h)_{h \in {1,\ldots,H}}$ is the allocation vector with strata sample sizes, and population parameters $A_0,\, A_h > 0,\, h = 1,\ldots,H$, do not depend on the $x_h,\, h = 1,\ldots,H$. The allocation problem mentioned, can be further complemented by imposing lower or upper bounds on sample sizes is strata.
Among stratified estimators and stratified sampling designs that jointly give rise to a variance of the above form, is the so called stratified $\pi$ estimator of the population total with stratified simple random sampling without replacement design, which is one of the most basic and commonly used stratified sampling designs. This case yields $A_0 = \sum_{h = 1}^H N_h S_h^2$, $A_h = N_h S_h,\, h = 1,\ldots,H$, where $S_h$ denotes stratum standard deviation of study variable and $N_h$ is the stratum size (see e.g. @sarndal, Result 3.7.2, p.103).
A minor modification of the classical optimum sample allocation problem leads to the minimum cost allocation. This problem lies in the determination of a vector of strata sample sizes that minimizes total cost of the survey, under assumed fixed level of the stratified $\pi$ estimator's variance. As in the case of the classical optimum allocation, the problem of minimum cost allocation can be complemented by imposing upper bounds on sample sizes in strata.
Package stratallo provides two user functions:
opt()optcost()
that solve sample allocation problems briefly characterized above as well as the following helpers functions:
var_st()var_st_tsi()asummary()ran_round()round_oric().
Functions var_st() and var_st_tsi() compute a value of the variance $V_{st}$.
The var_st_tsi() is a simple wrapper of var_st() that is dedicated for the
case when $A_0 = \sum_{h = 1}^H N_h S_h^2$ and $A_h = N_h S_h,\, h = 1,\ldots,H$.
asummary() creates a data.frame object with summary of the allocation.
Functions ran_round() and round_oric() are the rounding functions that can
be used to round non-integers allocations (see section Rounding, below).
The package comes with three predefined, artificial populations with 10, 507
and 969 strata. These are stored under pop10_mM, pop507 and pop969
objects, respectively.
Minimization of the variance with opt() function
The opt() function solves the following three problems of the optimum sample
allocation, formulated in the language of mathematical optimization. User of
opt() can choose whether the solution computed will be for
Problem 1, Problem 2 or Problem 3. This is achieved with the proper
use of m and M arguments of the function. Also, if required, the inequality
constraints can be removed from the optimization problem. For more details, see
the help page for opt() function.
Problem 1 (one-sided upper bounds)
Given numbers $n > 0,\, A_h > 0,\, M_h > 0$, such that $M_h \leq N_h,\, h = 1,\ldots,H$, and $n \leq \sum_{h=1}^H M_h$, \begin{align} \underset{\x \in \R_+^H}{\texteq{minimize ~\,}} & \quad f(\x) = \sum_{h=1}^H \tfrac{A_h^2}{x_h} \ \texteq{subject ~ to} & \quad \sum_{h=1}^H x_h = n \ & \quad x_h \leq M_h, \quad{h = 1,\ldots,H,} \end{align} where $\x = (x_h)_{h \in {1,\ldots,H}}$.
There are four different algorithms available to use for Problem 1, RNA (default), SGA, SGAPLUS, COMA. All these algorithms, except SGAPLUS, are described in detail in @wesolowski2021. The SGAPLUS is defined in @wojciak2019 as Sequential Allocation (version 1) algorithm.
Examples
library(stratallo)
Define example population.
N <- c(3000, 4000, 5000, 2000) # Strata sizes. S <- c(48, 79, 76, 16) # Standard deviations of a study variable in strata. A <- N * S n <- 190 # Total sample size.
Tschuprow-Neyman allocation (no inequality constraints).
xopt <- opt(n = n, A = A) xopt sum(xopt) == n # Variance of the st. estimator that corresponds to the optimum allocation. var_st_tsi(xopt, N, S)
One-sided upper bounds.
M <- c(100, 90, 70, 80) # Upper bounds imposed on the sample sizes in strata. all(M <= N) n <= sum(M) # Solution to Problem 1. xopt <- opt(n = n, A = A, M = M) xopt sum(xopt) == n all(xopt <= M) # Does not violate upper-bounds constraints. # Variance of the st. estimator that corresponds to the optimum allocation. var_st_tsi(xopt, N, S)
Problem 2 (one-sided lower bounds)
Given numbers $n,\, A_h > 0,\, m_h > 0$, such that $m_h \leq N_h,\, h = 1,\ldots,H$, and $n \geq \sum_{h=1}^H m_h$, \begin{align} \underset{\x \in \R_+^H}{\texteq{minimize ~\,}} & \quad f(\x) = \sum_{h=1}^H \tfrac{A_h^2}{x_h} \ \texteq{subject ~ to} & \quad \sum_{h=1}^H x_h = n \ & \quad x_h \geq m_h, \quad{h = 1,\ldots,H,} \end{align} where $\x = (x_h)_{h \in {1,\ldots,H}}$.
The optimization Problem 2 is solved by the LRNA that in principle is based on the RNA and it is introduced in @wojciak2023.
Examples
m <- c(50, 120, 1, 2) # Lower bounds imposed on the sample sizes in strata. n >= sum(m) # Solution to Problem 2. xopt <- opt(n = n, A = A, m = m) xopt sum(xopt) == n all(xopt >= m) # Does not violate lower-bounds constraints. # Variance of the st. estimator that corresponds to the optimum allocation. var_st_tsi(xopt, N, S)
Problem 3 (box constraints)
Given numbers $n,\, A_h > 0,\, m_h > 0,\, M_h > 0$, such that $m_h < M_h \leq N_h,\, h = 1,\ldots,H$, and $\sum_{h=1}^H m_h \leq n \leq \sum_{h=1}^H M_h$, \begin{align} \underset{\x \in \R_+^H}{\texteq{minimize ~\,}} & \quad f(\x) = \sum_{h=1}^H \tfrac{A_h^2}{x_h} \ \texteq{subject ~ to} & \quad \sum_{h=1}^H x_h = n \ & \quad m_h \leq x_h \leq M_h, \quad{h = 1,\ldots,H,} \end{align} where $\x = (x_h)_{h \in {1,\ldots,H}}$.
The optimization Problem 3 is solved by the RNABOX which is a new algorithm proposed by the authors of this package and it will be published soon.
Examples
m <- c(100, 90, 500, 50) # Lower bounds imposed on sample sizes in strata. M <- c(300, 400, 800, 90) # Upper bounds imposed on sample sizes in strata. n <- 1284 n >= sum(m) && n <= sum(M) # Optimum allocation under box constraints. xopt <- opt(n = n, A = A, m = m, M = M) xopt sum(xopt) == n all(xopt >= m & xopt <= M) # Does not violate any lower or upper bounds constraints. # Variance of the st. estimator that corresponds to the optimum allocation. var_st_tsi(xopt, N, S)
Minimization of the total cost with optcost() function
The optcost() function solves the following minimum total cost allocation
problem, formulated in the language of mathematical optimization.
Problem 4
Given numbers $A_h > 0,\, c_h > 0,\, M_h > 0$, such that $M_h \leq N_h,\, h = 1,\ldots,H$, and $V \geq \sum_{h=1}^H \tfrac{A_h^2}{M_h} - A_0$, \begin{align} \underset{\x \in \R_+^H}{\texteq{minimize ~\,}} & \quad c(\x) = \sum_{h=1}^H c_h x_h \ \texteq{subject ~ to} & \quad \sum_{h=1}^H \tfrac{A_h^2}{x_h} - A_0 = V \ & \quad x_h \leq M_h, \quad{h = 1,\ldots,H,} \end{align} where $\x = (x_h)_{h \in {1,\ldots,H}}$.
The algorithm that solves Problem 4 is based on the LRNA and it is described in @wojciak2023.
Examples
A <- c(3000, 4000, 5000, 2000) A0 <- 70000 unit_costs <- c(0.5, 0.6, 0.6, 0.3) # c_h, h = 1,...4. M <- c(100, 90, 70, 80) V <- 1e6 # Variance constraint. V >= sum(A^2 / M) - A0 xopt <- optcost(V = V, A = A, A0 = A0, M = M, unit_costs = unit_costs) xopt sum(A^2 / xopt) - A0 == V all(xopt <= M)
Rounding
stratallo comes with 2 functions: ran_round() and round_oric() that can
be used to round non-integer allocations.
Examples
m <- c(100, 90, 500, 50) M <- c(300, 400, 800, 90) n <- 1284 # Optimum, non-integer allocation under box constraints. xopt <- opt(n = n, A = A, m = m, M = M) xopt xopt_int <- round_oric(xopt) xopt_int
Installation
You can install the released version of stratallo package from CRAN with:
install.packages("stratallo")
Note on finite precision arithmetic
Consider the following example
N <- c(3000, 4000, 5000, 2000) S <- c(48, 79, 76, 17) a <- N * S n <- 190 xopt <- opt(n = n, A = A) # which after simplification is (n / sum(a)) * a xopt
and note that
sum(xopt) == n
which results from the fact that
options(digits = 22) sum(xopt) sum((n / sum(A)) * A) == n # mathematically, it should be TRUE!
References
Try the stratallo 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.