R Package ‘stratifyR’

Type Package

Title Optimal Stratification of Univariate Populations

Version 2.0-2

Date 2026-10-14

Author Karuna G. Reddy, M.G.M. Khan

Maintainer Karuna G. Reddy <karuna.reddy@auckland.ac.nz>

Description This implements the stratification of univariate populations under stratified sampling designs using the method of Khan et al. (2002) https://doi.org/10.1177/0008068320020518, Khan et al. (2008) https://www150.statcan.gc.ca/n1/pub/12-001-x/2008002/article/10761-eng.pdf and Khan et al. (2015) https://doi.org/10.1080/02664763.2015.1018674. It determines the Optimum Strata Boundaries (OSB) and Optimum Sample Sizes (OSS) for the study variable, \(y\), using the best-fit frequency distribution of a survey variable (if data is available) or a hypothetical distribution (if data is not available). The method formulates the problem of determining the OSB as mathematical programming problem which is solved by using a dynamic programming technique. Version 2.0 adds three further solvers (COBYLA, GLOBAL and EXACT), described in Section 7. If a dataset of the population is available to the surveyor, the method estimates its best-fit distribution and determines the OSB and OSS under Neyman allocation directly. When the dataset is not available, stratification is made based on the assumption that the values of the study variable, \(y\), are available as hypothetical realizations of proxy values of \(y\) from recent surveys. Thus, it requires certain distributional assumptions about the study variable. At present, it handles stratification for the populations where the study variable follows a continuous distribution, namely, Pareto, Triangular, Right-triangular, Weibull, Gamma, Exponential, Uniform, Normal, Log-normal and Cauchy distributions. In this vignette, the two major functionalities in the package are applied to a number of real and simulated populations and to some hypothetical populations.

License GPL (>= 3)

LazyData yes

NeedsCompilation no

Depends R (>= 4.1.0)

Imports stats, utils, graphics, grDevices, MASS, fitdistrplus, nloptr, actuar, mc2d

Repository CRAN

Using stratifyR through the interactive app

The quickest and easiest way to use stratifyR is through its built-in interactive application. The app gives you the complete workflow — load your own data or pick a built-in dataset, compute the Optimum Strata Boundaries (OSB) and Optimum Sample Sizes (OSS) with the DP, COBYLA, GLOBAL or EXACT solvers, compare design efficiencies, explore the boundaries visually, and download the results — entirely through a point-and-click interface, with no R programming required. For most users, and especially for anyone who wants results quickly, this is the recommended way to work with the package.

To launch it, install the app’s helper packages once (if you do not already have them) and then call stratifyRApp():

# One-time install of the app's helper packages (skip any you already have)
install.packages(c("shiny", "bslib", "DT"))

# Optional extras: interactive 3-D / slider plots and Excel (.xlsx) import
install.packages(c("plotly", "readxl"))

# Launch the app
library(stratifyR)
stratifyRApp()

The application opens in your web browser (or as a standalone desktop-style window). The remainder of this vignette documents the underlying R functions, for readers who prefer to script their analyses or embed stratifyR in a larger workflow.

1 Introduction

The main aim of stratification is to produce estimators with a small variance when a population characteristic \((y)\) is under study. A simple method that can be used to create strata for this population, if \(y\) itself is the stratification variable. The ideal situation is that the distribution of such a study variable is known and the OSB can be determined by placing boundaries on the range of this distribution at suitable cut-points. This problem of determining the OSB, when both the estimation and stratification variables are the same, was first discussed by Dalenius (1950). He provided equations for the determination of stratum boundaries that minimize the variance of population estimates under optimal allocation. Dalenius (1957) futher proposed a solution to the problem by taking equal intervals of the cumulative square root of frequency scale of the stratification variable.

One of the many kinds of stratification methods that has been proposed in the literature is due to Bühler and Deutler (1975). They formulated the problem of determining OSB as an optimization problem and developed a computational technique to solve the problem by using Dynamic Programming (DP). A good review of this method can be found in Khan et al. (2008). The procedure is also applied by Khan et al. (2002), Khan et al. (2003), Khan et al. (2005), Khan et al. (2008), Khan et al. (2009), Khan, Rao, et al. (2015), Khan, Reddy, et al. (2015) for determining OSB for many different distributions. With the known frequency function of the study variable, they considered the problem of finding OSB as an equivalent problem of determining Optimum Strata Width (OSW), which is formulated as a Mathematical Programming Problem (MPP) and solved by using DP technique. They applied the technique to several univariate populations where the study variables followed different probability distributions. In this package, a univariate stratification technique is developed, which is based on the probability distribution assumed by the stratification variable as discussed by the authors above.

As implemented in this package, there are many advantages of using the method. The real advantage of the stratifyR package is that when the dataset of the study variable is not available, which may occur in practice, the package is still able to construct OSB and OSS based on the distributional assumption on the data. Moreover, the population survey can be a costly and time-consuming affair, hence, this approach also has the advantage of determining OSB for the study variable that is not available prior to conducting the survey. This, of course, requires certain distributional or parametric assumptions on the study variables, which can easily be obtained from recent or past surveys.

Other advantages of the method are that it leads to substantial gains in the precision of the estimates over other available methods. Results reveal that the variances get smaller with increasing number of strata \((L)\), and they get much smaller at a much faster rate than other available methods. Once the OSB have been determined, the Optimum Sample Sizes (OSS) can be easily calculated for each stratum using Neyman allocation (Neyman (1934)). The algorithm may be a little slow, however, it does provide very efficient results which is probably the most important objective of our survey estimation efforts.

2 The Package in a Glance

The stratifyR package implements the DP technique (from various literature by Khan et al) as a stratification procedure for univariate populations, when the stratification variable follows a continuous probability distribution, namely: uniform, triangular, right-triangular, pareto, exponential, normal, log-normal, cauchy, weibull and gamma. The package is able to determine the OSB and OSS from real data and as well as hypothetical situations when the dataset is not available. The latter requires some assumptions about the distribution and its initial value, range, parameter values, fixed sample size, etc.

There are two major functions which basically solve the two types of stratification problems: strata.data() which carries out univariate stratification for those univariate populations where dataset is available and strata.distr() which performs stratification when dataset is not available prior to conducting the survey.

In the former case, data on the study variable, number of strata \((h)\), fixed sample size \((n)\) and population size \((N)\) are used as the input arguments to the strata.data() function in the package. In the latter case, strata.distr() function is called which requires the distribution to be assumed, its parameters, the inital value and the estimated range of the distribution, fixed sample \((n)\) and population sizes \((N)\). When executed, both the functions output the OSB and OSS, amongst other quantities such as stratum weight \((W_{h})\), stratum variance \((S_{h}^{2})\), stratum objective function values \((W_{h}S_{h})\), stratum sample sizes \((n_{h})\), stratum population sizes \((N_{h})\) and stratum sampling fraction \((f_{h})\).

The following sections show the general formulation of the problem of stratification, the DP solution procedure, the concept of Neyman allocation as the method of determining sample size and the overview of package functionalities. The package is then applied to numerous survey populations with real and simulated data to illustrate its application.

3 General Formulation of the Univariate Stratification Problem

Let the target population of the variable under study be stratified into \(L\) strata where the estimation of the mean of the study variable \((y)\) is of interest. If a simple random sample of size \(n_{h}\) is to be drawn from \(h^{th}\) stratum with sample mean \(\bar{y}_{h}\), then the stratified sample mean, \(\bar{y}_{st}\), is given by \[\begin{eqnarray} \bar{y}_{st}=\sum_{h=1}^{L}W_{h}\bar{y}_{h}, \tag{1} \end{eqnarray}\] where \(W_{h}\) (stratum weight) is the proportion of the population contained in the \(h^{th}\) stratum.

When the finite population correction factors are ignored, under the Neyman (1934) allocation, the variance of \(\bar{y}_{st}\) is given by \[\begin{eqnarray} Var(\bar{y}_{st})=\dfrac{\left(\sum_{h=1}^{L}W_{h}S_{h}\right)^{2}}{n}, \tag{2} \end{eqnarray}\] where \(S_{h}^{2}\) is the stratum variance for the study variable in the \(h^{th}\) (where \(h=1, 2, ..., L\)) stratum and \(n\) is the preassigned total sample size.

Let \(f(y);\;a\leq y\leq b\) be the frequency function of the study variable, \(y\), on which OSB are to be constructed. If the population mean of this study variable is estimated under Neyman allocation, then the problem of determining OSB is to cut up the range, \(d=b-a\), at \((L-1)\) intermediate points \(a=y_{0} \leq y_{1} \leq y_{2} \leq, ..., \leq y_{L-1} \leq y_{L}=b\,\) such that (2) is minimum. The lower and upper bounds of the study variable are denoted by \(a\) and \(b\) respectively and the cut-points \(y_{0}, y_{1}, y_{2}, ..., y_{L-1}\) are the OSB.

For a fixed sample size \(n\), minimizing the expression of the right hand side of equation (2) is equivalent to minimizing \[\begin{eqnarray} \sum_{h=1}^{L} W_{h}S_{h} \tag{3} \end{eqnarray}\]

If \(f(y)\) is known and integrable, the stratum weight \((W_{h})\), stratum variance \((S_{h}^{2})\) and stratum mean \((\mu_{h})\) can be obtained as a function of the boundary points \(y_{h}\) and \(y_{h-1}\) by using the following expressions: \[\begin{eqnarray}\label{stratumweight} W_{h} = \int_{y_{h-1}}^{y_{h}} f(y)dy \tag{4} \end{eqnarray}\]

\[\begin{eqnarray} S_{h}^{2}= \dfrac{1}{W_{h}}\int_{y_{h-1}}^{y_{h}}y^{2}f(y)dy - \mu_{h}^{2} \tag{5} \end{eqnarray}\]

\[\begin{eqnarray} \textrm{where}\;\;\;\mu_{h}=\dfrac{1}{W_{h}}\int_{y_{h-1}}^{y_{h}}yf(y)dy \tag{6} \end{eqnarray}\]

where \((y_{h-1}, y_{h})\) are the boundaries of \(h^{th}\) stratum.

Thus, the objective function in (3) could be expressed as a function of boundary points \(y_{h}\) and \(y_{h-1}\) only. We further define \[\begin{eqnarray} l_{h}=y_{h}-y_{h-1}; \; h=1,2, ..., L \tag{7} \end{eqnarray}\]

where \(l_{h}\geq 0\) denotes the range or width of the \(h^{th}\) stratum. Then, the range of the distribution, \(d = b - a\), is expressed as a function of stratum width as: \[\begin{eqnarray}\label{215} \sum_{h=1}^{L}l_{h}=\sum_{h=1}^{L}(y_{h}-y_{h-1}) = b-a = y_{L} - y_{0} = d \tag{8} \end{eqnarray}\] The \(h^{th}\) stratification point \(y_{h};\;h=1,2,..., L\,\) is then expressed as \(y_{h}=y_{h-1}+l_{h}\) and from (8), the problem can be treated as an equivalent problem of determining Optimum Strata Widths (OSW), \(l_{1}, l_{2}, ..., l_{L}\). Due to the special nature of functions, the problem may be treated as a function of \(l_{h}\) alone and can be expressed as: \[\begin{eqnarray}\label{genMPP} &\textrm{Minimize}& \;\;\; \sum_{h=1}^{L}\phi_{h}(l_{h}),\nonumber \\ &\textrm{subject to}&\;\;\;\; \sum_{h=1}^{L}l_{h} = d,\nonumber \\ &\textrm{and}& \;\;\;\; l_{h} \geq 0; \;\; h = 1,2,...,L. \tag{9} \end{eqnarray}\]

4 Dynamic Programming Technique as a Solution Procedure

Many multistage decision problems can be formulated as a Mathematical Programming Problem (MPP). The Dynamic Programming technique, developed by Bellman, R E (1957), is a computational method which is well suited for solving MPPs that may be treated as a multistage decision problem. The DP technique determines the optimum solution of a multistage problem by decomposing it into stages, each stage comprising of a single stage. The advantage of the decomposition is that the optimization process at each stage involves one variable only, which simplifies the computational task by dealing with all variables simultaneously.

The solution to an MPP is achieved in a sequential manner starting from one stage problem, moving onto a two stage problem, to a three stage problem and so on until finally all stages are included. The solution for \(n\) stages is obtained by adding the \(n^{th}\) stage to the solution of \(n - 1\) stages.

The basic concept of DP technique is contained in the principle of optimality proclaimed by Bellman, R E (1957), which implies that given the initial state of a system, an optimal policy for the subsequent stages does not depend upon the policy adopted at the preceding stages. It determines the optimum solution of a multi-variable problem by decomposing it into stages, each stage compromising a single variable sub-problem. A dynamic programming model is basically a recursive equation which links the different stages of the problem in a manner which guarantees that each stage’s optimal feasible solution is also optimum and feasible for the entire problem.

Consider the following sub-problem of (9) for first \(k(<L)\) strata:

\[\begin{eqnarray} &\text{Minimize}& \;\;\;\; \sum_{h=1}^{k}\phi_{h}(l_{h}),\nonumber \\ &\text{subject to}&\;\;\;\;\; \sum_{h=1}^{k}l_{h} = d_{k},\nonumber \\ &\text{and}& \;\;\;\;\; l_{h} \geq 0; \;\; h = 1,2,...,k.\tag{10} \end{eqnarray}\]

where \(d_{k}< d\) is the total width available for division into strata or the state value at stage \(k\). Note that \(d_{k} = d\) for \(k = L\).

The transformation functions are given by:

\[\begin{eqnarray*} d_{k}\;\;&=&\;\;l_{1}+l_{2}+...+l_{k},\\ d_{k-1}\;\;&=&\;\;l_{1}+l_{2}+...+l_{k-1}\;=\;d_{k}-l_{k},\\ d_{k-2}\;\;&=&\;\;l_{1}+l_{2}+...+l_{k-2}\;=\;d_{k-1}-l_{k-1},\\ &\vdots&\;\;\;\;\;\;\;\;\;\;\;\;\vdots\\ d_{2}\;\;&=&\;\;l_{1}+l_{2}\;=\;d_{3}-l_{3},\\ d_{1}\;\;&=&\;\;l_{1}\;=\;d_{2}-l_{2}. \end{eqnarray*}\]

Let \(\Phi_{k}(d_{k})\) denote the minimum value of the objective function of MPP (10), that is,

\[\begin{equation} \Phi_{k}(d_{k}) = \text{min}\left[ \sum_{h=1}^{k}\phi_{h}(l_{h})\bigg|\sum_{h=1}^{k}l_{h}=d_{k}, \,\textrm{and}\,l_{h}\;\geq 0;\,h=1,2,...,k \;\textrm{and}\; 1\leq k \leq L \right].\nonumber \end{equation}\]

With the above definition of \(\Phi_{k}(d_{k})\), (10) is equivalent to finding \(\Phi_{L}(d)\) recursively by finding \(\Phi_{k}(d_{k})\) for \(k = 1, 2, ..., L\) and \(0 \leq d_{k} \leq d.\)

