Improved Inference for Respondent-Driven Sampling Data with Application to HIV Prevalence Estimation

Researchers and government officials are often interested in characteristics of human populations for which there are no practicable sampling frames for direct sampling. In some such hidden populations, members are connected through social networks. A common approach is to collect a sample using a variant of link-tracing (Thompson, 2006a, b), such as a snowball sample (Goodman, 1961), where subsequent sample members are selected based on their relations with previously sampled individuals. When the initial sample is not a probability sample, this approach does not result in a probability sample. However, most available alternatives, (Muhib et al., 2001; Peterson et al., 2008; Watters and Biernacki, 1989) also fail to produce a probability sample of the population.

Respondent-Driven Sampling (RDS, introduced by Heckathorn 1997, 2002, see also Salganik and Heckathorn, 2004, Volz and Heckathorn 2008) is a recently introduced variant of link-tracing sampling which increases the ease of sampling and produces samples which, it is argued, approach probability samples as sampling progresses. The lack of satisfactory alternatives has spurned a strong demand for RDS (Johnston et al., 2008; Lansky et al., 2007).

RDS presents two main innovations: a sampling design and a corresponding approach to estimation. The former is highly effective in many settings (Lansky et al., 2007); the respondent-driven design relies on the respondents at each wave to select the next wave by distributing uniquely identified coupons to others in the target population, who can choose to return the coupons to enroll in the study. Thus, the sampling exploits the network of social relations while also avoiding the confidentiality concerns associated with recording the names of contacts.

The key innovation for estimation is that through many waves of sampling, the dependence of the final sample on the initial sample is reduced, allowing researchers more confidence in making approximate probability statements about the resulting samples. This insight allows for statistical inference in settings where the initial sample is typically selected by a convenience mechanism. Although current inference is likely superior to the alternative non-probability methods, existing methods are sensitive to deviations from many assumptions (Goel and Salganik, 2009; Gile and Handcock, 2010; Neely, 2009). This paper offers a modification of the existing theoretical formulation of respondent-driven sampling, and corresponding inference to address what we find to be a serious conceptual weakness of existing work: the known inaccuracy of the with-replacement approximation to the sampling process.

In the next section, we begin by introducing current RDS estimation, particularly the estimator introduced by Volz and Heckathorn (2008), and illustrate the sensitivity of this estimator to the with-replacement sampling assumption in cases of substantial sample fractions. In Sections 3 and 4, we then introduce a new model for RDS sampling based on successive sampling, and introduce a new estimator based on that model. In Section 5, we use a simulation study to illustrate the superior performance of the new estimator. We also include sensitivity analyses concerning inaccurate estimation of the size of the target population and the characteristics of the initial sample. In Section 6, we apply our estimator to data collected in three populations of drug users and men who have sex with men. We conclude with a broader discussion of the method and its limitations.

The basic ideas underlying estimation from RDS data are clever and important. They allow for something like valid statistical inference, in a sampling setting where the target population cannot be effectively reached using a traditional sampling frame. The original article, (Heckathorn, 1997) made very strong assumptions about the sampling procedure so as to assume that the sample proportions were representative of the population proportions. Salganik and Heckathorn (2004) introduced a Markov chain argument for population mixing, and proposed an estimator based on equating the number of cross-relations between pairs of sub-populations of interest, based on the referral patterns of each group. This estimator is currently in wide use, and is implemented in the standard RDS analysis software (Volz et al., 2007). Volz and Heckathorn (2008) connect RDS estimation to mainstream survey sampling through the use of a generalized Horvitz-Thompson estimator form. This estimator relies on the estimation of the inclusion probabilities of the sampled units, $\pi_{i}$. Based on an argument for treating the sample as independent draws from the stationary distribution of a random walk on the nodes of an undirected graph, Volz and Heckathorn (2008) approximate the sampling probabilities as proportional to nodal degree, $d_{i}$, or number of incident edges in the graph. This estimator avoids the problem of potentially unknown population size $N$ by using the generalized Horvitz-Thompson form, normalizing by an estimator of $N$ as follows:

