As the complexity and cost of biological experiments has grown considerably in recent years, partly due to technological advances (high throughput technologies and more), there is an increasing need to design experiments that maximize the information content of the collected sample. For most standard statistical analyses, where the aim is to estimate some population parameter, maximizing information translates into minimizing the variance associated with a parameter’s estimate. In many situations, researchers observe multiple outcomes for each unit in the sample and wish to make inferences on a parameter of the underlying population’s joint distribution, routinely this is done via estimating the population mean vector. It is often the case that some or all of the individual components of this response vector are costly, risky (complications due to biopsy), or even destructive (requiring animal sacrifice). In such cases it may be desirable, for monetary or ethical reasons, to extract information from each unit that is sampled, without taking the exact measurement of the response of interest for each unit.

The most common approach for data collection method for making inference about population parameter is simple random sample (SRS) from a population. Even though each subject selected by SRS has an equal chance of being selected from a population to ensure the representativeness of a population, there is no guarantee that the selected sample will truly represent the population. However, the only guarantee one can have is that if the sampling process is being repeated over and over again, then the average of the attribute of interest for multiple SRS would provide the good estimator of the population value of the attribute. Ranked set sampling (RSS) (McIntyre, 1952) is a type of sampling scheme which allows researchers to use information from each unit in the sample, without taking every unit’s exact measurement. The overall goal of the RSS is to obtain the sample from a population that is more likely to span the full range of the values in the population to have a more representative sample than the SRS of similar sample size. Traditionally, RSS can be used provided there is a reliable ranking mechanism available, which should be cheaper or safer than exact measurement, for the response of interest. The ranked but unmeasured units provide increased information over SRS of the same size improving parameter inference. The additional information provided by ranking is due to the fact that aspects of population structure are encoded through the order statistics. Knowledge of observations’ order statistic and exact measurement improve inference since ranked units target different population attributes, unlike the identically distributed unit from a SRS. This has been shown in many works to translate into improvement in parameter inference compared to simple random samples of the same size.

In many situations, the outcome of interest is correlated with some auxiliary variable which may be easier to measure than the outcome of interest. For instance weight may be correlated with fasting blood glucose and may be easily obtained whereas some lab measurement would be necessary for blood glucose measurement. The application of RSS has appeared in series of papers. See for example, Chen (1999), Demir and Çıngı (2000), Huang et al. (2016), Jabrah et al. (2017), Kaur et al. (1996), Samawi and Al-Sagheer (2001).

An outline of the paper is as follows. In Section 2 we introduce the necessary notation and prove that mean estimation is unbiased with a smaller variance for RSS as compared to SRS. In addition, in Section 2, we also derived the limiting distribution of Hotelling’s statistics (Q) as well as the multivariate regression estimator using RSS. In Section 3, we perform a simulation study to compare the performance of RSS to SRS in terms of estimation as well as hypothesis testing. In Section 4, we apply the method on a real data set in the context of public health. We give concluding remarks and future directions for the method in Section 5.

2.1. Ranked set sampling procedure

In this section, we will briefly describe how a ranked set sample may be collected in this section for a univariate random variable. To select the RSS of size n based on the auxiliary variable (X), the following steps should be performed.

• Select the SRS of size r, from a population based on the auxiliary variable (X). r is referred as the set size which is typically between 2 and 5 although any size is possible. However, sizes larger than 5 may become impractical (Takahasi and Wakimoto, 1968).

• Order the auxiliary variable and choose the minimum of (X₍₁₎). Measure the multivariate outcome of interest Y₍₁₎.

• Select the SRS of size r and order it based on the auxiliary variable again. Choose the second minimum (X₍₂₎) and measure the multivariate outcome of interest Y_[2].

• Repeat this process until the X_(r) and Y_(r) of r^th independent SRS are obtained.

• The entire process of obtaining (X₍₁₎, X₍₂₎, …, X_(r)) and (Y₍₁₎, Y₍₂₎, …, Y_(r)) is called a cycle.

• Repeat m independent cycles to obtain a RSS of size n = rm.