We can write:

\[\begin{equation} \Phi_{k}(d_{k}) = \text{min}\left[\phi_{k}(l_{k})+ \sum_{h=1}^{k-1}\phi_{h}(l_{h})\bigg|\sum_{h=1}^{k-1}l_{h}=d_{k}-l_{k}, \;\textrm{and}\;l_{h}\;\geq\;0;\;h=1,2,...,k\right].\nonumber \end{equation}\]

For a fixed value of \(l_{k}\); \(0 \leq l_{k} \leq d_{k}\),

\[\begin{equation} \Phi_{k}(d_{k}) = \phi_{k}(l_{k})+ \text{min}\left[ \sum_{h=1}^{k-1}\phi_{h}(l_{h})\bigg|\sum_{h=1}^{k-1}l_{h}=d_{k}-l_{k}, \;\textrm{and}\;l_{h}\;\geq 0;\;h=1,2,...k-1\;\textrm{and}\;\; 1\leq k \leq L\right].\nonumber \end{equation}\]

Using the Bellman’s principle of optimality, a forward recursive equation can be written as:

\[\begin{equation} \Phi_{k}(d_{k}) = {\text{min} \atop 0 \leq l_{k} \leq d_{k}}\left[\phi_{k}(l_{k}) + \Phi_{k-1}(d_{k}-l_{k})\right],\;\;k\;\geq\;2.\tag{11} \end{equation}\]

For the first stage, that is, for \(k=1\):

\[\begin{equation} \Phi_{1}(d_{1}) = \phi_{1}(d_{1})\; \Longrightarrow\;\;l_{1}^{*}=d_{1},\tag{12} \end{equation}\]

where \(l_{1}^{*}=d_{1}\) is the optimum width of the first stratum. The relations (11) and (12) are solved recursively for each \(k=1, 2, ..., L\,\) and \(0 \leq d_{k} \leq d\), and \(\Phi_{L}(d)\) is obtained. From \(\Phi_{L}(d)\) the optimum width of \(L^{th}\) stratum, \(l_{L}^{*}\), is obtained. From \(\Phi_{L-1}(d - l_{L}^{*})\) the optimum width of \((L - 1)^{th}\) stratum, \(l_{L-1}^{*}\), is obtained and so on until \(l_{1}^{*}\) is obtained.

The above algorithm of the Dynamic Programming solution procedure to solve the MPP given in (9) is summarized with the following steps:

1. Start at \(k=1\). Set \(\Phi_{0}(d_{0})=0\).
2. Calculate \(\Phi_{1}(d_{1})\), the minimum value of RHS of \((12)\) for \(l_{1}=d_{1},\; 0\leq l_{1}\leq d_{1}\), and \(0\leq d_{1}\leq d\).
3. Record \(\Phi_{1}(d_{1})\) and \(l_{1}\).
4. For \(k=2\), express the state variable as \(d_{k-1}=d_{k}-l_{k}\).
5. Set \(\Phi_{k}(d_{k})=0\) if \(l_{k}>d_{k}\), where \(0\leq d_{k}\leq d\).
6. Calculate \(\Phi_{k}(d_{k})\), the minimum value of RHS of \((11)\) for \(l_{k};\;0\leq l_{k}\leq d_{k}\).
7. Record \(\Phi_{k}(d_{k})\) and \(l_{k}\).
8. For \(k\geq3, ..., L\), go to Step 4.
9. At \(k=L, \; \Phi_{L}(d)\) is obtained and hence the optimum value \(l_{L}^{*}\) of \(l_{L}\) is obtained.
10. At \(k=L-1\), using the backward calculation for \(d_{L-1}=d-l_{L}^{*}\), read the value of \(\Phi_{L-1}(d_{L-1})\) and hence the optimum value \(l_{L-1}^{*}\) of \(l_{L-1}\).
11. Repeat Step \((10)\) until the optimum value \(l_{1}^{*}\) of \(l_{1}\) is obtained from \(\Phi_{1}(d_{1})\).

5 Optimum Sample Sizes Using Neyman Allocation

When OSB \((y_{h}, y_{h-1})\) have been determined, the Optimum Sample Sizes (OSS) \(n_{h}; h=1,2,...,L\,\) that minimizes the variance of the estimate can easily be computed. The sample size \(n_{h}\) are obtained for a fixed total sample of size \(n\) under the Neyman allocation for \(h=1,2,...,L\,\) and given as follows:

\[\begin{equation}\label{Ney_alloc} n_{h} = n\,\frac{W_{h}S_{h}}{\sum_{h=1}^{L}W_{h}S_{h}} \tag{13} \end{equation}\]

where \(W_{h}\) and \(S_{h}\) are derived in terms of the optimum boundary points \((y_{h}, y_{h-1})\).

In Neyman allocation, the total sample size is proportional to the stratum size multiplied by the standard deviation of the stratum. If the variances are specified correctly, Neyman allocation will give an estimator with smaller variance compared to proportional allocation (Lohr (2009)).

In equation (13), it is also worth noting that the OSB \((y_{h}, y_{h-1})\) are so obtained that \(n_{h}\) must satisfy the restrictions:

\[\begin{equation}\label{res} 1\leq n_{h}\leq N_{h}, \tag{14} \end{equation}\]

where \(N_{h}=NW_{h}\). Thus, the restriction (14) indicates that the \(h^ {th}\) stratum must form with at least one unit and also avoid the problem of over-sampling.

6 Overview of Package Functionalities

For the numerical illustrations and application of the package, some of the real datasets such as ‘sugarcane’ of Khan, Reddy, et al. (2015), ‘anaemia’ of Reddy et al. (2014); ‘hies’ and ‘math’ data are provided in the stratifyR. The ‘quakes’ and ‘Boston’ data provided in the datasets package in R are also used for illustration purposes. The stratifyR package is also tested on some published and commonly-used datasets such as ‘UScities’ and ‘UScolleges’ data from Cochran (1961), ‘Debtors’ data of Gunning and Horgan (2004), ‘REV84’ variable for ‘Swedish municipalities’ data from Särdnal et al. (1992) and ‘MRTS’ variable from ‘Statistics Canada Monthly Retail Trade Survey’ together with ‘SHS’ data collected in ‘Statistics Canada Survey of Household Spending’. For those distributions where real data is not found in literature, data is simulated to demonstrate the application of the package in this documentaion.

For a user, there are two different routes available in the package and these are basically dependent on the type of stratification problem to be solved. It could either be a data-based (i.e., when the dataset of the stratification variable is available) or a distribution-based (i.e., when dataset is not available) stratification problem. Whether stratification is based on data or not, the idea is that the problem is formulated as an MPP using the estimated (with available data) or assumed (with unavailable data) distribution of the data set. There are numerous functions created in the package, however, there are only a few major functions that are used by the two different types of problems being studied in univariate stratification.

If it is a data-based problem, the function used is strata.data() and the user has to provide as input arguments: the data, the number of strata (\(L\)) and the fixed sample size (\(n\)). For the distribution-based problem, the function used is strata.distr() and the user has to provide the name of the assumed distribution, number of strata (\(L\)), possible range of data (distance), fixed sample size (\(n\)) and the population size (\(N\)). The following two subsections delve a little deeper into the workings surrounding the two functions: strata.data() and strata.distr().

6.1 The Function strata.data()

This function computes the OSB, OSS, and other important quantities from univariate survey populations by employing the methodology proposed by Khan et al. (2002) Khan et al. (2003), Khan et al. (2005), Khan et al. (2008), Khan et al. (2009), Nand and Khan (2009), Khan, Rao, et al. (2015), Khan, Reddy, et al. (2015), Reddy et al. (2018) and Reddy and Khan (2019). Their method uses the estimated distribution of the data and formulates the problem of determining OSB as a Mathematical Programming Problem which is an optimization problem that is solved by the DP technque. The OSB obtained are optimal solutions that are used to calculate the OSS under Neyman allocation. The usage and arguments are as follows:

    strata.data(data, h, n, cost=FALSE, ch=NULL,
                method=c("dp","cobyla","global","exact"),
                n_starts=20L, max_iter=2000L, tol=1e-9, verbose=FALSE)

    data     - A vector: data containing every unit of the survey population
    h        - A numeric: number of strata to be sampled. The default is 2
    n        - A numeric: fixed total sample size
    cost     - A logical: stratum cost. Default cost=FALSE.
    ch       - A numeric: denotes a vector of stratum costs. Default ch=NULL.
    method   - A character: optimisation method. One of "dp" (default), "cobyla",
               "global" or "exact" (see Section 7).
    n_starts - An integer: number of random restarts for COBYLA. Default 20.
    max_iter - An integer: maximum function evaluations per COBYLA run. Default 2000.
    tol      - A numeric: convergence tolerance for COBYLA. Default 1e-9.
    verbose  - A logical: if TRUE, prints per-start diagnostics. Default FALSE.
    

The steps below describe the default DP solver. The other solvers replace steps 4 and 5 (see Section 7).

To show the application of the strata.data() function, an example of the command used and its output from the package is given below. The problem uses the ‘mag’ variable from the ‘quakes’ data (with a population of \(N=1000\)) available from the datasets package in R. To construct a 2-strata solution with a fixed sample size of \(n=300\), we use the following codes:

data(quakes)
head(quakes)
#>      lat   long depth mag stations
#> 1 -20.42 181.62   562 4.8       41
#> 2 -20.62 181.03   650 4.2       15
#> 3 -26.00 184.10    42 5.4       43
#> 4 -17.97 181.66   626 4.1       19
#> 5 -20.42 181.96   649 4.0       11
#> 6 -19.68 184.31   195 4.0       12
mag <- quakes$mag
length(mag)
#> [1] 1000
hist(mag) #to see the distribution

res <- strata.data(mag, h = 2, n=300) # a 2-strata solution
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 300 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 1000</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [4.0000, 6.4000]  (d = 2.4000)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> lnorm</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter | Estimate|
#> |:---------|--------:|
#> |meanlog   | 1.526810|
#> |sdlog     | 0.085036|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	4.68122	0.58	0.03	0.109	140	585	0.24
2	6.4	0.42	0.09	0.124	160	415	0.39
Total		1.00	0.12	0.233	300	1000	0.30

The boundaries are computed at full precision and printed to six significant figures. Other quantities are rounded for display, so the individual stratum values in the tables may not add up exactly to the totals.

6.2 The Function strata.distr()

This function is also used to compute the OSB, OSS, and other important quantities from univariate survey populations by employing the methodology proposed by Khan et al given earlier. Its algorithm is quite similar to that of the strata.data(), however, its functionality is applied to the case where the dataset of the population is not available and the distributonal assumptions of the study variable are strictly required. Another caveat for such distribution-based stratification is that the distr.alloc() function uses the probability density functions of the assumed distributions and integration rules presented by equations (4)-(6) to calculate the stratum sample sizes. It must be noted that this function works on ideal distributions that assumes the parameters chosen by the user. The usage and arguments are as follows:

    strata.distr(h, initval = NULL, dist = NULL,
                 distr = c("pareto", "triangle", "rtriangle", "weibull", "gamma",
                 "exp", "unif","norm", "lnorm", "cauchy"), params = c(shape=0,
                 scale=0, rate=0, gamma=0, location=0, mean=0, sd=0, meanlog=0,
                 sdlog=0, min=0, max=0, mode=0), n, N, cost=FALSE, ch=NULL,
                 method=c("dp","cobyla","global"),
                 n_starts=20L, max_iter=2000L, tol=1e-9, verbose=FALSE)

    h        - A numeric: number of strata to be sampled
    initval  - A numeric: initial value of the assumed distribution
    dist     - A numeric: distance or range of the assumed distribution
    distr    - A character: the assumed distribution of the hypothetical population
    params   - A list: parameters of the assumed distribution
    n        - A numeric: fixed total sample size
    N        - A numeric: fixed population size
    cost     - A logical: stratum cost. Default cost=FALSE.
    ch       - A numeric: denotes a vector of stratum costs. Default ch=NULL.
    method   - A character: optimisation method. One of "dp" (default), "cobyla" or "global".
    n_starts - An integer: number of random restarts for COBYLA. Default 20.
    max_iter - An integer: maximum function evaluations per COBYLA run. Default 2000.
    tol      - A numeric: convergence tolerance for COBYLA. Default 1e-9.
    verbose  - A logical: if TRUE, prints per-start diagnostics. Default FALSE.
    

The algorithm for strata.distr() is quite similar to the strata.data() for the contruction of OSB. It is only at the sample allocation (OSS) stage that the two are different. strata.distr() is the main function and once the OSB have been computed, this function uses the distr.alloc() function which uses the numerical integration rules for different distibutions (i.e., the equations (4)-(6)) to compute the OSS.

To show the application of the strata.distr() function, let us construct a 2-strata solution assuming that the dataset of the stratification variable is not available but its distribution and estimated parameters are. Let us consider the ‘depth’ variable from the ‘quakes’ dataset from the datasets package, which has a Triangular distribution with parameters \(min=39.99998, max=680, mode=39.99999\). It starts at an initial value of \(40\) and has a distance (range) of \(640\) with a fixed sample size of \(n=300\) from a population of \(N=1000\) seismic events. Thus, we use the following commands:

data(quakes) 
depth <- quakes$depth
hist(depth) #see distribution

min(depth); max(depth); d=max(depth)-min(depth);d #min, max and range of data 
#> [1] 40
#> [1] 680
#> [1] 640
# the 2-strata solution is
res <- strata.distr(h=2, initval=40, dist=640, distr = "triangle",
             params = c(min=39.99998, max=680, mode=39.99999), n=300, N=1000)
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 300 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 1000</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [40.0000, 680.0000]  (d = 640.0000)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> triangle</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter |  Estimate|
#> |:---------|---------:|
#> |min       |  39.99998|
#> |max       | 680.00000|
#> |mode      |  39.99999|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	266.723	0.58	4217.46	37.862	145	583	0.25
2	680	0.42	9488.76	40.619	155	417	0.37
Total		1.00	13706.22	78.481	300	1000	0.30

7 Alternative Optimisation Methods: COBYLA, GLOBAL and EXACT

7.1 Motivation

The Dynamic Programming solver described in Section 4 is the foundation of stratifyR. It searches a grid of 1,000 candidate positions and is globally optimal on that grid for the distribution it is given. Two things limit it in practice. First, the search grows quickly with the number of strata \(h\), so it can take a minute or more at \(h = 5\) or \(6\). Second, in strata.data() the DP optimises the fitted distribution, not the data themselves. When the fitted family describes the data well the two problems have almost the same solution, but when the fit is poor (for example in a very heavy tail) the DP boundaries can be noticeably worse for the actual population.

stratifyR 2.0 therefore offers four solvers through the method argument:

• method = "dp" (default): the original DP grid search on the fitted distribution;
• method = "cobyla": a fast multi-start derivative-free local optimiser;
• method = "global": a two-phase global search (DIRECT-L followed by BOBYQA refinement);
• method = "exact" (strata.data() only): DP run directly on the sorted data, which returns the best possible stratification of the population.