$\displaystyle\hat{\mu}_{\rm VH}=\frac{\sum_{i=1}^{N}\frac{{\bf S}_{i}{\bf z}_{i}}{{\bf d}_{i}}}{\sum_{i=1}^{N}\frac{{\bf S}_{i}}{{\bf d}_{i}}},$

(1)

where ${\bf S}_{i}=1$ indicates that the $i^{th}$ unit has been selected for sampling, and ${\bf S}_{i}=0$ indicates it has not been selected, and ${\bf z}_{i}$ represents the variable of interest measured on the $i^{th}$ unit.

We refer to this approach as the Volz-Heckathorn (VH) estimator. Gile and Handcock (2010) illustrate that this estimator consistently out-performs the Salganik-Heckathorn estimator, and we believe it is the most principled estimator currently available for RDS. For these reasons, we use the VH estimator as the standard for comparison in this paper.

Many critical assumptions required to justify this estimator are explored in Gile and Handcock (2010) and listed in Table 2. In this paper, we focus on the reliance on a with-replacement sampling model. We present an estimator based on an alternative sampling model reflecting the without-replacement nature of the sampling process.

The Volz-Heckathorn estimator requires many waves of sampling to justify its reliance on a stationary distribution. In practice the number of waves is small (almost always fewer than 20, and often 5 or fewer). Also, we wish to consider a without-replacement process for which stationarity does not apply. It is therefore instructive to consider a special case when stationarity is not necessary.

Consider a graph formed in the following manner: Begin with $N$ vertices, designated by indices $1:N$. Assign to each vertex a number of edge-ends according to arbitrary fixed degree distribution $\mathbb{N}=\mathbb{N}_{1},\mathbb{N}_{2},\ldots\mathbb{N}_{K}$, where $\mathbb{N}_{j}$ is the number of nodes of degree $j$, $1$ and $K$ are the minimum and maximum degrees, respectively, and subject to the constraint that twice the number of edges, $2E=\sum_{k\in 1\ldots K}k\mathbb{N}_{k}$, is even. Now select pairs of edge-ends completely at random and assign an edge connecting each pair. This procedure results in a variant of the so-called configuration model for networks, a popular null model for networks, especially in the physics literature (Molloy and Reed, 1995). Note that this formulation does allow loops, or links to oneself as well as multiple edges between the same pair of vertices, although the rate of these events decreases with increasing population size for fixed maximum degree $K$, such that several authors have suggested they are negligible for $K<(\bar{{\bf d}}N)^{1/2}$ or $K<N^{1/2}$, where $\bar{{\bf d}}$ is the mean degree of the network (Chung and Lu, 2002; Burda and Krzywicki, 2003; Boguna et al., 2004; Catanzaro et al., 2005; Foster et al., 2007).

Now consider a random walk ${\bf G}^{*}$ on a set of vertices with degrees given by ${\bf d}_{1},{\bf d}_{2},\ldots,{\bf d}_{N}$, such that $\sum_{i}\mathbb{I}({\bf d}_{i}=k)=\mathbb{N}_{k}$, where $\mathbb{I}(A)$ is the indicator function on $A$. Volz and Heckathorn (2008) consider the stationary distribution of this random walk for a fixed graph. Instead, we consider the transition probabilities of the corresponding walk over the distribution of all networks of fixed degree distribution constructed as above, in which the $j^{th}$ node visited, ${\bf G}^{*}_{j}$, is selected from the distribution of possible edges from node ${\bf G}^{*}_{j-1}$. The transition probabilities are then given by:

