Unbiased estimation strategies for respondent driven sampling

Abstract

In this paper, we focus on respondent-driven sampling (RDS), which is a valuable survey methodology to estimate the size and the characteristics of hidden or hard-to-measure population groups. The RDS methodology makes it possible to gather information on these populations by exploiting the relationships between their components. However, RDS suffers from the lack of an estimation methodology that is sufficiently robust to accommodate the varying conditions under which it is applied. In this paper, we address the estimation problem of the RDS methodology and, by approaching it as a particular indirect sampling technique, we propose three unbiased estimation methods as possible solutions.

1. Introduction

In this paper, we focus on respondent-driven sampling (RDS), which is a valuable survey methodology for both national and international organisations to estimate the size and characteristics of hidden (e.g., homeless people, undocumented immigrants) or hard-to-measure population groups (e.g., minorities, indigenous people).

The principle of “leaving no one behind” is at the heart of the 2030 Agenda and a key requirement for many Sustainable Development Goals (SDG) indicators is to be available for the most vulnerable and marginalised population groups. Nevertheless, halfway through the implementation of the 2030 Agenda, most SDG indicators are still not available at the needed level of disaggregation to monitor the socioeconomic conditions of hidden and hard-to-count population groups. As a result, it is neither possible to produce reliable structural data on the needed disaggregation dimensions nor to monitor the developments of emerging phenomena that need to be approached with targeted evidence-based policy interventions.

The RDS methodology makes it possible to gather information on these populations by exploiting the relationships between their components. Moreover, the effectiveness of the RDS can be further increased by employing an integrated approach in which the RDS is used in conjunction with other information sources, such as administrative or geographical data.

The RDS is a network-based sampling technique [1, 2] that was developed first by Eckathorn [3]. RDS has been the favourite survey method for sampling populations that are difficult to reach due to the potential of a viable sampling technique with reasonable inferential approaches. As a result, since its establishment, it has been employed in countless investigations of these populations across many countries [4]. RDS starts with a small sample of individuals (“seeds”) with which the researchers are familiar. Each participant is then given a small number of coupons with unique identifiers to distribute to their contacts in the target population, enrolling them in the study and growing the sample size until the sample includes the desired number of respondents. The RDS process stops either when, in the selection process, we encounter only units already identified in the previous steps or at a predetermined data-collection step (e.g., the fifth step). Picture 1.1 below illustrates an example of a network sampling process articulated into three steps. The blue lines are the links observed in the sample. Up to and including step 2, participants g, b, c, a, 3, and 4 are kept in the sample. Participants d, 1, 2, 5, and 6 are not observed.

Figure 1.

Example of network sampling process ${}^{*}$ . ${}^{*}$ The thin grey lines represent links existing in the population but not observed in the sample.

We may view RDS as a specific extension of the extensively used group of convenience sampling methods known as “snowball sampling networks” [5] which are frequently employed as a last option when a traditional sample frame is not available [6]. Compared to those of more traditional snowball sampling, RDS offers two key benefits. First, respondents receive few coupons. This enables statistical inference to be more appropriately defined and makes it more plausible to approximate the final sample as a probability sample. Second, asking respondents to pass coupons to their contacts in a potentially stigmatised community reduces potential confidentiality issues. Due to this innovation, RDS is a compelling method for gathering data from marginalized and difficult-to-reach populations.

However, the RDS methodology suffers from the lack of an estimation methodology that is sufficiently robust to accommodate the varying conditions under which it is applied. Although it is quite robust for estimating mean and proportion values [7], the accuracy of the total estimates depends on several features including the nature of the network connecting the individuals in the population and the lack of a rigorous sampling approach to select the sampling units.

In this paper, we address the estimation problem and by approaching the RDS methodology as a particular indirect sampling technique [8], we propose three unbiased estimation methods as possible solutions. In particular, the first method assumes a random sampling of the initial individuals. In contrast, the second method, which considers purposive sample selection, creates a nonbiased estimation if the initial sample of respondents falls into all the clusters of networks that characterise the population of interest. Finally, leveraging the generalised capture-recapture estimation approach [9], we propose an estimator that accounts for the noncoverage of two independent indirect samplings.

The paper is organised as follows. In Section 2, we summarise the traditional methodology of the RDS methodology, illustrating the data collection technique and the Volz and Heckathorn estimator [10], which has been very successful in practical applications due to its lack of computational complexity. Section 3 introduces the basic symbology, and we show how the RDS can be seen as a particular specification of indirect sampling in which each survey wave represents the indirect basis for the subsequent RDS phase. In Section 4, we expand the sampling aspects in the RDS. Section 5 introduces the three estimators. Section 6 concludes the work, and we begin to outline a strategy to overcome information gaps for SDG indicators for hard-to-reach populations, focusing on indigenous peoples.

2. Data collection and estimation in the classical RDS approach

RDS is frequently carried out by using techniques suggested by Salganik et al. [11] and outlined in protocols such as those proposed by White et al. [4] and Johnston [12].

A preliminary sample of typically 2 to 10 seeds is chosen. Aiming to represent all the key socioeconomic subpopulations that researchers anticipate may exist in the target population, seeds are selected to be as varied as possible. The rationale of this derives from the fact that each subpopulation may represent a separate network (or a cluster) of target individuals. If we select a seed in a given subpopulation, we can explore the network of related individuals. In contrast, if we do not select any individual in a subpopulation, the specific cluster of individuals cannot be observed in the RDS process. Therefore, picking up in the initial sample all possible distinct networks increases the possibility of constructing unbiased inferences on the target population.