DP remains the default, so code written for version 1.x gives the same kind of result as before.

7.2 The COBYLA Method (method = "cobyla")

COBYLA (Constrained Optimisation BY Linear Approximations, Powell (1994)) is a derivative-free method that minimises a function subject to inequality constraints by building linear models of the objective and constraints. stratifyR treats the \(h-1\) interior boundaries as continuous variables. Adjacent boundaries must be at least one thousandth of the range apart, as in the problem formulation, and every stratum must hold at least two units.

The solver is run from several starting points (argument n_starts, default 20): equally spaced boundaries, equal-probability quantiles, and random restarts. From each start both COBYLA and BOBYQA are run and the better result is kept. The best result over all starts is returned. The random restarts do not change the user’s random number stream.

• For strata.distr(), COBYLA minimises the same analytical objective as the DP, but over continuous boundaries rather than a grid.
• For strata.data(), COBYLA minimises the empirical objective \(\sum_h W_h S_h\) computed from the data. This is exactly the quantity reported as WhShTot.

Advantages: fast, usually about a second even for six strata, so it is well suited to trying many values of \(h\) or many variables. On data it usually lands on the same partition as method = "exact", and you can check this in any particular case by running both.

Limitations: it is a local method, so a global optimum is not guaranteed. The summary() output reports a Cauchy-Schwarz optimality gap, KKT residuals and a convergence flag to help judge the solution.

7.3 The GLOBAL Method (method = "global")

Phase 1, DIRECT-L. DIRECT-L (Dividing RECTangles, locally biased; Gablonsky and Kelley (2001)) splits the search box into smaller boxes and samples their centres, refining boxes that are either promising or still large and unexplored. It is deterministic and uses a budget of \(\max(5000, 5000(h-1))\) evaluations.

Phase 2, BOBYQA. The best DIRECT-L point is polished with the BOBYQA local optimiser. If neither phase finds a feasible design (one with at least two units in every stratum), the solver falls back to multi-start COBYLA.

Advantages: deterministic and reproducible; on smooth problems and moderate \(h\) it usually reaches the same solution as COBYLA.

Limitations: when the data cover a very wide range, the fixed evaluation budget may not be enough to find narrow strata near the bottom of a skewed distribution, especially for larger \(h\). If GLOBAL and COBYLA disagree, trust the lower WhShTot, or check with method = "exact".

7.4 The EXACT Method (method = "exact")

When you have every value of the stratification variable, sorting them shows that a set of strata is just a way of cutting the sorted list into \(h\) consecutive pieces. method = "exact" tries every such cut, using the dynamic programming idea of Section 4 so that this stays fast, and returns the cut with the smallest WhShTot. In other words, no other choice of boundaries for these data can give a smaller Neyman cost. Each stratum keeps at least two units, and units with the same value are always kept together.

In practice:

• It needs the data, so it is available in strata.data() only. With strata.distr() there is nothing to sort, and the DP is already optimal for the assumed distribution.
• It is quick for typical survey frames: a few thousand units take a second or two, and around 14,000 units with \(h = 6\) take well under half a minute.
• For very large frames (more than 20,000 distinct values) it searches 20,000 cut points spread evenly through the sorted data, which is close to exact. You can raise the limit, at the cost of time and memory, with options(stratifyR.exact_max_cand = 50000). The element exact_full of the result tells you whether every cut point was searched.
• The boundary it reports is the largest value in each stratum, so a unit equal to a boundary belongs to the lower stratum.

7.5 Which method should I use?

A few questions settle it:

• Do you have the data, or only an assumed distribution? With only a distribution, use strata.distr() with the default "dp", or "cobyla" if you want results instantly.
• Is the design for this very frame? Then "exact" gives the best boundaries for it, usually in seconds.
• Will the boundaries be reused, for example in later rounds of the survey or on an updated frame? The default "dp" places the boundaries for the fitted distribution rather than for every detail of one dataset, so they tend to change less from one frame to the next. It is worth comparing its WhShTot with "exact": a small difference means the fitted family describes your data well.
• Are you exploring, for example many values of \(h\) or several variables? Use "cobyla" for speed, then confirm the final choice with "exact" or "dp".

7.6 Application: Comparing Methods for strata.distr

The following example uses a Weibull distribution with the parameters fitted to the Iron variable of the ‘anaemia’ data, and compares the three solvers available in strata.distr() for a 3-strata solution.

# Weibull(shape = 2.1446, scale = 13.7907), as fitted to the Iron variable
data(anaemia)
Iron <- anaemia$Iron

# DP (reference solution)
res_dp_d <- strata.distr(h=3, initval=1.5, dist=33.2,
              distr="weibull", params=c(shape=2.1446, scale=13.7907),
              n=350, N=724, method="dp")
#> The program is running, it'll take some time!
summary(res_dp_d)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 3 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 350 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 724</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [1.5000, 34.7000]  (d = 33.2000)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> weibull</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter | Estimate|
#> |:---------|--------:|
#> |shape     |   2.1446|
#> |scale     |  13.7907|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	9.66013	0.37	4.53	0.775	116	270	0.43
2	16.3069	0.39	3.55	0.733	110	281	0.39
3	34.7	0.24	11.94	0.822	124	173	0.72
Total		1.00	20.02	2.330	350	724	0.48

# COBYLA
res_cob_d <- strata.distr(h=3, initval=1.5, dist=33.2,
               distr="weibull", params=c(shape=2.1446, scale=13.7907),
               n=350, N=724, method="cobyla")

# GLOBAL (DIRECT-L + BOBYQA)
res_glob_d <- strata.distr(h=3, initval=1.5, dist=33.2,
                distr="weibull", params=c(shape=2.1446, scale=13.7907),
                n=350, N=724, method="global")

# Direct comparison of the three methods
cat("DP     OSB:", round(res_dp_d$OSB, 3),   "  WhShTot:", res_dp_d$WhShTot, "\n")
#> DP     OSB: 9.66 16.307 34.7   WhShTot: 2.33
cat("COBYLA OSB:", round(res_cob_d$OSB, 3),  "  WhShTot:", res_cob_d$WhShTot, "\n")
#> COBYLA OSB: 9.66 16.307 34.7   WhShTot: 2.33
cat("GLOBAL OSB:", round(res_glob_d$OSB, 3), "  WhShTot:", res_glob_d$WhShTot, "\n")
#> GLOBAL OSB: 9.66 16.307 34.7   WhShTot: 2.33

All three solvers minimise the same smooth objective here, and they agree to within the DP grid spacing. COBYLA and GLOBAL take a fraction of a second, while the DP takes much longer.

7.7 Application: Comparing Methods for strata.data

The next example uses the Iron data directly and compares all four solvers for a 2-strata solution.

data(anaemia)
Iron <- anaemia$Iron

res_dp_s    <- strata.data(Iron, h=2, n=350, method="dp")
#> The program is running, it'll take some time!
res_cob_s   <- strata.data(Iron, h=2, n=350, method="cobyla")
res_glob_s  <- strata.data(Iron, h=2, n=350, method="global")
res_exact_s <- strata.data(Iron, h=2, n=350, method="exact")
summary(res_exact_s)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 350 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 724</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [1.5000, 34.7000]  (d = 33.2000)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> weibull</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter |  Estimate|
#> |:---------|---------:|
#> |shape     |  2.145048|
#> |scale     | 13.791617|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	13.2	0.62	9.63	1.938	194	452	0.43
2	34.7	0.38	17.35	1.565	156	272	0.57
Total		1.00	26.99	3.503	350	724	0.48

# Direct comparison
cat("DP     OSB:", round(res_dp_s$OSB, 3),    "  WhShTot:", res_dp_s$WhShTot, "\n")
#> DP     OSB: 12.739 34.7   WhShTot: 3.511366
cat("COBYLA OSB:", round(res_cob_s$OSB, 3),   "  WhShTot:", res_cob_s$WhShTot, "\n")
#> COBYLA OSB: 13.22 34.7   WhShTot: 3.502752
cat("GLOBAL OSB:", round(res_glob_s$OSB, 3),  "  WhShTot:", res_glob_s$WhShTot, "\n")
#> GLOBAL OSB: 13.282 34.7   WhShTot: 3.502752
cat("EXACT  OSB:", round(res_exact_s$OSB, 3), "  WhShTot:", res_exact_s$WhShTot, "\n")
#> EXACT  OSB: 13.2 34.7   WhShTot: 3.502752

COBYLA, GLOBAL and EXACT reach the same minimum cost. Their printed boundaries can differ slightly, because any boundary that lies between the same two consecutive data values gives the same strata; EXACT reports the largest value in the lower stratum. The DP cost is slightly higher, because the DP optimises the fitted Weibull model rather than the data. The difference here is small (about 0.25%), which shows that the Weibull describes these data well.

8 Working with the Results

This section walks through the tools that help you turn a set of boundaries into a survey design: reading the output, choosing the number of strata, comparing the design with simpler ones, allowing for different costs per stratum, and plotting. The examples use household expenditure from the hies data (\(N = 3{,}566\) households) with a sample of \(n = 400\).

data(hies)
expend <- hies$Expenditure
summary(expend)
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#>    991.2   6707.5  10200.4  13933.7  16449.1 136539.1

8.1 What is in a result

Every call to strata.data() or strata.distr() returns an object of class "strata". Printing it gives a short overview, and summary() gives the full table. The pieces you will use most often are:

Element	Meaning
OSB	the optimum strata boundaries; the last one is the top of the range
Nh, nh	the number of population units and of sampled units in each stratum
Wh, Vh	stratum weights and variances
WhShTot	the Neyman cost \(\sum_h W_h S_h\); smaller means a more precise design
method, distr	the solver used and the fitted (or assumed) distribution
converged, optimality_gap, kkt_residuals	checks on the solution, mainly for "cobyla" and "global"

res <- strata.data(expend, h = 4, n = 400, method = "exact")
res
#> -- stratifyR 2.0 -- Optimal Stratification Results ------------------
#>   Distribution : lnorm
#>   Strata (H)   : 4
#>   Sample size  : n = 400,  N = 3566
#>   Method       : exact
#>   Converged    : yes
#>   Optimality gap (C-S): 0.0000%
#> 
#>   OSB  : 1.008e+04 | 1.964e+04 | 4.179e+04 | 1.365e+05
#>   nh   : 116 | 95 | 105 | 84
#>   Wh   : 0.5000 | 0.3200 | 0.1500 | 0.0300
#>   WhSh : 1010.71300 | 827.02900 | 915.22400 | 730.17800
#>   V*   : 30330.719148
#> ----------------------------------------------------------------------
res$OSB
#> [1]  10083.26  19642.56  41791.24 136539.06
#> attr(,"converged")
#> [1] TRUE
#> attr(,"optimality_gap")
#> [1] 0
#> attr(,"kkt_residuals")
#> [1] NA NA NA
#> attr(,"exact_full")
#> [1] TRUE
#> attr(,"n_candidates")
#> [1] 3559
res$nh
#> [1] 116  95 105  84
sum(res$nh)    # the stratum sample sizes always add up to n
#> [1] 400

A unit whose value equals a boundary belongs to the stratum below it. The sample sizes are whole numbers that add up to \(n\) and never exceed the number of units in a stratum. If Neyman allocation asks for more units than a stratum has, that stratum is fully enumerated and the rest of the sample is shared among the other strata.

8.2 How many strata?

Each extra stratum lowers the Neyman cost, but by less each time, and every stratum adds fieldwork and reporting effort. A quick way to choose \(h\) is to compute WhShTot for a range of values with the fast solver and look for the point where the gains level off.

h_values <- 2:7
cost_h <- sapply(h_values, function(h)
  strata.data(expend, h = h, n = 400, method = "cobyla")$WhShTot)
round(cost_h, 1)
#> [1] 6687.3 4598.9 3483.1 2776.5 2325.9 2015.0
round(100 * (1 - cost_h[-1] / cost_h[-length(cost_h)]), 1)   # % drop at each step
#> [1] 31.2 24.3 20.3 16.2 13.4
plot(h_values, cost_h, type = "b", pch = 19,
     xlab = "Number of strata (h)", ylab = "WhShTot",
     main = "Neyman cost by number of strata")

For these data the cost falls by about 31% from two to three strata and by about 24% from three to four, after which each extra stratum gains less. Four strata is a sensible choice here. The Strata count tab of the app does the same calculation and marks the elbow for you.

8.3 Comparing with simpler designs

compare_designs() takes a result and compares, for the same sample size, the variance of the estimated mean under simple random sampling (SRS), stratified sampling with proportional allocation, and stratified sampling with Neyman allocation.

compare_designs(res)
#> 
#>   -- Design Efficiency Comparison ------------------------------
#>      H = 4 strata   |   n = 400   |   S^2 = 1.60793e+08
#> 
#>   Design                  V(ybar)              SE(ybar)    DEFF
#>   -------------------------------------------------------------
#>   SRS  (baseline)    4.0198e+05   634.021686   1.000
#>   Proportional alloc 5.9051e+04   243.003447   0.1469
#>   Neyman (optimal)   3.0331e+04   174.157168   0.0755
#>   -------------------------------------------------------------
#>   Neyman vs SRS       :  13.25x more efficient
#>   Neyman vs Propnl    :  1.95x more efficient
#>   Equivalent SRS n    :  5302  (vs n = 400 stratified)
#>   Sample-size saving  :  92.5%

The design effect (DEFF) is the variance relative to SRS. The equivalent SRS sample size is how many households a simple random sample would need to match the precision of the stratified design; the difference is the saving that stratification buys.

8.4 Different costs per stratum

If some strata are more expensive to survey, for example because large households take longer to interview, give the cost per sampled unit in each stratum with cost = TRUE and ch. The solvers then minimise \(\sum_h W_h S_h \sqrt{c_h}\), so both the boundaries and the allocation change: cheaper strata receive more of the sample.

ch <- c(1, 1, 2, 4)    # relative cost of one interview in each stratum
res_cost <- strata.data(expend, h = 4, n = 400, method = "exact",
                        cost = TRUE, ch = ch)
rbind(equal_costs = res$nh, with_costs = res_cost$nh)
#>             [,1] [,2] [,3] [,4]
#> equal_costs  116   95  105   84
#> with_costs   179  157   53   11
c(equal_costs = sum(res$nh * ch), with_costs = sum(res_cost$nh * ch))   # field cost
#> equal_costs  with_costs 
#>         757         486

With the same 400 interviews, the cost-aware design moves the sample towards the cheaper lower strata and needs noticeably less field effort. Its variance is a little higher than the Neyman design’s, because Neyman allocation is the most precise choice for a fixed \(n\); the cost-aware design is the better choice when the budget, rather than the sample size, is fixed.

8.5 Plots

plot() shows the population with the boundaries marked. The default type = "2d" is a histogram with the fitted density overlaid. With the plotly package installed, type = "3d" and type = "interactive" give interactive views that you can rotate and explore in the RStudio viewer or a browser.

plot(res, data = expend)

plot(res, type = "3d")            # needs plotly
plot(res, type = "interactive")   # needs plotly