$\displaystyle P({\bf G}^{*}_{j}={\bf g}^{*}_{j}|{\bf G}^{*}_{1},{\bf G}^{*}_{2},\ldots{\bf G}^{*}_{j-1}={\bf g}^{*}_{1},{\bf g}^{*}_{2},\ldots{\bf g}^{*}_{j-1})=\left\{\begin{array}[]{cl}\frac{{\bf d}_{{\bf g}^{*}_{j}}}{2E-1}&{\bf g}^{*}_{j}\neq{\bf g}^{*}_{j-1}\\ \\ \frac{{\bf d}_{{\bf g}^{*}_{j}}-1}{2E-1}&{\bf g}^{*}_{j}={\bf g}^{*}_{j-1},\end{array}\right.$

(5)

where this probability is taken over the space of all possible configuration model graphs of given degree distribution, as well as over the steps of the random walk.

This procedure results in stationary distribution proportional to ${\bf d}_{i}$, and, in fact, selection probabilities at each step very nearly proportional to ${\bf d}_{i}$. Thus, for a network structure given by such a configuration model, the Volz-Heckathorn estimator constitutes a Hansen and Hurwitz (1943) estimator, without further requirement of sufficient waves for convergence.

Now consider the corresponding self-avoiding random walk ${\bf G}$. In this procedure,

$\displaystyle P({\bf G}_{j}={\bf g}_{j}|{\bf G}_{1},{\bf G}_{2},\ldots{\bf G}_{j-1}={\bf g}_{1},{\bf g}_{2},\ldots{\bf g}_{j-1})=\left\{\begin{array}[]{cl}\frac{{\bf d}_{{\bf g}_{j}}}{2E-\sum_{i=1}^{j-1}{\bf d}_{{\bf g}_{i}}}&{\bf g}_{j}\notin{\bf g}_{1}\ldots{\bf g}_{j-1}\\ 0&{\bf g}_{j}\in{\bf g}_{1}\ldots{\bf g}_{j-1}.\end{array}\right.$

(8)

This sampling procedure is mathematically equivalent to successive sampling (SS) or probability proportional to size without replacement sampling (PPSWOR), dating back to Yates and Grundy (1953), and typically defined by the following sampling process:

• •

  Begin with a population of $N$ units, denoted by indices $1\ldots N$ with varying sizes represented by ${\bf d}_{1},{\bf d}_{2},\ldots{\bf d}_{N}$, with $\sum_{i=1}^{N}{\bf d}_{i}=2E$, for total edges $E$.

• •

  Sample the first unit ${\bf G}_{1}$ from the full population $\{1\ldots N\}$ with probability proportional to size ${\bf d}_{i}$.

• •

  Select each subsequent unit with probability proportional to size from among the remaining units, such that conditional sampling probabilities are given by (8).

In the survey sampling literature, mostly based on the work of Raj (1956) and Murthy (1957), this sampling design is referred to as probability proportional to size without replacement (PPSWOR), and typically used in instances where the desired probability proportional to size (PPS) design is not feasible. In such cases, the sizes of population units are all known, the sampling design is implemented, and the main interest is in estimating the population total or mean of a variable measured on sampled units. The analytical properties of this procedure are quite difficult, as suggested by more recent work including Rao et al. (1991) and Kochar and Korwar (2001). In fact, even the marginal unit sampling probabilities are not available in closed form.

In the geological discovery literature, successive sampling is not a purposive design, but an approximation to a non-designed sampling process. Important work in this area includes Andreatta and Kaufman (1986), Nair and Wang (1989), and Bickel et al. (1992). These authors address the case of oil field discovery, in which successive fields are discovered with probabilities proportional to some measure of their size, typically taken to be their volume. This literature does not assume the full population of sizes to be known, and typically takes some function of the sizes, such as the sum of the sizes of undiscovered reserves, as the object of inference. Our application of successive sampling also assumes the unobserved sizes to be unknown, however as in Raj (1956) and Murthy (1957), our object is to use these sizes to estimate the population characteristics on another variable measured on all sampled units. To do so, we both develop a novel algorithm and leverage more recent work by Fattorini (2006).

Under successive sampling, for population of sizes given by $\mathbb{N}$ and sample size $n$, there is a function $f_{{\mbox{\boldmath$\pi$}}}(k;n,\mathbb{N})$ mapping the size $k$ of a unit to its sampling probability ${\mbox{\boldmath$\pi$}}_{i}$. Our proposed estimator is based on estimated sampling weights, which are based on the $\pi$ given by the successive sampling procedure, applied to nodal sampling units and “sizes” given by nodal degrees. There are three key challenges for this approach. The mapping depends on first, the known population size $N$ and, second, the degree distribution, $\mathbb{N}$, neither of which is known in the general RDS case. Finally, given $\mathbb{N}$ and the sample size $n$, the mapping is not explicit. The lack of explicit mapping has been addressed by Fattorini (2006), who suggests estimating the mapping by simulation. For known $\mathbb{N}$ and given $n$, he simulates the successive sampling procedure, then estimates the inclusion probability ${\mbox{\boldmath$\pi$}}_{i}$ associated with unit $i$ by:

$\displaystyle\tilde{{\mbox{\boldmath$\pi$}}}_{i}=\frac{{\bf U}_{i}+1}{M+1},$

(9)

where ${\bf U}_{i}$ is the number of times unit $i$ is sampled in the $M$ trials. He proposes using these estimated probabilities in the standard Horvitz-Thompson estimator:

$\displaystyle\sum_{j:{\bf S}_{j}=1}\frac{{\bf z}_{j}}{\tilde{{\mbox{\boldmath$\pi$}}}_{j}}.$

(10)

In most of this paper, we assume the population size, $N$, is known. We evaluate the sensitivity of our results to that assumption in Section 5.2, and in the examples in Section 6.

Given the known population size, we present a novel approach to estimating the degree distribution $\mathbb{N}$ jointly with inclusion probabilities when degrees are only observed for sampled nodes. This procedure can be applied beyond the RDS setting whenever the population distribution of unit sizes corresponding to a sample collected through successive sampling is unknown.

Our approach iteratively estimates the population distribution $\mathbb{N}$ and the mapping $f_{{\mbox{\boldmath$\pi$}}}(k;n,\mathbb{N}):k\to{\mbox{\boldmath$\pi$}}$. This approach relies on two key points:

• •

  For known population of degrees $\mathbb{N}$, the mapping $f_{{\mbox{\boldmath$\pi$}}}(k;n,\mathbb{N})$ can be estimated by simulation in a form similar to (9).

• •

  For known mapping $f_{{\mbox{\boldmath$\pi$}}}(k;n,\mathbb{N})$, the number (or proportion) of the population of degree $k$ can be estimated using a form similar to (10).

We leverage these two points to propose the following procedure for the estimation of population mean $\mu$ in the case of nodal degrees observed only for sampled units.

Let $\mathop{\rm E\,}\nolimits[\cdot;n,\mathbb{N}]$ denote expectation with respect to a sample of size $n$ sampled by successive sampling from a population with degree counts $\mathbb{N}=\{\mathbb{N}_{1},\mathbb{N}_{2},\ldots,\mathbb{N}_{K}\}.$ Then:

$\displaystyle\mathop{\rm E\,}\nolimits[{\bf V}_{k};n,\mathbb{N}]=\mathbb{N}_{k}f_{{\mbox{\boldmath$\pi$}}}(k;n,\mathbb{N})~~~~~~~~~k=1,\ldots,K$

(11)

where ${\bf V}_{k}$ is the random variable representing the number of sample units with degree $k$, $f_{{\mbox{\boldmath$\pi$}}}(k;n,\mathbb{N})=\mathop{\rm E\,}\nolimits[{\bf S}_{j}:{\bf d}_{j}=k;n,\mathbb{N}]$ is the (common) inclusion probability of a node of degree $k$, and $K$ is the maximum degree, $K<N$. This suggests first order moment equations for the unknown true $\mathbb{N}:$

$\displaystyle\mathop{\rm E\,}\nolimits[{\bf V}_{k};n,\mathbb{N}]={\bf v}_{k}~~~~~~~~~k=1,\ldots,K$

(12)

where ${\bf v}_{k}$ is the observed number of sample units with degree $k$.

The algorithm is then:

1. 1.

   Initial estimate:

   $\displaystyle f_{{\mbox{\boldmath$\pi$}}}(k;n,\mathbb{N}^{0})=\frac{k}{N}\sum_{l=1}^{K}\frac{v_{l}}{l},$

   (13)

   that is $f_{{\mbox{\boldmath$\pi$}}}(k;n,\mathbb{N}^{0})$ proportional to $k$.

2. 2.

   For $i=1,\ldots,r$ iterate the following steps:

   1. (a)

      Estimate the population distribution of degrees:

      $\displaystyle{\mathbb{N}_{k}^{i}=N\cdot\frac{\frac{{\bf v}_{k}}{f_{{\mbox{\boldmath$\pi$}}}(k;n,\mathbb{N}^{i-1})}}{\sum_{l=1}^{K}\frac{{\bf v}_{l}}{f_{{\mbox{\boldmath$\pi$}}}(l;n,\mathbb{N}^{i-1})}}}~~~~~~~~~~~k=1,\ldots,K$

      (14)

      where ${\bf v}_{k}$ is the observed number of sample units with degree $k$.

   2. (b)

      Estimate inclusion probabilities:

      1. i.

         Simulate $M$ successive sampling samples of size $n$ from a population with composition $\mathbb{N}^{i}.$

      2. ii.

         Estimate the inclusion probabilities:

         $\displaystyle f_{{\mbox{\boldmath$\pi$}}}(k;n,\mathbb{N}^{i})=\frac{\mathop{\rm E\,}\nolimits[{\bf V}_{k};n,\mathbb{N}^{i}]}{\mathbb{N}_{k}^{i}}\approx\frac{{\bf U}_{k}+1}{M\cdot\mathbb{N}_{k}^{i}+1},$

         (15)

         where ${\bf U}_{k}$ is the total number of observed units of size $k$ in the $M$ simulations.

3. 3.

   Estimate the population distribution of degrees and the corresponding inclusion probabilities via, respectively: ${\hat{\mathbb{N}}}=\mathbb{N}^{r}$ and ${\hat{\pi}}(\cdot)=f_{{\mbox{\boldmath$\pi$}}}(\cdot;n,\mathbb{N}^{r}).$

4. 4.

   Use the resulting mapping ${\hat{\pi}}$ to estimate $\mu$ via the generalized Horvitz-Thompson estimator:

   $\displaystyle\hat{\mu}_{SS}=\frac{\sum_{j=1}^{N}\frac{{\bf S}_{j}{\bf z}_{j}}{{\hat{\pi}}({\bf d}_{j})}}{\sum_{j=1}^{N}\frac{{\bf S}_{j}}{{\hat{\pi}}({\bf d}_{j})}}.$

   (16)

For computational efficiency, most simulations in this paper were conducted with $M=500$ and $r=3$, with good results. In general, we recommend at least $M=2000$ and $r=3$, and we have used these parameters for the simulations of the standard error procedure in the supplemental materials, in the application to real data, and in the extension to $N>1000$ in the discussion. Estimation time scales with sample size, population size, and $M$. In our simulations, with $N=1000$ and $M=500$, estimates require about $1.5$ seconds on a personal computer, increasing to about $6$ seconds when $M=2000$. In practice, these parameters can be adjusted for desired precision in the solution to (12). Higher values of $M$ are particularly helpful for more dispersed degree distributions and larger population sizes.

Simulations have shown the results of this procedure to be at least as good as those provided by a first order asymptotic approximation provided by Andreatta and Kaufman (1986). This approach is novel and its theoretical properties are not well understood. It is appealing is to consider the algorithm as a variant of an EM algorithm, such as an ECM algorithm (Meng and Rubin, 1993), with the unobserved part of the degree sequence, $d_{n+1},\ldots d_{N}$ as the latent variable. Unfortunately, in the design-based frame, these values also fully determine the unknown parameters $N_{k}=\sum_{i=1}^{N}{\mathbb{I}}(d_{i}=k)$, resulting in a degenerate likelihood form.

Our simulation study is designed to highlight the treatment of without-replacement sampling in the case of a large sample fraction. To increase the realism of the study, we chose parameters to match the characteristics of the pilot data from the CDC surveillance program (Abdul-Quader et al., 2006) wherever possible. The general procedure was as follows:

• •

  1000 networks are simulated under each test condition.

• •

  An RDS sample is simulated from each sampled network.

• •

  RDS estimators are computed from each sample.

Because the CDC’s surveillance system aims for a sample size of 500, and many RDS studies approach exhaustion of their populations of interest (Johnston, 2009; CDC, 2008), we fix all sample sizes at 500. We also consider a mean degree of 7, close to the mean of the pilot data from the CDC study (Abdul-Quader et al., 2006).

We assign a discoverable class to each member of the simulated population. In reference to studies designed to estimate the prevalence of infectious disease, we refer to this characteristic as “infection status,” assigning the “infected” status (${\bf z}_{i}=1$) to 20% of simulated population members in each simulation. Note that $\hat{\mu}_{SS}$ could also be applied to a continuous or categorical variable.

We consider networks sampled from models from the Exponential-family Random Graph Model (ERGM) class (Snijders et al., 2006). Here the relations ${\bf y}$ are represented as a realization of the random variable ${\bf Y}$ with distribution:

$\displaystyle P_{\eta}({\bf Y}={\bf y}|{\bf x})=\exp\{{\mbox{\boldmath$\eta$}}{\cdot}{\bf g}({\bf y},{\bf x})-\kappa({\mbox{\boldmath$\eta$}},{\bf x})\}\quad\quad{\bf y}\in{\cal Y},$

(19)

where ${\bf x}$ are nodal or dyadic covariates, ${\bf g}({\bf y},{\bf x})$ is a $p$-vector of network statistics, ${\mbox{\boldmath$\eta$}}\in\mathbb{R}^{p}$ is the parameter vector, ${\cal Y}$ is the set of all possible undirected graphs, and $\exp\{\kappa({\mbox{\boldmath$\eta$}},{\bf x})\}=\sum_{{\bf u}\in{\cal Y}}\exp\{{\mbox{\boldmath$\eta$}}{\cdot}{\bf g}({\bf u},{\bf x})\}$ is the normalizing constant (Barndorff-Nielsen, 1978). The structure of the networks represented is determined by the choice of ${\bf g}({\bf y},{\bf x})$.

In this study, we vary the structure of the networks in three ways: First, we consider populations of sizes (i.e. numbers of nodes) 1000, 835, 715, 625, 555, and 525, such that samples of size 500 constitute about 50%, 60%, 70%, 80%, 90%, and 95% of the target population. We focus on this range of sample fractions because it highlights settings in which we might expect to see finite population biases such as those $\hat{\mu}_{SS}$ is intended to address. Another critical feature is the activity ratio $w$, equal to the ratio of the mean degree of infected nodes to the mean degree of uninfected nodes,

$\displaystyle w=\frac{\sum_{i=1}^{N}d_{i}z_{i}}{\sum_{i=1}^{N}z_{i}}\frac{\sum_{i=1}^{N}(1-z_{i})}{\sum_{i=1}^{N}d_{i}(1-z_{i})}.$

(20)

This measures the differential tendency for the groups to be socially connected in the population. We consider $w\in\{0.5,0.8,1,1.1,1.4,1.8,2.5,3\}$.

We also induce homophily on infection status in these simulations, parameterized as the relative probability of an edge between two infected nodes, and an edge between an infected and an uninfected node. This is an intuitive parameterization, also used in Gile and Handcock (2010), but differs from that used in other analyses of RDS. Except in Section 5.3, the edge probability between the two infected nodes is fixed at five times that of the mixed dyad. For $w=1$, this, along with the 20% infected, implies that an edge between two uninfected nodes is twice as likely as an edge in a mixed dyad. We consider several other levels of homophily, and in Section 5.3 we show that this feature has important implications for the bias induced by biased initial samples, however for the case of initial samples at random with respect to ${\bf z}$, as in most of the simulations presented here, bias was not affected by level of homophily, although increasing homophily does increase the variance of both $\hat{\mu}_{SS}$ and $\hat{\mu}_{\rm VH}$, and is therefore important to consider in standard error estimation.

These features are represented in the ERGM by choosing network statistics to represent the mean degree, the activity ratio $w$, and homophily (based on using the “infected”  status as a nodal covariate). These three values were specified using the “nodemix” statnet model term, which includes three parameters corresponding to the three cells of the mixing matrix on infection. The parameter $\eta$ was chosen so the expected values of the statistics were equal to the values given above (van Duijn et al., 2009). Samples from the resulting models were taken using the statnet R package (Handcock et al., 2003, 2008).

The RDS sampling mechanism is again designed to mimic that of the CDC’s pilot study. Ten initial sample nodes were chosen for each sample, selected sequentially with probability proportional to degree, (i.e. by successive sampling). In the sensitivity analysis, we also consider initial sample selection regimes dependent on ${\bf z}$. Subsequent sample waves were selected without-replacement by sampling up to two nodes at random from among the un-sampled alters of each sampled node. Exactly two alters were sampled whenever possible. This process typically resulted in the sampling of four complete waves and part of a fifth wave, stopping when a sample size of 500 was attained.

We augment our basic results with two sub-studies evaluating the sensitivity of the estimator $\hat{\mu}_{SS}$ to assumptions not required for $\hat{\mu}_{\rm VH}$: the accuracy of the assumed population size $\hat{N}$ and the dependence on the initial sample.

The key insight of this paper is the recognition that the true mapping $f_{{\mbox{\boldmath$\pi$}}}(k;n,\mathbb{N})$ from nodal degree ${\bf d}_{i}$ to inclusion probability, ${\mbox{\boldmath$\pi$}}_{i}$ under Respondent-Driven Sampling is better approximated by successive sampling than by the linear mapping assumed by Volz and Heckathorn (2008). In addition, we introduce a novel approach to estimating the unit size (or degree) distribution in the population, based on the population size and sizes of observed units under successive sampling. Combining this insight and estimation strategy, we present a new estimator for population proportions based on an RDS sample.

The contribution of this new estimator is illustrated in Figure 6. In cases with no correlation between degree and quantity of interest ($w=0$), and initial sample selected at random with respect to infection status, the naive sample mean, $\hat{\mu}$ is an unbiased estimator of the population mean, as are existing estimators such as $\hat{\mu}_{\rm VH}$ and the proposed estimator $\hat{\mu}_{SS}$. This is illustrated in the first three columns of the figure.

However, when the variable of interest is related to nodal degree, classes with higher degrees are over-represented in the population mean, resulting in the positive bias in the sample mean $\hat{\mu}$ in the second three columns of Figure 6. This bias shrinks as the sample fraction increases. The existing estimator $\hat{\mu}_{\rm VH}$ adjusts reasonably well for this effect when the sample fraction is small, however for larger sample fractions, this estimator over-compensates for the effect of degree distribution, resulting in bias opposite that of $\hat{\mu}$. The contribution of the proposed estimator, $\hat{\mu}_{SS}$ is that it correctly adjusts for the joint effects of varying degree distributions and large sample fractions.

The last three columns of Figure 6 illustrate a shortcoming of all three of these estimators. In this case, all initial samples, or seeds, are selected from among the infected nodes. This results in increased positive bias in all three of these estimators. All three estimators are also subject to other sources of bias discussed in Gile and Handcock (2010), including bias induced by the systematically biased passing of RDS coupons.

Because RDS is in wide usage, and often in cases where the sample fraction is large, this new estimator may improve estimation in many contexts. Its applicability is limited, however, by the requirement that the population size $N$ is known. Figures 3 and 4 illustrate that inaccurate estimates of $N$ can introduce bias into $\hat{\mu}_{SS}$, although this new estimator still out-performs $\hat{\mu}_{\rm VH}$ with the level of inaccuracy considered in this study.

While most of our simulation study has focused on sample fractions consistent with the finite population effects that $\hat{\mu}_{SS}$ is intended to address, it is of interest to note that $\hat{\mu}_{SS}$ is nearly identical to $\hat{\mu}_{\rm VH}$ for small sample fractions. We have replicated much of the simulation study here with population size $N=10,000$. We find no significant differences between $\hat{\mu}_{SS}$ and $\hat{\mu}_{\rm VH}$ in these simulations, over a full range of values of $w$. When using inaccurate estimates of $N$ as in (21), we found only one case of significant difference between $\hat{\mu}_{\rm VH}$ and $\hat{\mu}_{SS}$, corresponding to a dramatic under-estimate of population size, $\hat{N_{s}}=5250,N=10,000$, and for high activity ratio, $w=3$. In this case, $\hat{\mu}_{SS}$ had bias $0.0139$ and $\hat{\mu}_{\rm VH}$ had bias $0.0105$. The difference in mean squared error was not significant. Thus, the proposed estimator is helpful in correcting for finite population biases when they exist, and sensitive to the estimated population size in these cases. In cases where finite population effects are not present, the proposed estimator nearly coincides with existing estimator $\hat{\mu}_{\rm VH}$, and, unless the inaccurate population size estimate is small enough to transition to the region of finite population sensitivity, $\hat{\mu}_{SS}$ is insensitive to population size estimates in these cases.

Known population size $N$ and random mixing are key assumptions of the estimator $\hat{\mu}_{SS}$. Additional assumptions required by $\hat{\mu}_{SS}$ are listed in Table 2. The assumptions with strikethrough are necessary for $\hat{\mu}_{\rm VH}$ but not for $\hat{\mu}_{SS}$, and the italic text indicates additional assumptions required by $\hat{\mu}_{SS}$ but not $\hat{\mu}_{\rm VH}$. In these terms, the key contribution of $\hat{\mu}_{SS}$ is to remove the dependence on a known inaccurate random walk model for estimating sampling probabilities. In return, $\hat{\mu}_{SS}$ sacrifices the theoretical robustness to the initial sample promised by the Markov chain model, although this robustness was not truly present in $\hat{\mu}_{\rm VH}$ either, and the proposed estimator performs no worse than the former in this respect (Table 1). More critically, $\hat{\mu}_{SS}$ relies on the assumption of known population size. The estimator $\hat{\mu}_{SS}$ is also sensitive to the degree of non-random mixing, or homophily in the network, although no more sensitive than $\hat{\mu}_{\rm VH}$ (Table 1), as well as to other assumptions such as accurate self-reported degree and an undirected network. Gile and Handcock (2010) illustrate the sensitivity of $\hat{\mu}_{\rm VH}$ to deviations from some other assumptions in Table 2. We do not expect $\hat{\mu}_{SS}$ to be more robust to such deviations than $\hat{\mu}_{\rm VH}$.

Our application to HIV prevalence estimation in Section 6 illustrates several uses for $\hat{\mu}_{SS}$. First, when the population size is known, $\hat{\mu}_{SS}$ provides a better estimate of population proportions than the other available methods. This can also be done when only a range of population sizes is known. In particular, in the case of MSM, officials are often willing to assume that the population ranges between 1% and 3%. When no information is available on population size, $\hat{\mu}_{SS}$ can be used to perform a sensitivity analysis. Finally, other information from the sampling process, such as the exhaustion of the population available for sampling, may suggest that the sample fraction is large. Additional information may also be gathered by asking respondents about their exposure to others who have been sampled.

Intuition suggests that in the case of a large sample fraction, an RDS estimator should negotiate between the infinite-population assumption of the Volz-Heckathorn estimator and the full-population sample assumption of the naive sample mean. In this paper, we provide an estimator based on a successive sampling model that appropriately negotiates between these two extremes. Whenever the hidden population size is known, it will be preferable to use this estimator. When the hidden population size is not known, this estimator still provides a helpful diagnostic check on the other available estimators.

Beyond the estimator itself, this paper contributes a new theoretical framework for understanding the sampling process in respondent-driven sampling. Previous understandings have relied on a with-replacement or infinite population assumption, an assumption known to be inaccurate, and critically so in some populations. The introduction of a sampling model appropriately accommodating the true without-replacement nature of the sampling process opens new possibilities for future research on respondent-driven sampling. We intend to make code available for these procedures in the R package RDS on CRAN.