The enumerators should include community opinion leaders in the initial seeds. Hence, their acceptance and support of the survey method may likely inspire widespread involvement from other target population members. This buy-in is crucial in target populations that are unlawful or stigmatised, especially if the population has any prior exposure to risky research practices. Following an interview, the seeds are given some coupons, each with a unique identification number, to spread to other population members. This number was used to reconcile the practical need to prevent the early termination of the sample trees with the inferential aim of limiting the branching of the sample. Members of the population who receive coupons visit a study centre where they are directly interviewed or given an interview appointment. Three coupons are likewise supplied to subsequent replies; this process continues until the sample size is approximately reached and the coupons are tapered or discontinued. Participants are paid for their time spent taking the survey. Additionally, for each successful recruiting, participants receive rewards. In the survey, the number of target population contacts of each respondent must be measured. This is typically done by asking questions that narrow the recruit’s references to the precise definition of the target population. Interviewers must also verify membership in the target population. Researchers must also assure participants do not participate in the survey more than once. Study staff are familiar enough with the target population in many settings to notice repeat participation attempts. In other cases, repetition is prevented by collecting nonidentifiable but unique information about participants, as in Johnston [13].

The RDS methodology can be applied alternatively to the entire population of individuals or by considering only the subpopulation at risk of belonging to the target population. For instance, if the target population coincides with forced labour people, we may observe people working in sugarcane.

Figure 2.

Example of target population divided in disjoint clusters.

Contrary to what was previously believed, Eckathorn [3] used a Markov modelling of the peer recruitment process to demonstrate that bias from the convenience sample of beginning participants from which the sample started gradually diminished as the sample increased wave by wave. By using the model, they demonstrated how the sample approached an equilibrium independent of the beginning location or independent of the convenience sample of seeds from which it started as it expanded wave by wave. The conclusion was that this sampling technique may become reliable if there were enough waves, meaning that any seed selection can eventually yield the same equilibrium sample composition. However, the researchers did not show how an unbiased estimate can be derived.

Eckathorn[14] introduced the first RDS population estimator based on the essential idea that in RDS, relationships tend to be reciprocal. This implies that if person A knows person B, then B knows A.

The estimator bases its validity on the principle of network balance between population subgroups. Up to a constant factor, the volume of network connections to and from each group can be approximated. For each pair of groups, this results in a set of balance equations that may be used to solve for the relative size of each group. Volz and Eckathorn [10] proposed a slightly biased estimator. In the following we call this estimator the VH estimator. The VH estimator is based on the following hypotheses [15]: 1. The network size is large compared to that of the realised sampling, including the initial seeds and the respondents recruited by the RDS process. 2. Homophily is weak enough, where homophily is the tendency for nodes to preferentially form network contacts with others like themselves. 3. Reciprocity of contacts. 4. With-replacement sampling. 5. Enough sample waves. 6. Accurate measurement. 7. The recruitment in the subsequent waves of the RDS process is random.

The first three hypotheses relate to the nature of the contact network, while Hypotheses 4 to 7 relate to sampling. Hypothesis 4 is the most critical and may introduce a certain level of bias, as the sampling process adopted assumes a link between the persons recruited.

Focusing on the first three assumptions, let us consider the example below in Fig. 2, where we assume that people of the target population belong to three disjoint clusters. If the starting sample includes only persons in one group, the traditional RDS can estimate only the total number of persons in that cluster suffering from a substantial undercoverage problem. Since people of the target population are often grouped geographically, observing each cluster’s units in the starting sample is appropriate.

The VH estimator has had great application success in the practice of real investigation due to its great computational ease. Subsequently, other estimators have been proposed in the literature (see among others [16, 17]) each of which overcomes some of the limitations associated with the assumptions made with the VH estimator. These estimators have higher computational complexity, and introduce some modelling assumptions on the cluster variables or the nature of the contact network.

To illustrate the classical VH estimator, let us consider the case where the total of a characteristic

$Y={\displaystyle \sum _{k\in U}}{y}_k$(1)

in the population of interest $U$ (of $N$ units) is to be estimated, where $y_k$ is the value of $y$ for unit $k\in U$ .

Let $S$ be the sample of different units at the end of the RDS process. Let ${\alpha }_k$ be the number of times unit $k$ is observed in $S$ , where

$m={\displaystyle \sum _{k\in S}}{a}_k.$

The probability of selection of unit $k$ is supposed to be proportional to their contacts

$p_k=L_k/N\overline{L},$

where

$L_k={\displaystyle \sum _{j\in U}}{\lambda }_{j,k},\overline{L}={\displaystyle \frac{L}{N}},\text{and}L={\displaystyle \sum _{k\in U}}{L}_k$

where ${\lambda }_{j,k}$ is the link (0,1) variable between individuals $j$ and $k$ . Let

${\widehat{p}}_k=L_k/N\widehat{\overline{L}}$

be the Hansen Hurwitz (HH) [18] estimate of $p_k$ , where

$\widehat{\overline{L}}=\left({\displaystyle \sum _{\mathrm{\ell }\in S}}{\displaystyle \frac{a_{\mathrm{\ell }}{L}_{\mathrm{\ell }}}{m{p}_{\mathrm{\ell }}}}\right)/\left({\displaystyle \sum _{\mathrm{\ell }\in S}}{a}_{\mathrm{\ell }}/m{p}_{\mathrm{\ell }}\right)$

is the HH estimate ratio of the average number of contacts in the population. The numerator on the right-hand side of the previous equality is given by:

${\displaystyle \sum _{\mathrm{\ell }\in S}}{\displaystyle \frac{a_{\mathrm{\ell }}{L}_{\mathrm{\ell }}}{m{p}_{\mathrm{\ell }}}}={\displaystyle \sum _{\mathrm{\ell }\in S}}{\displaystyle \frac{a_{\mathrm{\ell }}{L}_{\mathrm{\ell }}}{m{L}_{\mathrm{\ell }}}}L$ $={\displaystyle \frac{1}{m}}L{\displaystyle \sum _{\mathrm{\ell }\in S}}{a}_{\mathrm{\ell }}=L.$

Figure 3.

