Optimum strata boundaries and sample sizes in health surveys using auxiliary variables

Using convenient stratification criteria such as geographical regions or other natural conditions like age, gender, etc., is not beneficial in order to maximize the precision of the estimates of variables of interest. Thus, one has to look for an efficient stratification design to divide the whole population into homogeneous strata that achieves higher precision in the estimation. In this paper, a procedure for determining Optimum Stratum Boundaries (OSB) and Optimum Sample Sizes (OSS) for each stratum of a variable of interest in health surveys is developed. The determination of OSB and OSS based on the study variable is not feasible in practice since the study variable is not available prior to the survey. Since many variables in health surveys are generally skewed, the proposed technique considers the readily-available auxiliary variables to determine the OSB and OSS. This stratification problem is formulated into a Mathematical Programming Problem (MPP) that seeks minimization of the variance of the estimated population parameter under Neyman allocation. It is then solved for the OSB by using a dynamic programming (DP) technique. A numerical example with a real data set of a population, aiming to estimate the Haemoglobin content in women in a national Iron Deficiency Anaemia survey, is presented to illustrate the procedure developed in this paper. Upon comparisons with other methods available in literature, results reveal that the proposed approach yields a substantial gain in efficiency over the other methods. A simulation study also reveals similar results.


Introduction
Stratified random sampling is an important sampling technique utilized in estimating the prevalence of diseases such as diabetes, anaemia, obesity hypertension, and smoking. In stratified sampling, the sampling-frame is divided into a number (say, L) of non-overlapping groups or strata in such a way that the strata constructed are internally homogeneous with respect to the variable (or main variable) under study, because that maximizes the precision of the estimator of the parameter of interest concerning the study variable, e.g. its mean [1]. An advantage of stratified sampling design is that when a stratum is homogeneous, the measurements of a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 the study variable (y) vary little from each other and the precise estimate of y can be obtained from a small sample in that stratum. Thus, combining these estimates from all L strata, the design produces a gain in the precision of estimate of the variable in the whole population [1]. However, in most practical situations, especially in health surveys, it is difficult to construct such optimum strata, and hence, more often the health surveyors stratify the population in most convenient manners such as the use of geographical regions (e.g. North, Central, South, etc.), administrative regions (e.g. provinces, districts, etc.) or other natural criteria (e.g. urbanrural, sex, age etc.). Moreover, the stratification by convenience manner is not always a reasonable criterion as the strata so obtained may not be internally homogeneous with respect to a variable of interest. Thus, one has to look for the Optimum Stratum Boundaries (OSB) that maximizes the precision of the estimators.
The problem of determining OSB for a variable, when its frequency distribution is known, is well-known in the sampling literature. The basic consideration involved in determining OSB is that the strata should be as internally homogenous as possible, that is, in order to achieve maximum precision, the stratum variances should be as small as possible [1,2]. When a single variable is of interest and the stratification is made based on this study variable, then an ideal situation is that the distribution of the study variable is known and the OSB can be determined by cutting the range of its distribution at suitable points. This problem of determining the OSB, when both the estimation and stratification variables are the same, was first discussed by Dalenius [3]. He presented a set of minimal equations which are usually difficult to solve for OSB because of their implicit nature. Hence, subsequently the attempts for determining approximately optimum stratum boundaries have been made by several authors [4][5][6][7][8][9].
Many authors have also attempted to determine the global OSB. For example, Unnithan [10] proposed an iterative method that requires a suitable initial solution. For a skewed population where a certainty stratum (some specific units are included in the sample where extremely large units are isolated so that they do not influence sampling variability) is necessary. Lavallée [11] proposed an algorithm to construct stratum boundaries for a power allocated (applying an exponential value q, where 0 < q < 1, to the stratum population value under Neyman Allocation to allow for a sufficient spread of the sample allocation) stratified sample. Later on, Hidiroglou [12] presented a more general form of the algorithm. After reviewing Lavallée and Hidiroglou's algorithm, a modified algorithm that incorporated the different relationships between the stratification and study variables was proposed [13,14].
There are several other algorithms available in the literature, for example, Niemiro [15] proposed a random search method and the simplex method [16] was used to present a new method of stratification [17]. Later on, Kozak [18] presented a modified random search algorithm. Gunning [19] proposed an alternative method to approximate stratification based on a geometric progression. This approach was compared with three other methods [8,11,20] which confirmed that the geometric progression method is more efficient [21]. The usefulness of Gunning and Horgan's geometric progression method was studied and it revealed that the geometric progression approach is less efficient than Lavallée and Hidiroglou's algorithm [22,23].
Another kind of stratification method that has been proposed in the literature is due to Khan et. al. [24][25][26][27][28][29][30]. When the distributions of the study variables were known, they formulated the problems of determining OSB as optimization problems, which are solved by developing computational techniques Dynamic Programming (DP). The DP technique was first proposed by Bühler & Deutler [31], which was also used for determining the OSB which would divide the population domain of two stratification variables into distinct subsets such that the precision of the variables of interest is maximized [11,32].
Numerous research have also been undertaken whereby auxiliary variable(s), which can be historical data, are used to improve the precision of the estimates of study variable y. When the frequency distribution of the auxiliary variable, x, is known, several approximation methods of determining OSB using the auxiliary variable(s) have been suggested and discussed by many authors [1,9,[33][34][35][36][37][38][39][40][41][42][43].
In this paper, a procedure for determining the OSB and sample size for each stratum of a variable of interest in health surveys is developed. The determination of the OSB and sample sizes, based directly on the survey variable (y), is not feasible in practice since the it is unavailable prior to conducting the survey. Hence, optimum stratification is made based on multiple auxiliary variables (x 1 , x 2 , . . ., x p ) that are readily available in health surveys. It shall be assumed that the population values of the study variable y are available as realizations of a stochastic background variable or at least can be realized as proxy values of y from previous or other recent surveys and y holds a regression model with the auxiliary variable(s) [2,14,30,[44][45][46]. Moreover, often y is highly correlated with x such that the regression of y upon x has homoscedastic errors. In situations like this, stratification can be achieved using the auxiliary variable (s). The application of the proposed methodology will be demonstrated with empirical investigations using real and simulated datasets. This proposed research deals with the problem of stratification for a study variable using the many auxiliary variables that are found in a multivariate survey. In health surveys, these auxiliary variables normally characterize positively skewed distributions that are families of the Gaussian distribution such as Weibull, Gamma, Log-normal, etc. Thus, this research investigates if the proposed parametric-based mathematical programming approach for determining the OSB yields a gain in efficiency over other methods that are well-known in literature. This research also tries to find out if the proposed method works for skewed distributions such as the Weibull or Gamma when both linear and nonlinear regression models are used in the MPP formulation of the stratification problem.
The problem of determining OSB is redefined as the problem of determining Optimum Strata Widths (OSW) and is formulated as a Mathematical Programming Problem (MPP) that seeks minimization of the variance of the estimated population parameter. Since the formulated MPP can be viewed as a multistage decision problem, it is solved using a DP technique. These OSB are then used to compute the sample size for each stratum under Neyman allocation. A numerical example with a real data set of skewed population, where the auxiliary variables follow Weibull distributions, is presented to illustrate the proposed procedure. The results are compared with the Dalenius & Hodges' cum ffiffi ffi f p method [20], Gunning & Horgan's geometric method [19] and Lavallée & Hidiroglou's method [11].

The general formulation of the problem of OSB as an MPP
Let the population be stratified into a fixed L strata based on p auxiliary variables, x 1 , x 2 , . . ., x p , and the estimation of the mean of 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 " y h ; ðh ¼ 1; 2; . . . ; LÞ, then an unbiased stratified sample mean, " y st , is given by where W h = N h /N is the proportion of the population contained in the h th stratum for the study variable y, where N is the total number of units in the population and is assumed to be known while N h is the total unknown number of units in each stratum. Then the variance of " y st is given by The finite population correction factors in (2) could be ignored [8,9,20,36]. Thus, under the Neyman allocation [47], that is, (2) is given by where σ hy is the stratum standard deviation of y in h th stratum; h = 1, 2, . . ., L and n is the preassigned total sample size. Consider that the study variable has the regression model of the form: where λ(x 1 , x 2 , . . ., x p ) is a linear (or nonlinear) function of x i ;(i = 1, 2, . . ., p) and is an error term such that E(|x 1 , x 2 , . . ., x p ) = 0 and V(|x 1 , The parameters in λ are assumed to be known from a recent survey.
Assuming that λ and are uncorrelated [4], it follows that where s 2 hlðx 1 ;x 2 ;...;x p Þ denotes the variance of λ(x 1 , x 2 , . . ., x p ) in the h th stratum and s 2 h is the variance of in the h th stratum. Eq (6) assumes homoscedasticity, i.e., homogeneity of the variance of over the distribution of the predictors x i (i = 1, 2, . . ., p), given the stratum h.
Let f(x i ) be the estimated frequency functions of the auxiliary variables, x i (i = 1, 2, . . ., p), that are used for the stratification of the main variable. If the population mean of the study variable y is to be estimated over a range (a, b) under the allocation (3), then the problem of determining the strata boundaries of y is to cut up the range, (a, b) at (L − 1) intermediate points a = y 0 y 1 y 2 , . . ., y L−1 y L = b such that (4) is minimum. Since the study variable is not available at the design stage, the range (a, b) could either be the compromise range derived from all the auxiliary variables or an estimated range that best explains the study variable, possibly chosen from previous surveys.
For a fixed sample size n, minimizing the expression of the right hand side of (4) is equivalent to minimizing P L h¼1 W h s hy . Thus, from (6), the following is minimized: If f(x i ) are known and integrable frequency functions of the auxiliary variables, then for the given λ(x 1 , x 2 , . . ., x p ), the first term inside the square root function in (7) can be expressed as the functions of the boundary points (y h−1 , y h ) by finding the stratum weight W hx i , mean μ hx i and variance s 2 hx i of ith auxiliary variable x i using the following expressions: where i = 1, 2, . . ., p. The quantities computed by Eqs (8)-(10) may be different for different auxiliary variables since it depends on their best-fit frequency distributions, for example, Weibull, Gamma, or any other skewed distribution. If two or more auxiliary variables are characterized by the same distribution, the quantities in (8)-(10) may still be different because their estimated parameters would certainly be different.
Note that the second term in (7) are also obtained as a function of boundary points using the frequency function of the regression error. Thus, the objective function (7) could be expressed as a function of boundary points (y h−1 , y h ) only: Then, the problem of determination of OSB can be expressed as the following optimization problem to determine y 1 , y 2 , . . ., y L .
We further define l h = y h − y h−1 ;h = 1, 2, . . ., L, where l h ! 0 denotes the range or width of the h th stratum. From this, the range of the distribution of y, d = b − a, can be expressed as a function of the stratum width.
The h th stratification point y h ; h = 1, 2, . . ., L is then expressed as Adding (12) as a constraint, the problem (11) can be treated as an equivalent problem of determining optimum strata widths (OSW), l 1 , l 2 , . . ., l L , and is expressed as: Note that if y 0 is known, the first term, ϕ 1 (l 1 , y 0 ), in the objective function of the MPP (14) is a function of l 1 alone. Once l 1 is known, the second term ϕ 2 (l 2 , y 1 ) will become a function of l 2 alone and so on. Due to the special nature of functions, the MPP (14) may be treated as a function of l h alone and is expressed as: In real-world situations, the study variable is unknown at the design stage, hence, readilyavailable auxiliary variables can be used to create OSB. The proposed technique carries out optimization through the MPP (15) on the compromise range (d) derived from all auxiliary variables. The technique also assumes that the parameters of the regression model in (15) are known from a recent survey. The best-fit distributions of the auxiliary variables, x i , are used in the formulation of MPP (15).

The solution procedure using dynamic programming technique
The problem (15) is a multistage decision problem in which the objective function and the constraint are separable functions of l h , which allows us to use a DP technique [28]. Dynamic programming determines the optimum solution of a multi-variable problem by decomposing it into stages, each stage compromising a single variable subproblem. A DP model is basically a recursive equation based on Bellman's principle of optimality [48]. This recursive equation links the different stages of the problem in a manner which guarantees that each stage's optimal feasible solution is also optimal and feasible for the entire problem [49].
Consider the following subproblem of (15) for first k(<L) strata: and where d k < d is the total width available for division into k strata or the state value at stage k. Note that d k = d for k = L and the transformation functions are given by denote the minimum value of the objective function of (16), that is, and l h ! 0; h ¼ 1; 2; . . . ; k and 1 k L

" #
With the above definition of F k (d k ), the MPP (15) is equivalent to finding F L (d) recursively by finding F k (d k ) for k = 1, 2, . . ., L and 0 d k d. We can write: and l h ! 0; h ¼ 1; 2; . . . ; k " # For a fixed value of l k ; 0 l k d k , and l h ! 0; h ¼ 1; 2; . . . k À 1; 1 k L # Using the Bellman's principle of optimality, we write a forward recursive equation of the DP technique for k ! 2 as: For the first stage, that is, for k = 1: where l Ã 1 ¼ d 1 is the optimum width of the first stratum. The relations (17) and (18) are solved in a forward manner first for k = 1, 2, . . ., L to determine the optimum subproblem objective and then solved in a backward manner second to determine the OSB.
The application of the above solution procedure is summarized in Appendix A in order to determine the OSB for MPP (15).

Determination of optimum sample size
When OSB (y h , y h−1 ) are determined as discussed in Sections 2-3, the optimum sample size n h ; h = 1, 2, . . ., L that minimizes the variance of the estimate can easily be computed.
If the study variable holds the regression model (5) with the auxiliary variables across the strata, using (2) and (7), the sample size n h are obtained for a fixed total sample of size n under Neyman allocation [47] for h = 1, 2, . . ., L and given as follows: where W h , s 2 hlðx 1 ;x 2 ;...;x p Þ and s 2 h are derived in terms of the optimum boundary points (y h , y h−1 ). It is also worth noting that the OSB (y h , y h−1 ) through the MPP (15) are so obtained that n h must satisfy the restrictions: where N h = NW h . The restriction 1 n h is added to the formulation so that the h th stratum must form with at least a unit and the restriction n h N h is added to avoid the over sampling.

Determination of optimum number of strata
This is one of the first issues that need to be considered in an optimal stratification design, however, it can be dependent on the OSB and the allocation of sample units among the strata.
The goal of stratification is to make all strata as homogenous as possible, which implies that the more the number of strata, the more the homogeneity within a stratum. This results in a reduction in the total variance of " y st , that is, Var ð" y st Þ. However, an increase in the number of strata may involve extra cost and resources in planning and drawing the samples.
Problem of determining optimum number of strata was first discussed by Dalenius [3] who postulated that uniformly distributed variable, Var ð" y st Þ is inversely proportional to L 2 . Later, Cochran [50] investigated the effect of the number of strata on Var ð" y st Þ for some skewed distributed populations with Neyman allocation. He confirmed that this relationship holds for skewed distribution and the rate of reduction in Var ð" y st Þ is independent of skewness of the population. The results indicated that only a little reduction in variance is to be expected beyond L = 6 unless the correlation between the auxiliary information and the survey population is greater than 0.95.
To apply the above idea to the current situation of the optimal number of strata, assume that fpc is negligible and consider the distribution of the data to be approximately uniform, as done by Cochran [1]. Then the range of the distribution of values of y [a, b] is d = b − a, and hence the variance of the distribution is S 2 = d 2 /12. The variance of the sample mean for a simple random sample of size n can therefore be calculated as: Thus, if a simple case of creating L strata of equal size is considered, stratum variance would then be calculated as S 2 = d 2 /12L 2 . It follows from W h = 1/L and Eq (4), This reveals that variance of the sample mean is inversely proportional to the square of the number of strata, L. This however, does not consider the relationship between the auxiliary variables and the study population. It can be extended by considering a linear relationship given in Eq (5). As suggested by Cochran [50], in this case, using (5), it can be shown that This again shows that the variance is inversely related to the square of the number of strata. Applying (23), one can empirically study the effect of increasing the number of strata. To complete this analysis, a cost function that shows the relation of cost with L, for planning and executing a survey, is required. However, whatever the form of cost function, [1] showed that the increase in L beyond 6 will seldom be profitable. Thus, if the extra cost involved in planning and executing the survey, which is incurred due to an increase in the number of strata is not of much importance, a reasonable approach to determining the optimum number of strata may be discussed as follows: Compute Vð" y st Þ for L = 1, 2, ‥, k, where k is a possible value of the candidate L. Now Vð" y st Þ decreases as L increases and Vð" y st Þ is minimum when L = k. Therefore, a surveyor may choose the optimum number of strata at the point where an increase in L is not useful as it gives only a small decrease in Vð" y st Þ. The approach is illustrated in Fig 1, which is a hypothetical plot of Vð" y st Þ against L. One can choose the desired number of strata as the point at which the "elbow" in the curve becomes apparent. Clearly, this requires judgment on the part of the surveyors.

Construction of OSB with Weibull auxiliary variables
The Weibull distribution is a two or three-parameter family of continuous probability distributions. Because of its versatility in the fitting of a variety of distributions, it is one of the most widely used distributions in applied statistics, especially in survival analysis, mortality or failure analysis, reliability, engineering to model manufacturing and delivery times, in extreme value theory and weather forecasting. Due to its moderately skewed profile, it also characterizes well a wide range of health data, including health monitoring data, epidemiological data such as episode durations of depression and gene expressions data [51][52][53][54].
If all the auxiliary variables, x i ;i = 1, 2, . . ., p, approximately follow Weibull distribution on the interval [x i,0 , x i, L ], its three-parameter probability density function with a state space x i ! 0 is given by: where r i > 0 is the shape parameter, θ i > 0 is the scale parameter and γ i is the location parameter of the distribution of i th auxiliary variable. The shape parameter gives the Weibull distribution its flexibility. By changing the value of the shape parameter, the distribution can model a wide variety of data that follows the Exponential distribution, the Rayleigh distribution, the Normal distribution or even the approximate Log-normal distribution.

Scaling the auxiliary variables
While stratifying the study variable based on multiple auxiliary variables, the raw data in the form of different auxiliary variables are generally of different scales (eg., kg, mg, dollar, etc.). The values of one variable may be less or more spread out than other variables. With the auxiliary variables exhibiting different distributions, the range of data, minimum and maximum values for these auxiliary variables will certainly be different from each other. Hence, this may affect the convergence of the MPP (15) and hence its ability to determine the OSB accurately. A way to encounter this problem is to standardize each variable by subtracting its mean and dividing by its standard deviation.
Another method, which this paper uses, is a simple scaling procedure whereby every variable is divided by its maximum value. While maintaining the original distribution of the auxiliary variables, this scaling procedure results in the auxiliary variables getting closer to each other, which in turn, helps in reducing the overall range or the search space of the optimal solution. One must note that the solution procedure of dynamic programming technique is generally advisable and feasible for small sets of units (N 20) [55], hence, scaling is a necessary means for faster convergence of an optimal solution. The MPP (15), when solved, provides the OSB of the scaled study variable and the OSB for the original study variable can be obtained by the usual re-scaling procedure.

Estimating the regression model
To illustrate the estimation of the regression model in formulating the problem of determining OSB as an MPP for a population with more than one auxiliary variable, we use a health survey data on Anameia, which was obtained from the 2004 Fiji National Nutritional Survey conducted by the National Food and Nutrition Centre (Fiji) and funded by AusAID, UNICEF and Government of Fiji. The data included a micronutrient survey where blood samples were drawn from women of childbearing age and measurements were made to record levels of Haemoglobin, Iron and Folate amongst many other variables. Whilst only tabulations are publicly available from http://ghdx.healthdata.org/record/fiji-national-nutrition-survey-2004, data used for the purpose of applying the proposed method is accessible from http://repository.usp. ac.fj/id/eprint/10439 where the main aim is to estimate the Haemoglobin content in Fijian women. The whole data was fully anonymized before making them accessible. The data cannot be de-anonymized because there is no public datasets available to cross-reference.
The data has the following three characteristics for each woman: 1. Level of Haemoglobin 2. Level of Iron

Level of Folate
Suppose that a survey on Iron Deficiency Anaemia is to be conducted in a country, where Haemoglobin (y) is the variable of interest and is to be stratified. Then, the levels of Iron and Folate collected in this study may be the reasonable choice for the auxiliary variables, x 1 and x 2 . In this example, Haemoglobin is available to us but in reality the main variable might not be available prior to the survey. Thus, Haemoglobin will be used purely as an example for numerical illustrations and comparison purposes.
To estimate the Haemoglobin content (y) in women, a multiple regression model (given by Eq (4)) was fitted using scaled data for the survey mentioned above. It was observed that the data significantly fitted a linear regression model with Iron and Folate levels (p < 0.001)-the estimated parameters for these two predictors were also highly significant (p < 0.001).
The coefficient of determination R 2 = SSR/SST, with an Adjusted R-squared value of 12.54% was found to be one of the highest for the linear model when compared with the model summary of all the other non-linear models available in standard statistical packages. Thus, this model fits the data best and gives us no reason to consider an alternative model. There is a small positive linear relationship between Haemoglobin and Iron (r = 0.350, p < 2.2e − 16), and Haemoglobin and Folate (r = 0.161, p < 1.31e − 05). Therefore, the Haemoglobin content (y), Iron level (x 1 ) and Folate level (x 2 ) are fairly assumed to follow a linear regression model given in (5): This idea can be applied in the ideal situation where the main variable is not available. The beta weights of the regression model, initial and final values could be taken as guestimates from prior surveys.

Estimating the distribution of the auxiliary variables
To determine the distributions of the auxiliary variables, f(x 1 ) and f(x 2 ), relative frequency histograms for Iron and Folate are constructed. The two histograms presented in Figs 2 and 3 reveal that the distributions of both auxiliary variables are right-skewed and match 3P Weibull distribution with different parameters.
Using the Kolmogorov-Smirnov test for each of the two variables, the maximum difference (D) between the observed distribution and the Weibull distribution is found to be non-significant (all p-values >0.05). This also supports the fact that all variables follow 3P Weibull distributions, where parameters are obtained by the maximum likelihood estimate (MLE) method. Optimum strata boundaries and sample sizes

Formulation of the MPP with Weibull distribution
Considering that y has a linear regression on x i ; (i = 1, 2, . . ., p). Then, from (5), the function λ (x 1 , x 2 , . . ., x p ) is of the form: Assume that the model in (26) holds for all strata. Then, Let all the auxiliary variables, x i , follow 3P Weibull distribution (i.e., x i * W(r i , θ i ), γ i ) with density function given by (24). By using (8)-(10), the quantities W hx i , μ hx i , and s 2 hx i can be obtained as a function of boundary points (y h−1 , y h ). Using the substitution of y h = y h−1 +l h , they are presented as follows: μ hλ can be expressed as: Optimum strata boundaries and sample sizes Let Γ(r, x) and Q(r, s) denote the upper incomplete gamma function and the regularized incomplete gamma function, respectively, given by Qðr; xÞ ¼ 1 GðrÞ Then, using Eqs (28)- (31), μ hx i can be simplified to be Similarly, the quantity s 2 hx i is reduced to where W hx i and and m 2 hx i are given by Eqs (28) and (32) respectively. Since the auxiliary variables follow Weibull distributions, W h and s 2 hl in the first term of (7) are given by (28) and (33) respectively. Thus, for the i th auxiliary variable, Using (34), the formulated MPP given in (15) could be generalised and expressed as the following MPP in order to determine the OSB for the main variable: ; Subject to where d in Eq (35) is the estimated or hypothetical range of the main study variable, β i are the regression coefficients, θ i and r i are parameters of the 3P Weibull distributions of i th auxiliary variable, Γ(Á) is the upper incomplete gamma function and Q(Á) is the upper regularized incomplete gamma function. Whereas, the term W 2 h s 2 h can be computed when the distribution of is known. For the current model, since this error term is normally distributed, the distribution is given by: Then, following from (8)- (10), W h and σ h are obtained as: where erf ðy h Þ À erf ðy hÀ 1 Þ ¼ 2 ffiffi p p R y h y hÀ 1 expðÀ u 2 Þ du and h = 1, 2, . . ., L.

Numerical illustrations
In this section, numerical results are presented to illustrate the application of the proposed technique to a real and a simulated population. The OSB for the main variable are obtained and presented together with the values of the objective function ð0 h ðl h Þ ¼ P L h¼1 W h s h Þ for L = 2, 3, . . ., 6 for different regression models.

Real data
The real data, as discussed earlier in Section 6.2, has Haemoglobin as the study variable while Iron and Folate are auxiliary variables that follow Weibull distributions with their estimated parameters. Haemoglobin is being used here purely for comparison purposes, in reality, the main variable is not available. Using the recursive Eqs (17) and (18), the MPP (35) with d = 10.9 (range of main variable) is solved by executing a C++ computer program developed to implement the proposed DP technique. R codes were also developed for computing the quantities such as the initial value (x 0 ) of the distribution, regression coefficients (β i ), Weibull parameters (α, β, γ), range (d) of the distribution, etc. required for determining the OSB using the C++ program. Users can easily stratify a population by executing the C++ program for the given value of L, x 0 , d, n, etc. in an open source IDE such as DEV C++. The C++ program and R codes can be made available on request from the authors.
The results for the OSB (y h ) along with optimum sample sizes (n h ) and the values of the objective function ð P L h¼1 W h s h Þ are presented in Table 1 for the following regression models: Model 1 :

Simulated data
A skewed population with two auxiliary variables (x 1 and x 2 ) and the study variable (y), each of size N = 5000, were randomly generated using the R software. This data had a relatively weak linear relationship between y and x 1 (r = 0.014, p = 0.34), and a weak linear relationship as well between y and x 2 (r = 0.023, p = 0.11). The simulated data was different from the real data in the sense that it had a very low predictive power in its regression models (Adj. R 2 = 0.03%). The ANOVA results from multiple linear regression also indicated a non-statistically significant model fit (p = 0.175). For the simulated data, the OSB (y h ) along with optimum sample sizes (n h ) and ∑W h σ h values are presented in Table 2 for the following regression models: Various other investigations related to OSB, sample size and the performance of the proposed technique, in both real and simulated data, are carried out and discussed in the following section and subsections.

Results and discussion
Primarily, this paper involves the usage of multiple auxiliary variables in determining the OSB for the study variable. Investigations into the performance of the proposed method are also carried out to investigate some of the very pertinent issues such as:  5. Sensitivity of the proposed method with linear regression against nonlinear regression; 6. Consistency of the results obtained for real data with a simulated data set.
Thus, in the following subsections, comparative results are presented for the three models that are to create OSB for the main variable in real and simulated data. These are done to ascertain the effects of using a single auxiliary variable and multiple auxiliary variables in terms of the changes observed in the OSB, sample sizes and the ∑W h σ h values achieved for L = 2, . . ., 6. Together with results on the performance of the proposed method against other methods, results for 3P Gamma distribution and nonlinear regression are also presented.

Use of single and multiple auxiliary variables
Tables 1 and 2 present the OSB, sample sizes and ∑W h σ h values for real and simulated data respectively. For real data in Model 2, which uses Folate, ∑W h σ h values are the lowest in the three models while for simulated data, it is lowest in Model 4 which uses variable x 1 . The ∑W h σ h values using the other models (i.e., models 1 & 3 in real data and 5 & 6 in simulated data) are close to each other. It is seen from the results that the ∑W h σ h values of the main variables in both real and simulated data appear to be declining exponentially as L increases in all the models. It must also be noted that in Table 2, the OSB and OSS are equivalent in all three models. This is due to the fact that the simulated data set is quite large and it results in a very precise-fitting of the distribution, which leads to equivalent OSB in all three models.
The findings in both data are similar in the sense that a single auxiliary variable model performs either better or worse than the model with multiple auxiliary variables. In real data, Model 2 performs better than Models 1 and 3 while in simulated data, Model 4 performs better than Models 5 and 6. This may be due to the fact that both Model 2 in real data and Model 4 in simulated data have a much weaker correlation with the dependent variable (see Tables 3  and 4). The results in Tables 3 and 4 presents key statistics such as the correlations, measure of regression error (RSE) and goodness of fit (AIC) for all the models in both data. It appears that the model with the auxiliary variable(s) that has the lowest correlation and Adjusted R 2 and the highest RSE or AIC performs the best for the proposed method. Thus, the proposed method of stratification works best with uncorrelated auxiliary variable(s).

Comparison with other available methods
For the purpose of the comparison of the performance of the proposed method, the following univariate methods available in the literature are considered: 1. Cum ffiffi ffi f p method [20].
3. Lavallée and Hidiroglou (Kozak) method [11,18] The stratification package developed by [56] in the R statistical software is used to determine the OSB and sample sizes for the main study variable, Haemoglobin. The OSB are then used to compute the variance of the estimated mean (i.e., the values of the objective function or ∑W h σ h ) in each of the six models so that a comparative analysis could be carried out between the established methods and the proposed method. Note that comparisons are only possible here since the main variable is available to us in this example. The three methods above need the main variable to work out the OSB, however, the proposed method can work on auxiliary variable(s) to compute OSB for the main variables, with a few assumptions on the main variable.
The results, based on Models 1, 2 and 3 for real data, are given in Table 5, which presents the ∑W h σ h values of the estimate for Cum ffiffi ffi f p method, Geometric method, Lavallée and Hidiroglou's method and the proposed method with a fixed total sample of size n = 500, for L = 2, 3, . . ., 6. The efficiencies of the proposed DP method over the other three methods are also presented in the table.
Upon examination of these results, it is noted that when a single auxiliary variable (Model 1) is used to determine OSB, the proposed method performs considerably well over the three methods and the efficiency of these OSB increases by about 2% to 50% for L = 1, 2, . . ., 6. Model 2 also produces much more efficient OSB over other methods and the efficiency increases from about 1% to 49%, which is quite similar to Model 1. Model 3 also increases the efficiencies from about 2% to 50%, being almost exactly similar to Model 1. Thus, with the use of auxiliary variables, either single or both, the proposed method increases the precision of the estimate compared to other univariate methods. Table 6 provides the OSB and sample sizes using the other methods which can be compared with the results of the proposed method presented in Table 1.
For simulated data, Table 7 presents the ∑W h σ h values for the three methods along with the proposed method for the three different models (Models 4-6) together with the efficiencies of the proposed method over the others. The results generally support the similar findings obtained for real data. Compared to all other methods, the proposed method increases the precision ranging from about 11% to 63% in Model 4 and 21% to 133% in Model 6. For Model 5, the proposed method increases the precision ranging from about 26% to 127% against Cum ffiffi ffi f p method and 25% to 135% against L-H (Kozak) method. It doesn't perform so well against Geometric method. Table 8 provides the OSB and sample sizes using the other methods which can be compared with the results of the proposed method presented in Table 2.
When considering Weibull distribution cases, the sample allocations under the proposed method (which uses Neyman allocation given by (19)) are given in Tables 1 and 2 for real and simulated data respectively. In the method, the overall size of strata (N h ) as well as variability (S 2 h ) of the auxiliary variable(s) affects the stratum sample sizes (n h ), i.e., n h / W h S h . It is noticeable that for both real and simulated examples, the stratum samples sizes given by the proposed method is a bit different from the sample sizes given by other methods presented in Optimum strata boundaries and sample sizes Optimum strata boundaries and sample sizes Tables 6 and 8. This is because of the differences seen in the OSB, and hence the W h , between the methods. To substantiate the results, the method of bootstrap re-sampling is used to investigate the behaviour of the findings made earlier on the real dataset. A large number (n = 10,000) of independent re-samples are drawn with replacement from the population data. The re-samples are of the same size as the Anaemia data (N = 724), creating many variants of the original data. Since there are three variables in the Anaemia population, bootstrap re-sampling is done on individuals, which means three variables are randomly generated for each population. From the large number of bootstrap re-samples, results for only 5 randomly selected samples are presented for the sake of brevity. We consider all three models given by equations in (39).
For all five bootstrap samples, Tables 9-13 present the OSB, OSS (n h ) and variances ð P L h¼1 W h s h Þ for all three models are calculated using the proposed method. It is again observed that Model 2 has the lowest variance and this means that it is the best model to use out of the three. To further investigate why Model 2 is the best, Table 14 is drawn up. It is found out that results are consistent in all five bootstrap samples. Model 2 performs the best because it has a low correlation with the main variable together with a high RSE, a very low adjusted R 2 and the highest AIC amongst the three models. Thus, whether it is a single or multiple auxiliary variables (ie., all models studied herein) used in the formulation of the problem of stratification, the gains in efficiency of the proposed method over other established methods are substantial. These are given by Tables 15-19 where we see that the variances given by the proposed method are lower than the other methods. Hence, with bootstrap re-sampling procedure, it is seen that we obtain consistent findings to what was seen in the original Anaemia data. Optimum strata boundaries and sample sizes Optimum strata boundaries and sample sizes Optimum strata boundaries and sample sizes Optimum strata boundaries and sample sizes Optimum strata boundaries and sample sizes Optimum strata boundaries and sample sizes

Number of strata
To study the relationship between the number of strata and the ∑W h σ h value, an investigation is carried out for the real and simulated data using the six models. The ∑W h σ h are calculated using the proposed method and the results are presented for L = 2, 3, . . ., 20. These are presented in Figs 4 and 5 where the curves appear to be on top of each other and all of them decrease exponentially. After L = 7, where the "elbow" is found, the rate of decrease in the ∑W h σ h values from there onwards is not as big as what is seen from L = 2 to 7. For argument's sake, one might even be comfortable with L = 6 as the appropriate number of strata. The finding supports the investigation carried out by Cochran [1] that the number of strata to be Optimum strata boundaries and sample sizes constructed beyond six is not much useful in terms of the relative gain in efficiency or the reduction of ∑W h σ h . All six models in real and simulated data are very similar when it comes to the relative gain in efficiencies and one can easily pick out L = 7 where the "elbow" appears, indicating that the percentage reduction thereafter is not worth investing in for a sample survey since additional costs are involved with increase in the number of strata. Increasing the number of strata to more than 7 may not be a good trade-off for a little gain of precision in the estimates.

Using skewed distributions other than Weibull
The distribution of the auxiliary variable can vary depending on how well the data fits a particular skewed distribution based on the best MLE of its parameters. Weibull is selected in this paper due it's versatility in fitting skewed distributions, especially for health data. To probe into the performance of Weibull distribution against any other skewed distribution, both auxiliary variables in the real and simulated data are fitted with a 3P Gamma distributions because of its moderately skewed profile as well. Three different linear regression models are again used for the comparison of results. The associated MPP is formulated and solved using the DP technique.
The OSB, sample sizes and ∑W h σ h values are presented in Table 20 for the real data while Optimum strata boundaries and sample sizes reveal that fitting the data with Weibull distribution yields a much more efficient set of OSB compared to fitting the data with Gamma distribution. This holds true for both single or multiple auxiliary variables. Results for the simulated data in Tables 2 and 21 also reveal similar findings. This may be due to the fact that Weibull was a better fit than Gamma for the two auxiliary variables,

Linear versus nonlinear regression
As shown in (5), the proposed method can incorporate linear as well as nonlinear regression for construction of OSB. In the preceding sections, it has been discussed that linear regression performs well in real as well as simulated data. To investigate the sensitiveness of linear regression over nonlinear regression, a simple case of quadratic regression is fitted in this section. Consider that the study variables are to be stratified using a single auxiliary variable (e.g., Iron in real & and x 2 in simulated data). Then, λ(x) in (5) reduces to: Simulated : The ANOVA results for this quadratic regression reveals that the model is statistically significant (p-value < 0.001) for both real and simulated data. Using the procedures discussed in Sections 3-7, the OSB and sample sizes are determined. Table 22 presents the results along with the ∑W h σ h values for real and simulated data respectively. The results reveal that for both data, the ∑W h σ h values from linear regression (Model 3 from Table 1 and Model 6 from Table 2) are lower than non-linear regression model which means that linear regression performs better than the nonlinear regression.
To investigate this further, Table 23 presents some key statistical measures such as measure of regression error (RSE) and goodness of fit (AIC) with regards to how the model under nonlinear regression performs against the models under linear regression for both real and simulated data. The measures for linear regression are presented in Tables 3 and 4. They reveal that the results are consistent with the findings earlier in the paper-that the model with the lowest Adjusted R 2 and the highest RSE or AIC performs the best. Thus, linear regression model performs better than the nonlinear regression model. Optimum strata boundaries and sample sizes

Conclusion
Stratified random sampling is an efficient and widely used sampling technique in health surveys to estimate the prevalence of diseases and many other parameters. Often, the surveyors encounter two major difficulties prior to drawing the samples and these are: (i) constructing the optimum strata within which the units are as homogeneous as possible and (ii) determining the optimum sample size to be drawn from each stratum, so that the precision of the estimates of parameters of the study or target variables are maximized. In this paper, a parametric-based method is proposed to address these two problems, which can be used to estimate parameters with more precision. The optimum stratification based on the study variable is not feasible in practice since it is unknown prior to conducting the survey. Thus, the proposed technique uses auxiliary information in designing the sampling plan. This paper investigates how the usage of one or more auxiliary variables influence the OSB and hence the effect on the efficiency of the stratum boundaries by fitting a distribution of Weibull family that characterize many health variables. It also investigates the sensitivity of the OSB and the performance of the proposed method by fitting with other skewed distributions such as Gamma. Together with investigating the optimum number of strata, the proposed method also sees the sensitiveness of linear and nonlinear regression modelling techniques in implementing the proposed method.
The problem of finding the OSB is formulated as an MPP that seeks minimization of the variance of the estimated population parameter and solved using a DP technique. The solution procedure is implemented through a C++ computer program and an R script to facilitate the computation of the OSB through the C++ program. Both materials can be made available on request from the authors. After obtaining the OSB, they are then used to compute the optimum sample size for each stratum using Neyman allocation. Numerical examples using a real data set and a simulated data set are presented to illustrate the application, the sensitivity and the usefulness of the proposed technique. This paper also presents the results from cum ffiffi ffi f p method [20], geometric method [19] and the generalized Lavallée and Hidiroglou's method [11,18] for a comparative analysis. It can be concluded that in the construction of strata for health populations, usage of both single or multiple auxiliary variables leads to substantial gains in the precision of the estimates over other available methods. It was also established that using uncorrelated auxiliary variable (s) to determine OSB for the main variable leads to much more efficient results. It was also found out that when another skewed distribution such as Gamma is used to characterize the distribution of the auxiliary variables, it performed well but not quite as accurate as Weibull. Hence, the best-fit distribution should always be chosen for more accurate calculation of OSB. It was also found out that when linear regression was used in formulating the problem of stratification, it performed better than nonlinear regression. This simply depends on the data and one must always choose the best regression technique to represent the relationship between the variables. 4. For k = 2, express the state variable as d k−1 = d k − l k . 5. Set F k (d k ) = 0 if l k > d k , where 0 d k d. 6. Calculate F k (d k ), the minimum value of RHS of (17) for l k ;0 l k d k .
7. Record F k (d k ) and l k .
8. For k ! 3, . . ., L, go to Step 4. 9. At k = L, F 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 F 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 F 1 (d 1 ).