8.6 Large frames and the exact solver

For frames with more than 20,000 distinct values, method = "exact" searches a large, evenly spread set of cut points, and res$exact_full is FALSE. If you need the full search, raise the limit before the call, keeping in mind that time and memory grow quickly with it:

options(stratifyR.exact_max_cand = 50000)
res_big <- strata.data(big_frame, h = 5, n = 2000, method = "exact")
res_big$exact_full

9 Application to Numerous Survey Populations

As discussed earlier, the stratifyR package is able to handle ten continuous distributions that are quite commonly-used in real-life situations. This section presents a brief overview of these distributions and the application of the proposed method of stratification using real or simulated data which follows a particular distribution. Examples for hypothetical distributions are also presented. For the sake of brevity, the mathematical formulations of the problem of determining the OSB and the DP solution procedure are presented only for the Pareto Type II variable. For all other distributions, only the examples are presented to illustrate the application of the package.

9.1 Stratification for a Survey Variable with Pareto Type II Distribution

Let the study variable follow the Pareto Type II distribution on the domain of [\(0, +\infty\)], its two-parameter probability density function with a state space \(y\geq 0\) is given by: \[\begin{equation} f(y; s,a)=\dfrac{a\,s^{a}}{(y+s)^{a+1}}, \;\;\;\;\;a,s > 0 \tag{15} \end{equation}\]

where \(\alpha > 0\) is the shape parameter and \(s>0\) is the scale parameter of the distribution.

9.1.1 MPP Formulation for Pareto Type II Distribution

If the study variable \(y\) follows Pareto II (or Lomax) distribution (i.e., \(y\sim P(a, s)\)) with density function given by (15). By using equations (4)-(6), the formulated MPP given in (10) could be expressed as: \[\begin{eqnarray} \textrm{Minimize} \;\;\; \sum_{h=1}^{L} \textit{SQRT} && \Biggl\{as^{2a} \left[\dfrac{(y_{h-1}+l_{h}+s)^{a}-(y_{h-1}+s)^{a}}{(y_{h-1}+s)^{a}(y_{h-1}+l_{h}+s)^{a}}\right] \nonumber \\[10pt] && \times\left[\dfrac{(y_{h-1}+l_{h}+s)^{2-a}}{2-a} - \dfrac{2s(y_{h-1}+l_{h}+s)^{1-a}}{1-a} \right. \nonumber\\[10pt] && - \dfrac{s^2(y_{h-1}+l_{h}+s)^{-a}}{a} -\dfrac{(y_{h-1}+s)^{2-a}}{2-a}\nonumber\\[10pt] && + \left. \dfrac{2s(y_{h-1}+s)^{1-a}}{1-a} + \dfrac{s^2(y_{h-1}+s)^{-a}}{a}\right] \nonumber\\[10pt] && \times \dfrac{s^{2a}}{(1-a)^{2}}\left[\dfrac{a(y_{h-1}+l_{h})+s}{(y_{h-1}+l_{h}+s)^{a}} - \dfrac{ay_{h-1}+s}{(y_{h-1}+s)^{a}} \right]^{2}\Biggr\}\nonumber\\[15pt] \textrm{subject to} && \;\;\;\; \sum_{h=1}^{L}l_{h} = d,\nonumber \\ \textrm{and} && \;\;\;\; l_{h} \geq 0; \;\; h = 1,2,...,L. \tag{16} \end{eqnarray}\]

where \(d=y_{L}-y_{0}\), \(a\) and \(s\) are parameters of the Pareto Type II distribution.

9.1.2 DP Solution for Pareto Type II Distribution

To solve the MPP (16) using the DP technique as a solution procedure, we apply the algorithm, that is, the solution proceure using Dynamic Progrmming technique discussed earlier in Section 4. After substituting the quatity \(y_{h-1}=y_{0} + d_{h} - l_{h}\), the recurrence relations (11) and (12) are reduced to:

For the first stage, \(k=1\), at \(l_{1}^{*}=d_{1}\):

\[\begin{eqnarray} \Phi_{1}d_{1}&=& \textrm{Sqrt}\Biggl\{ as^{2a}\left[\dfrac{(d_{1}+y_{0}+s)^{a}-(y_{0}+s)^{a}}{(y_{0}+s)^{a}(d_{1}+y_{0}+s)^{a}}\right] \nonumber\\[10pt] &&\times\left[\dfrac{(d_{1}+y_{0}+s)^{2-a}}{2-a} - \dfrac{2s(d_{1}+y_{0}+s)^{1-a}}{1-a} \right. \nonumber\\[10pt] &&- \dfrac{s^2(d_{1}+y_{0}+s)^{-a}}{a} -\dfrac{(y_{0}+s)^{2-a}}{2-a}\nonumber\\[10pt] &&+\left. \dfrac{2s(y_{0}+s)^{1-a}}{1-a} + \dfrac{s^2(y_{0}+s)^{-a}}{a}\right] \nonumber\\[10pt] &&\times \dfrac{s^{2a}}{(1-a)^{2}}\left[\dfrac{a(d_{1}+y_{0})+s}{(d_{1}+y_{0}+s)^{a}} - \dfrac{ay_{0}+s}{(y_{0}+s)^{a}} \right]^{2}\Biggr\} \tag{17} \end{eqnarray}\]

And for the stages \(k\geq2\):

\[\begin{eqnarray} \Phi_{k}d_{k}&=&\min_{0\leq l_{k}\leq d_{k}}\Biggl\{ \textrm{Sqrt}\Biggl\{ as^{2a}\left[\dfrac{(d_{k}+y_{0}+s)^{a}-(d_{k}+l_{k}+y_{0}+s)^{a}}{(d_{k}+l_{k}+y_{0}+s)^{a}(d_{k}+y_{0}+s)^{a}}\right] \nonumber\\[10pt] &&\times\left[\dfrac{(d_{k}+y_{0}+s)^{2-a}}{2-a} - \dfrac{2s(d_{k}+y_{0}+s)^{1-a}}{1-a} \right. \nonumber\\[10pt] &&- \dfrac{s^2(d_{k}+y_{0}+s)^{-a}}{a} -\dfrac{(d_{k}+l_{k}+y_{0}+s)^{2-a}}{2-a}\nonumber\\[10pt] &&+\left. \dfrac{2s(d_{k}+l_{k}+y_{0}+s)^{1-a}}{1-a} + \dfrac{s^2(d_{k}+l_{k}+y_{0}+s)^{-a}}{a}\right] \nonumber\\[10pt] &&\times \dfrac{s^{2a}}{(1-a)^{2}}\left[\dfrac{a(d_{k}+y_{0})+s}{(d_{k}+y_{0}+s)^{a}}- \dfrac{a(d_{k}+l_{k}+y_{0})+s}{(d_{k}+l_{k}+y_{0}+s)^{a}} \right]^{2} \Biggr\} + \Phi_{k-1}(d_{k}-l_{k})\Biggr\}. \tag{18} \end{eqnarray}\]

The recurrence relations (17) and (18) are solved using the DP technique to determine the OSB.

9.1.3 A Numerical Example for Pareto Type II Distribution

A dataset for a univariate population of size \(N=5000\) with the study variable that follows Pareto Type II distribution (\(pareto\_data\)) was simulated using parameters \(shape=5\) and \(scale=8\) to demonstrate the application of the strata.data() function to determine the OSB and other quantites. The data exhibits a 2-parameter Pareto Type II distribution with the MLE estimates of the parameters as \(shape=5.026907\) and \(scale=8.191676\). The minimum and maximum values in the simulated data are \([y_{0}, y_{L}] = [0.0002193, 38.56871]\), which implies that \(d=38.56849\).

To construct the OSB (a 2-strata solution, i.e., \(h = 2\)) for the \(pareto\_data\) with a fixed total sample size of \(500\), we use the following codes:

set.seed(8235411)
pareto_data <- actuar::rpareto(5000, shape=5, scale=8)

dpareto <- actuar::dpareto
ppareto <- actuar::ppareto
qpareto <- actuar::qpareto

head(pareto_data)
#> [1] 1.8464786 0.9226109 1.1079785 9.7699864 4.5600912 0.0441124
hist(pareto_data, breaks=100)

min(pareto_data); max(pareto_data); d=max(pareto_data)-min(pareto_data);d
#> [1] 0.0002192696
#> [1] 38.56871
#> [1] 38.56849
fit <- fitdistrplus::fitdist(pareto_data, "pareto", start = list(shape = 1, scale = 500))
fit
#> Fitting of the distribution ' pareto ' by maximum likelihood 
#> Parameters:
#>       estimate Std. Error
#> shape 5.026907  0.4337687
#> scale 8.191676  0.8358576
res <- strata.data(pareto_data, h = 2, n=500) # a 2-strata solution
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 500 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 5000</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [0.0002, 38.5687]  (d = 38.5685)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> weibull</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter | Estimate|
#> |:---------|--------:|
#> |shape     | 0.875786|
#> |scale     | 1.895872|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	2.72604	0.76	0.54	0.560	205	3803	0.05
2	38.5687	0.24	11.34	0.806	295	1197	0.25
Total		1.00	11.88	1.366	500	5000	0.10

Similarly, in order to find the OSB and other quantities, we can apply the strata.distr() function to a Pareto Type II population. Let us assume from past knowledge that the study variable in the population follows Pareto Type II distribution with the given attributes such as the \(shape=5.05, scale=8.20, initial\,value=0.15\) and \(distance=38.55\). Then, if a sample of size \(n=500\) is drawn from the population of size \(N=5000\), we can execute the following command to obtain the results:

res <- strata.distr(h=2, initval=0.15, dist=38.55, distr = "pareto",
             params = c(shape=5.05, scale=8.20), n=500, N=5000)
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 500 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 5000</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [0.1500, 38.7000]  (d = 38.5500)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> pareto</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter | Estimate|
#> |:---------|--------:|
#> |shape     |     5.05|
#> |scale     |     8.20|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	3.20939	0.81	0.67	0.592	241	4057	0.06
2	38.7	0.19	11.42	0.637	259	943	0.27
Total		1.00	12.09	1.229	500	5000	0.10

9.2 Stratification for a Survey Variable with Triangular Distribution

Let the study variable follow Triangular distribution on the domain of [\(a, b\)], its three-parameter probability density function with a state space \(y\geq 0\) is given by: \[\begin{equation} f(y) = \begin{cases} \dfrac{2(y-a)}{(b-a)(c-a)}; & y \in [a, c]\\[10pt] \dfrac{2(b-y)}{(b-a)(b-c)}; & y \in (c,b] \\ \end{cases} \tag{19} \end{equation}\]

where \(a\) is the location parameter, \(b\) is the scale parameter and \(c\) is the shape parameter of the distribution.

Then, formulating the problem as an MPP and solving the recurrence relations as discussed above for Pareto II variate, the OSB are obtained.

9.2.1 DP Solution for Triangular Distribution

To solve the MPP formulated for Triangular distribution (19), we apply the algorithm using the DP technique discussed in Section 4. The recurrence relations used to determine the OSB are given by:

For the first stage, \(k=1\), at \(l_{1}^{*}=d_{1}\):

\[\begin{eqnarray} \Phi_{1}d_{1}&=& \textrm{Sqrt}\Biggl\{ \sum_{h=1}^{\lambda_{1}} \dfrac{d_{1}^{2}\sqrt{d_{1}^{2}+6(y_{0}-a)d_{1}+6(y_{0}-a)^{2}}}{3\sqrt{2}(b-a)(c-a)} + \sum_{h=\lambda_{2}}^{L}\dfrac{d_{1}^{2}\sqrt{6(b-y_{0})^{2}-6(b-y_{0})d_{1}+d_{1}^{2}}}{3\sqrt{2}(b-a)(b-c)}\Biggr\} \tag{20} \end{eqnarray}\]

And for the stages \(k\geq2\):

\[\begin{eqnarray} \Phi_{k}d_{k}&=&\min_{0\leq l_{k}\leq d_{k}}\Biggl\{ \textrm{Sqrt}\Biggl\{\sum_{h=1}^{\lambda_{1}} \dfrac{l_{k}^{2}\sqrt{l_{k}^{2}+6a_{k}l_{k}+6a_{k}^{2}}}{3\sqrt{2}(b-a)(c-a)} +\sum_{h=\lambda_{2}}^{L}\dfrac{l_{k}^{2}\sqrt{6b_{k}^{2}-6b_{k}l_{k}+l_{k}^{2}}}{3\sqrt{2}(b-a)(b-c)}\Biggr\} + \Phi_{k-1}(d_{k}-l_{k})\Biggr\} \tag{21} \end{eqnarray}\]

where \(a_{k}=d_{k}-l_{k}+y_{0}-a\) and \(b_{k}=b-d_{k}+l_{k}-y_{0}+a\).

Substituting the values of \(a\), \(b\), \(c\), \(y_{0}\) and \(d\), the OSW (\(l_{h}^{*}\)) and the OSB (\(y_{h}^{*}=y_{h-1}^{*}-l_{h}^{*}\)) are obtained by executing the strata.data() function.

9.2.2 A Numerical Example for Triangular Distribution

Data on Mathematics marks of first year students in a University in Fiji, with a size of \(N=487\) called ‘math’ data is used to demonstrate the application of the stratifyR package on a Triangular population. In this example, the variable `final_marks’ is used - it exhibits a 3-parameter Triangular distribution. The minimum and maximum values in the ‘math’ data are \([y_{0}, y_{L}] = [7, 97]\), which implies that \(d=90\). Fitting a Triangular distribution by fixing \(\min\) and \(\max\) just outside the sample bounds and estimating only the mode by maximum likelihood gives \(mode \approx 54\). (The triangular density is zero at its end points, so the likelihood cannot be evaluated if \(\min\) and \(\max\) equal the smallest and largest marks.) We therefore use \(min = 7\), \(max = 97\) and \(mode = 54\).

To construct the OSB (\(h = 2\)) for the `final_marks’ data with a fixed total sample size of \(150\), we use the following code:

data(math)
final_marks <- math$final_marks
hist(final_marks)

res <- strata.data(final_marks, h = 2, n=150) # a 2-strata solution
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 150 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 487</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [7.0000, 97.0000]  (d = 90.0000)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> weibull</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter |  Estimate|
#> |:---------|---------:|
#> |shape     |  3.371526|
#> |scale     | 59.909598|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	53.4373	0.49	114.72	5.278	74	240	0.31
2	97	0.51	116.01	5.463	76	247	0.31
Total		1.00	230.73	10.741	150	487	0.31

Similarly, in order to find the OSB and other quantities, we can apply the strata.distr() function to a Triangular population. Based on the assumption from past knowledge that the population follows Triangular distribution with the given attributes such as the initial value, distance, parameters, etc., we can execute the following command to obtain the results:

data(math)
final_marks <- math$final_marks
a <- min(final_marks); b <- max(final_marks)
a; b; d <- b - a; d
#> [1] 7
#> [1] 97
#> [1] 90