Example of construction of the totals $Y_{S_0,R}$ .

The denominator can be expressed as

${\displaystyle \sum _{\mathrm{\ell }\in S}}{\displaystyle \frac{a_{\mathrm{\ell }}}{m{p}_{\mathrm{\ell }}}}={\displaystyle \sum _{\mathrm{\ell }\in S}}{\displaystyle \frac{a_{\mathrm{\ell }}L}{m{L}_{\mathrm{\ell }}}}$

$\widehat{\overline{L}}$ and ${\widehat{p}}_k$ are, therefore, given by:

$\widehat{\overline{L}}={\displaystyle \frac{L}{{\sum }_{\mathrm{\ell }\in S}\frac{a_{\mathrm{\ell }}L}{m{L}_{\mathrm{\ell }}}}}={\displaystyle \frac{m}{{\sum }_{\mathrm{\ell }\in S}{a}_{\mathrm{\ell }}\frac{1}{L_{\mathrm{\ell }}}}},{\widehat{p}}_k=L_k/N{\displaystyle \frac{m}{{\sum }_{\mathrm{\ell }\in S}{a}_{\mathrm{\ell }}\frac{1}{L_{\mathrm{\ell }}}}}.$

The VH estimator of $Y$ is

${\widehat{Y}}_{\text{𝑉𝐻}}={\displaystyle \sum _{k\in S}}{\displaystyle \frac{a_k{y}_k}{m{\widehat{p}}_k}}={\displaystyle \sum _{k\in S}}{y}_k{w}_{\text{𝑉𝐻},k}$(2)

where

$w_{\text{𝑉𝐻},k}={\displaystyle \frac{a_k}{m{\widehat{p}}_k}}.$(3)

is the sample weight assigned to unit $k$ .

3. Totals of interests and a formalisation of the RDS as a particular case of indirect sampling

The RDS can be formalised as an indirect sampling scheme. In this type of sampling, there is an initial $U^A$ population of $N^A$ units from which the research starts, and a $U^B$ population of $N^B$ units that constitute the study’s target population.

In our case, $U^B\equiv U$ means that it coincides with the target population $U$ , which implies $N^B=N$ .

The specific unit $j$ of the initial population $U^A$ consists of the unit $j$ itself and all its contacts. Let $j$ and $k$ be the labels identifying the population units in $U^A$ and $U^B$ , respectively. Let ${\lambda }_{j,k}$ be the link variable between units $j\in U^A$ and $k\in U^B$ , where ${\lambda }_{j,k}=1$ if $j$ is directly linked to $k$ . We have ${\lambda }_{j,k}={\lambda }_{k,j}$ and ${\lambda }_{j,j}=1.$ Let

${\overline{y}}_j^A=-{\displaystyle \sum _{k\in U^B}}{\displaystyle \frac{{\lambda }_{j,k}}{L_k^B}}{y}_k$(4)

be the value of the characteristic of interest for unit $j\in U^A$ , where

$L_k^B={\displaystyle \sum _{j\in U^A}}{\lambda }_{j,k}.$

In the initial population $U^A$ , each unit in contact with unit $j$ contributes to the $y$ value of that unit in a weighted manner, where the reciprocal of the total number of contacts gives the weight.

The two populations $U^A$ and $U^B$ have the same numbers of units: $N^A=N^B=N$ . We have

$N^A={\displaystyle \sum _{j\in U^A}}{\displaystyle \sum _{k\in U^B}}{\displaystyle \frac{{\lambda }_{j,k}}{L_k^B}}$ $={\displaystyle \sum _{k\in U^B}}{\displaystyle \sum _{j\in U^A}}{\displaystyle \frac{{\lambda }_{j,k}}{L_k^B}}$(5) $={\displaystyle \sum _{k\in U^B}}1=N^B=N.$

Moreover, the target parameter may also be expressed as the sum of the ${\overline{y}}_j^A$ values over the population $U^A$ :

$Y={\displaystyle \sum _{j\in U^A}}{\overline{y}}_j^A.$(6)

Indeed, it is

$Y={\displaystyle \sum _{j\in U^A}}{\overline{y}}_j^A={\displaystyle \sum _{j\in U^A}}{\displaystyle \sum _{k\in U^B}}{\displaystyle \frac{{\lambda }_{j,k}}{L_k^B}}{y}_k={\displaystyle \sum _{k\in U^B}}{\displaystyle \sum _{j\in U^A}}{\displaystyle \frac{{\lambda }_{j,k}}{L_k^B}}{y}_k={\displaystyle \sum _{k\in U^B}}{y}_k.$(7)

In addition to the total $Y$ , another total that plays a crucial role in the RDS methodology is the aggregate, $Y_{S_{0,R}}$ . That is, the total of the variable $y$ where starting from the sample $S_0$ (which constitutes step 0 of the process), additional units are considered in subsequent steps through all their contacts. The total considers this aggregate after the step of this process. We can consider this as a search process on a graph.

Figure 3 illustrates this process starting from sample $S_0$ , which includes only unit $j$ .

To clarify how this total can be constructed, let us compute it step by step. To distinguish among the search processes of the different steps, let $j_r$ (with $j_r\in U^A$ and $r=0,1,\mathrm{\dots },R-1)$ be the subscript of the population unit involved in step $r$ of this search process on a graph.

At step 0, the total $Y_{S_0,0}$ is simply the sum of the values $y_j$ for the units included in the sample $S_0$ :

$Y_{S_0,0}={\displaystyle \sum _{j_0\in S_0}}{y}_j.$(8)

In step 1, we observe the links of units in $S_0$ . We compute $Y_{S_0,1}$ by adding to the total $Y_{S_0,0}$ the values of $y$ of the new units individuated in the search process starting from the units in $S_0$ . We can formalise this process as