Table 1 represents the structure of RSS. For more details about RSS (Jozani and Johnson, 2011; Kowalczyk, 2004; Patil et al., 1995; Takahasi and Futatsuya, 1998).

2.2. Multivariate naive estimator

Our population of interest is an univariate auxilliary variable X and a d dimensional multivariate outcome Y, with a covariance structure on the joint distribution of (Y, X) given by $\mathrm{\Sigma }=\left[\begin{array}{cc}{\mathrm{\Sigma }}_{11}& {\mathrm{\Sigma }}_{12}\\ {\mathrm{\Sigma }}_{21}& {\mathrm{\Sigma }}_{22}\end{array}\right]$. Assuming that we have collected m RSS cycles of set size r, where the ranking has been done on X, we denote the data (X_(i)k, Y_[i]k), i = 1, 2, …, r, K = 1, 2, …,m. Note that the subscript on X indicates that ranking has been done on X and the subscript on Y indicates that ranking on X may result in imperfect ranking on elements of Y. The naive estimator is defined as ${\widehat{\mathit{\mu }}}_{y\text{RSS}}=(1/rm){\sum }_{k=1}^m{\sum }_{i=1}^r{\mathit{Y}}_{[i]k}$. It is straightforward to show this is an unbiased estimator of the mean of Y. Since ${\sum }_{i=1}^r{f}_{X_{(i)}}(x)=r{f}_{X_{(i)}}(x)$ (Dell and Clutter, 1972)

Similarly for the variance (Dell and Clutter, 1972), by defining μ_[i] = EY_[i] we have

It is clear that ${\sum }_{i=1}^r({\mathit{\mu }}_{[i]}-\mathit{\mu }){({\mathit{\mu }}_{[i]}-\mathit{\mu })}^{\top }$ is positive semi-definite since ∀u ∈ ℝ^d we have u^⊤(μ_[i]−μ)(μ_[i]− μ)^⊤u ≥ 0. Then ∀u ∈ ℝ^d u^⊤(Var( μ̂_ySRS) – Var( μ̂_yRSS))u ≥ 0, or equivalently Var( μ̂_ySRS) ≽ Var( μ̂_yRSS). Under the additional assumption that r ≥ d and X is correlated with each component of Y we have strict inequality.

2.3. Multivariate regression estimator

Regression estimators are used to increase precision in mean estimation by incorporating information in an auxilliary variable. In this case, we assume a linear regression of Y on X

where X and ε are independent and ε is a mean zero residual vector with covariance ∑_ε. Then the regression equation with corresponding data from RSS is

It is worth noting that typically the mean of X, μ_x, is unknown. However, since the auxilliary variable X may be much cheaper to measure one may use the r²m units collected from the first stage of sampling to estimate this quantity as ${\overline{\mu }}_x=(1/r^{2}m){\sum }_k^m{\sum }_i^r{\sum }_j^r{X}_{ijk}$.

Then the regression estimator for the mean of the response is given by

where

It is straightforward to show that μ̄_x and μ̂_x are unbiased estimates of μ_x using similar arguments as in the previous section. When (2.3) holds conditional expectation implies that Eβ̂ = β and EȲ_reg = μ_y, so that the regression estimator based on RSS is unbiased. Also

Since E_Y (Ȳ_reg|X) = μ_y + β(μ̄_x – μ̂_x) the second term above is $(1/r^{2}m){\sum }_{i=1}^r{\sigma }_{X(i)}^2\mathit{\beta }{\mathit{\beta }}^{\top }$. For the first term Cov( μ̂_yRSS, β(μ̄_x – μ̂_x)|X) = 0 so that

2.4. Testing for H₀ : μ = μ₀

Theorem 2

Let {${\mathit{Y}}_{[i]k}^{(t)}$}, i = 1, 2, … r, k = 1, 2, …m and t = 1, 2 are two RSS samples fromand. Let