# Estimate the mode by maximum likelihood, with min and max fixed just outside
# the sample range (the density is zero at min and max, so they cannot equal
# sample values). The likelihood of a triangular mode is maximised at one of
# the data values, so we simply evaluate it at each distinct value. This is
# exact and gives the same answer on every platform.
lo <- a - 0.01 * d; hi <- b + 0.01 * d
loglik <- function(m) sum(log(mc2d::dtriang(final_marks, min = lo, mode = m, max = hi)))
cand <- sort(unique(final_marks))
mode_hat <- cand[which.max(sapply(cand, loglik))]
mode_hat
#> [1] 54

# 2-strata solution
res <- strata.distr(h=2, initval=7, dist=90, distr = "triangle",
      params = c(min=7, max=97, mode=54), n=150, N=352)
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 150 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 352</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [7.0000, 97.0000]  (d = 90.0000)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> triangle</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter | Estimate|
#> |:---------|--------:|
#> |min       |        7|
#> |max       |       97|
#> |mode      |       54|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	53.0008	0.5	117.56	5.424	77	176	0.44
2	97	0.5	107.45	5.180	73	176	0.41
Total		1.0	225.01	10.604	150	352	0.43

9.3 Stratification for a Survey Variable with Right-Triangular Distribution

If a study variable follows the Right-Triangular distribution on the domain of [\(a, b\)], its two-parameter probability density function is given by:

\[\begin{equation} f(y) = \begin{cases} \dfrac{2(b-y)}{(b-a)^2}; & y \in [a, b]\\ 0; & otherwise \\ \end{cases} \tag{22} \end{equation}\]

where \(a\) is the location parameter and \(b\) is the scale parameter of the distribution.

Note that the Right-Triangular distribution is a special case of the Triangular distribution discussed in Section 9.2 where the parameters \(a=c\), i.e., minimium value is equal to the mode. Thus, the density function (22) is a two-parameter distribution.

9.3.1 DP Solution for Right-Triangular Distribution

To solve the MPP formulated for Right-Triangular distribution (22), we apply the algorithm using the DP technique discussed in Section 4. The recurrence relations used to determine the OSB are given by:

For the first stage, \(k=1\), at \(l_{1}^{*}=d_{1}\):

\[\begin{eqnarray} \Phi_{1}d_{1}&=& \textrm{Sqrt}\Biggl\{ \left[\dfrac{d_{1}(2(y_{0}-a)-d_{1})}{(b-a)^{2}}\right]^{2}\times\left[\dfrac{d_{1}^{2}(d_{1}^{2}-6(y_{0}-a)d_{1}+6a_{h}^{2})}{18(2(y_{0}-a)-d_{1})^{2}}\right]\Biggr\} \tag{23} \end{eqnarray}\]

And for the stages \(k\geq2\):

\[\begin{eqnarray} \Phi_{k}d_{k}&=&\min_{0\leq l_{k}\leq d_{k}}\Biggl\{ \textrm{Sqrt}\Biggl\{\left[\dfrac{l_{k}(2a_{k}-l_k)}{(b-a)^{2}}\right]^{2} \times\left[\dfrac{l_{k}^{2}(l_{k}^{2}-6a_{k}l_{k}+6a_{k}^{2})}{18(2a_{k}-l_{k})^{2}}\right]\Biggr\} + \Phi_{k-1}(d_{k}-l_{k})\Biggr\} \tag{24} \end{eqnarray}\]

where \(a_{k}=b-d_{k}+l_{k}-y_{0}\).

Substituting the values of \(a\), \(b\), \(y_{0}\) and \(d\), the OSW (\(l_{h}^{*}\)) and the OSB (\(y_{h}^{*}=y_{h-1}^{*}-l_{h}^{*}\)) are obtained by executing the strata.data() function.

9.3.2 A Numerical Example for Right-Triangular Distribution

A data following Right-Triangular distribution, of size \(N=5000\), was simulated to demonstrate the application of the stratifyR package. The simulated data takes in three parameters where \(a=c\) to indicate that it’s a Right-Triangular distribution. Upon fitting the data, its best-fit distribution exhibits a three-parameter Triangular distribution with the parameters \(min = 1.987776\), \(max=7.935599\) and \(mode=2.026685\). Note that because the \(min\) and \(mode\) parameters are very close to each and not exactly the same, it is treated as a Triangular distribution even thought the data simulated was for a Right-Triangular distribution. The minimum and maximum values in the simulated data are \([y_{0}, y_{L}] = [2.000052, 7.83871]\) with \(d=5.838658\).

To construct the OSB (\(h = 2\)) for the Right-Triangular data with a fixed total sample size of \(500\), we use the following code:

#Generate RT data
set.seed(12546)
data <- rtriangle(n=1000, a=2, b=8, c=2) #right-triangular since a=c
hist(data)

res <- strata.data(data, h = 2, n=500) # a 2-strata solution
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 500 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 1000</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [2.0001, 7.8387]  (d = 5.8387)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> lnorm</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter | Estimate|
#> |:---------|--------:|
#> |meanlog   | 1.327465|
#> |sdlog     | 0.353496|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	4.17145	0.58	0.38	0.359	246	582	0.42
2	7.83871	0.42	0.78	0.369	254	418	0.61
Total		1.00	1.16	0.727	500	1000	0.50

We see in the above example that the simulated data is a Right-Triangular distribution, however, when data is fitted using MLE method in the package, it turns out that it best-fits Triangular distribution because the \(min\) is not exactly equal to \(mode\).

Similarly, in order to find the OSB and other quantities, we can apply the strata.distr() function to a Right-Triangular population. Based on the assumption from past knowledge that the population follows Right-Triangular distribution with the given attributes such as the initial value, distance, parameters, etc., we can execute the following command to obtain the results:

res <- strata.distr(h=2, initval=1.007202, dist=0.992781, distr = "rtriangle",
         params = c(min=2, max=10, mode=2), n=500, N=1000)
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 500 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 1000</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [1.0072, 2.0000]  (d = 0.9928)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> rtriangle</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter | Estimate|
#> |:---------|--------:|
#> |min       |        2|
#> |max       |       10|
#> |mode      |        2|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	1.49635	-0.13	0.02	0.019	-130	-130	1.00
2	1.99998	1.13	0.02	0.019	630	1130	0.56
Total		1.00	0.04	0.038	500	1000	0.50

The results show that this fits a two-paramter Right-Triangular distribution. Do note that in the above command, one has to specify all three parameters \(min=2, max=10\) and \(mode=2\), even for a Right-Triangular distribution, where \(min=mode\).

9.4 Stratification for a Survey Variable with Weibull Distribution

If the study variable follows the Weibull distribution on the interval \([y_{0}, y_{L}]\), its two-parameter probability density function with a state space \(y\geq 0\) is given by: \[\begin{equation} f(y; \theta, r)=\dfrac{r}{\theta}\left(\dfrac{y}{\theta}\right)^{r-1}e^{-\left(y/\theta\right)^{r}}, \;\;\;\;\;y\geq 0 \tag{25} \end{equation}\]

where \(r > 0\) is the shape parameter and \(\theta > 0\) is the scale parameter of the distribution.

9.4.1 DP Solution for Weibull Distribution

To solve the MPP formulated for Weibull distribution (25), we apply the algorithm using the DP technique discussed in Section 4. The recurrence relations involve the Upper Regularized Incomplete Gamma function \(Q(s, z)\), defined as: \[\begin{equation} Q(s,z)=\dfrac{1}{\Gamma(s)}\int_{z}^{\infty}t^{s-1}e^{-t}\,dt,\;\;\;\;\;s,z>0,\tag{25a} \end{equation}\] which in R is computed as . The recurrence relations used to determine the OSB are given by:

For the first stage, \(k=1\), at \(l_{1}^{*}=d_{1}\):

\[\begin{eqnarray} \Phi_{1}d_{1}&=& \textrm{Sqrt}\Biggl\{ \theta^{2}\,\Gamma\left(1+\dfrac{2}{r} \right)\left[e^{-\left(\frac{y_{0}}{\theta}\right)^{r}}-e^{-\left(\frac{d_{1}+y_{0}}{\theta}\right)^{r}}\right] \nonumber \\[10pt] &&\times \left[Q\left(1+\dfrac{2}{r},\left(\dfrac{y_{0}}{\theta}\right)^{r}\right)- Q\left(1+\dfrac{2}{r},\left(\dfrac{d_{1}+y_{0}}{\theta}\right)^{r}\right)\right]\nonumber\\[10pt] &&-\theta^{2}\left[\Gamma\left(1+\dfrac{1}{r}\right)\left [Q\left(1+\dfrac{1}{r},\left(\dfrac{y_{0}}{\theta}\right)^{r}\right)\right.\right.\nonumber\\[10pt] &&-Q\left.\left.\left(1+\dfrac{1}{r},\left(\dfrac{d_{1}+y_{0}}{\theta}\right)^{r}\right)\right]\right]^{2}\Biggr\} \tag{26} \end{eqnarray}\]

And for the stages \(k\geq2\):

\[\begin{eqnarray} \Phi_{k}d_{k}&=&\min_{0\leq l_{k}\leq d_{k}}\Biggl\{ \textrm{Sqrt}\Biggl\{ \theta^{2}\,\Gamma\left(1+\dfrac{2}{r} \right)\left[e^{-\left(\frac{d_{k}-l_{k}+y_{0}}{\theta}\right)^{r}}-e^{-\left(\frac{d_{k}+y_{0}}{\theta}\right)^{r}}\right]\nonumber\\[10pt] &&\times \left[Q\left(1+\dfrac{2}{r},\left(\dfrac{d_{k}-l_{k}+y_{0}}{\theta}\right)^{r}\right)- Q\left(1+\dfrac{2}{r},\left(\dfrac{d_{k}+y_{0}}{\theta}\right)^{r}\right)\right]\nonumber\\[15pt] &&-\theta^{2}\left[\Gamma\left(1+\dfrac{1}{r}\right)\left [Q\left(1+\dfrac{1}{r},\left(\dfrac{d_{k}-l_{k}+y_{0}}{\theta}\right)^{r}\right)\right.\right.\nonumber\\[15pt] &&-Q\left.\left.\left(1+\dfrac{1}{r},\left(\dfrac{d_{k}+y_{0}}{\theta}\right)^{r}\right)\right]\right]^{2}\Biggr\}+\Phi_{k-1}(d_{k}-l_{k})\Biggr\}. \tag{27} \end{eqnarray}\]

The recurrence relations (26) and (27) are solved using the DP technique to determine the OSB.

9.4.2 A Numerical Example for Weibull Distribution

A health data of size \(N=724\), called ‘anaemia’ data (Reddy et al. (2014)), is used to demonstrate the application of the stratifyR package on Weibull population. The ‘anaemia’ data comes from the National Nutritional Survey on the ``Micronutrient Status of Women in Fiji” and has many variables such as level of Iron, Folate, Zinc, etc. In this example, the variable Iron is used since it exhibits a 2-parameter Weibull distribution with the shape and scale parameters as \(r = 2.144586\) and \(\theta = 13.790744\) respectively. The minimum and maximum values are \([y_{0}, y_{L}] = [1.5, 34.7]\), which implies that \(d=33.2\).

To construct the OSB (\(h = 2\)) for the Iron data with a fixed total sample size of \(500\), we use the following codes:

data(anaemia) #using the anaemia data
Iron <- anaemia$Iron
hist(Iron)

min_value <- min(Iron) # Find the minimum value
Iron[Iron == min_value] <- -0.001 # Replace the minimum value with -0.001
res <- strata.data(Iron, h = 2, n=500) # a 2-strata solution
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 500 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 724</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [-0.0010, 34.7000]  (d = 34.7010)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> norm</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter |  Estimate|
#> |:---------|---------:|
#> |mean      | 12.194472|
#> |sd        |  6.026547|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	12.5777	0.58	8.83	1.724	244	420	0.58
2	34.7	0.42	18.49	1.805	256	304	0.84
Total		1.00	27.32	3.530	500	724	0.69

Similarly, in order to find the OSB and other quantities, we can apply the strata.distr() function to a Weibull population. Based on the assumption from past knowledge that the population follows Weibull distribution with the given attributes such as the initial value, distance, parameters, etc., we can execute the following command to obtain the results:

res <- strata.distr(h=2, initval=2.9, dist=55.9, distr = "weibull",
       params = c(shape=2.144586, scale=13.790744), n=500, N=5000)
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 500 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 5000</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [2.9000, 58.8000]  (d = 55.9000)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> weibull</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter |  Estimate|
#> |:---------|---------:|
#> |shape     |  2.144586|
#> |scale     | 13.790744|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	13.0489	0.59	7.49	1.516	239	2943	0.08
2	58.8	0.41	16.28	1.660	261	2057	0.13
Total		1.00	23.77	3.176	500	5000	0.10

9.5 Stratification for a Survey Variable with Gamma Distribution

If the study variable follows the Gamma distribution (i.e., \(y\sim \Gamma(r, \theta)\)) on the interval \([y_{0}, y_{L}]\), it has the following two-parameter probability density function:

\[\begin{equation} f(y;r,\theta)=\dfrac{1}{\theta^{r}\Gamma(r)}\,y^{r-1}e^{-\dfrac{y}{\theta}},\;\;\;\;\;\;\;y>0;\;\;r,\theta>0, \tag{28} \end{equation}\]

where is a shape parameter and \(\theta\) is the scale parameter and \(\Gamma(r)\) is a Gamma function defined by

\[\begin{equation}\label{222} \Gamma(r)=\int_{0}^{\infty}t^{r-1}e^{-t}\,dt,\;\;\;\;\;\; r > 0. \tag{29} \end{equation}\]

The function in equation \((\ref{222})\) is also defined by an upper incomplete gamma function \(\Gamma(r,x)\) and a lower incomplete gamma function \(\gamma(r,x)\), respectively, as follows:

\[\begin{eqnarray} \Gamma(r,y)=\int_{y}^{\infty}t^{r-1}e^{-t}\,dt;\label{223} \tag{30}\\ \gamma(r,y)=\int_{0}^{y}t^{r-1}e^{-t}\,dt.\label{224} \tag{31} \end{eqnarray}\]

There also exist regularized/normalised incomplete Gamma functions which give a value restricted between \(0\) and \(1\) and can be stated as:

\[\begin{eqnarray} Q(r,y)=\dfrac{1}{\Gamma(r)}\int_{y}^{\infty}t^{r-1}e^{-t}\,dt,\;\;\;\;\;\;r,y>0;\;\;\Gamma(r)\neq0;\label{225} \tag{32}\\ P(r,y)=\dfrac{1}{\Gamma(r)}\int_{0}^{y}t^{r-1}e^{-t}\,dt,\;\;\;\;\;\;r,y>0;\;\;\Gamma(r)\neq0,\label{226} \tag{33} \end{eqnarray}\]

