Procedural steps
The optimization of the sampling design starts by making the sampling frame available, defining the target estimates of the survey and establishing the precision constraints on them. Then, a choice has to be made with regard to what stratification variables to choose among those available in the frame, on the basis of an analysis of the correlation existing between the two sets of variables (stratification and target). When the chosen stratification variables are both categorical and continuous, in order to make them homogeneuous the continuous ones should be categorized (by using for instance a clustering k-means algorithm). Then, the optimization step making use of the ‘atomic’ method can be executed. If, conversely, the stratification variables are all of the continuous type the optimization step can be directly executed by making use of the method ‘continuous.’ It is also possible to perform both kinds of optimization, compare the results and choose the more convenient.
Prior to the optimization based on the use of the genetic algorithm, it is advisable to run a different quick optimization task based on the use of the k-means algorithm, with a twofold purpose:
- give a hint on a suitable number of final strata;
- get an initial ‘good’ solution to be given as a ‘suggestion’ to the genetic algorithm in order to speed its convergence to the final solution.
In the optimization step it is possilble to indicate the set of sampling units that must be selected in any case (‘take-all’ strata). After the optimization, it is possible to evaluate the quality of the solution by simulating the selection of a high number of samples from the frame, and calculating sampling variance and bias for all the target variables. It is also possible to ‘adjust’ the sample size of the optimized solution on the basis of the available budget: if a higher size is allowable, sampling rates in strata are increased proportionally until the new total sample size is reached; the opposite is done in case we are obliged to reduce the sample size.
Finally, we proceed with the selection of the sample.
In the following, each step will be illustrated starting from a real sampling frame, the one that comes with the R package sampling (the dataframe swissmunicipalities).
Preparation of inputs required by the optimization step
Frame
For sake of simplicity, let us consider a subset of the swissmunicipalities dataset:
library(SamplingStrata)
data(swissmunicipalities)
swissmun <- swissmunicipalities[swissmunicipalities$REG < 4,
c("REG","COM","Nom","HApoly",
"Surfacesbois","Surfacescult",
"Airbat","POPTOT")]
head(swissmun)
#> REG COM Nom HApoly Surfacesbois Surfacescult Airbat POPTOT
#> 2 1 6621 Geneve 1593 67 31 773 177964
#> 3 3 2701 Basel 2391 97 93 1023 166558
#> 4 2 351 Bern 5162 1726 1041 1070 128634
#> 5 1 5586 Lausanne 4136 1635 714 856 124914
#> 9 2 371 Biel (BE) 2123 976 196 463 48655
#> 10 2 942 Thun 2158 425 694 523 40377In order to limit the processing time we have selected only the first 3 regions and only the variables of interest for our example. Each row in this dataset contains information on a swiss municipalities, identified by COM and Nom, and belonging to one of three selected regions (REG).
Suppose we want to plan a sampling survey where the target estimates Ys are the totals of woods area (Surfacesbois) and buildings area (Airbat) for each of the 3 regions (domains of interest). Suppose also that in each municipalities are always updated the values of total area (HApoly) and total population (POPTOT). From the correlation matrix:
cor(swissmun[,c(4:8)])
#> HApoly Surfacesbois Surfacescult Airbat POPTOT
#> HApoly 1.00000000 0.76920101 0.3398758 0.2602289 0.09940795
#> Surfacesbois 0.76920101 1.00000000 0.4492188 0.2968918 0.09881328
#> Surfacescult 0.33987579 0.44921881 1.0000000 0.3230496 0.11451750
#> Airbat 0.26022890 0.29689183 0.3230496 1.0000000 0.86896631
#> POPTOT 0.09940795 0.09881328 0.1145175 0.8689663 1.00000000we see that the correlations between Surfacesbois and HApoly, from one side, and between Airbat and POPTOT from the other side, are high (respectively 0.77 and 0.87), so we decide that both HApoly and POPTOT play the role of stratification variables Xs in our frame.
In a first moment we decide to treat the stratification variables as categorical, so we have to categorize them. A suitable way to do so, is to apply a k-means clustering method (see Hartigan and Wong 1979) by using the function var.bin:
swissmun$HApoly.cat <- var.bin(swissmun$HApoly,15)
table(swissmun$HApoly.cat)
#>
#> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
#> 311 345 281 214 158 126 111 84 66 32 29 32 15 11 8
swissmun$POPTOT.cat <- var.bin(swissmun$POPTOT,15)
table(swissmun$POPTOT.cat)
#>
#> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
#> 452 388 266 193 124 109 79 59 42 33 35 28 5 6 4We can now define the frame dataframe in the format required by SamplingStrata. Function buildFrameDF permits to organize data in a suitable mode for next processing:
frame1 <- buildFrameDF(df = swissmun,
id = "COM",
X = c("POPTOT.cat","HApoly.cat"),
Y = c("Airbat","Surfacesbois"),
domainvalue = "REG")
head(frame1)
#> id X1 X2 Y1 Y2 domainvalue
#> 1 6621 15 8 773 67 1
#> 2 2701 15 9 1023 97 3
#> 3 351 15 12 1070 1726 2
#> 4 5586 15 11 856 1635 1
#> 5 371 14 9 463 976 2
#> 6 942 14 9 523 425 2Strata
This dataframe is not explicitly required, as it is automatically produced from the frame dataframe by the optimStrata function. Norwithstanding, it is worth while to produce it in order to analyse the initial stratification of the frame, and what could be the associated sample size without optimization. The function buildStrataDF is the one to produce the strata dataframe:
strata1 <- buildStrataDF(frame1, progress=F)
#>
#> Computations are being done on population data
#>
#> Number of strata: 350
#> ... of which with only one unit: 130
head(strata1)
#> STRATO N M1 M2 S1 S2 COST CENS DOM1 X1 X2
#> 1*1 1*1 64 5.359375 34.07812 2.846968 24.54161 1 0 1 1 1
#> 1*10 1*10 3 17.333333 858.66667 13.299958 264.73928 1 0 1 1 10
#> 1*11 1*11 3 9.666667 842.00000 1.699673 349.82567 1 0 1 1 11
#> 1*12 1*12 1 9.000000 1100.00000 0.000000 0.00000 1 0 1 1 12
#> 1*13 1*13 5 12.200000 1085.40000 3.187475 537.57924 1 0 1 1 13
#> 1*15 1*15 1 10.000000 387.00000 0.000000 0.00000 1 0 1 1 15Each row in this dataframe report information related to a given stratum (obtained by cross-classifying each unit with the values of the X variables) regarding:
- the identifier of the stratum (named ‘strato’), concatenation of the values of the X variables;
- the values of the m auxiliary variables (named from X1 to Xm) corresponding to those in the frame;
- the total number of units in the population (named ‘N’);
- a flag (named ‘cens’) indicating if the stratum is to be censused (=1) or sampled (=0);
- a variable indicating the cost of interviewing per unit in the stratum (named ‘cost’);
- for each target variable y, its mean and standard deviation, named respectively ‘Mi’ and ‘Si’);
- the value of the domain of interest to which the stratum belongs (‘DOM1’).
Precision constraints
The errors dataframe contains the accuracy constraints that are set on target estimates. This means to define a maximum coefficient of variation for each target variable and for each domain value. Each row of this frame is related to accuracy constraints in a particular subdomain of interest, identified by the domainvalue value. In our case, we have chosen to define the following constraints:
ndom <- length(unique(swissmun$REG))
cv <- as.data.frame(list(DOM=rep("DOM1",ndom),
CV1=rep(0.10,ndom),
CV2=rep(0.10,ndom),
domainvalue=c(1:ndom) ))
cv
#> DOM CV1 CV2 domainvalue
#> 1 DOM1 0.1 0.1 1
#> 2 DOM1 0.1 0.1 2
#> 3 DOM1 0.1 0.1 3This example reports accuracy constraints (maximum CV allowable equal to 10%) on variables Y1 and Y2 that are the same for all the 3 different subdomains (Swiss regions) of domain level DOM1. Of course we can differentiate the precision constraints region by region. It is important to underline that the values of ‘domainvalue’ are the same than those in the frame dataframe, and correspond to the values of variable ‘DOM1’ in the strata dataframe.
We check that the dataframes that we have defined so far are correct:
checkInput(errors = checkInput(errors = cv,
strata = strata1,
sampframe = frame1))
#>
#> Input data have been checked and are compliant with requirementsFor instance, this function controls that the number of auxiliary variables is the same in the frame and in the strata dataframes; that the number of target variables indicated in the frame dataframe is the same than the number of means and standard deviations in the strata dataframe, and the same than the number of coefficient of variations indicated in the errors dataframe.
So far so good. Now we want to determine the total sample size, and related allocation, under the given strata, using the function bethel (Bethel (1989)):
allocation <- bethel(strata1,cv[1,])
sum(allocation)
#> [1] 570This is the total sample size (570) required to satisfy precision constraints under the current stratification, before the optimization.
Optimization
The function optimStrata is the one performing the optimization step. Actually, this function is a ‘wrapper’ that calls three different optimization functions:
- optimizeStrata (method = atomic, if stratification variables are categorical, or reduced to categorical);
- optimizeStrata2 (method = continuous, if stratification variables are continuous);
- optimizeStrataSpatial (method = spatial, if stratification variables are continuous and there is spatial correlation among units in a geo-coded sampling frame).
For continuity reasons, these functions are still available to be used standalone, and in some situations it may be useful to use them by a direct call (see related help for details).
Here we report the most important parameters related to the three methods (for the others see the help):
| Parameter | Description |
|---|---|
| framesamp | The name of the dataframe containing the information related to the |
| sampling frame. | |
| framecens | The name of the dataframe containing the units to be selected in any case. |
| It has same structure than “framesamp” dataframe. | |
| nStrata | The number of final optimized strata to be obtained. |
| model | In case the Y variables are not directly observed, but are estimated by |
| means of other explicative variables, in order to compute the anticipated | |
| variance, information on models are given by a dataframe “model” with as | |
| many rows as the target variables. Default is NULL. | |
| errors | The dataframe containing the precision levels expressed in |
| terms of maximum allowable coefficient of variation on a given estimate. | |
| minnumstr | Minimum number of units that must be allocated in each stratum |
| (default is 2). | |
| iter | Maximum number of iterations (= generations) of the genetic algorithm |
| (default is 50). | |
| pops | The dimension of each generations in terms of individuals (default is 20). |
| suggestions | Optional parameter for genetic algorithm that indicates a suggested solution |
| to be introduced in the initial population. The most convenient is the | |
| one found by the function “KmeanSolution” (default is NULL). |
Method ‘atomic’
As a first run we execute the optimization step using the method atomic (required as the stratification variables are of the categorical type).
set.seed(1234)
solution1 <- optimStrata(method = "atomic",
errors = cv,
nStrata = rep(10,ndom),
framesamp = frame1,
iter = 50,
pops = 10)
# Number of strata: 350
# ... of which with only one unit: 130
# *** Starting parallel optimization for 3 domains using 3 cores
# |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=07s
#
# *** Sample size : 306
# *** Number of strata : 22
# ---------------------------
The execution of optimStrata produces the solution of 3 different optimization problems, one for each domain. The graphs illustrate the convergence of the solution to the final one starting from the initial one (i.e. the one related to the atomic strata). Along the x-axis are reported the executed iterations, from 1 to the maximum, while on the y-axis are reported the size of the sample required to satisfy precision constraints. The upper (red) line represent the average sample size in each iteration, while the lower (black) line represents the best solution found until the i-th iteration.
We can calculate (analytically) the expected CVs by executing the function:
expected_CV(solution1$aggr_strata)
#> cv(Y1) cv(Y2)
#> DOM1 0.0990294 0.0998385
#> DOM2 0.1000000 0.1000000
#> DOM3 0.0989679 0.0952105The obtained total size of the sample required to satisfy precision constraint is much lower than the one obtained with the simple application of the Bethel algorithm to the initial atomic stratification, but maybe not yet satisfactory.
In order to explore other solutions we may want that each unit in the sampling frame be considered as an atomic stratum, and let to the optimization step to aggregate them on the basis of the values of the Y variables. In any case, as we have to indicate at least one X variable, we can use to this purpose a simple progressive number:
swissmun$progr <- c(1:nrow(swissmun))
frame2 <- buildFrameDF(df = swissmun,
id = "COM",
X = "progr",
Y = c("Airbat","Surfacesbois"),
domainvalue = "REG")
head(frame2)
#> id X1 Y1 Y2 domainvalue
#> 1 6621 1 773 67 1
#> 2 2701 2 1023 97 3
#> 3 351 3 1070 1726 2
#> 4 5586 4 856 1635 1
#> 5 371 5 463 976 2
#> 6 942 6 523 425 2We can use this approach because the number of units in the frame is small: it would not be possible to consider each unit as a stratum in case of real population registers or even with business registers.
Even so, the processing of the 1,823 strata may be slow.
In order to speed up the convergence towards the optimal solution, an initial one can be given as a “suggestion” to ‘optimizeStrata’ function. The function KmeansSolution produces this initial solution by clustering atomic strata considering the values of the means of all the target variables Y. For any given number of clusters, the correspondent aggregation of atomic strata is considered as input to the function bethel. The number of clusters for which the value of the sample size necessary to fulfil precision constraints is the minimum one, is retained as the optimal one. Also, the optimal number of clusters is determined inside each domain. It is possible to indicate a maximum number of aggregated strata to be obtained by using the parameter maxcluster:
strata2 <- buildStrataDF(frame2, progress=F)
initial_solution2 <- KmeansSolution(strata = strata2,
errors = cv,
maxclusters = 10) 
The overall solution is obtained by concatenating optimal clusters obtained in domains. The result is a dataframe with two columns: the first indicates the clusters, the second the domains. On the basis of these, we can calculate the most convenient number of final strata for each domain:
nstrata2 <- tapply(initial_solution2$suggestions,
initial_solution2$domainvalue,
FUN=function(x) length(unique(x)))
nstrata2
#> 1 2 3
#> 10 8 10and we can supply both the initial_solution and the number of strata nstrata as inputs to the optimization step:
set.seed(1234)
solution2 <- optimStrata(method = "atomic",
errors = cv,
framesamp = frame2,
iter = 50,
pops = 10,
nStrata = nstrata2,
suggestions = initial_solution2)
# Number of strata: 1823
# ... of which with only one unit: 1823
# *** Starting parallel optimization for 3 domains using 3 cores
# |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=08s
#
# *** Sample size : 164
# *** Number of strata : 28
# ---------------------------
Notice that the obtained solution in this run in terms of sample size is significantly better than in the previous one.
outstrata2 <- solution2$aggr_strata
expected_CV(outstrata2)
#> cv(Y1) cv(Y2)
#> DOM1 0.0994007 0.0301118
#> DOM2 0.0995359 0.0236793
#> DOM3 0.0985103 0.0206549Method ‘continuous’
The last thing to do is to test also the continuous method.
First, we have to redefine the frame dataframe in this way:
frame3 <- buildFrameDF(df = swissmun,
id = "COM",
X = c("POPTOT","HApoly"),
Y = c("Airbat","Surfacesbois"),
domainvalue = "REG")
head(frame3)
#> id X1 X2 Y1 Y2 domainvalue
#> 1 6621 177964 1593 773 67 1
#> 2 2701 166558 2391 1023 97 3
#> 3 351 128634 5162 1070 1726 2
#> 4 5586 124914 4136 856 1635 1
#> 5 371 48655 2123 463 976 2
#> 6 942 40377 2158 523 425 2Also in this case we want to produce an initial solution:
set.seed(1234)
init_sol3 <- KmeansSolution2(frame=frame3,
errors=cv,
maxclusters = 10) 
nstrata3 <- tapply(init_sol3$suggestions,
init_sol3$domainvalue,
FUN=function(x) length(unique(x)))
nstrata3
# 1 2 3
# 10 10 9 initial_solution3 <- prepareSuggestion(init_sol3,frame3,nstrata3)Note that this time we call a different function (KmeansSolution2) that requires, instead of the strata dataframe, directly the frame dataframe. Moreover, we need an intermediate step to prepare the suggestion for the optimization, that is the execution of the function prepareSuggestion.
We are now able to perform the optimization step with the continuous method:
set.seed(1234)
solution3 <- optimStrata(method = "continuous",
errors = cv,
framesamp = frame3,
iter = 50,
pops = 10,
nStrata = nstrata3,
suggestions = initial_solution3)
# *** Starting parallel optimization for 3 domains using 3 cores
# |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=09s
#
# *** Sample size : 93
# *** Number of strata : 26
# ---------------------------
This solution requires a total sample size that is by far the best among those we produced, so we decide to select this one.