Then, $\mathit{Q}=\{(r_1{m}_1·r_2{m}_2)/(r_1{m}_1+r_2{m}_2)\}{({\overline{\mathit{Y}}}_{rss}^{(1)}-{\overline{\mathit{Y}}}_{rss}^{(2)})}^{\top }{S}_{rss}^{-1}({\overline{\mathit{Y}}}_{rss}^{(1)}-{\overline{\mathit{Y}}}_{rss}^{(2)})$, for large samples, has the limiting distribution as χ²with d degrees of freedom under H₀ : μ⁽¹⁾ = μ⁽²⁾.

For small to moderate samples, for SRS, under H₀ the Q statistics is distributed as {(N – 1) d}/(N – p) F_d,N–d (Seber, 2009). As explicit distribution of Q statistics is not known, for small or moderate size of RSS samples, we recommend performing hypothesis testing by Bootstrap method. Resampling method for RSS was proposed by (Chen et al., 2004; Modarres et al., 2006). They suggest a natural method to obtain bootstrap samples from each row (within cycle) of a RSS.

In this section, we conducted the simulation study to estimate the multivariate outcome mean and the performance of the hypothesis testing by RSS scheme. We also studied the performance of testing hypothesis of equality of multivariate outcome means for two groups. For estimation of α of testing H_o : μ = μ₀ vs. H_a : μ ≠ μ₀, we considered four multivariate outcomes Y_i (i = 1, 2, 3, 4) with μ = [0.3, 0.3, 0.3, 0.2], variances as ${\sigma }_1^2={\sigma }_2^2={\sigma }_3^2={\sigma }_4^2=4$ and covariances as σ₁₂ = 2.39, σ₁₃ = 1.59, σ₁₄ = 2.83, σ₂₃ = 3.19, σ₂₄ = 1.18, and σ₃₄ = 2.24. The auxiliary covariate (X) was simulated with mean 0 and variance ${\sigma }_x^2=1$. For this simulation study, we considered unstructured covariance among multivariate outcome Y_i as shown below. Moreover, we used autoregressive covariance structure between auxiliary variable X and Y_i with correlation parameter ρ.

The RSS for of X and Y_i were simulated from multivariate normal with mean μ and above variance covariance matrix by following the steps as described in Section 2.1. For comparisons of estimation of α for SRS and RSS, different sample sizes (n = rm) were evaluated by varying the ρ, set size and cycle size. This entire process was repeated 2,000 times. For details of the parameter values, referred to Table 2. Table 2 results demonstrate that we can achieve nominal value for α by using RSS with moderate to large samples, however, for smaller sample bootstrap RSS sampling can achieve nominal value for α.

For estimation of the power of testing H_o : μ = 0 vs. H_a : μ ≠ 0, similar simulation settings were considered as described above except with μ = [0.6, 0.6, 0.6, 0.4]. In addition to that bootstrap power was also calculated by taking 1,000 bootstrap samples for each simulated RSS. Furthermore, MSE of SRS, MSE of RSS and the multivariate naive estimator efficiency were calculated. Table 3 reports the simulation results for estimating the power of testing hypothesis under various simulation settings. We can also report that the power of the test increases as the set size increases with RSS, however, for testing hypothesis RSS gives more power than SRS. As expected, Table 3 also shows that RSS provides more efficient estimates of the multivariate naive estimator in terms of smaller MSEs.

Furthermore, the performance of testing hypothesis of equality of multivariate outcome means for two groups, we simulated two groups with multivariate outcome (Y_i ) (i = 1, 2, 3, 4) with means for the first group μ₁ = [0.3, 0.3, 0.3, 0.2] and mean for the second group μ₂ = [0.6, 0.6, 0.6, 0.4] with similar covariance matrix of Y as described above (Cov(X, Y_i)). Table 4 represents the results of estimation of power of the testing hypothesis H_o : μ₁ = μ₂ vs. H_a : μ₁ ≠ μ₂ with various parameter values of ρ, set size and cycle sizes. Overall, from Table 4, we can conclude that RSS is more powerful for testing hypothesis of equality of multivariate outcome means for two groups compared to SRS.