where \(Q(r,y)\) denotes the Upper Regularized Incomplete Gamma function while \(P(r,y)\) denotes the Lower Regularized Incomplete Gamma function (Abramowitz and Stegun (1972), Chapter 6). Note that \(Q(r,y)=1-P(r,y)\). In R, \(Q(r,y)\) is evaluated as , which is the regularized form. It is important to note that \(Q(r,y)\) is distinct from the unregularized upper incomplete gamma function \(\Gamma(r,y)\) (equation 30); they are related by \(Q(r,y)=\Gamma(r,y)/\Gamma(r)\). The DP recurrence relations for both Weibull (equations 26–27) and Gamma (equations 34–35) use the regularized \(Q(r,y)\), and stratifyR 2.0 implements these correctly using .

9.5.1 DP Solution for Gamma Distribution

To solve the MPP formulated for Gamma distribution (28), we apply the algorithm using the DP technique discussed in Section 4. The recurrence relations used to determine the OSB are given by:

For the first stage, \(k=1\), at \(l_{1}^{*}=d_{1}\):

\[\begin{eqnarray} \Phi_{1}d_{1}&=& \textrm{Sqrt}\Biggl\{ \theta^{2}\,r(r+1)\left[ Q\left(r,\dfrac{y_{0}}{\theta}\right)-Q\left(r,\dfrac{d_{1}+y_{0}}{\theta}\right)\right]\nonumber\\[10pt] &&\times\left[Q\left(r+2,\dfrac{y_{0}}{\theta}\right)-Q\left(r+2,\dfrac{d_{1}+y_{0}}{\theta}\right) \right]\nonumber\\[10pt] &&-\theta^{2}r^{2}\left[ Q\left(r+1,\dfrac{y_{0}}{\theta}\right)-Q\left(r+1,\dfrac{d_{1}+y_{0}}{\theta}\right)\right]^{2}\Biggr\} \tag{34} \end{eqnarray}\]

And for the stages \(k\geq2\):

\[\begin{eqnarray} \Phi_{k}d_{k}&=&\min_{0\leq l_{k}\leq d_{k}}\Biggl\{ \textrm{Sqrt}\Biggl\{ \theta^{2}\,r(r+1)\left[ Q\left(r,\dfrac{d_{k}-l_{k}+y_{0}}{\theta}\right.\right)\nonumber\\[10pt] &&-\left.Q\left(r,\dfrac{d_{k}+y_{0}}{\theta}\right)\right] \times\left[ Q\left(r+2,\dfrac{d_{k}-l_{k}+y_{0}}{\theta}\right.\right)\nonumber\\[10pt] &&-\left.Q\left(r+2,\dfrac{d_{k}+y_{0}}{\theta}\right) \right]-\theta^{2}r^{2}\times\left[ Q\left(r+1,\dfrac{d_{k}-l_{k}+y_{0}}{\theta}\right)\right.\nonumber\\[10pt] &&-\left.Q\left(r+1,\dfrac{d_{k}+y_{0}}{\theta}\right)\right]^{2}\Biggr\}+\Phi_{k-1}(d_{k}-l_{k})\Biggr\}. \tag{35} \end{eqnarray}\]

The recurrence relations (34) and (35) are solved using the DP technique to determine the OSB.

9.5.2 A Numerical Example for Gamma Distribution

Again, the health data of size \(N=724\), derived from the ‘National Nutritional Survey’ on the ‘Micronutrient Status of Women in Fiji’ is used to demonstrate the application of the stratifyR package on Gamma population. In this example, the variable Folate is used since it exhibits a 2-parameter Gamma distribution with the shape and scale parameters as \(r = 6.9922\) and \(\theta = 2.5785\) respectively. The minimum and maximum values are \([y_{0}, y_{L}] = [4.9, 45.4]\), which implies that \(d=40.5\).

To construct the OSB (\(h = 2\)) for the Folate data with a fixed total sample size of \(500\), we use the following codes:

data(anaemia)
Folate <- anaemia$Folate
hist(Folate)

min_value <- min(Folate) # Find the minimum value
Folate[Folate == min_value] <- -0.001 # Replace the minimum value with -0.001
res <- strata.data(Folate, h = 2, n=500) # a 2-strata solution
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 500 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 724</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [-0.0010, 45.4000]  (d = 45.4010)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> norm</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter |  Estimate|
#> |:---------|---------:|
#> |mean      | 18.017124|
#> |sd        |  6.985322|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	18.1314	0.56	9.85	1.752	210	404	0.52
2	45.4	0.44	29.95	2.419	290	320	0.91
Total		1.00	39.80	4.170	500	724	0.69

Similarly, in order to find the OSB and other quantities, we can apply the strata.distr() function to a Gamma population. Based on the assumption from past knowledge that the population follows Gamma distribution with the given attributes such as the initial value, distance, parameters, etc., we can execute the following command to obtain the results:

res <- strata.distr(h=2, initval=0.5, dist=50, distr = "gamma",
       params = c(shape=3.835768, rate=0.340328), n=500, N=12000)
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 500 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 12000</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [0.5000, 50.5000]  (d = 50.0000)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> gamma</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter | Estimate|
#> |:---------|--------:|
#> |shape     | 3.835768|
#> |rate      | 0.340328|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	12.1776	0.63	6.95	1.653	247	7522	0.03
2	50.5	0.37	20.50	1.689	253	4478	0.06
Total		1.00	27.45	3.342	500	12000	0.04

9.6 Stratification for a Survey Variable with Exponential Distribution

If the study variable follows the Exponential distribution on the interval \([y_{0}, y_{L}]\), its one-parameter probability density function is given by:

\[\begin{equation} f(y; \lambda)= \begin{cases} \lambda\mathrm{e}^{-\lambda y}; & y > 0\\ \;\;\;\;0; & elsewhere \\ \end{cases} \tag{36} \end{equation}\]

where \(\lambda\) is the continuous rate parameter (or the inverse scale parameter).

9.6.1 DP Solution for Exponential Distribution

To solve the MPP formulated for Exponential distribution (36), we apply the algorithm using the DP technique discussed in Section 4. The recurrence relations used to determine the OSB are given by:

For the first stage, \(k=1\), at \(l_{1}^{*}=d_{1}\):

\[\begin{eqnarray} \Phi_{1}d_{1}&=& \textrm{Sqrt}\Biggl\{ \left(e^{-\lambda y_{0}}\right)^{2}\left[\dfrac{1}{\lambda^2}\left(1-e^{-\lambda l_{k}}\right)^{2}-l_{k}^{2}e^{-\lambda l_{k}}\right]\Biggr\} \tag{37} \end{eqnarray}\]

And for the stages \(k\geq2\):

\[\begin{eqnarray} \Phi_{k}d_{k}&=&\min_{0\leq l_{k}\leq d_{k}}\Biggl\{ \textrm{Sqrt}\Biggl\{ \left(e^{-\lambda (d_{k}-l_{k}+y_{0})}\right)^{2}\left[\dfrac{1}{\lambda^2}\left(1-e^{-\lambda l_{k}}\right)^{2}-l_{k}^{2}e^{-\lambda l_{k}}\right]\Biggr\} + \Phi_{k-1}(d_{k}-l_{k})\Biggr\} \tag{38} \end{eqnarray}\]

Substituting the values of \(\lambda\), \(y_{0}\) and \(d\), the OSW (\(l_{h}^{*}\)) and the OSB (\(y_{h}^{*}=y_{h-1}^{*}-l_{h}^{*}\)) are obtained by executing the strata.data() function.

9.6.2 A Numerical Example for Exponential Distribution

A data following Exponential distribution (herein called \(Exp\) data), of size \(N=10000\) was simulated to demonstrate the application of the stratifyR package on an Exponential population. The data exhibits an Exponential distribution with the parameter \(rate = 1.359205\). The minimum and maximum values in the simulated data are \([y_{0}, y_{L}] = [5.747904e-05, 8.016871]\), which implies that \(d=8.016814\).

To construct the OSB (\(h = 2\)) for the data that follows Exponential distribution with a fixed total sample size of \(500\), we use the following codes:

set.seed(28951)
data <- rexp(5000, rate = 1.36)
hist(data)

res <- strata.data(data, h = 2, n=500) # a 2-strata solution
#> The program is running, it'll take some time!
#> Warning in mledist(data, distname, start = arg_startfix$start.arg, fix.arg =
#> arg_startfix$fix.arg, : The BFGS method cannot be used with bounds without
#> provided the gradient. The method is changed to L-BFGS-B.
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 500 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 5000</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [0.0001, 8.0169]  (d = 8.0168)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> exp</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter | Estimate|
#> |:---------|--------:|
#> |rate      | 1.359205|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	0.928033	0.72	0.07	0.187	235	3596	0.07
2	8.01687	0.28	0.57	0.212	265	1404	0.19
Total		1.00	0.64	0.399	500	5000	0.10

Similarly, in order to find the OSB and other quantities, we can apply the strata.distr() function to a Exponential population. Based on the assumption from past knowledge that the population follows Exponential distribution with the given attributes such as the initial value, distance, parameters, etc., we can execute the following command to obtain the results:

set.seed(28951)
data <- rexp(5000, rate = 1.36)
min(data); max(data); d=max(data)-min(data);d
#> [1] 5.747904e-05
#> [1] 8.016871
#> [1] 8.016814
fit <- fitdist(data, distr="exp", method="mle")
fit
#> Fitting of the distribution ' exp ' by maximum likelihood 
#> Parameters:
#>      estimate Std. Error
#> rate 1.359205 0.01922205
res <- strata.distr(h=2, initval=5.748e-05, dist=8.017, distr = "exp", 
             params = c(rate=1.36), n=500, N=5000) #a 2-strata solution
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 500 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 5000</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [0.0001, 8.0171]  (d = 8.0170)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> exp</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter | Estimate|
#> |:---------|--------:|
#> |rate      |     1.36|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	0.927493	0.72	0.07	0.185	235	3584	0.07
2	8.01706	0.28	0.54	0.208	265	1416	0.19
Total		1.00	0.60	0.392	500	5000	0.10

9.7 Stratification for a Survey Variable with Uniform Distribution

If the study variable follows the Uniform distribution on the interval \([y_{0}, y_{L}]\), its two-parameter probability density function is given by:

\[\begin{equation} f(y; a,b)= \begin{cases} \dfrac{1}{b-a}; & y > 0\\ \;\;\;\;0; & otherwise \\ \end{cases} \tag{39} \end{equation}\]

where \(a\) and \(b\) are the continuous boundary parameters, i.e., \(a\) and \(b\) are the minimum and maximum values respectively.

9.7.1 DP Solution for Uniform Distribution

To solve the MPP formulated for Uniform distribution (39), we apply the algorithm using the DP technique discussed in Section 4. The recurrence relations used to determine the OSB are given by:

For the first stage, \(k=1\), at \(l_{1}^{*}=d_{1}\):

\[\begin{eqnarray} \Phi_{1}d_{1}&=& \textrm{Sqrt}\Biggl\{ \sum\limits_{h=1}^{L} \dfrac{d_{1}^{2}}{2\sqrt{3}(b-a)}\Biggr\} \tag{40} \end{eqnarray}\]

And for the stages \(k\geq2\):

\[\begin{eqnarray} \Phi_{k}d_{k}&=&\min_{0\leq l_{k}\leq d_{k}}\Biggl\{ \textrm{Sqrt}\Biggl\{ \sum\limits_{h=1}^{L} \dfrac{l_{k}^{2}}{2\sqrt{3}(b-a)}\Biggr\}+ \Phi_{k-1}(d_{k}-l_{k})\Biggr\} \tag{41} \end{eqnarray}\]

Substituting the values of \(a\), \(b\), \(y_{0}\) and \(d\), the OSW (\(l_{h}^{*}\)) and the OSB (\(y_{h}^{*}=y_{h-1}^{*}-l_{h}^{*}\)) are obtained by executing the strata.data() function.

9.7.2 A Numerical Example for Uniform Distribution

A data following Uniform distribution of size \(N=5000\) was simulated to demonstrate the application of the stratifyR package on a Uniform population. The data for Uniform distribution is simulated with the parameters \(min = 2\) and \(max=15\). When fitted, the minimum and maximum values in the simulated data are \([y_{0}, y_{L}] = [2.006522, 14.99764]\), which implies that \(d=12.99112\).

To construct the OSB (\(h = 2\)) for the Uniformly-distributed data with a fixed total sample size of \(450\), we use the following codes:

set.seed(15669)
data <- runif(5000, min = 2, max = 15)
hist(data)

res <- strata.data(data, h = 2, n=450) # a 2-strata solution
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 450 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 5000</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [2.0065, 14.9976]  (d = 12.9911)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> unif</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter |  Estimate|
#> |:---------|---------:|
#> |min       |  2.006522|
#> |max       | 14.997645|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	8.50216	0.49	3.55	0.919	220	2437	0.09
2	14.9976	0.51	3.52	0.962	230	2563	0.09
Total		1.00	7.07	1.880	450	5000	0.09

Note that the above results indicate that the best-fit distribution is Triangular and not Uniform. This is because when stratifyR package assesses the data to ascertain the best-fit distribution, it is not able to calculate its AIC value for Uniform distribution. Hence, the distribution that gives the lowest AIC is taken to be the best-fit distribution. Since AIC for Uniform distribution is not available, the next best-fit distribution (Triangular) is be chosen. The results, you would find, are still quite accurate!

Similarly, we can apply the strata.distr() function for the Uniform distribution where the arguments can be assumed from past knowledge. To construct the OSB for \(h = 2\) for a hypothetical variable that follows Uniform distribution (with parameters \(min=3\) and \(max=15\)) with a fixed total sample size of \(450\) from a population of \(5000\), the following command is used:

# For a hypothetical uniform distribution, it does give a result
res <- strata.distr(h=2, initval=3, dist=12, distr = "unif",
                 params = c(min=3, max=15), n=450, N=5000)
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 450 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 5000</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [3.0000, 15.0000]  (d = 12.0000)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> unif</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter | Estimate|
#> |:---------|--------:|
#> |min       |        3|
#> |max       |       15|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	9	0.5	3	0.866	225	2500	0.09
2	15	0.5	3	0.866	225	2500	0.09
Total		1.0	6	1.732	450	5000	0.09

9.8 Stratification for a Survey Variable with Normal Distribution

If the study variable follows the Normal distribution on the interval \([y_{0}, y_{L}]\), it has the following two-parameter probability density function: \[\begin{equation} f(y;\mu, \sigma)=\dfrac{1}{\sigma\sqrt{2\pi}}\,exp\left\{-\dfrac{1}{2}\left(\dfrac{y-\mu}{\sigma}\right)^{2}\right\},\;\;\;\;\;\ -\infty < y < \infty \tag{42} \end{equation}\]

where \(\sigma > 0\) is a scale parameter and \(\mu\) is the location parameter.

The following definitions of error function are worth noting since they are needed to simplify the integrations used to derive the stratum weight, mean and variance due to normal distribution. \[\begin{eqnarray} erf(z) = \dfrac{2}{\sqrt{\pi}}\int_{0}^{z}exp\left\{-y^{2}\right\}\,dy \tag{43} \end{eqnarray}\] It can also be written as \[\begin{eqnarray}\label{erf} \dfrac{1}{\sqrt{2\pi}}\int_{0}^{z}exp\left\{-\dfrac{1}{2}y^{2}\right\}\,dy = \dfrac{1}{2}\,erf\left(\dfrac{z}{\sqrt{2}}\right) \tag{44} \end{eqnarray}\]

