Row–Column Sampling Design Using Auxiliary Ranking Variables

In this paper, we present a new class of row–column sampling design when a sampling region is presented by non-overlapping quadrats of an rp-by-p row–column grid. Quadrats are rank-ordered on a relevant auxiliary variable, first in the rp rows of the grid and then in the so formed p ranked columns. This ranking enhances the precision of the sampling using ideas from standard rank set sampling methods. The sampling design falls in the traditional sampling framework where the sample selection probabilities are independent of the variables of interest. Selection probabilities are governed by a spatial design and the auxiliary information to facilitate better sampling coverage over the sampling region. The paper constructs unbiased estimators for population mean, total and provides variance estimates for them. It also shows when the new sampling design performs better than its competitors. The new sampling design is applied to a forest and two agricultural field samples to estimate the population means. Supplementary materials accompanying this paper appear online.


INTRODUCTION
In agricultural, ecological or environmental sampling, population units at a geographical location are usually arranged in a number of quadrats that naturally or by design fall into a grid structure. In this grid structure, a subset of grid cells is selected under a plan produced on the basis of well-specified sampling methodology. The sampling methodology signifies the underlying probability structure (or sampling variation) of samples generated from the population with certain selection probabilities. These selection probabilities make it possible to draw inferences about overall population units in a finite population setting.
The choice of sampling methodology depends on the objective of the study, the nature of the underlying population, prior knowledge of the sampling region and spatial features of the grid structure. If the sampling region has some sort of spatial pattern, a sampling design that provides a good spatial coverage over the region may provide more information than a sampling plan providing a poor coverage (Bhave and Sadhwani 2021). A spatial pattern in a region may be a function of environmental factors, or chemical, physical and microbiological patterns of the soil map. Often, prior information is available for soil structure from previous experiments, yield maps and soil analysis of the sampling region. On a small scale, variation can also arise from a design difference within agricultural equipment-for example, variation caused by the tine geometry of a seeder. Environmental factors can be assessed by ground elevation, wind pattern and vegetation density of the sampling region or surrounding areas. This paper uses all available prior knowledge on grid structure to construct a sampling design to produce a better spatial pattern coverage for a population under study.
The prior knowledge of the sampling region could be subjective, such as in the form of visual inspection of the sampling region, or could be in the form of an auxiliary variable highly correlated with the variable of interest. The auxiliary variable may be obtained by flying a drone equipped with image processing software over the sampling region. The images can be processed to obtain a rough predictor of the variable of interest on each grid cell over the sampling region. Readers are referred to Kasprzak (2021) for further information regarding image processing in rectangular grids. Prior information regardless of its origin is used to determine relative position (rank) of within-row and/or within-column sampling units prior to selection of samples. For example, ranks, 1, 2, . . . , p, of p quadrats in a row allow within-row units to be stratified based on their relative positions (not based on their actual unknown values ) with respect to each other. Ranks in rows and columns then lead to a two-way stratification (row-column structure) of the sampling region. We introduce a sampling method to provide a better spatial coverage than a simple random sampling method over this two-way stratification. In this paper, sampling (rank) coverage, representative sample or spatial pattern coverage are used interchangeably to describe a representative sample that consists of observations across the ranks of the auxiliary ranking variable.
To highlight the significance of the use of relative position information, we consider a sugarcane disease caused by the fungus Pachymetra chaunorhiza (Chau). This disease is associated with yield loss up to 40% in Australia. The fungus invades individual roots and breaks down the internal root tissues. The variable of interest (y) is the amount of disease (log-DNA copies/g sample) that sugarcane plants have in each grid cell. The sampling region consists of 49 grid cells organized into 7 by 7 square plots of sugarcane plants. The values of y are presented in panel (b) of Table 1. The population mean and standard deviation are 15.87 and 15.56 log-DNA copies/g, respectively. In this population, we consider the estimation of the population mean. Before selecting the sample units, we create a two-way stratification based on the relative position (rank) of within-row and within-column units of the square. Ranks can be obtained from a proper ranking procedure that does not require actual measurement of the units, such as brightness of leaf colors, which can be obtained by flying a drone equipped with image processing software. The quadrats in rows and columns then can be ranked independently using the brightness or color intensity of the leaves. Since the row and column rankings are independent, the ranking in the entire square may not have a joint perfect ranking even if the rankings in rows and columns are perfect. On the other hand, it provides a reasonable stratification in both row and column spaces. Once the twoway stratification based on ranks is completed, we select a better representative sample than a Latin square sample by exploiting the ranking structure of the quadrats in the sampling region. The detail of the sampling procedure is provided in Sect. 2.
The row-column structure is similar to a Latin square sampling design, where the population units are organized in a p by p square of grid cells based on two stratification criteria, row and column stratification. A Latin square sample consists of the selection of p cells such that exactly one unit from each of the square's rows and columns is selected. To our knowledge, Vickery et al. (1949) is the first publication on Latin square sampling, where the average content of a mineral in plant leaves is estimated in a controlled environment. Since then there have been many contributions on the topic. Further details and relevant references about Latin square samples and related sampling designs can be found in Borkowski (2003), Cochran (1972), Kish (1965, Munholland and Borkowski (1996), Raj (1968) and Yates (1960).
The sampling design in this paper differs from existing designs in the cited references. It has potential benefits for constructing a better representative sample by exploring the heterogeneity of the sampling region. The new design, unlike the Latin square sampling which relies only on the grid structure of the region, uses both grid structure and the row and column ranks of the quadrats. Since ranks of quadrats are more closely related to the values of the variable than the positions of the quadrats in a Latin square, they provide a better row-column stratification and lead to a better representative sample. This paper presents a procedure relevant to the Latin hypercube sampling with ancillary variables by Minasny and McBratney (2006). Their procedure requires a knowledge of probability distributions of ancillary random variables X 1 and X 2 , which would be a strong assumption to meet, whereas our procedure relies on relative closeness of units in rows and columns. Hence, our procedure would be less restrictive for practical application.
Use of ranking information in sampling is not new. It has been successfully used in ranked set sampling (RSS) to reduce the sampling cost by stratifying sampling units based on their relative positions in a field experiment (McIntyre 2005). An RSS of set size H and sample size n is constructed in three basic steps.
Step 1 selects n H grid cells at random from the population in the sampling region, and partitions them into n sets each of size H . Step 2 ranks the H grid cells in each of those n sets using within-set ranking information and selects cells for the next step.
Step 3 measures the actual value of the variable of interest in the h-th ordered grid cell (ranked based on prior information) in the n h sets, Y i [h] , i = 1, . . . , n h , where H h=1 n h = n, and n h is the number of observations from the h-th ranking group. The final set of the n measured cells is called a ranked set sample. Details of RSS design can be found in Deshpande et al. (2006), McIntyre (2005, Murray et al. (2000), Ozturk (2017Ozturk ( , 2019, Kravchuk (2021, 2022), Ridout (2003), Takahasi and Wakimoto (1968), Wolfe (2004Wolfe ( , 2012 and the references contained therein. RSS design uses relative closeness in a comparison set of a one dimensional space, but ignores the closeness in row and column structure of a rectangular grid. This paper combines the stratification principle in RSS with the stratification principle in Latin square sampling to create a two-way structure for the sampling region. Stratification of rows uses similarity (or dissimilarity) based on ranks of within-column units (or adjacent units within the column) while stratification of columns uses the relative closeness of within-row units with respect to the variable of interest in the sampling region. The sampling method provides a reasonably good representative sample in both row and column directions of the sampling region. Section 2 introduces the design structure and its properties in a finite population setting. Section 3 constructs an estimator for the variance of the sample mean of the new sampling design. Section 4 provides empirical results to compare the proposed sampling design with its competitors. Section 5 investigates the design effects of the new sampling methodology for three examples from field sampling. Section 6 provides some concluding remarks. All proofs are provided in the appendix.

THE DESIGN STRUCTURE AND ITS PROPERTIES
In addition to the fungus population described in the introduction, we provide two more populations to motivate the ranking structure of a sampling region. Munholland and Borkowski (1996) considered estimating the total abundance of trees in the Rocky Mountain area in the USA. The region was subdivided into a square of 64 grid cells, which has a distinct elevation gradient with higher elevations occurring north to south through its center. The abundance of trees appears to be influenced by the change of elevation in the region. The data set is presented in panel (a) of Table 1 for easy reference. We use this data set as our population and estimate the mean of tree abundance in this population. The true values of the mean and standard deviation of abundance of trees are 15.92 and 7.38, respectively. This population naturally falls into a grid structure. Even though Latin square sampling uses this grid structure to create a better representative sample than a simple random sample, it fails to accommodate the use of available auxiliary information, such as vegetation cover of the sampling region. Total abundance of trees in each grid cell is highly correlated with the vegetation cover. The assessment of vegetation cover may be performed using smartphone based digital technologies in agricultural research (Heinonen and Mattila 2021). If there is no strong decreasing or increasing trend in row and column directions, one can first rank within-row cells and then independently rank the within-column cells using the vegetation cover from the images of aerial photographs. The grid cells are then rearranged based on their row and column ranks to create a two-way stratification based on row and column rank structure. The new sampling methodology selects the sample units from this new ranking structure.
The third example is again related to sugarcane plants. It involves another disease, caused by the nematode Pratylenchus zeae (Zeae). The nematode causes mechanical damage and discrete light-brown lesions in the root cortex; the lesions become necrotic, darken and coalesce as the nematode population increases. The disease eventually kills the plant or causes a significant amount of yield loss. In this example, the variable of interest is the amount of disease due to nematodes (log-DNA copies/g sample). The population values of y for 49 grid cells are given in panel (c) of Table 1. The mean and standard deviation of this population are 277 and 212.84 log-DNA copies/g, respectively. In this case, bright green leaf colors (or freshness) could be associated with healthy plants and faded green color could be related to the severity of the disease (or large values of y). Flying a drone equipped with an online camera could provide the images of all 49 grid cells. These images lead to a two-way stratification based on row and column ranking similar to the one in the abundance population. In Fig. 1, we also provide heat map plots of the sampling regions of these three populations to show the spatial patterns in row and column directions. The abundance data shows a strong quadratic trend in rows and relatively strong decreasing trend in columns. The fungus data does not show a particularly strong trend in either direction, but a somewhat increasing trend in columns. The nematode data shows a weak quadratic pattern in columns, and little trend in rows. We show in Sect. 5 that weak spatial patterns in fungus and nematode populations can be improved using one-or two-way stratification using within-row and or within-column ranks.
In these three examples, ranking population units does not require their manual measurement. In general, ranking can be performed based on visual inspection, a yield map, soil analysis of the sampling region or an auxiliary variable obtained by flying a drone equipped with image processing software. Hence, ranking cost will be considered as minimal and have not been included.
We now provide a general setup for the sampling methodology. We consider a finite population of sampling units u i, j , i = 1 . . . , r p, j = 1, . . . , p, arranged in a grid of dimension r p rows and p columns, such that N = r p 2 . The variable of interest measured on unit u i, j is denoted by y i, j . The population values of y will be denoted with y = y 1,1 , . . . , y r p, p . Prior to sampling, we rank the population units in each row, column or both row and column independently and arrange them based on their rank order as is explained in the three examples above. For example, if we consider ranking within-row units only, we arrange the unit having rank 1 as the first unit in the row, the unit having rank 2 as the next, etc. Rank order of units does not necessarily replace the neighboring units. Rank ordering is constructed using relative closeness based on perceived values of y on the experimental units. The perceived value describes the value of y based on the perception of the ranker using all available ranking information. Hence, units with rank 1 and p could be located next to each other in physical location, but in the two-way grid layout they could be replaced as the first and last in a row. After completing the within-row ranking process and arranging the units by their within-row ranks to form rank columns, if there is no strong increasing or decreasing trend for within-column units, we perform a similar ranking process for within-column units independent of the within-row ranking process and re-arrange this time the units by their within-column ranks to form rank rows. This creates a two-way stratification based on within-row and within-column ranking of the units. We note that row and column rankings are performed independently, Hence, joint ranking may not be in perfect order even if the row and column rankings are perfect.
We superimpose a row-column treatment allocation design on row-column ordered rectangular grid cells. Latin letters are used to label the treatment levels in a row-column design in a design of experiment, but these letters are used here to sample population units rather than for a treatment allocation as in the design of experiment. The basic structure of a row-column sampling design is similar to a row-column design in a design of experiment.
To have a better sampling coverage, we require that each letter appears in each row and column exactly one and r (r > 0) times, respectively. The choices of row-column structures of Latin letters will be denoted with L l , l = 1, . . . , N . In order to account for variation due to different row-column structures, we permute independently the rows and columns of the matrix of Latin letters. For p = 5 and r = 2, the rectangular grid with a particular row-column structure (L 1 ) is represented in Table 2. In Table 2, the two 5 by 5 Latin square letters are stacked together for the sake of illustration. This need not be the case in practice since we independently permute the rows and columns of Latin letters. Sampling units need not be arranged in a rectangular grid. If the sampling frame is arranged in a linear region, such as a riparian region in a river system, a conceptual r p × p rectangular grid can be defined with each row corresponding to the runs of p consecutive units and columns representing the ranks of p units in each run. The rectangular grid here corresponds to r p runs along a linear region. In this case, the row-column sampling design not only selects a unit in each run, but also a unit(s) from every rank sequence position.
To simplify the notation, we rewrite the population units in a vector form v and denotes the transpose of the vector. The entire set of population units is represented in The vector v is constructed in such a way that its first p entries are the population units in the first row, the second p entries are the population units in the second row of the ordered grid population, etc.
In this population structure, a row-column sampling procedure with sampling design L l involves selecting m (1 < m ≤ p) letters from the collection of p letters at random and identifying corresponding population units as a row-column sample (RC S l ). For example, with design L 1 in Table 2, we have p = 5 , m = 2 and r = 2. If a particular random selection of letters from {A, B, C, D, E} yields B and D, then our RC S 1 would contain the units In general RC S l contains rmp observations, r p of which are selected using the first letter, the next r p units are selected using the second letter, etc. For our example in Table 2, when p = 5, m = 2 and r = 2, RC S l sample could be defined as follows where s i, j is the j-th ( j = 1, . . . , m) selected sample unit in row i of the row-column grid population. We note that the RC S l is a balanced sample since there are m and mr observations from each row and column ranking groups, respectively. The RC S l structure helps to explore the two-dimensional space. Ranking between row units explores the heterogeneity among rows while ranking between column units explores the heterogeneity among columns through ranking of experimental units prior to sampling. We do not assume that the ranking process is exact. It could be subjective and in error. The ties in ranking can be resolved at random. Hence, it covers a wide range of sampling designs. For example if r = 1 and the ranking process is completely random in both row and column spaces, the RC S l sample is equivalent to a simple random sample of size mp. Any improvement in ranking process yields a different design and provides higher design efficiencies. The design efficiency with respect to Latin square sample depends on monotonic trends in row and column spaces and the quality of the ranking process. If the quality of ranking is strong enough to improve the existing trend in row and column spaces, the RC S sample improves the L S sample. Otherwise, it may be less efficient than a Latin square sample.
In this finite population structure, we are interested in estimating the population mean or total. Since N is known, we only consider the estimation of population mean. The population total can be estimated by multiplying the sample mean by N .
Theorem 1 In a row-column sample, for a given row-column structure of Latin letters L l , inclusion probabilities are given as follows where N * L is the number of total unique row-column designs, and I kt (l) = 1 if v k and v t share the same letters in RC S l , and zero otherwise.
We now introduce our estimator of the population mean where Z k = 1 if the unit v k is in RC S l and zero otherwise. In this formulation, Z k is a random variable, and y k is a constant. This estimator can be written in a slightly different formȳ Note that Y s j is a random variable here since s j is random. We note that the estimators for both the RC S l and S RS designs are simply the sample means since the inclusion probabilities are always m/ p. For comparison purposes, we also consider a simple random sample (SRS) using the sampling without replacement selection procedure from the same population The marginal and joint inclusion probabilities for S RS are given by π * k = rmp r p 2 and π * k,t = (rmp−1)rmp (r p 2 −1)r p 2 . An estimator for the population mean using the S RS sample is given bȳ where y * s j is the value of y on unit v * s j .

Theorem 2
The estimatorsȳ RC S l andȳ S RS are unbiased for the population mean μ = N k=1 y k /N . The conditional variance ofȳ RC S l , for a given row-column design L l , is given by σ 2 The variance ofȳ S RS is given by We note that σ 2 RC S l is a conditional variance for a given design L l . Since the estimator y RC S l is unbiased for τ for l = 1, . . . , N * L , and each design is equally likely, the uncondi-tional variance of the RC S design is given by We now investigate the conditions under which the RC S estimator is more efficient than a simple random sample estimator. Both sampling designs have the same marginal inclusion probabilities (π k = π * k ), but they differ in their joint inclusion probabilities (π kt (l) = π * kt ). It follows that σ 2 RC S l ≤ σ 2 S RS if and only if In the above inequality, we partition the summations with respect to I kt (l) = 1 and I kt (l) = 0 and organize the terms to write Note that the summation on the left has r p 2 (r p − 1)/2 terms, while the summation on the right has r p 3 ( p − 1)/2 terms. Hence, if all y values are positive, we expect that this inequality will hold for most populations we encounter in practice. Note that using the conditional expectation we can write where I s 1 s 2 (l) = 1 if sample units v s 1 and v s 2 share the same letter in design L l and zero otherwise. Using these equalities in the inequality (2), we conclude σ 2 RC S l ≤ σ 2 S RS if and only if In inequality (3), the right-hand side has a larger number of joint expectations from different ranking groups than the left-hand side. Since order statistics are positively correlated, joint expectations from different ranking groups are usually larger than the joint expectation from the same ranking groups. Hence, this inequality indicates that the proposed RC S design is more efficient than a simple random sampling design if the ranking groups are well separated from each other. Since the largest separation happens under perfect ranking, we conjecture that the largest efficiency improvement is achieved under perfect ranking.

VARIANCE ESTIMATE
All joint inclusion probabilities are positive. Hence, the variance estimator of the Horvitz-Thompson estimator is unbiased. Using the Horvitz-Thompson variance estimator, an unbiased estimator of σ 2 RC S l is given bŷ In a similar fashion, an unbiased estimator for σ 2 S RS is given bŷ We note that even thoughσ 2 RC S l andσ 2 S RS are unbiased estimators for σ 2 RC S l and σ 2 S RS they are unstable estimators. If the sample units s j and s r do not share the same letters in RC S l and have large y values, the estimator in Eq. (4) frequently produces negative values for some sample selections. For this reason, we introduce another estimator for σ 2 RC S l based on bootstrap methods.
The RC S is a probability sampling that considers the spatial structure in the grid population. A bootstrap resampling procedure should capture this spatial structure of the sampling region. For a given RC S sample RC S l = Y s j , j = 1, . . . , r pm , we construct the bootstrap resample RC S b l = Y b s j , j = 1, . . . , r pm by selecting an RC S sample from the population grid using design L l for b = 1, . . . , B. All empty units in the bootstrap resample RC S b l are then filled with the nearest neighbors in RC S l sample. The bootstrap variance estimate is then given byσ 2 We use the following algorithm to compute the bootstrap variance estimate.

Algorithm 1. Bootstrap variance estimate and confidence interval:
• Input: RC S l = Y s j , j = 1, . . . , r pm and design L l .
• For b = 1, . . . , B do • Permute the rows and columns in Latin square design L l , • Select m letters from the pool of letters in L l , say, letters A and C for m = 2, • Select all units corresponding to these m letters, RC S b l = Y b s j , j = 1, . . . , r pm , • Replace each unit Y b s j with its nearest neighbor in RC S l , The neighborhood of a unit in a grid cell is defined in terms of its Euclidean distance from another point. Table 3 presents a bootstrap resample RC S b l constructed from an RC S l using the fungus population described in Sect. 5. In this table, the underlined entries represent the RC S l that corresponds to the randomly selected two Latin letters. The bootstrap resample units are identified by a bold face font that corresponds to randomly selected Latin letters. For example in Table 3, the entry in row 1 and column 4 is a unit in the bootstrap resample. Since this point is not observed in the original sample (not underlined), it has to be filled with its nearest neighbors from the underlined units. The closest underlined units are 5 and 4 in cells (1,3) and (2,4), respectively. So the bootstrap resample for this point becomes the average of 5 and 4. The bootstrap resample for cell (1, 3) is 5 since it is already in theRC S l The other bootstrap resample points are computed in a similar fashion.
An approximate (1 − α)100-percent bootstrap percentile confidence interval for the population mean μ can be obtained from Algorithm 1. We can also construct a (1 − α)100percent confidence interval using a normal approximation where t a,d f is the a-th upper quantile of a t-distribution with degrees of freedom d f = mr p − 1.

DESIGN EFFECT
In this section, we investigate the design effect of the RC S design. Since RC S considers two-way stratification, we look at six different populations having different spatial structures. The values of the variable y in these populations are generated with different mean functions (μ i j ) in row and column space. For notational convenience, we define i 1 = i − 1 and j 1 = j − 1. The mean functions for population P t , t = 1, . . . , 6, are given below where sin and cos are the sine and cosine functions. In these six populations, all means are shifted to make the smallest cell mean 10. The population values y of each cell is then constructed by adding a random number between 0 and 5 generated by a uniform distribution. For a fair comparison among the six populations, all population units are scaled by dividing each value by the population standard deviation. The heat maps of these populations are presented in Fig. 2. The population P1 does not show any spatial pattern. The populations P2 and P3 have decreasing trends in row and column directions, respectively. The populations P4, P5 and P6 use mean functions for the corrugated inclined plane, peak function and bird functions in Robertson et al. (2018). These six populations cover a wide range of spatial patterns that we might see in practice. We performed a simulation study to investigate the design effect of RCS design with respect to L S sampling design using these six populations. Simulation parameters are selected to be m = 2, r = 1 and p = 30. For row and column ranking, we used the Dell and Clutter (1972) model. This model creates a ranking variable U correlated with a response variable y, ρ = corr (U, y). In the simulation study, the values of ρ are selected to be ρ = 0, 0.2, . . . , 1. The details of this ranking model can be found in Ozturk (2019) and Ozturk and Kravchuk (2021). The simulation size is taken to be 5000. For one-way stratification based on ranks, in each replication of the simulation, before selecting an RC S sample only the within-row (within-column) units are ranked independently using the Dell and Clutter model with a given value of ρ. No ranking is performed for within-column (within-row) units. Rows and columns of design L l are permuted so that each replication uses a different Latin square design. Hence, the simulation study accounts for the variation due to within-row ranking error and selection of row-column design L l , l = 1, . . . , N * L . For two-way stratification based on ranks, in each iteration, within-row units in the population are first ranked independently in each row using the Dell and Cutter model, and then within-column units in each column are ranked independently. Again rows and columns of design L l are permuted in each iteration. Figure 3 presents the design effect of the RC S sample with respect to the L S sample for this simulation setting. The design effect is defined as the ratio of variances whereσ 2 RC S is the average of σ 2 RC S l in equation (1) from the 5000 simulation replications. The values of D e f f greater than one indicate that the RC S design is superior to the L S design.   Figure 3 provides the D e f f curves of the six populations for within-row (Row), withincolumn (Col.) and row-column (RowCol.) rankings when m = 2, r = 1, and p = 30. The design effects for each ranking category are plotted against the ranking correlation ρ = 0, 0.2, . . . , 0.8. The design effects are much larger under the perfect ranking (ρ = 1) than under the imperfect ranking (ρ < 1). Hence, they are reported separately in Table 4 to illustrate the detailed feature of D e f f in an appropriate plotting scale when ρ < 1. The heat maps in Fig. 2 indicate that the population P1 does not show any spatial trend either in the row or column spaces. Hence, ranking in one of the row, column or row-column directions yields higher design effect for the RC S sample. For the population P1, all three curves lie above 1, indicating that the RC S sample is more efficient than the L S sample. The populations P2 and P3 have a strong increasing and decreasing trend only in one of the row and column spaces. They show no trend (uniform) in the other direction. For example, the within-row units in the population P2 have increasing order, while the within-column units are uniformly distributed. For the population P2, D e f f values are larger than one for the within-column units orderings for all ρ values. The row and row-column orderings perform poorly for the population P2. This intuitively makes sense since there is already a strong increasing trend in the within-row units, and any ranking in this direction with ρ < 1 will destroy this trend and reduce the degree of stratification of the within-row units. Similar observations can be made in the population P3, where there is a strong negative trend in the column direction. Hence, any ranking of the within-column units with (ρ < 1) will reduce the design effect of the RC S sample with respect to the L S sample.
The populations P4, P5 and P6 have strong spatial patterns, but they do not necessarily show a strong increasing or decreasing trend in either the row or column space. Hence, ranking procedures in the RC S sample create strong increasing trends in the within-row and column units, and lead to a better two-way stratification than a Latin square sample. This yields a higher design effect than a LS sample for a moderate ranking correlation ρ. Overall patterns in these populations indicate that the effect of the RC S design is more pronounced for populations lacking a strong trend in row and/or column space. In this case, the ranking procedure exploits the relative closeness (ranks) of sampling units in the grid structure to create a strong two-way stratification, and provide a better spatial coverage in the sampling region. Table 4 presents the design effects of the RC S sample under a perfect ranking (ρ = 1). In this case, the RC S sample outperforms the L S sample regardless of the spatial patterns. The reason for this is that even if the grid structure of the population has a strong linear trend, a perfect ranking does not destroy this trend. Hence, it always improves the L S sample. The improvement is substantial. The smallest D e f f value in Table 4 is 4.883 for the population P5 and row-column ranking, which indicates that the RC S sample is approximately 5 times more efficient than an L S sample.
We performed another simulation study to investigate the properties of the bootstrap variance estimates of RC S sample means. In this part of the simulation, we only considered population P6. Simulation parameters are selected as m = 1, r = 1. Simulation size was set to be 50000 replications. Again population units are ranked by using Dell and Clutter model with ranking correlation coefficient ρ = 0, 0.2, 0.4, 0.6, 0.8, 0.9, and 1. Figure 4 illustrates the box plots of the bootstrap variance estimates of the RC S mean estimator for different values of ρ. The first, second and third rows present the box plots for the within-row, within-column and row-column units ordering, respectively. The true value of the variance (variance of the 50000 RC S sample means) is marked with a triangle on each box plot. The box plots indicate that bootstrap variance estimates are skewed right and provide under estimation for the true variance. The amount of bias decreases for large values of ρ. When ρ = 1 the true variance is at around the median lines in the box plots. Under estimation of the variance of the estimator is a common problem in a spatially balanced sampling designs in the literature (Stevens and Olsen 2003;Grafström and Schelin 2014;Robertson et al. 2018).
We also computed the coverage probabilities of the approximate 95% bootstrap and normal approximated confidence intervals of the population mean using the RC S sample means. Table 5 indicates that coverage probabilities are lower than the nominal coverage probability 0.95 for both bootstrap and normal (t-distribution) approximated confidence intervals. On the other hand, the confidence intervals based on normal approximation provide better coverage probabilities than the bootstrap confidence intervals, specially for large ρ values. For example, if ρ ≥ 0.8 the coverage probabilities are all greater than 0.90, and closer to 0.95 for ρ = 0.9 and ρ = 1.

EXAMPLES
In this section, we performed another simulation study using the RC S and L S sampling designs to estimate the population mean μ of the tree abundance, fungus and nematode populations that we introduced in Sects. 1 and 2. For the simulation parameters, we selected m = 1, r = 1, and ranking correlations ρ = 0.1, 0.2, . . . , 1.0. As before, the Dell and Clutter model is used to rank the within-row and -column units. Simulation size is taken to be 5000. Simulation results are presented in Fig. 5. The three curves in the panels represent the design effect of the RC S sample with respect to the L S sample for row, column and row-column ordered population units.
For the tree abundance population, the RC S sample is less efficient than the L S sample when ρ < 0.8. Note that the heat map of this population in Fig. 2 indicates that there is a strong trend in both row and column directions. Hence, use of the imperfect ranking weakens the strong spatial structure and reduces the efficiency of the RC S sample. Since the trend in the column direction is not as strong as in the row direction, the dashed red curve corresponding to the within-column units ordering lies above the other two curves corresponding to the within-row units and row-column units orderings. Under perfect ranking, the RC S sample is superior to the L S sample in both one-way and two-way ranking stratification. In Fig. 2, while the heat maps are indicating some increasing and decreasing trends in the column direction, they show a small increasing and decreasing trend in the row direction for the fungus and nematode populations, respectively. As expected, the within-row and rowcolumn units ranking (solid black and dotted green curves) out performs within-column unit ranking. These findings are consistent with our simulation results in Section 4, where we showed that RC S samples outperform L S samples if there is no strong increasing or decreasing trend in row and column directions.

CONCLUDING REMARKS
A Latin square sample ensures a good coverage for the sampling region, but it does not provide a variance estimate for the estimator unless a replicated sample is generated. This article shows that when there is a spatial trend among the units that is not strongly increasing or decreasing in row and column directions in the sampling region, latinized row-column ranked sampling design is generally more efficient than a simple random sample and a Latin square sample of the same sample size. A row-column sample creates a two-way stratification in rows and columns. The cells in a rectangular grid in the sampling region are reorganized based on the relative closeness (ranks) of units within rows and/or columns. The sampling method provides a better sampling coverage by re-indexing physical grid quadrats and ensuring that at least one observation is selected from each rank-row and rank-column. The ranking can be done on the variable of interest and/or one or more than one auxiliary variables. In this paper, we consider a ranking procedure with a single auxiliary variable. We have constructed point estimators, for the population mean and total and their variances. We have also constructed an approximate (1 − α)100% confidence interval for the population mean.
We have investigated the design effect of the new sampling method for a wide ranges of spatial population structures. The ranking structure in the RC S sampling design is superior to the usual L S sampling if there is no strong increasing or decreasing trend in the row and/ or column space of the grid structure. For reasonable ranking information, such as when the correlation between the ranking variable and the variable of interest is greater than 0.4, the proposed sampling design performs better than Latin square and simple random sampling designs of the same sample size as long as the sampling region does not show a very strong increasing or decreasing trend in the rows and columns of the grid structure.
Two confidence intervals are suggested for the population mean. The first interval is the bootstrap percentile confidence interval. The second interval is constructed using the bootstrap variance estimator of the sample mean with a normal approximation. The coverage probabilities for bootstrap confidence intervals are smaller than the nominal coverage probability of 0.95. They are also smaller than the coverage probabilities of the confidence intervals based on normal approximation. The coverage probabilities of the confidence intervals based on normal approximation approach to the nominal coverage probability of 0.95 for larger values of ρ.
In some sampling regions, the size of the grid cells may be large and it may not be easy to quantify the value of the variable of interest y. In this case, we may consider a two-stage sampling design. The stage-one sample would select an RC S sample from the primary population of a rectangular grid. The stage-two sample would select a sample from each of the grid cells in the stage-one sample. Development of this sampling design will be a future project.
Proof of Theorem 2. We first consider the expected value ofȳ RC S l E(ȳ RC S l ) = 1 p 2 N k=1 E(Z k )y k π k = 1 p 2 N k=1 y k = τ.
For the variance, we write This completes the proof. The proof of σ 2 S RS is similar. Hence, it is omitted.