We also conducted a simulation study to show that the multivariate regression estimator for RSS is more efficient than SRS. We considered multivariate outcomes Y with mean μ = (0.3, 0.3, 0.3, 0.2) and the variance-covariance matrix (Cov(X, Y_i))) as described above in this section. We also simulated correlated auxiliary covariate (X) with mean 0 and variance 1. Table 5 shows that for various parameter settings, the RSS is more efficient than SRS in estimating multivariate regression estimator.

In this section, we illustrate the efficient ranked set sampling method via ranking on baseline covariate to estimate the multivariate outcome mean, investigate the performance of the hypothesis testing for two groups and estimation of multivariate regression estimator by using the China Health and Nutrition Survey (CHNS) for year 2009. The CHNS is the only large-scale household based survey in China (Yan et al., 2003). As a part of the survey, anthropometry were collected on 10,242 children and adults aged ≥ 7 in year 2009 along with other demographic information. Only 9,986 individuals agreed to provide the fasting blood samples which were evaluated for many biomarkers of diabetes and cardio-metabolic risk factors. For illustration purposes, we focused on the variables such as age of the individuals as our ranking auxiliary variable, and cardio-metabolic biomarkers, for example, Apolipoprotein A, Total cholesterol and Hemoglobin A1c. We treated the survey data as a population and selected the range of RSS (N = set * cycle) as shown in Table 6 by ranking on the baseline covariate age. SRS of similar size N was also selected from CHNS data to evaluate the performance of the hypothesis testing and the efficiency of the sampling procedure compared to RSS in estimating the multivariate outcome mean. The correlations (ρ) between age and biomarkers Apolipoprotein A, Total cholesterol and Hemoglobin A1c are 0.12, 0.32, and 0.22 respectively. The mean for Apolipoprotein A, Total cholesterol and Hemoglobin A1c are 1.14 (g/L), 4.78 mmol/L and 5.67 mmol/L respectively, and for comparison purposes, they can be treated as the true parameters. Table 7 represents the power comparison of RSS with SRS for multivariate means of males and females. Table 7 represents that we can achieve more power with RSS compared to SRS with similar sample sizes. Table 8 shows the results for multivariate regression estimation for biomarker data. We also took 1,000 samples of SRS and RSS of sample size 80 (set = 4 and cycle = 20) and plotted the confidence regions as shown in Figure 1. From Figure 1, we can see that the confidence region for SRS (blue nets) lies completely outside of the confidence region of RSS (red).

In statistics, it is important to have a sampling method which is cost effective. RSS is one the important method which can be used to have a more efficient multivariate mean estimator compared to most commonly used method of SRS. The samples taken by using RSS method are more representative samples due to its inherent structure imposed by ranking based on easy-to-available covariates. In this paper, we demonstrated that the RSS is more efficient in estimating the multivariate mean as well as in hypothesis testing for one and two independent samples. Simulation studies for the performance of hypothesis testing showed that the RSS is more powerful compared to SRS. In general, in estimation of the population mean, RSS improves the precision relative to SRS with the same sample size, n. This is true even if the correlation between the auxiliary variable X and multivariate outcome Y is moderate to high (±0.4 to ±0.8). However, when the correlation between X and Y is very low (such as ± 0.001), RSS is equivalent to SRS and the ranking is not better than random. In practice, the key issue is whether the increase in precision is sufficient to justify the increased costs associated with the ranking process. In contrast, when the correlation between X and Y is very high (±0.9 or higher), the precision in estimating the population mean will be very high as this will improve the ranking of X on Y (Ridout, 2003).

Missing data is a very common problem in all most every research and can have a very significant impact on the inferences drawn from the collected data such as biased estimation of population parameters and loss of statistical power (Little and Rubin, 2014). The valid statistical analysis which has appropriate missing data mechanisms assumptions (missing completely at random, missing at random, or missing not at random) should be performed in SRS and in RSS. There is an extensive literature available on how to deal with missing data for RSS in auxiliary variable X and univariate response Y (Bouza-Herrera, 2013). However, handling the missing data in multivariate Y with monotone or arbitrary missing pattern is still the active area of research.