9.8.1 DP Solution for Normal Distribution

To solve the MPP formulated for Normal distribution (42), we apply the algorithm using the DP technique discussed in Section 4. The recurrence relations used to determine the OSB are given by:

For the first stage, \(k=1\), at \(l_{1}^{*}=d_{1}\):

\[\begin{eqnarray} \Phi_{1}d_{1}&=& \textrm{Sqrt}\Biggl\{ \dfrac{\sigma^{2}}{2\sqrt{2\pi}} \left[erf\left(\dfrac{d_{1}+y_{0}-\mu}{\sigma\sqrt{2}}\right) - erf\left(\dfrac{y_{0}-\mu}{\sigma\sqrt{2}}\right)\right]\nonumber \\[10pt] && \times\left[\left(\dfrac{y_{0}-\mu}{\sigma}\right)exp\left(-\left(\dfrac{y_{0}-\mu}{\sigma\sqrt{2}}\right)^{2}\right)- \left(\dfrac{d_{1}+y_{0}-\mu}{\sigma}\right)exp\left(-\left(\dfrac{d_{1}+y_{0}-\mu}{\sigma\sqrt{2}}\right)^{2}\right)\right]\nonumber\\[10pt] && + \dfrac{\sigma^{2}}{4}\left[erf\left(\dfrac{d_{1}+y_{0}-\mu}{\sigma\sqrt{2}}\right) - erf\left(\dfrac{y_{0}-\mu}{\sigma\sqrt{2}}\right)\right]^{2}\nonumber\\[10pt] && - \dfrac{\sigma^{2}}{2\pi}\left[exp\left(-\left(\dfrac{y_{0}-\mu}{\sigma\sqrt{2}}\right)^{2}\right)- exp\left(-\left(\dfrac{d_{1}+y_{0}-\mu}{\sigma\sqrt{2}}\right)^{2}\right)\right]^{2}\Biggr\} \tag{45} \end{eqnarray}\]

And for the stages \(k\geq2\):

\[\begin{eqnarray} \Phi_{k}d_{k}&=&\min_{0\leq l_{k}\leq d_{k}}\Biggl\{ \textrm{Sqrt}\Biggl\{ \dfrac{\sigma^{2}}{2\sqrt{2\pi}} \left[ erf\left(\dfrac{d_{k}+y_{0}-\mu}{\sigma\sqrt{2}}\right) - erf\left(\dfrac{(d_{k}-l_{k}+y_{0})-\mu}{\sigma\sqrt{2}}\right) \right]\nonumber \\[10pt] && \times\left[\left(\dfrac{(d_{k}-l_{k}+y_{0})-\mu}{\sigma}\right)exp\left(-\left(\dfrac{(d_{k}-l_{k}+y_{0})-\mu}{\sigma\sqrt{2}}\right)^{2}\right)\right. \nonumber\\[10pt] && - \left. \left(\dfrac{d_{k}+y_{0}-\mu}{\sigma}\right)exp\left(-\left(\dfrac{d_{k}+y_{0}-\mu}{\sigma\sqrt{2}}\right)^{2}\right)\right]\nonumber\\[10pt] && + \dfrac{\sigma^{2}}{4}\left[erf\left(\dfrac{d_{k}+y_{0}-\mu}{\sigma\sqrt{2}}\right) - erf\left(\dfrac{(d_{k}-l_{k}+y_{0})-\mu}{\sigma\sqrt{2}}\right)\right]^{2}\nonumber\\[10pt] && - \dfrac{\sigma^{2}}{2\pi}\left[exp\left(-\left(\dfrac{(d_{k}-l_{k}+y_{0})-\mu}{\sigma\sqrt{2}}\right)^{2}\right)-exp\left(-\left(\dfrac{d_{k}+y_{0}-\mu}{\sigma\sqrt{2}}\right)^{2}\right)\right]^{2}\Biggr\} + \Phi_{k-1}(d_{k}-l_{k})\Biggr\} \tag{46} \end{eqnarray}\]

Substituting the values of \(\mu\), \(\sigma\), \(y_{0}\) and \(d\), the OSW (\(l_{h}^{*}\)) and the OSB (\(y_{h}^{*}=y_{h-1}^{*}-l_{h}^{*}\)) are obtained by executing the strata.data() function.

9.8.2 A Numerical Example for Normal Distribution

A data following Normal distribution (herein called \(Norm\) data), of size \(N=5000\) was simulated to demonstrate the application of the stratifyR package on a Normal population. The data exhibits an Normal distribution with the parameters \(mean = 16.010776\) and \(sd=1.662357\). The minimum and maximum values in the simulated data are \([y_{0}, y_{L}] = [9.923816, 22.51267]\), which implies that \(d=10.62118\).

To construct the OSB for \(h = 2\) using the \(Norm\) data with a fixed total sample size of \(580\), the command below can be used:

set.seed(89821)
data <- rnorm(5000, mean = 16, sd = 1.65)
hist(data)

res <- strata.data(data, h = 2, n=500) #construct a 2-strata solution
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 500 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 5000</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [9.9238, 22.5127]  (d = 12.5889)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> norm</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter |  Estimate|
#> |:---------|---------:|
#> |mean      | 16.010776|
#> |sd        |  1.662357|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	16.0115	0.5	0.98	0.493	244	2489	0.1
2	22.5127	0.5	1.05	0.515	256	2511	0.1
Total		1.0	2.03	1.009	500	5000	0.1

Similarly, in order to find the OSB and other quantities, we can apply the strata.distr() function to a Normal population. Based on the assumption from past knowledge that the population follows Normal distribution with the given attributes such as the initial value, distance, parameters, etc., we can execute the following command to obtain the results:

set.seed(89821)
data <- rnorm(5000, mean = 16, sd = 1.65)
min(data); max(data); d=max(data)-min(data);d
#> [1] 9.923816
#> [1] 22.51267
#> [1] 12.58885
fit <- fitdist(data, distr="norm", method="mle")
fit
#> Fitting of the distribution ' norm ' by maximum likelihood 
#> Parameters:
#>       estimate Std. Error
#> mean 16.010776 0.02350928
#> sd    1.662357 0.01662355
res <- strata.distr(h=2, initval=9.923816, dist=12.58885, distr = "norm",
             params = c(mean=16.010776, sd=1.662357), n=500, N=5000)
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 500 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 5000</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [9.9238, 22.5127]  (d = 12.5889)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> norm</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter |  Estimate|
#> |:---------|---------:|
#> |mean      | 16.010776|
#> |sd        |  1.662357|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	16.0115	0.5	1	0.5	250	2501	0.1
2	22.5127	0.5	1	0.5	250	2499	0.1
Total		1.0	2	1.0	500	5000	0.1

9.9 Stratification for a Survey Variable with Log-Normal Distribution

If the study variable follows the Log-normal distribution on the interval \([y_{0}, y_{L}]\), it has the following two-parameter probability density function:

\[\begin{equation} f(y;\mu,\sigma)=\dfrac{1}{y\sigma\sqrt{2\pi}}\,exp\left\{-\dfrac{1}{2}\left(\dfrac{ln(y)-\mu}{\sigma}\right)^{2}\right\},\;\;\;\;\;\ y>0 \tag{47} \end{equation}\]

where \(\sigma > 0\) is a scale parameter and \(\mu\) is the location parameter.

9.9.1 DP Solution for Log-Normal Distribution

To solve the MPP formulated for Log-Normal distribution (47), we apply the algorithm using the DP technique discussed in Section 4. The recurrence relations used to determine the OSB are given by:

For the first stage, \(k=1\), at \(l_{1}^{*}=d_{1}\):

\[\begin{eqnarray} \Phi_{1}d_{1}&=& \textrm{Sqrt}\Biggl\{\frac{1}{4}exp\left(2\mu+2\sigma^{2}\right)\left[ erf\left(\frac{ln(d_{1}+y_{0})-\mu-2\sigma^{2}}{\sigma\sqrt{2}}\right)\right.\nonumber \\[10pt] && \left. - erf\left(\frac{ln(y_{0})-\mu-2\sigma^{2}}{\sigma\sqrt{2}}\right) \right] \left[ erf\left(\frac{ln(d_{1}+y_{0})-\mu}{\sigma\sqrt{2}}\right) - erf\left(\frac{ln(y_{0})-\mu}{\sigma\sqrt{2}}\right)\right]\nonumber \\[10pt] && -\frac{1}{4}exp\left(2\mu+\sigma^{2}\right)\left[erf\left(\frac{ln(d_{1}+y_{0})-\mu-\sigma^{2}}{\sigma\sqrt{2}}\right) - erf\left(\frac{ln(y_{0})-\mu-\sigma^{2}}{\sigma\sqrt{2}}\right) \right]^{2}\Biggr\} \tag{48} \end{eqnarray}\]

And for the stages \(k\geq2\):

\[\begin{eqnarray} \Phi_{k}d_{k}&=&\min_{0\leq l_{k}\leq d_{k}}\Biggl\{ \textrm{Sqrt}\Biggl\{ \frac{1}{4}exp\left(2\mu+2\sigma^{2}\right)\left[erf\left(\frac{ln(d_{k}+y_{0})-\mu-2\sigma^{2}}{\sigma\sqrt{2}}\right)\right.\nonumber \\[10pt] && - \left. erf\left(\frac{ln(d_{k}-l_{k}+y_{0})-\mu-2\sigma^{2}}{\sigma\sqrt{2}}\right) \right] \left[ erf\left(\frac{ln(d_{k}+y_{0})-\mu}{\sigma\sqrt{2}}\right)\right.\nonumber \\[10pt] && -\left. erf\left(\frac{ln(d_{k}-l_{k}+y_{0})-\mu}{\sigma\sqrt{2}}\right) \right]-\frac{1}{4}exp\left(2\mu+\sigma^{2}\right)\nonumber \\[10pt] &&\times \left[ erf\left(\frac{ln(d_{k}+y_{0})-\mu-\sigma^{2}}{\sigma\sqrt{2}}\right) - erf\left(\frac{ln(d_{k}-l_{k}+y_{0})-\mu-\sigma^{2}}{\sigma\sqrt{2}}\right) \right]^{2}\Biggr\} + \Phi_{k-1}(d_{k}-l_{k})\Biggr\}. \tag{49} \end{eqnarray}\]

Substituting the values of \(\mu\), \(\sigma\), \(y_{0}\) and \(d\), the OSW (\(l_{h}^{*}\)) and the OSB (\(y_{h}^{*}=y_{h-1}^{*}-l_{h}^{*}\)) are obtained by executing the strata.data() function.

9.9.2 A Numerical Example for Log-Normal Distribution

The ‘hies’ data of size \(N=3566\) is used to demonstrate the application of the stratifyR package on Log-normal population. The ‘hies’ data comes from the HIES survey conducted in Fiji in the year 2010. The data contains only two aspects of the survey, namely Income and Expenditure. In this example, the variable Expenditure is used since it exhibits a 2-parameter Log-normal distribution with the shape and scale parameters as \(meanlog = 9.2804934\) and \(sdlog = 0.6917842\) respectively. The minimum and maximum values are \([y_{0}, y_{L}] = [991.24, 136539.1]\), which implies that \(d=135547.8\).

To construct the OSB (\(h = 2\)) for the Iron data with a fixed total sample size of \(500\), we use the following codes:

data(hies)
Expenditure <- hies$Expenditure
head(Expenditure);length(Expenditure)
#> [1] 10880.58  3633.28  3882.84  4862.67  9237.19  9637.24
#> [1] 3566
hist(Expenditure)

min(Expenditure); max(Expenditure); d=max(Expenditure)-min(Expenditure);d
#> [1] 991.24
#> [1] 136539.1
#> [1] 135547.8
fit <- fitdist(Expenditure, distr="lnorm", method="mle")
fit
#> Fitting of the distribution ' lnorm ' by maximum likelihood 
#> Parameters:
#>          estimate  Std. Error
#> meanlog 9.2804934 0.011584571
#> sdlog   0.6917842 0.008191452
res <- strata.data(Expenditure, h = 2, n=500) 
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 500 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 3566</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [991.2400, 136539.0600]  (d = 135547.8200)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> lnorm</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter | Estimate|
#> |:---------|--------:|
#> |meanlog   | 9.280493|
#> |sdlog     | 0.691784|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	17206.7	0.77	14022946	2886.776	214	2749	0.08
2	136539	0.23	285033008	3868.016	286	817	0.35
Total		1.00	299055954	6754.792	500	3566	0.14

Similarly, in order to find the OSB and other quantities, we can apply the strata.distr() function to a Log-normal population. Based on the assumption from past knowledge that the population follows Log-normal distribution with the given attributes such as the initial value, distance, parameters, etc., we can execute the following command to obtain the results:

res <- strata.distr(h=2, initval=10, dist=188, distr = "lnorm",
             params = c(meanlog=3.23, sdlog=0.65), n=500, N=1588)
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 500 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 1588</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [10.0000, 198.0000]  (d = 188.0000)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> lnorm</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter | Estimate|
#> |:---------|--------:|
#> |meanlog   |     3.23|
#> |sdlog     |     0.65|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	39.6614	0.76	63.79	5.423	247	1200	0.21
2	198	0.24	516.96	5.535	253	388	0.65
Total		1.00	580.75	10.958	500	1588	0.31

9.10 Stratification for a Survey Variable with Cauchy Distribution

If the study variable follows the Cauchy distribution on the interval \([y_{0}, y_{L}]\), its two-parameter probability density function is given by:

\[\begin{eqnarray} f(y; \mu, \sigma)=\dfrac{1}{{\pi}\sigma\left[1+\left(\dfrac{y-\mu}{\sigma}\right)^2\right]} \;\;\;\;-\infty < y < +\infty \tag{50} \end{eqnarray}\]

where \(\mu\) is the location parameter which specifies the location of the peak of the distribution and \(\sigma\) is the scale parameter, which specifies the half-width at half maximum of the distribution.

9.10.1 DP Solution for Cauchy Distribution

To solve the MPP formulated for Cauchy distribution (50), we apply the algorithm using the DP technique discussed in Section 4. The recurrence relations used to determine the OSB are given by:

For the first stage, \(k=1\), at \(l_{1}^{*}=d_{1}\):

