--- title: "Optimum Sample Allocation in Stratified Sampling with `stratallo`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Optimum Sample Allocation in Stratified Sampling with `stratallo`} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} csl: "asa_modified_url.csl" references: - id: neyman type: article-journal title: "On the Two Different Aspects of the Representative Method: The Method of Stratified Sampling and the Method of Purposive Selection" author: - family: Neyman given: Jerzy container-title: Journal of the Royal Statistical Society volume: 97 number: 4 page: 558-606 issued: year: 1934 - id: Tschuprow type: article-journal title: On the Mathematical Expectation of the Moments of Frequency Distributions in the Case of Correlated Observations author: - family: Tschuprow given: Alexander Alexandrovich container-title: Metron volume: 2 number: 4 page: 461-493, 636-680 issued: year: 1923 - id: wesolowski2021 type: article-journal title: Optimality of the Recursive Neyman Allocation author: - family: Wesołowski given: Jacek - family: Wieczorkowski given: Robert - family: Wójciak given: Wojciech container-title: "Journal of Survey Statistics and Methodology" DOI: 10.1093/jssam/smab018 URL: https://arxiv.org/abs/2105.14486 issued: year: 2021 - id: wojciak2023 type: article title: Another Solution of Some Optimum Allocation Problem author: - family: Wójciak given: Wojciech container-title: Statistics in Transition new series volume: 24 number: 5 issued: year: 2023 URL: https://arxiv.org/abs/2204.04035 - id: wojciak2019 type: thesis title: Optimal Allocation in Stratified Sampling Schemes author: - family: Wójciak given: Wojciech container-title: MSc Thesis URL: 'http://home.elka.pw.edu.pl/~wwojciak/msc_optimal_allocation.pdf' publisher: Warsaw University of Technology issued: year: 2019 - id: sarndal type: book title: Model Assisted Survey Sampling author: - family: Särndal given: Carl-Erik - family: Swensson given: Bengt - family: Wretman given: Jan publisher: Springer issued: year: 1993 --- ```{r, include = FALSE} 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 ```{r load_package} library(stratallo) ``` Define example population. ```{r pop} 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). ```{r opt_Neyman} 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. ```{r opt_M} 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 ```{r opt_m} 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 ```{r opt_box} 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 ```{r optcost} 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 ```{r rounding} 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](https://CRAN.R-project.org) with: ``` r install.packages("stratallo") ``` ## Note on finite precision arithmetic Consider the following example ```{r finit_prec1} 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 ```{r finit_prec2} sum(xopt) == n ``` which results from the fact that ```{r finit_prec3} options(digits = 22) sum(xopt) sum((n / sum(A)) * A) == n # mathematically, it should be TRUE! ``` ## References