$Y_{S_0,1}={\displaystyle \sum _{j_0\in S_0}}{\displaystyle \sum _{k\in U^B}}{y}_k{\displaystyle \frac{{\lambda }_{j_0,k}}{L_k^{S_0}}},$(9)

where $L_k^{S_0}={\sum }_{j_0\in S_0}{\lambda }_{j_0,k}$ .

To clarify the notation, we note that the $j_0$ unit is also one of the $k$ units. We also note that $L_k^{S_0}$ is the total number of links to unit k (identified in step 1 of the process) that can be computed starting from the units in $S_0$ .

We can reverse the order of the sums and formulate as

$Y_{S_0,1}={\displaystyle \sum _{k\in U^B}}{y}_k{\displaystyle \sum _{j_0\in S_0}}{\displaystyle \frac{{\lambda }_{j_0,k}}{L_k^{S_0}}}$(10)

where $L_k^{S_0}={\sum }_{j_0\in S_0}{\lambda }_{j_0,k}$ .

We note that Eq. (10) avoids the multiple counting of a unit linked to different elements of the initial sample, as illustrated in Fig. 4, where unit $b$ is connected to both units $j$ and a of $S_0$ . Indeed, it is

${\displaystyle \sum _{j_0\in S_0}}{\displaystyle \frac{{\lambda }_{j_0,b}}{L_b^{S_0}}}={\displaystyle \frac{{\lambda }_{a,b}}{L_b^{S_0}}}+{\displaystyle \frac{{\lambda }_{b,b}}{L_b^{S_0}}}+{\displaystyle \frac{{\lambda }_{j,b}}{L_b^{S_0}}}$ $={\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{3}}=1.$

Figure 4.

Example of Steps 0 and 1.

At step 2, we have

$Y_{S_0,2}={\displaystyle \sum _{j_0\in S_0}}{\displaystyle \sum _{j_1\in U^A}}{\displaystyle \sum _{k\in U^B}}{y}_k{\displaystyle \frac{{\lambda }_{j_0,j_1}}{L_{j_1}^{S_0}}}{\displaystyle \frac{{\lambda }_{j_1,k}}{L_k^B}}.$(11)

We note that the last summation is over $U^B$ , i.e., the target population. The middle summation is on $U^A$ , i.e., on the initial population. The first summation is limited to the initial sample from which the research starts. This kind of organisation of the order of summations also appears in the following formulas. The last summation is always on $U^B$ , and the first is on $S_0$ . In contrast, the intermediate summations are always on the starting population $U^A$ .

Reversing the order of the summations, we have

$Y_{S_0,2}={\displaystyle \sum _{k\in U^B}}{y}_k{\displaystyle \sum _{j_0\in S_0}}{\displaystyle \sum _{j_1\in U^A}}{\displaystyle \frac{{\lambda }_{j_0,j_1}}{L_{j_1}^{S_0}}}{\displaystyle \frac{{\lambda }_{j_1,k}}{L_k^B}}.$(12)

Let us note that in this case, the unit $k\in U^B$ is counted only once, avoiding the multiple counting of a unit linked to different elements of the initial sample and its links. This is illustrated in Appendix 1.

$\mathrm{\dots }$

Continuing recursively the above process, at step $R$ we have

$Y_{S_0,R}={\displaystyle \sum _{j_0\in S_0}}{\displaystyle \sum _{j_1\in U^A}}\mathrm{\dots }{\displaystyle \sum _{j_{R-1}\in U^A}}{\displaystyle \sum _{k\in U^B}}{y}_k{\displaystyle \frac{{\lambda }_{j_0,j_1}}{L_{j_1}^{S_0}}}{\displaystyle \frac{{\lambda }_{j_1,j_2}}{L_{j_2}^B}}\times \mathrm{\dots }\times {\displaystyle \frac{{\lambda }_{j_{R-2},j_{R-1}}}{L_{j_{R-1}}^B}}{\displaystyle \frac{{\lambda }_{j_{R-1},k}}{L_k^B}}={\displaystyle \sum _{k\in U^B}}{y}_k{\displaystyle \sum _{j_0\in S_0}}{\displaystyle \sum _{j_1\in U^A}}\mathrm{\dots }{\displaystyle \sum _{j_{R-1}\in U^A}}{\displaystyle \frac{{\lambda }_{j_0,j_1}}{L_{j_1}^{S_0}}}{\displaystyle \frac{{\lambda }_{j_1,j_2}}{L_{j_2}^B}}\times \mathrm{\cdots }\times {\displaystyle \frac{{\lambda }_{j_{R-2},j_{R-1}}}{L_{j_{R-1}}^B}}{\displaystyle \frac{{\lambda }_{j_{R-1},k}}{L_k^B}}.$(13)

Even in this case, we avoid the multiple counting of a unit linked to different elements collected in the $R$ steps of the RDS process.

Based on the above expressions, we note the following.

Each step of network sampling can be formalised as an indirect sampling mechanism.

In a given step of the RDS process, the participants from which the search starts constitute the source list $U^A$ , and their links are the target population $U^B$ .

In the subsequent step, the people found in the target population $U^B$ become the people of the initial population $U^A$ , from which a new search starts.

This switch in the role of the sample participants, from the target population $U^B$ to the initial population $U^A$ of the next step, occurs at each wave of the RDS search chain.

Remark 1. A path in graph theory is a finite or infinite sequence of edges that joins a sequence of vertices. We note that if $R$ is greater than the maximum of the minimum paths between any pair of nodes in each cluster of the units of $S_0$ , then the following implications follow.

$Y_{S_0,R}$ represents the total of the $y$ variable related to all the groups to which the elements of $S_0$ belong.

if $S_0$ does not include all clusters characterising the population of interest, then $Y_{S_0,R}<Y$ .

We obtain these important outcomes illustrated in remarks 2 and 3 below by interpreting the result of Remark 1 in an alternative way.

Remark 2. $Y_{S_0,R}=Y$ if