\[\begin{eqnarray} \Phi_{1}d_{1}&=& \textrm{Sqrt}\Biggl\{\dfrac{1}{\pi^{2}}\left[\tan^{-1}\left(\dfrac{d_{1}+y_{0}-\mu}{\sigma}\right)-\tan^{-1}\left(\dfrac{y_{0}-\mu}{\sigma}\right)\right]\nonumber\\[10pt] &&\times\left[\mu\sigma\ln\left(\left(d_{1}+y_{0}-\mu\right)^2+\sigma^2\right)+\left(\mu^2-\sigma^2\right)\tan^{-1}\left(\dfrac{d_{1}+y_{0}-\mu}{\sigma}\right)\right.\nonumber\\[10pt] &&+\sigma\left(d_{1}+y_{0}\right)-\mu\sigma\ln\left(\left(y_{0}-\mu\right)^2+\sigma^2\right)\nonumber\\[10pt] &&-\left.\left(\mu^2-\sigma^2\right)\tan^{-1}\left(\dfrac{y_{0}-\mu}{\sigma}\right)-\sigma y_{0}\right]\nonumber\\[10pt] &&-\left[\dfrac{1}{4\pi^{2}}\left[\sigma\ln\left(\left(d_{1}+y_{0}-\mu\right)^2+\sigma^2\right)+2\mu\tan^{-1}\left(\dfrac{d_{1}+y_{0}-\mu}{\sigma}\right)\right.\right.\nonumber\\[10pt] &&-\left.\left.\sigma\ln\left(\left(y_{0}-\mu\right)^2+\sigma^2\right)-2\mu\tan^{-1}\left(\dfrac{y_{0}-\mu}{\sigma}\right)\right]\right]^{2}\Biggr\} \tag{51} \end{eqnarray}\]

And for the stages \(k\geq2\):

\[\begin{eqnarray} \Phi_{k}d_{k}&=&\min_{0\leq l_{k}\leq d_{k}}\Biggl\{ \textrm{Sqrt}\Biggl\{ \dfrac{1}{\pi^{2}}\left[\tan^{-1}\left(\dfrac{d_{k}+y_{0}-\mu}{\sigma}\right)-\tan^{-1}\left(\dfrac{d_{k}-l_{k}+y_{0}-\mu}{\sigma}\right)\right]\nonumber\\[10pt] &&\times\left[\mu\sigma\ln\left(\left(d_{k}+y_{0}-\mu\right)^2+\sigma^2\right)+\left(\mu^2-\sigma^2\right)\tan^{-1}\left(\dfrac{d_{k}+y_{0}-\mu}{\sigma}\right)\right.\nonumber\\[10pt] &&+\sigma\left(d_{k}+y_{0}\right)-\mu\sigma\ln\left(\left(d_{k}-l_{k}+y_{0}-\mu\right)^2+\sigma^2\right)\nonumber\\[10pt] &&-\left.\left(\mu^2-\sigma^2\right)\tan^{-1}\left(\dfrac{d_{k}-l_{k}+y_{0}-\mu}{\sigma}\right)-\sigma y_{h-1}\right]\nonumber\\[10pt] &&-\left[\dfrac{1}{4\pi^{2}}\left[\sigma\ln\left(\left(d_{k}+y_{0}-\mu\right)^2+\sigma^2\right)+2\mu\tan^{-1}\left(\dfrac{d_{k}+y_{0}-\mu}{\sigma}\right)\right.\right.\nonumber\\[10pt] &&-\left.\left.\sigma\ln\left(\left(d_{k}-l_{k}+y_{0}-\mu\right)^2+\sigma^2\right)-2\mu\tan^{-1}\left(\dfrac{d_{k}-l_{k}+y_{0}-\mu}{\sigma}\right)\right]\right]^{2} \Biggr\} + \Phi_{k-1}(d_{k}-l_{k})\Biggr\}. \tag{52} \end{eqnarray}\]

Substituting the values of \(\mu\), \(\sigma\), \(y_{0}\) and \(d\), the OSW (\(l_{h}^{*}\)) and the OSB (\(y_{h}^{*}=y_{h-1}^{*}-l_{h}^{*}\)) are obtained by executing the strata.data() function.

9.10.2 A Numerical Example for Cauchy Distribution

The Boston data from the MASS package is used to demonstrate the application for Cauchy population. In this example, the variable ‘black’ is used since it exhibits a 2-parameter Cauchy distribution with the location and scale parameters as \(location = 393.864307\) and \(scale = 4.710457\) respectively. The minimum and maximum values are \([y_{0}, y_{L}] = [0.32, 396.9]\), which implies that \(d=396.58\).

data(Boston) #Housing Values in Suburbs of Boston
black = Boston$black
hist(black)

min(black); max(black); d=max(black)-min(black);d
#> [1] 0.32
#> [1] 396.9
#> [1] 396.58
fit <- fitdist(black, distr="cauchy", method="mle")
fit
#> Fitting of the distribution ' cauchy ' by maximum likelihood 
#> Parameters:
#>            estimate Std. Error
#> location 393.863075  0.3060843
#> scale      4.709346  0.3395534
res <- strata.data(black, h = 2, n=500)
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 500 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 506</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [0.3200, 396.9000]  (d = 396.5800)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> cauchy</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter |   Estimate|
#> |:---------|----------:|
#> |location  | 393.864327|
#> |scale     |   4.708422|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	364.295	0.2	18335.38	27.028	101	101	1.00
2	396.9	0.8	58.14	6.103	399	405	0.99
Total		1.0	18393.52	33.131	500	506	0.99

When the data include negative values, as they can for a Cauchy variable, the families that are defined only for positive values cannot be fitted. The package skips them quietly and chooses among the remaining families by AIC.

Similarly, in order to find the OSB and other quantities, we can apply the strata.distr() function to a Cauchy population. Based on the assumption from past knowledge that the population follows Cauchy distribution with the given attributes such as the initial value, distance, parameters, etc., we can execute the following command to obtain the results:

#for a cauchy distribution with initial value of x0=-1, d=2 and 
#location and scale parameters 0 and 1 respectively
res <- strata.distr(h=2, initval=-1, dist=2, distr = "cauchy",
             params = c(location=0, scale=1), n=500, N=5000)
#> The program is running, it'll take some time!
summary(res)
#> <div style='margin-bottom:0.5rem'><span style='color:#0b7285;font-weight:700'>Stratification summary</span></div>
#> <div><span style='color:#0b7285;font-weight:600'>Strata (h):</span> 2 &nbsp; <span style='color:#0b7285;font-weight:600'>n:</span> 500 &nbsp; <span style='color:#0b7285;font-weight:600'>N:</span> 5000</div>
#> <div><span style='color:#0b7285;font-weight:600'>Data range:</span> [-1.0000, 1.0000]  (d = 2.0000)</div>
#> <div><span style='color:#0b7285;font-weight:600'>Best-fit distribution:</span> cauchy</div>
#> 
#> 
#> 
#> Table: Parameter estimate(s)
#> 
#> |Parameter | Estimate|
#> |:---------|--------:|
#> |location  |        0|
#> |scale     |        1|

Per-stratum results and totals

Strata	OSB	Wh	Vh	WhSh	nh	Nh	fh
1	0	0.5	0.08	0.07	250	2500	0.1
2	1	0.5	0.08	0.07	250	2500	0.1
Total		1.0	0.16	0.14	500	5000	0.1

10 Summary of Method Comparison Across All Supported Distributions

To compare the four solvers on the examples of Section 9, we ran strata.data() with \(h = 3\) on the data used in each distribution section. Table 1 shows the Neyman cost \(\sum W_h S_h\) (WhShTot) for each solver, and the gap of the DP cost above the exact optimum. The values were computed with stratifyR 2.0-2; the code to reproduce them is given after the table.

Table 1: Neyman cost (WhShTot) of the four solvers with h = 3 on the examples of Section 9. ‘AIC family’ is the distribution chosen by AIC, which the DP optimises. ‘DP gap’ is the DP cost above the exact optimum. Lower WhShTot means a more efficient design.

Section	Data	N	AIC family	DP	COBYLA	GLOBAL	EXACT	DP gap (%)
Pareto	Simulated	5000	weibull	0.93820	0.92710	0.92710	0.92710	1.19
Triangular	math final marks	487	weibull	7.11720	7.09160	7.09160	7.09160	0.36
Right-triangular	Simulated	1000	lnorm	0.49590	0.48850	0.48850	0.48850	1.50
Weibull	anaemia Iron	724	norm	2.46380	2.44340	2.44340	2.44340	0.84
Gamma	anaemia Folate	724	norm	2.93570	2.89850	2.89850	2.89850	1.28
Exponential	Simulated	5000	exp	0.27360	0.27341	0.27341	0.27341	0.07
Uniform	Simulated	5000	unif	1.24710	1.24700	1.24700	1.24700	0.01
Normal	Simulated	5000	norm	0.71698	0.71685	0.71685	0.71685	0.02
Log-normal	hies Expenditure	3566	lnorm	4678.50000	4598.90000	4598.90000	4598.90000	1.73
Cauchy	Boston black	506	cauchy	20.67200	14.92100	14.92100	14.92100	38.50

Three points stand out.

1. COBYLA, GLOBAL and EXACT agree everywhere. On all ten examples, the two fast solvers found the exact optimum. They took between 0.2 and 2.5 seconds, and EXACT took between 0.03 and 1.3 seconds.

2. DP is within 2% of the optimum when the fitted model is reasonable. The gap is below 0.1% for the exponential, uniform and normal data, where the chosen family matches the data well, and between 0.4% and 1.7% for the skewed examples. The DP took 20 to 50 seconds per run at \(h = 3\).

3. The DP depends on the fitted model. For the Boston black variable, the AIC chooses the Cauchy, but the data are strongly skewed to the left and the fitted Cauchy describes them poorly. The DP boundaries, which are optimal for that model, cost 38.5% more than the exact optimum for the actual data. Fitting a distribution is what makes strata.distr() possible when no data are available, but when the data are at hand, it is worth checking the DP result against method = "exact".

Note that the family chosen by AIC can differ from the section name. For example, the Iron and Folate data (with the minimum replaced as in Section 9) are fitted best by the normal. The family chosen, and therefore the DP result, can also differ slightly between platforms when a maximum-likelihood fit is close to failing. The COBYLA, GLOBAL and EXACT results do not depend on the fit.

In summary, "exact" gives the best design for a given population, "cobyla" gives the same answer in almost all cases and is the fastest, and "dp" remains the default because it optimises a smooth fitted model and is the method of the original package.

The table can be reproduced with the following code (it takes a few minutes because of the DP runs):

library(MASS); library(triangle); library(actuar)
pops <- list()
set.seed(8235411); pops$Pareto <- actuar::rpareto(5000, shape=5, scale=8)
data(math); pops$Triangular <- math$final_marks
set.seed(12546); pops$`Right-triangular` <- rtriangle(n=1000, a=2, b=8, c=2)
data(anaemia)
x <- anaemia$Iron;   x[x == min(x)] <- -0.001; pops$Weibull <- x
x <- anaemia$Folate; x[x == min(x)] <- -0.001; pops$Gamma   <- x
set.seed(28951); pops$Exponential <- rexp(5000, rate=1.36)
set.seed(15669); pops$Uniform <- runif(5000, min=2, max=15)
set.seed(89821); pops$Normal <- rnorm(5000, mean=16, sd=1.65)
data(hies); pops$`Log-normal` <- hies$Expenditure
pops$Cauchy <- Boston$black

tab <- t(sapply(pops, function(d) {
  n <- min(500, floor(length(d) / 2))
  sapply(c("dp", "cobyla", "global", "exact"),
         function(m) strata.data(d, h=3, n=n, method=m)$WhShTot)
}))
tab

References

Abramowitz, M, and I A Stegun. 1972. “Handbook of Mathematical Functions Dover Publications.” New York, 361.

Bellman, R E. 1957. Dynamic Programming. Princeton University Press.

Bühler, W, and T Deutler. 1975. “Optimal stratification and grouping by dynamic programming.” Metrika 22 (1): 161–75. https://doi.org/10.1007/BF01899725.

Cochran, W G. 1961. “Comparison of methods for determining stratum boundaries.” Bulletin of the International Statistical Institute 38 (2): 345–58.

Dalenius, T. 1950. “The Problem of Optimum Stratification.” Scandinavian Actuarial Journal 1950 (3-4): 203–13.

Dalenius, T. 1957. Sampling in Sweden: contributions to the methods and theories of sample survey practice. Almqvist; Wiksell.

Gablonsky, J M, and C T Kelley. 2001. “A Locally-Biased Form of the DIRECT Algorithm.” Journal of Global Optimization 21 (1): 27–37.

Gunning, P, and J M Horgan. 2004. “A new algorithm for the construction of stratum boundaries in skewed populations.” Survey Methodology 30 (2): 159–66.

Khan, E A, M G M Khan, and M J Ahsan. 2002. “Optimum stratification: a mathematical programming approach.” Calcutta Statistical Association Bulletin 52: 323–33.

Khan, M G M, N Ahmad, and Sabiha Khan. 2009. “Determining the Optimum Stratum Boundaries Using Mathematical Programming.” Journal of Mathematical Modelling and Algorithms 8 (4): 1–15. https://doi.org/10.1007/s10852-009-9115-3.

Khan, M G M, E A Khan, and M J Ahsan. 2003. “Theory & Methods: An Optimal Multivariate Stratified Sampling Design Using Dynamic Programming.” Australian & New Zealand Journal of Statistics 45 (1): 107–13.

Khan, M G M, N Nand, and N Ahmad. 2008. “Determining the Optimum Strata Boundary Points Using Dynamic Programming.” Survey Methodology 34 (2): 205–14.

Khan, M G M, D K Rao, A H Ansari, and M J Ahsan. 2015. “Determining Optimum Strata Boundaries and Sample Sizes for Skewed Population with Log-Normal Distribution.” Communications in Statistics-Simulation and Computation 44 (5): 1364–87.

Khan, M G M, K G Reddy, and D K Rao. 2015. “Designing Stratified Sampling in Economic and Business Surveys.” Journal of Applied Statistics 42 (10): 2080–99.

Khan, M G M, N Sehar, and M J Ahsan. 2005. “Optimum stratification for exponential study variable under Neyman allocation.” Journal of the Indian Society of Agricultural Statistics 59 (2): 146–50.

Lohr, S. 2009. Sampling: Design and Analysis. Nelson Education.

Nand, N, and M G M Khan. 2009. “Optimum Stratification for Cauchy and Power Type Study Variable.” Journal of Applied Statistical Science 16 (4): 453.

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, 558–625.

Powell, M J D. 1994. A Direct Search Optimization Method That Models the Objective and Constraint Functions by Linear Interpolation. (Dordrecht), 51–67.

Reddy, K G, M G M Khan, and D K Rao. 2014. “Computing optimal stratum boundaries for multivariate surveys.” Computer Science and Engineering (APWC on CSE), 2014 Asia-Pacific World Congress on, 1–7.

Reddy, Karuna Garan, Mohammad GM Khan, and Sabiha Khan. 2018. “Optimum Strata Boundaries and Sample Sizes in Health Surveys Using Auxiliary Variables.” PloS One 13 (4): e0194787.

Reddy, Karuna G, and MGM Khan. 2019. “Optimal Stratification in Stratified Designs Using Weibull-Distributed Auxiliary Information.” Communications in Statistics-Theory and Methods 48 (12): 3136–52.

Särdnal, C E, B Swensson, and J H Wretman. 1992. Model Assisted Survey Sampling.