$R$ is greater than the maximum of the minimum paths between any pair of nodes in each cluster of the units of $S_0$ .

The $S_0$ sample must include people from all clusters with people from the target population.

Remark 3. We consider the case of two interconnected units in which they know each other. However, these units may belong to two separate clusters if the RDS search rules call for stopping the search if a connection is identified with a person living in a geographic location distinct from that of the original contact. In this case, to ensure the equality $Y_{S_0,R}=Y$ , in addition to the two conditions in Remark 2, there is the additional condition that the $S_0$ sample must include people from all geographic locations with people in the target population.

4. Sampling the RDS research chain

Unlike in the previous section, in the RDS process, not all links of a unit kept in the process are observed, but only a randomly selected sample of them is observed. The RDS process starts at step $r=0$ , with sample $S_0$ (which may be randomly selected or not), and in subsequent steps $r=1,2,\mathrm{\dots },R$ , we form samples $S_1\subset S_2\subset \mathrm{\dots }\subset S_r\subset \mathrm{\dots }\subset S_R$ , each incorporating the sample from the previous step. To illustrate the formation of the generic sample $S_{r+1}$ , we introduce additional symbology below. Let us consider the $j_r\in S_r$ unit and denote by

$L_{j_r}^{S_r}={\displaystyle \sum _{k\in U^B}}{\lambda }_{j_r,k}{\delta }_k\left(S_r\right)={\displaystyle \sum _{k\in S_r}}{\lambda }_{j_r,k}$

the total number of contacts of the unit that have been selected in sample $S_r$ , where ${\delta }_k(A)=1$ if unit $k$ belongs to set $A$ and ${\delta }_k(A)=0$ otherwise. For the same unit $j_r$ , let

$L_{j_r}^{C_r}=L_{j_r}^B-L_{j_r}^{S_r}$

be the number of contacts that have not been selected in sample $S_r$ , and can be selected in sample $S_{r+1}$ . For each unit $j_r$ included in the $S_r$ sample, we select, independently of the other units included in $S_r,m_{j_r}+1$ units for the $S_{r+1}$ sample where $m_{j_r}=\text{𝑀𝑖𝑛}(m,L_{j_r}^{C_r})$ , being $m$ is a fixed number (e.g., $m=2$ or 3) that remains unchanged in the different steps of the RDS. Of these $m_{j_r}+1$ units, one is the $j_r$ unit itself, and the other $m_{j_r}$ units are selected with a simple random sampling without replacement (SRSWOR) out from the $L_{j_r}^{C_r}$ units. The conditional probability that unit $j_{r+1}$ is selected in sample $S_{r+1}$ , given $j_r\in S_r$

${\tau }_{j_{r+1}\mid j_r\in S_r}=\{\begin{array}{cc}1\hfill & \text{if}({\lambda }_{j_r,j_{r+1}}=1)\text{and}\hfill \\ & \left[{\delta }_{j_{r+1}}\left(S_r\right)=1\right]\hfill \\ {\displaystyle \frac{m_{j_r}}{L_{j_r}^{C_r}}}\hfill & \text{if}({\lambda }_{j_r,j_{r+1}}=1)\text{and}\hfill \\ & \left[{\delta }_{j_{r+1}}\left(S_r\right)=0\right]\hfill \\ 0\hfill & \text{if}{\lambda }_{j_r,j_{r+1}}=0\hfill \end{array}.$

Figure 5.

Example of the formation of sample $S_1$ from sample $S_0$ with $m=$ 2.

Remark 4 on the feasibility of the selection process. The illustrated selection process avoids the dead-end loops typical of graph sampling. However, to make it feasible, it is essential to know the $L_{j_r}^{C_r}$ quantity, obtained as difference of two quantities, $L_{j_r}^B$ and $L_{j_r}^{S_r}$ . The value of $L_{j_r}^B$ can be requested directly from respondent $j_r$ . Remembering that the relationships explored in RDS have the character of reciprocity, $L_j^B$ corresponds to the total number of people who unit $j_r$ knows and can point to in turn. Operationally, the $L_{jr}^{S_r}$ quantity can be obtained in alternative ways. Suppose nonidentifiable but unique information about contacts of units included in the $S_r$ sample [13] is available in the data-collection APP used by the interviewer. In that case, a specific software application can be launched that identifies units not included in $S_r$ and proceeds to select $m_{jr}$ units to be included in sample $S_{r+1}$ randomly. Alternatively, the same software application can be run by the study centre (see Section 2 above) that supports the survey operations, and the results can be reported and provided in real-time to the interviewer who makes the $S_{r+1}$ sample selection. Depending on the specific conditions of the survey, other feasible operational mechanisms can be defined.

Remark 5 on the research chain for a subpopulation. We consider the case of constructing the RDS sampling search chain only on the units of a subpopulation, for example, only on the people of a particular class-age. We denote by $x_j$ a dichotomous variable that takes value 1 if unit $j$ belongs to the subpopulation and takes value 0 otherwise. In this case, the link variables are defined as

${\lambda }_{\left(x\right)j,k}={\lambda }_{j,k}{x}_j{x}_k.$

The values $L_k^B$ are modified accordingly.

5. Estimators

Next, we present three estimators. The first assumes a random selection of the $S_0$ sample, and the second adopts the traditional RDS methodology while considering a non-probabilistic $S_0$ sample. The third estimator is developed under a capture-recapture approach [12] while allowing for the smoothing of the coverage problems that may affect both of the first two estimators.

5.1 $S_0$ selected with a random sampling

Let us suppose a random sample $S_0$ of fixed size $n_0$ is selected from $U^A$ without replacement and with inclusion probabilities ${\pi }_{j0}$ , where

${\pi }_{j_0}>0\text{for}{j}_0=1,\mathrm{\dots },N^A\text{and}$ ${\displaystyle \sum _{j_0\in U^A}}{\pi }_{j_0}=n_0.$(14)

To facilitate the understanding of the calculation method, we construct the estimator step by step. In each step, we obtain a sampling unbiased estimate of the total $Y$ . However, as the steps of the RDS process progress, the estimate is based on a more significant number of observations.

At step 0, we have the classical HT estimator:

${\widehat{Y}}_0={\displaystyle \sum _{j_0\in S_0}}{y}_{j_0}{\displaystyle \frac{1}{\pi j_0}}={\displaystyle \sum _{j_0\in S_0}}{y}_{j_0}{w}_{j_0},$(15)

where $w_{j_0}=1/{\pi }_{j_0}$ is the sampling weight.

In Step 1, we have

${\widehat{Y}}_1={\displaystyle \sum _{j_0\in S_0}}{\displaystyle \sum _{k\in S_1}}{y}_k{\displaystyle \frac{{\lambda }_{j_0,k}}{L_k^B}}{\displaystyle \frac{1}{{\pi }_{j_0}}}{\displaystyle \frac{1}{{\tau }_{k|j_0\in S_0}}}.$(16)

As illustrated in Appendix 2, denoting $E(\cdot )$ , the sampling expectation operator, we have $E({\widehat{Y}}_1)=Y$ , meaning that ${\widehat{Y}}_1$ is a sampling unbiased estimate of $Y$ . Reversing the order of sums, we can express ${\widehat{Y}}_1$ in the classical weighted form as

${\widehat{Y}}_1={\displaystyle \sum _{k\in S_1}}{y}_k{w}_{k_1}$(17)

where

$w_{k_1}={\displaystyle \sum _{j_0\in S_0}}{\displaystyle \frac{{\lambda }_{j_0,k}}{L_k^B}}{\displaystyle \frac{1}{{\pi }_{j_0}}}{\displaystyle \frac{1}{{\tau }_{k|j_0\in S_0}}}.$

Estimator ${\widehat{Y}}_1$ can also be formulated referring to Eq. (6) as:

${\widehat{Y}}_1={\displaystyle \sum _{j_0\in S_0}}{\displaystyle \frac{1}{{\pi }_{j_0}}}{\widehat{\overline{y}}}_{j_0}^A$(18)

where ${\widehat{\overline{y}}}_{j_0}^A={\sum }_{k\in S_1}{y}_k\frac{{\lambda }_{j_0,k}}{L_k^B}\frac{1}{{\tau }_{k\mid j_0\in S_0}}$ is the unbiased estimate of ${\overline{y}}_{j_0}^A$

In Step 2, taking Eqs (16), (17), and (18) of step 1 as a reference, the unbiased estimator ${\widehat{Y}}_2$ can be expressed according to these three alternative ways

${\widehat{Y}}_2={\displaystyle \sum _{j_0\in S_0}}{\displaystyle \sum _{j_1\in S_1}}{\displaystyle \sum _{k\in S_2}}{y}_k{\displaystyle \frac{{\lambda }_{j_0,j_1}}{L_{j_1}^B}}{\displaystyle \frac{{\lambda }_{j_1,k}}{L_k^B}}{\displaystyle \frac{1}{{\pi }_{j_0}}}{\displaystyle \frac{1}{{\tau }_{j_1\mid j_0\in S_0}}}{\displaystyle \frac{1}{{\tau }_{k\mid j_1\in S_1}}},$ ${\widehat{Y}}_2={\displaystyle \sum _{k\in S_2}}{y}_k{w}_{k_2},$ ${\widehat{Y}}_2={\displaystyle \sum _{j_0\in S_0}}{\displaystyle \sum _{j_1\in S_1}}{\displaystyle \frac{1}{{\pi }_{j_0}}}{\displaystyle \frac{{\lambda }_{j_0,j_1}}{L_{j_1}^B}}{\displaystyle \frac{1}{{\tau }_{j_1\mid j_0\in S_0}}}{\widehat{\overline{y}}}_{j_1}^A,$

where

$w_{k_2}={\displaystyle \sum _{j_0\in S_0}}{\displaystyle \sum _{j_1\in S_1}}{\displaystyle \frac{{\lambda }_{j_0,j_1}}{L_{j_1}^B}}{\displaystyle \frac{{\lambda }_{j_1,k}}{L_k^B}}{\displaystyle \frac{1}{{\pi }_{j_0}}}{\displaystyle \frac{1}{{\tau }_{j_1\mid j_0\in S_0}}}{\displaystyle \frac{1}{{\tau }_{k\mid j_1\in S_1}}},$

and ${\widehat{\overline{y}}}_{j_1}^A={\sum }_{k\in S_2}{y}_k\frac{{\lambda }_{j_1,k}}{L_k^B}\frac{1}{{\tau }_{k\mid j_1\in S_1}}$ is the unbiased estimate of ${\overline{y}}_{j_1}^A$ . We have $E\left({\widehat{Y}}_2\right)=Y$ .

Recursively using the previous procedure, at the ultimate step $R$ of the RDS process, we have:

${\widehat{Y}}_R={\displaystyle \sum _{j_0\in S_0}}{\displaystyle \sum _{j_1\in S_1}}\mathrm{\dots }{\displaystyle \sum _{j_{R-1}\in S_{R-1}}}{\displaystyle \sum _{k\in S_R}}{y}_k\times {\displaystyle \frac{{\lambda }_{j_0,j_1}}{L_{j_1}^B}}\times \mathrm{\dots }\times {\displaystyle \frac{{\lambda }_{j_{R-1},k}}{L_k^B}}{\displaystyle \frac{1}{{\pi }_{j_0}}}{\displaystyle \frac{1}{{\tau }_{j_1\mid j_0\in S_0}}}\times \mathrm{\dots }\times {\displaystyle \frac{1}{{\tau }_{k\mid j_{R-1}\in S_{R-1}}}},$ ${\widehat{Y}}_R={\displaystyle \sum _{k\in S_R}}{y}_k{w}_{k_R,}$ ${\widehat{Y}}_R={\displaystyle \sum _{j_0\in S_0}}{\displaystyle \sum _{j_1\in S_1}}\mathrm{\dots }{\displaystyle {\sum }_{j_{R-1}\in S_{R-1}}}{\displaystyle \frac{{\lambda }_{j_0,j_1}}{L_{j_1}^B}}\times \mathrm{\dots }\times {\displaystyle \frac{{\lambda }_{j_{R-1},k}}{L_k^B}}{\displaystyle \frac{1}{{\pi }_{j_0}}}{\displaystyle \frac{1}{{\tau }_{j_1\mid j_0\in S_0}}}\times \mathrm{\dots }\times {\displaystyle \frac{1}{{\tau }_{j_{R-1}\mid j_{R-2}\in S_{R-2}}}}{\widehat{\overline{y}}}_{j_{R-1}}^A,$

where

$w_{k_R}={\displaystyle \sum _{j_0\in S_0}}{\displaystyle \sum _{j_1\in S_1}}\mathrm{\dots }{\displaystyle \sum _{j_{R-1}\in S_{R-1}}}{\displaystyle \frac{{\lambda }_{j_0,j_1}}{L_{j_1}^B}}\times \mathrm{\dots }\times {\displaystyle \frac{{\lambda }_{j_{R-1},k}}{L_k^B}}{\displaystyle \frac{1}{{\pi }_{j_0}}}{\displaystyle \frac{1}{{\tau }_{j_1\mid j_0\in S_0}}}\times \mathrm{\dots }\times {\displaystyle \frac{1}{{\tau }_{k\mid j_{R-1}\in S_{R-1}}}}$

and

${\widehat{\overline{y}}}_{j_{R-1}}^A={\displaystyle \sum _{k\in S_{R-1}}}{y}_k{\displaystyle \frac{{\lambda }_{j_{R-1},k}}{L_k^B}}{\displaystyle \frac{1}{{\tau }_{k\mid j_{R-1}\in S_{R-1}}}}.$

In Appendix 2, we see $E\left({\widehat{Y}}_R\right)=Y$

Remark 6 on estimating the sampling variance. We can approximate the RDS design with multistage sampling with replacement, where each step may be viewed as a specific sampling stage, and the replacement refers to a single unit. In that way, we may derive an estimate of the sampling variance [19][(Formula 11.35)];

$v\left({\widehat{Y}}_R\right)={\displaystyle \frac{1}{n_0\left(n_0-1\right)}}{\displaystyle \sum _{j_0\in S_0}}\left({\displaystyle \frac{1}{z_{j_0}}}{\widehat{Y}}_{j_0}-{\widehat{Y}}_R\right)$

where

${\widehat{Y}}_{j_0}={\displaystyle \sum _{j_1\in S_1}}\mathrm{\dots }{\displaystyle \sum _{j_{R-1}\in S_{R-1}}}{\displaystyle \sum _{k\in S_R}}{y}_k{\displaystyle \frac{{\lambda }_{j_0,j_1}}{L_{j_1}^B}}\times \mathrm{\dots }\times {\displaystyle \frac{{\lambda }_{j_{R-1},k}}{L_k^B}}{\displaystyle \frac{1}{{\pi }_{j_0}}}{\displaystyle \frac{1}{{\tau }_{j_1\mid j_0\in S_0}}}\times \mathrm{\dots }\times {\displaystyle \frac{1}{{\tau }_{k\mid j_{R-1}\in S_{R-1}}}},$

and $z_{j_0}={\pi }_{j_0}/n_0$ .

Remark 7 on the estimator for a subpopulation. As illustrated in Remark 5, the link variables are defined as ${\lambda }_{(x)j,k}=$ ${\lambda }_{j,k}{x}_j{x}_k$ , and the variables $L_k^B$ are modified accordingly. Moreover, the target variable $y_k$ is modified as $y_{(x)k}=y_k{x}_k$ .

Remark 8 on type of estimator. Considering the above expressions, we can see how the ${\widehat{Y}}_R$ estimator can be seen as a particular case of the generalised weight share method (GWSM) estimator.

Remark 9 on the starting sampling. The sampling design should maximise the number of observed individuals of the target population in the sample $S_0$ by adopting proper choices in the first and ultimate stages (or phases) of the sampling process. First-stage selection tends to oversample the areas where the researchers have some a priori information of a high concentration of the target population. Final-stage sampling tends to oversample the target people by modelling the inclusion probabilities on variables predictive of the phenomenon available in the sampling frames adopted to select the final-stage units.

5.2 $S_0$ selected with a nonrandom sampling

If the sample $S_0$ is selected nonrandomly, we can obtain a nonbiased estimate only of the total $Y_{S_0,R}$ . We illustrate this case by referring only to step $R$ and the weighted form of the estimator. An unbiased (see Appendix 2) estimator of the total $Y_{S_0,R}$ is

${\widehat{Y}}_{S_0,R}={\displaystyle \sum _{k\in S_R}}{y}_k{w}_{k_R,S_0}$

where

$w_{k_R,S_0}={\displaystyle \sum _{j_0\in S_0}}{\displaystyle \sum _{j_1\in S_1}}\mathrm{\dots }{\displaystyle \sum _{j_{R-1}\in S_{R-1}}}{\displaystyle \frac{{\lambda }_{j_0,j_1}}{L_{j_1}^{S_1}}}\times \mathrm{\cdots }\times {\displaystyle \frac{{\lambda }_{j_{R-1},k}}{L_k^B}}{\displaystyle \frac{1}{{\tau }_{j_1\mid j_0\in S_0}}}\times \mathrm{\dots }\times {\displaystyle \frac{1}{{\tau }_{k\mid j_{R-1}\in S_{R-1}}}}.$

We note that the formulation of ${\widehat{Y}}_{S_0,R}$ is similar to that of the ${\widehat{Y}}_R$ estimator, except that in ${\widehat{Y}}_{S_0,R}$ , the first weighting factor $\left(1/{\pi }_{j_0}\right)$ of ${\widehat{Y}}_R$ is equal to 1 , and $L_{j_1}^{S_0}$ replaces the $L_{j_1}^B$ factor.

Remark 10 on type of estimator. We can straightforwardly see how the ${\widehat{Y}}_{S_0,R}$ estimator can be seen as a particular case of the GWSM estimator.

Remark 11 on the starting sampling. To ensure that estimator ${\widehat{Y}}_{S_0,R}$ is an unbiased estimate of the total $Y$ , i.e., that condition $E({\widehat{Y}}_{S_0,R})=Y$ is met, the initial sample $S_0$ should respect the three conditions illustrated in remarks 2 and 3.

5.3 Generalised capture-recapture estimator

Even if the $S_0$ sample is randomly selected, the first estimator ${\widehat{Y}}_R$ may be biased. Indeed, undercoverage may occur if respondents do not trust the interviewers and tend to hide their status.

Likewise, if the $S_0$ sample is nonrandomly chosen, even the second estimator ${\widehat{Y}}_{S_0,R}$ can undercover the total $Y$ if the following conditions are not met: (i) $R$ is greater than the maximum of the minimum paths between any pair of nodes in each cluster of the units of $S_0$ , and (ii) $S_0$ does not include all clusters characterising the population of interest.

The generalised capture-recapture (CReG) estimator allows us to overcome the abovementioned undercoverage by leveraging a capture-recapture perspective. A comprehensive treatment of this estimator and how it mitigates undercoverage deserves much more space in this article than can be devoted. The interested reader can undoubtedly look to the extensive work reported in [9].

Let us consider two independent surveys based on the RDS methodology, and we suppose they are articulated in $R$ steps. The first starts from an initial random sampling, while the second starts from a nonrandom sample. Furthermore, we suppose that the two sample selections are independent. The CReG estimator of $Y$ may be expressed as

${\widehat{Y}}_{\mathit{\text{CReG}}}={\displaystyle \frac{{\widehat{Y}}_R{\widehat{Y}}_{S_0,R}}{{\widehat{Y}}_{\text{𝑖𝑛𝑡𝑒𝑟𝑠𝑒𝑐𝑡},R}}}$

where

${\widehat{Y}}_{\text{𝑖𝑛𝑡𝑒𝑟𝑠𝑒𝑐𝑡},R}={\displaystyle \sum _{k\in S_{R,\mathit{\text{intersect}}}}}{w}_{k_R}{w}_{k_R,S_0}{y}_k$

where $S_{R,\mathit{\text{intersect}}}$ is the intersection sample between the two samples that are generated from random and nonrandom sampling after $R$ steps.

A useful approach is applying the estimator CReG on the two nonrandom starting samples but with a different mechanism of undercoverage of the two respondent groups.

6. Conclusions

The disaggregation of data for SDG indicators on hidden or hard-to-count population groups presents several critical issues that are difficult to overcome to produce reliable official statistics at the national and international levels. In this context, it is impossible to estimate the characteristics of these groups through models as in other situations.

Considering, for example, indigenous populations, data availability varies widely from country to country. Few countries provide up-to-date and high-quality data. Many other countries have only data that are scattered over time. Or the data they provide is not supported by a sufficiently robust methodology, both for precisely identifying the subpopulation group of interest and for the sampling technique adopted.

Given the current context of producing official statistics on the subject, it is unrealistic that this lack of data on indigenous peoples is going to improve soon.

Therefore, it becomes necessary to define and apply an implementation program that can improve this situation relatively quickly.

This implementation program should be based on the following pillars:

The first pillar is to develop a valuable data collection strategy for sample surveys that is capable of measuring the number of people belonging to the indigenous population in a given area. This strategy needs to cover various aspects, including the formulation of the questionnaire for identifying persons belonging to the indigenous population and the characteristics of special sampling techniques that can maximise the efficiency of surveys aimed at obtaining data on these hidden or hard-to-measure population groups. Specifically, the data collection strategy consists of technical manuals, open software modules, ad hoc training materials, etc. In short, anything that enables or helps conduct surveys or specific survey modules to estimate the size of indigenous populations. Regarding the questionnaire, the data collection strategy should develop a set of standard questions on the key characteristics of the specific population of interest and not adopt generic questions on whether specific individuals belong to an indigenous population group.

The second pillar is to adapt the data collection strategy to ongoing survey programs. For example, the indigenous module may be applied to a large-scale national survey conducted by a national statistics office. Another example can be to promote the application of the indigenous sampling module to international surveys, such as the World Bank’s LSMS survey. Regarding the sampling aspects, it is helpful to consider an overall sampling strategy to maximise the number of observed individuals of the target population in the sample by combining the first and final sampling stages. First-stage methods should tend to oversample the areas where the researchers have some a priori information on the geographical concentration of the target population. Final-stage sampling should oversample the target population by adopting strategies, such as respondent-driven sampling (RDS), that leverage the existing hidden relationships among the individuals of the target populations.

The third pillar is adopting estimation methods that allow unbiased estimates of phenomena of interest in target populations. In this paper, we have considered the RDS, widely used for observing hidden or rare populations but which lacks an estimation methodology that is sufficiently robust to accommodate the varying conditions under which it is applied. We have proposed three unbiased estimators. The first assumes a random selection of the starting sample, and the second considers a nonprobabilistic starting sample. The third estimator is developed under a capture-recapture approach and allows for smoothing of the coverage problems that may affect the estimators. We have studied their sampling expectation and have indicated the conditions that may be fulfilled to guarantee their unbiasedness concerning the target totals.