arXiv:1509.04704v2 [stat.ME] 28 Aug 2016

Central limit theorems for network driven
sampling
Xiao Li
School of Mathematical Sciences
Peking University
Karl Rohe
Department of Statistics
University of Wisconsin-Madison
Abstract
Respondent-Driven Sampling is a popular technique for sampling hidden populations. This paper models Respondent-Driven Sampling as a Markov process
indexed by a tree. Our main results show that the Volz-Heckathorn estimator is
asymptotically normal below a critical threshold. The key technical difficulties
stem from (i) the dependence between samples and (ii) the tree structure which
characterizes the dependence. The theorems allow the growth rate of the tree to
exceed one and suggest that this growth rate should not be too large. To illustrate
the usefulness of these results beyond their obvious use, an example shows that in
certain cases the sample average is preferable to inverse probability weighting. We
provide a test statistic to distinguish between these two cases.

1

Introduction

Classical sampling requires a sampling frame, a list of individuals in the target population with a method to contact each individual (e.g. a phone number). For many populations, constructing a sampling frame is infeasible. Network driven sampling enables
researchers to access populations of people, webpages, and proteins that are otherwise
difficult to reach. These techniques go by many names: web crawling, RespondentDriven Sampling, breadth-first search, snowball sampling, co-immunoprecipitation,
and chromatin immunoprecipitation. In each application, the only way to reach the
population of interest is by asking participants to refer friends.
Respondent-Driven Sampling (RDS) serves as a motivating example for this paper. The Centers for Disease Control, the World Health Organization, and the Joint
United Nations Programme on HIV/AIDS have invested in RDS to reach marginalized
and hard-to-reach populations [Heckathorn, 1997, WHO, 2013]. Each individual i in
the population has a corresponding feature yi (e.g. yi ∈ {0, 1} and yi = 1 if i is
HIV+). Using only the sampled individuals, we wish to make inferences about the average value of yi across the entire population, denoted as µ (e.g. the proportion of the
1

population that is HIV+). Extensive previous statistical research has proposed various
estimators of µ which are approximately unbiased based upon various types of models for an RDS sample [Salganik and Heckathorn, 2004, Volz and Heckathorn, 2008,
Gile, 2011]. We note that in the papers cited above (except [Gile, 2011]), RDS is assumed to sample with replacement. Previous research has also explored the variance of
these estimators [Goel and Salganik, 2009, Rohe, 2015]. This paper studies the asymptotic distribution of statistics related to these estimators.
Results on asymptotic distributions for RDS are useful for two obvious reasons.
First, they allow us to construct asymptotic confidence intervals for µ. Second, they
provide essential tools to test various statistical hypotheses. The only central limit
theorem associated considered in the RDS literature studied the case when the tree
indexed process reduces to a Markov chain [Goel and Salganik, 2009]; this presumes
that each individual refers exactly one person. Previous research suggests that the
number of referrals from each individual is fundamental in determining the variance of
common estimators [Rohe, 2015]. This paper establishes two central limit theorems in
settings which allow for multiple referrals.
The main results apply to both the sample average and the Volz-Heckathorn estimator, which is an approximation of the inverse probability weighted estimator (cf Remark 1). Because the inverse probability weighted (IPW) estimator and its extensions
are asymptotically unbiased, these estimators are often preferred to the sample average.
However, sometimes survey weights are not needed and they only introduce additional
variance to the estimator [Bollen et al., 2016]. This issue is particularly salient when
sampling weights are highly heterogeneous, as is often the case in RDS. Proposition
3 shows that if the outcomes yi are uncorrelated with the sampling weights, then the
sample average is unbiased. Theorem 3 extends this result to RDS to show that the IPW
estimator can have a larger variance than the sample average. Taken together, these results imply that the sample average can have a lower mean squared error (MSE) than
the IPW estimator. Section 4 introduces an estimator of the bias of the sample average.
The main results provide a path to test the null hypothesis that the bias is zero. This
can be used to select between the sample average and the IPW estimator. Section 6.2
studies this routine with the AddHealth social network.

2

Notation

Following [Goel and Salganik, 2009] and [Rohe, 2015], the results below model the
network sampling mechanism as a tree indexed Markov process on a graph. There
are many assumptions in this model which are incorrect in practice. However, like
the i.i.d assumption, it allows for tractable calculations. In the simulations, we show
that the theory derived from this model provides a good approximation for a more
realistic sampling model. [Lu et al., 2012] studies the sensitivities of the estimators to
this model.
Let G = (V, E) be a finite, undirected, and simple graph with vertex set V =
{1, ..., N } and edge set E. V contains the individuals in the population and E describes
how they are related to one another. As discussed in the introduction, y : V → R is
a fixed real-valued function on the state space V ; these are the node features that are
2

PN
measured on the sampled nodes. The target of RDS is to estimate µ = N −1 i=1 y(i).
If each sampled node referred exactly one friend, then the Markov sampling procedure would be a Markov chain. Several classical central limit theorems exist for this
model; see [Jones et al., 2004] for a review. The results herein allow for each sampled
node to refer more than one node. This is a Markov process indexed not by a chain,
but rather by a tree. Denote the referral tree as T. Where the node set of G indexes the
population, the node set of T indexes the samples. That is, we observe a subset of the
individuals in G with the sample {Xτ }τ ∈T ⊂ V . An edge (σ, τ ) in the referral tree
denotes that sampled individual Xσ referred individual Xτ into the sample. Mathematically, T is a rooted tree–a connected graph with n nodes, no cycles, and a vertex 0
which indexes the seed node. To simplify notation, σ ∈ T is used synonymously with
σ belonging to the vertex set of T.
For each non-root node τ ∈ T, denote p(τ ) ∈ T as the parent of τ (i.e. the node
one step closer to the root). This paper presumes that {Xτ }τ ∈T is a tree-indexed random walk on G, which was a model introduced by [Benjamini and Peres, 1994]. This
model generalizes a Markov chain on G; each transition Xp(τ ) → Xτ is an independent and identically distributed Markov transition with transition matrix P . Following
[Benjamini and Peres, 1994], we will call this process a (T, P )-walk on G. Unless
stated otherwise, it will be presumed throughout that the root node of the random walk
X0 is initialized from the equilibrium distribution of P . It follows that Xσ has distribution π for all σ ∈ T.
Unless stated otherwise, the results in this paper allow for the transition matrix
P to be constructed from a weighted graph G. Let wij be the weight of the edge
(i, j) ∈ E; if (i, j) 6∈ E, define wij = 0. If the graph is unweighted,
then let wij = 1
P
for all (i, j) ∈ E. Define the degree of node i as deg(i) = j wij . If the graph is
unweighted, then deg(i) is the number of connections to node i. Throughout this paper,
the graph is undirected. So, wij = wji for all pairs i, j. Given that {Xp(τ ) = i}, the
probability of {Xτ = j} is proportional to wij ;

wij
.
P Xτ = j|Xp(τ ) = i =
deg(i)
We use the term simple random walk for the Markov chain constructed on the unweighted graph (i.e. wi,j ∈ {0, 1} for all i, j). The simple random walk presumes that
each participant selects a friend uniformly and independently at random from their list
of friends.
In order to estimate µ, we observe y(Xτ ) for all τ ∈ T. Because G is undirected, P
is reversible and has stationary distribution π with πi ∝ deg(i) for all i ∈ G; this fact
is helpful for creating an asymptotically unbiased estimator for µ, particularly under
the simple random walk assumption [Volz and Heckathorn, 2008].
PN
Remark 1. In general, the quantity of interest µ = N −1 i=1 y(i) is not equal to
Eπ (y). As such, the sample average of y(Xτ )’s is a biased estimator for µ. With
inverse probability weighting, define a new function y 0 (i) = y(i)(N πi )−1 and the
respective estimator
1 X y(Xσ )
1X 0
µ̂IP W =
y (Xσ ) =
,
n
n
N πXσ
σ∈T

σ∈T

3

where n = |T| is the sample size. Then, Eπ (µ̂IP W ) = Eπ (y 0 ) = µ. As such, the
sample average of the y 0 (Xτ )’s is an unbiased estimator of µ. Unfortunately, the
values πi are unknown. In practice, RDS participants are asked various questions to
measure how many friends they have in G. Under the simple random walk assumption,
πi is proportional to the number of friends of i; this result also requires that the edges
in G are undirected, something that will be presumed throughout the paper. Therefore
the Volz-Heckathorn estimator
µ̂V H =

X y(Xσ )/deg(Xσ )
P
τ ∈T 1/deg(Xτ )

σ∈T

is in essence a Hájek estimator based upon deg(i) [Volz and Heckathorn, 2008]. Under the simple random walk assumption, this estimator provides an asymptotically unbiased estimator of µ.
For each node τ ∈ T, let |τ | be the distance of the node from the root; this is
also called the “wave” of τ . For every pair of node σ, τ ∈ T, define d(σ, τ ) to be
the distance between σ and τ on T (as a graph). For each non-leaf node σ ∈ T, let
η(σ) be the number of offspring of σ. A tree is said to be an m-tree of height h if
η(σ) = m for all σ ∈ T with |σ| < h and η(σ) = 0 for all |σ| = h. Here, both
m and h are a natural numbers (i.e. m, h ∈ N). T is said to be Galton-Watson if
η(σ) are i.i.d random variables in N. While the theorems below only study 2-trees; the
computational experiments in Section 6.1 suggest that the conclusions of the analytical
results are highly robust to replacing the 2-tree with a Galton-Watson tree.
[Levin et al., 2009] serves as this paper’s key P
reference for Markov processes. FolN
lowing the notation in that text, define Eπ (y) = i=1 πi y(i) and varπ (y) = Eπ (y −
Eπ (y))2 for the function y.
There are two primary concerns about the model and estimator used in the main
results below. First, the Markov model allows for resampling. Second, the results
below only apply to m-trees, not more general trees. The simulations in Section 6.1
suggest that the analytic results continue to hold under a more realistic setting that
addresses both of these concerns.

3

Main Results

The threshold m < λ−2
2 was previously identified in [Rohe, 2015] as being a critical
threshold for the design effect of network driven sampling; beyond this threshold, the
variance of the standard estimator does not decay at the standard rate. In other words,
var( p

1

X

|σ ∈ T : |σ| ≤ h| σ∈T:|σ|≤h

y(Xσ )) → ∞

as h → ∞. As such, using the traditional scaling, no central limit theorem holds above
the critical threshold. Because of this, the theorems focus on the case m < λ−2
2 . When
m > λ−2
,
the
simulations
in
Section
6.1
suggest
that
the
central
limit
theorem
does
2
not hold for any scaling.
4

Theorem 1 is a central limit theorem for an estimator constructed from the treeindexed Markov chain. The theorem holds for any function y, any reversible transition
matrix with second largest eigenvalue satisfying |λ2 | =
6 1, and any m < λ−2
2 .
Theorem 1. Suppose that P is a reversible transition matrix with respect to the equilibrium distribution π, and that the eigenvalues of P are 1 = λ1 > |λ2 | ≥ ... ≥ |λN |.
Without loss of generality, suppose that Eπ (y) = 0. Define
1
Yi = √
mi

X

y(Xτ ).

τ ∈T:|τ |=i

If T is an m-tree with m < λ2−2 , then
h

1 X
√
Yi → N (0, σ02 )
h i=1
√
√
in distribution, where σ02 = varπ (( mP − I)−1 y) − varπ (P ( mP − I)−1 y).
The sequence of random variables considered in Theorem 1 are not exactly sample
averages, but a reweighted form of sample average. Samples in the same wave are
equally weighted, while samples from different waves are not. The following theorem
provides theoretical guarantee on the distribution of sample average for a specific class
of transition matrix and node feature.
Theorem 2. Let TPbe a 2−tree. Without loss of generality, suppose that Eπ (y) = 0.
Define µ̂h = √1 h σ∈T,|σ|≤h y(Xσ ). Suppose that
2

(c1) E(µ̂2k+1
) = 0 for all h, k ∈ N;
h
(c2) for any function f on V satisfying Eπ f = 0, ||P f ||∞ ≤ |λ2 |||f ||∞ ;
(c3)|λ2 | < √12 ;
then
µ̂h → N (0, σ02 )
in distribution for some σ02 .
Remark 2. Condition (c1) is a technical condition on the symmetry of µ̂h that is necessary in the proof. The following proposition provides a sufficient condition for (c1).
Proposition 1. Suppose that y is symmetric, i.e. for any i ∈ V there exists j such
that y(j) = −y(i). If p(u, v) = P (y(Xσ ) = v|y(Xp(σ) ) = u) is well-defined and
p(u, v) = p(−u, −v) for all u, v ∈ y(V ), then condition (c1) is satisfied.
Proof. Under the conditions of the proposition, the distribution of µ̂h is symmetric
with respect to 0. Thus E(µ̂2k+1
) = 0 for all h, k ∈ N.
h
Conditions (c2)-(c3) can be substituted by the following condition (c2’):
(c2’) There exists c < √12 such that for any function f on V satisfying Eπ f = 0,
||P f ||∞ ≤ c||f ||∞ .

5

Condition (c2’) is weaker than (c2) and (c3) combined, but is stronger than (c3)
alone. To see this, let f be the eigenfunction of the second eigenvalue, and it follows
that |λ2 | < √12 . It can be easily seen that one necessary condition for (c2’) is that
X
j

1
|Pij − πj | < √
2

for all i ∈ V . In other words, all the rows of P must be close to π. As previously
discussed, condition (c3) is actually a necessary condition for the central limit theorem
[Rohe, 2015], in the sense that the variance of µ̂h tends to infinity if |λ2 | ≥ √12 .
For clarity in the exposition of the theorem and the proof, we have only proved
the theorem for the 2-tree. We believe that similar results are likely to hold for more
general m-trees.

3.1

Extension to the Volz-Heckathorn estimator

When P is restricted to be the transition matrix of the simple random walk on G, the
following corollary shows that Theorem 2 can be extended to the Volz-Heckathorn
estimator [Volz andP
Heckathorn, 2008].
Denote d¯ = N1 i∈V deg(i) as the average node degree. Following Remark 1, the
¯
IPW estimator contains 1/(N πi ) which is equal to d/deg(i).
The Volz-Heckathorn es¯
timator first estimates d with the harmonic mean of the observed degrees. Because this
harmonic mean converges to d¯ in probability, the following corollary applies Slutsky’s
Theorem to give a central limit theorem for the Volz-Heckathorn estimator.
Corollary 1. Let T be a 2-tree. Suppose in particular that P is the transition matrix of
the simple random walk on G. Define a new node feature y 0 (i) = y(i)/deg(i). Without
loss of generality, suppose that Eπ y 0 = 0 (this is not equivalent to Eπ y = 0). Define
1
µ̂h,V H = µ̂h dˆ¯ = √
2h
where

X

ˆ¯
y 0 (Xσ )d,

σ∈T,|σ|≤h

2h+1 − 1
.
σ∈T,|σ|≤h 1/deg(Xτ )

¯= P
dˆ

If the new node feature y 0 and the transition matrix P satisfy conditions (c1)-(c3) in
Theorem 2, then
2
µ̂h,V H → N (0, σ0,V
H)
2
in distribution for some σ0,V
H.

3.2

Illustrating the conditions with a blockmodel

Consider G as coming from a blockmodel with two blocks [Lorrain and White, 1971].
In this blockmodel, each node i is given a label z(i) ∈ {1, 2} and every edge weight
wi,j = Bz(i),z(j) for some symmetric 2 × 2 matrix B. Suppose that yi = yj if z(i) =
6

z(j). This model was previously studied in [Goel and Salganik, 2009] and it serves as
an approximation to the Stochastic Blockmodel.
Given the structural equivalence of nodes within the same block, it is sufficient to
study the transition matrix between blocks, P ∈ R2×2 . If P is a symmetric matrix
with entries p11 = p22 = p and p12 = p21 = 1 − p for some value p, then it can be
easily verified that condition (c1) is satisfied. Moreover, if 2p−1 < √12 , then conditions
(c2) and (c3) are also satisfied. Our theorem asserts that the Volz-Heckathorn estimator
converges to the normal distribution in this model.
More generally, suppose that the nodes in the blockmodel for G are equally balanced between 2K blocks with node features {y1 , −y1 , . . . , yK , −yK } and that the
1−p
I(u 6= v). We can verify that all the
transition matrix p(u, v) = pI(u = v) + 2K−1
1
1
1
conditions are satisfied as long as 2K < p < 2K + 2√
.
2

4

Comparing the variance of inverse probability weighting to the bias of the sample average

An estimator with a small mean square error (i.e. E(µ̂ − µ)2 ) has a small bias and
a small variance. It is generally known that inverse probability weighting provides
an unbiased estimator of µ. However, survey weights can also drastically inflate the
variance of the estimator. This matter has been heavily studied by survey statisticians and substantial literature have devoted to the methodologies and issues regarding the use of sampling weights; see [Pfeffermann, 1996], [Biemer and Christ, 2007],
[Valliant et al., 2013], and [Bollen et al., 2016] for a review. To determine whether one
should use sampling weights in RDS, this section gives a test statistic for the null hypothesis that the sample average is unbiased. The results in the previous section suggest
a confidence region for this test statistic.
Denote n = |T| as the number of samples. The next results compare µ̂IP W to the
sample average
1X
µ̂ =
y(Xτ ).
n
τ ∈T

Proposition 2 and Theorem 3 highlight the dangers of inverse probability weighting
by showing that it can increase the variance. Proposition 2 studies the simplified case
where the samples are i.i.d from the stationary distribution. Following the proposition,
Theorem 3 studies the more relevant setting of the (T, P )-walk on G. To simplify the
statements of Proposition 2 and Theorem 3 and their proofs, the node features y(i)
are presumed to be random variables. This condition could be removed with further
technical conditions on the moments of y : V → R and its relationship to π.
Proposition 2. Suppose that X1 , ..., Xn are sampled independently from the stationary distribution π and that y(1), ..., y(N ) are N uncorrelated and identically distributed random variables with finite second moment µ2 . Let C1 = max1≤i≤N N πi
PN
and var(π) = N1 i=1 πi2 − N12 , then
var(µ̂IP W ) − var(µ̂) ≥ µ2 N (
7

N
− 1)var(π).
nC12

(1)

Thus, as long as N > C12 n, which can be easily satisfied in practice,
var(µ̂IP W ) > var(µ̂).
Theorem 3. Suppose that {Xτ : τ ∈ T} is a sample from the (T, P )-walk on G, and
that y(1), ..., y(N ) are N uncorrelated and identically distributed random variables
with finite second moment µ2 . Assume that there exist constants C1 , C2 and C3 (not
the same constants as in Proposition 2) such that C1 N ≤ di ≤ C2 N for all i and
N 2 var(π) > C3 . Then
var(µ̂IP W ) − var(µ̂) ≥ µ2 (

C2
N 2 C12
var(π) −
),
nC22
N C12

(2)

and there exists C independent of n such that var(µ̂IP W ) > var(µ̂) as long as N >
Cn.
Proposition 2 shows that as var(π) increases, the difference between the variance
of µ̂IP W and µ̂ becomes larger. Similarly, Theorem 3 also shows that the difference
between the variance of µ̂IP W and µ̂ increases as var(π) increases, given that the
relative upper and lower bound of d (i.e. C1 and C2 ) remain fixed. Recall that when
P is the simple random walk, the probabilities πi are proportional to node degree. An
extensive literature (e.g. [Strogatz, 2001, Clauset et al., 2009]) has found that empirical
networks have highly heterogeneous node degrees. As such, Equation 1 shows that the
variance of µ̂IP W can be dramatically greater than the variance of µ̂. Moreover, both
Proposition 2 and Theorem 3 suggest that the variance difference var(µ̂IP W )−var(µ̂)
can be considerable if we sample only a small proportion of the whole population. This
problem is particularly salient when the population is large.
The bias-variance tradeoff presents a dilemma between inverse probability weighting and the sample average. This bias can be estimated. For every i ∈ V , define a new
node feature,
1
d¯
y 0 (i) = y(i)(1 −
) = y(i)(1 −
).
(3)
N πi
deg(i)
Proposition 3. The mean of the new node feature satisfies
Eπ (y 0 ) = Eπ (y) − µ,
which is the true bias of the sample average. Therefore, under the null hypothesis,
H0 : Eπ (y 0 ) = 0,
theP
sample average is an unbiased estimator. If π and y are uncorrelated (i.e.
1
i∈V y(i)), then H0 is satisfied and the sample average is unbiased.
N

P

i∈V

πi y(i) =

Under the conditions of Proposition 2, π and y satisfy the condition of Proposition
3 in expectation (i.e. they are uncorrelated in expectation). Both of these conditions
imply that the outcome is unrelated to the sampling weight (in some way). Under
such conditions, both estimators are unbiased. If it is also true that var(π) is large,
then µ̂ has a smaller variance. In this scenario, µ̂ is preferable to µ̂IP W in terms of
8

MSE. However, if the sample average is biased, then we must compare the bias and
variance of the two estimators. Asymptotically, the variance of both estimators will
vanish, while the bias stays constant. So, for sufficiently large sample size, one should
use µ̂IP W . For smaller sample sizes, the bias of the sample average could be small
(relative to the difference in the variances). In such settings, there will be a crossover
point, a sample size at which µ̂IP W becomes preferable to µ̂. To distinguish between
these two cases, we want to test the null hypothesis that the sample average is unbiased,
i.e. Eπ (y 0 ) = 0. Or, more generally, we want to provide a confidence region for the
bias of the sample average.
If d¯ is unknown, as is generally the case, we can estimate d¯ by the harmonic mean
[Salganik and Heckathorn, 2004]
n
.
dˆ¯ = P
1/deg(X
τ)
τ ∈T
¯ in Equation 3 yields the new node feature based on the VolzSubstituting d¯ for dˆ
Heckathorn estimator [Volz and Heckathorn, 2008]
yV0 H (i) = y(i)(1 −

dˆ¯
).
deg(i)

Similarly, define
X
1X
dˆ¯
d = 1
yV0 H (Xσ ) =
y(Xσ ) 1 −
bias
n
n
deg(Xτ )
σ∈T

!
= µ̂ − µ̂V H ,

(4)

σ∈T

d = µ̂ − µ̂V H is an asymptotically unbiased estimator for the bias of µ̂. It
then bias
serves as a test statistic for the null hypothesis H0 : Eπ (y 0 ) = 0.
d converges to the normal distribution. The
The theorems above suggest that bias
rejection region is then
σˆ0
d > 1.96 √
W = {|bias|
}.
n
where σˆ2 is an estimate of the variance.
0

5

Estimating the variance

d
For some node feature ỹ (e.g. HIV status y or the y 0 in Equation (3) that motivate bias),
2
let µ̃ denote the sample average. Denote σµ̃ as V arT,P (µ̃), where the subscript T, P
denotes that the data is collected via a (T, P )-walk on G. This subsection studies a
simple plug-in estimator for σµ̃2 .
The following function is essential to expressing σ [Rohe, 2015].
Definition 1. Select two nodes I, J uniformly at random from the tree T. Define the
random variable D = d(I, J) to be the graph distance in T between I and J. Define
G as the probability generating function for D,
G(z) = E(z D ).
9

In practice, T is observed. So, the function G can be computed. In many studies
there are multiple seed nodes. In these cases, we suggest computing d(I, J) on a tree
which has an artificial root node that connects to all of the seeds; this root node could
be imagined as an individual that is responsible for finding the seed nodes. In this tree,
two different seed nodes would be distance 2 apart.
Denote the autocorrelation at lag 1 of ỹ(Xτ ) as
R=

Cov(ỹ(Xp(τ ) ), ỹ(Xτ ))
.
varπ (ỹ)

Both Cov(ỹ(Xp(τ ) ), ỹ(Xτ )) and varπ can be estimated with plug-in quantities. Because the data has been sampled proportional to π, the plug-in quantity for varπ should
not explicitly adjust for π.
1X
\
(ỹ(Xτ ) − µ̂V H )2 .
var
π (ỹ) =
n
τ ∈T

Similarly,
Cov(ỹ(X\
p(τ ) ), ỹ(Xτ )) =

X
1
(ỹ(Xp(τ ) ) − µ̂V H )(ỹ(Xτ ) − µ̂V H ),
n−1
τ ∈T\0

where {T \ 0} contains all nodes except the root node 0 (because p(0) does not exist).
Using these plug-in quantities, define R̂. Then, the estimator is
\
σ̂µ̃2 = G(R̂)var
π (ỹ).
A popular bootstrap technique for estimating σ̂µ̃2 resamples y(Xτ ) as a Markov
process (i.e. in addition to Xτ being a Markov process, the bootstrap procedure also
assumes that y(Xτ ) is Markov) [Salganik, 2006]. This model is akin to the blockmodel with two blocks in Section 3.2. The following assumption is weaker than this
assumption.
Assumption 1: ỹ(i) = µ + σf (i), where µ, σ ∈ R and f : V → R is an eigenfunction of P with varπ (f ) = 1.
Proposition 4. Under Assumption 1,
σµ̃2 = G(R)varπ (ỹ).
While Assumption 1 is weaker than the previous assumption in [Salganik, 2006],
the next proposition highlights the danger of this assumption. It uses a different assumption which is a rather weak assumption.
Assumption 2: G is convex on [λmin , 1], where λmin is the smallest eigenvalue of
P.
Because G is a probability generating function, it is always convex on [0, 1]. As
such, we on need to be worried about negative√values. Recall, that the central limit
theorems above only hold when |λmin | < 1/ 2 ≈ .7 (the smallest possible value
for λmin is −1). Some simulated trees given in the appendix suggest that if G is not
10

convex,
√ it often fails in the neighborhood of −1. As such, the assumption that |λmin | <
1/ 2 ≈ .7 is likely to imply Assumption 2. In practice, one observes the referral tree
T. Thus, one can compute the second derivative of G. Eigenvalues of P close to
negative one arise in antithetic sampling, where adjacent samples are dissimilar. For
example, if the population in G was heterosexuals and edges in G represent sexual
contact, then men would only refer women and vice versa. In this case, λmin would be
exactly −1. While easily imagined, such settings are not current practice for RDS. As
such, large an negative values are uncommon; λmin is likely close to zero.
The following proposition follows from an application of Jensen’s inequality.
Proposition 5. Under Assumption 2,
σµ̃2 ≥ G(R)varπ (ỹ).
Because Assumption 2 is not very restrictive, the inequality in Proposition 5 highlights the danger in breaking Assumption 1 (and thus the Markov model in [Salganik, 2006]);
breaking Assumption 1 will lead to σ̂µ̃2 underestimating the variance.

5.1

A bias adjusted estimator for µ

The test above also allows us to derive a more robust estimator of µ. Define
n

X
d = 1
y 0 (Xi ).
bias
n i=1
d = µ̂ − µ̂IP W . Using the hypothesis test to choose between the sample averThen bias
age and inverse probability weighting is akin to hard thresholding the bias adjustment.
Define

d if reject H0
bias
d
biast =
0
otherwise.
d t (BA for bias adjustment). This
The final estimator of µ is then µ̂BA = µ̂ − bias
estimator is explored in Section 6.2.

6
6.1

Numerical results
Simulation

In this section we illustrate the theoretical results on simulated data. The simulations
are performed on networks simulated from the Stochastic Blockmodel [Holland et al., 1983].
The four colors in Figure 1 correspond to four different networks that are simulated
from four different Stochastic Blockmodels. Each of the four networks has N =5,000
nodes, equally balanced between group zero and group one. The probability of a
connection between two nodes in different blocks is r and the probability of connection between two nodes in the same block is p. To control the eigenvalues of the

11

5000 × 5000 transition matrix, consider the transition matrix between classes given by
P = E(D)−1 E(A). The second eigenvalue of P is [Rohe et al., 2011]
λ2 (P) =

p−r
p+r

where expectations are under the Stochastic Blockmodel. In our simulation, the second
eigenvalue of the actual transition matrix is typically very close to λ2 (P). We take
p + r = 0.01 in all four Stochastic Blockmodels so that the average degree is about 25.
As such, λ2 (P) is actually controlled by p − r.
For each of the four networks we carry out four different sampling designs. Let
T be either a 2−tree or a Galton-Watson tree with E(η(σ)) = 2. For the GaltonWatson tree, the distribution of η(σ) is uniform on {1, 2, 3}. For each T, we consider
both with replacement sampling (i.e. the (T, P )-walk on G) and without replacement
sampling (i.e. referrals are sampled uniformly from the friends that have not yet been
sampled). Note that the conditions of Theorem 2 may be violated when either the
Galton-Watson tree or without-replacement sampling is used. We take the first 8 waves
of T as our sample. As such, the sample size is roughly N/10. For each social network
and
design, we repeat the sampling process 2000 times and compute µ̂ =
Psampling
n
1
y(X
)
i for each sample. The Quantile-Quantile (Q-Q) plot of µ̂ is shown in the
i=1
n
left panel of Figure 1; note that the QQ plot centers and scales each distribution to have
mean zero and standard deviation one. In addition, we repeat the above simulation for
the Volz-Heckathorn estimator, and the QQ plot of µ̂V H is shown in the right panel of
Figure 1.
√It is clear from Figure 1 that there are two patterns of distribution: when λ2 <
1/ m ≈ 0.7, i.e. λ2 =0.5 or 0.6, the Q-Q plots of µ̂√
and µ̂V H approximately lie on the
line y = x for all sampling design; when λ2 > 1/ m ≈ 0.7, i.e. λ2 =0.8 or 0.9, the
Q-Q plot of µ̂ and µ̂V H departs from the line y = x. Taken together, Figure 1 suggests
that the distribution of µ̂ and µ̂V H converges to Gaussian distribution if and only if
m < λ−2
2 . Actually, the right panel of Figure 1 implies that there are two modes in
the asymptotic distribution of µ̂ and µ̂V H when m > λ−2
2 . The relationship between
the expectation of the offspring distribution and the second eigenvalue of the social
network determines the asymptotic distribution of RDS estimators, regardless of the
node feature, the particular structure of the tree or the way we handle replacement.

6.2

Analysis of Adolescent Health data

To illustrate Theorem 1 with the test statistic in Section 3, this section presents simulation results that use the Comm15 friendship network from the National Longitudinal
Survey of Adolescent Health. This simulation compares the MSE of µ̂ and µ̂IP W for
two different node features y. When y is correlated with π, then µ̂IP W has a smaller
MSE. When y is weakly correlated with π, then µ̂ has a smaller MSE. For settings in
which µ̂IP W clearly outperforms µ̂, the test statistic from Section 3 rejects the null
hypothesis that the sample average is unbiased (i.e. H0 : Eπ µ̂ = µ).
In the Comm15 network, N = 1089 students from two sister schools were asked to
list up to 10 friends; these friends can be inside or outside of the school. The students
also supplied information including their gender, grade and race. The analysis below
12

Figure 1: Q-Q plots of the sample average (left panel) and the Volz-Heckathorn estimator (right panel) for different social network and sampling designs. For each scenario we draw 2000 network driven samples of size
√ ≈ 500 from a network containing 5,000 nodes.
Here
the
threshold
for
λ
is
1/
2 ≈ 0.707. For the two settings
2
√
with λ2 <
1/
2,
the
distributions
appear
normal.
However,
for the two settings with
√
λ2 > 1/ 2, the distributions do not appear normal. Across all values of λ2 , there is
no apparent difference between the four different designs (i.e. replacement sampling
vs without replacement sampling and 2-tree vs Galton-Watson tree).
studies two node covariates: gender (0/1) and total nominations of friends (integer
between 0 and 10). Before the simulation, the network was symmetrized (i.e. consider
the new adjacency matrix Ã = A + AT ), yielding a network with average node degree
d¯ = 8.06. Because the students were only allowed to list up to 10 friends, the standard
deviation of the degrees is 4.7. This is drastically smaller than typical social networks.
However, even in this setting, the variance of π is sufficiently large to illustrate the
advantages of the sample average.
For both gender and the number of nominations, Table 1 displays (i) the correlation
between π and y, (ii) the bias of the sample average, and (iii) the crossover point. Recall
that the crossover point is the sample size at which the µ̂IP W has a smaller MSE than
µ̂; this calculation is based upon the simulations described below. The table shows that
gender is weakly correlated with π. As such, the sample average has a small bias and
the crossover point is large. Contrast this with the number of total nominations, which
is highly correlated with π. This makes the sample average clearly biased. Because of
this, it has a small crossover point. These two examples illustrate a range of possibilities
in terms of cor(π, y).
Before estimating the crossover points shown in Table 1, we first study the hypothesis test H0 : E µ̂ = µ for both gender and total nominations. To provide a
benchmark, this simulation compares RDS to independent sampling. Let P be the
transition matrix of the simple random walk on the network. The second largest eigenvalue of P is λ2 = 0.93. Let T be a sample from the Galton-Watson process with
0
Eη(σ) = 1.1 < λ−2
2 . For a node covariate y (gender or nominations), let y be the

13

Gender
Total nominations

cor(π, y)
0.11
0.54

Bias of sample average
0.037
1.1

Crossover point (sample size)
> 500
20

Table 1: Gender has a weak correlation with the sample weights. As such, µ̂ has a
small bias. The crossover point shows that if the samples were drawn independently
from the distribution π, then µ̂ has a smaller MSE when n < 500. Total nominations
has a larger correlation and thus a larger bias. When the sample size exceeds twenty,
µ̂IP W outperforms the sample average.
node feature defined in Eq.4 in Section 3. Recall that Eπ (y 0 ) is the bias of the sample
average. We consider the following methods of generating samples Y1 , . . . , Yn and
computing or estimating the variance σ 2 .
1. Y1 , ..., Yn is an independent and identically distributed sample from the normal distribution N (Eπ (y 0 ), varπ (y 0 )) and σ̂ 2 = varπ (y 0 ). Here the variance
is known.
2. Yi = y 0 (Xi ) for all i, where X1 , ..., Xn is an independent and identically distributed sample from the equilibrium distribution of P , and σ̂ 2 = varπ (y 0 ). Here
the variance is known.
3. Yi = y 0 (Xi ) for all i, where X1 , ..., Xn is a sample from the (T, P )-walk P
on G,
n
and σ̂ 2 is the true variance that is only computable in a simulation, var( √1n i=1 Yi ).
4. Yi = y 0 (Xi ) for all i, where X1 , ..., Xn is a sample from the (T, P )-walk on G,
0
\
and σ̂ 2 = var
π (y )G(R̂), the estimator discussed in Section 5.
For a sample of size n, let the rejection region be
n
1 X
Yi | > 1.96},
W = {| √
nσ̂ i=1

with Yi and σ̂ to be defined above. Here, the null hypothesis is rejected at α = 0.05.
For each scenario, Figure 2 plots the power of our test pr(W ) as a function of sample
size. The power under scenario (1) can be calculated exactly and serves as a benchmark
(black line). The power under scenario (2)-(4) is calculated by taking 1000 independent
samples and counting the number of samples that fall in W (red, blue, yellow lines
respectively).
Because scenario (4) underestimates the true variance, this technique is conservative in rejecting H0 and adopting µ̂IP W . For gender, none of the scenarios quickly
reject the null hypothesis. Compare this to the number of total nominations. Here, H0
is rejected even for small sample sizes.
The final figure plots the mean square error of µ̂, µ̂IP W , and µ̂BA ; this last estimator is the bias adjusted estimator from Section 4. This simulation uses scenario (4),
the most realistic of the previous scenarios. After drawing the sample, compute the
following (for both gender and total nominations)
14

1.0

1.0

0.4

Power

0.6

0.8

0.8
0.6
0.4
0.0

0.2

0.2

Power

0

100

200

300

400

500

0

Sample size n

100

200

300

400

500

Sample size n

(a) Gender

(b) Total nominations

Figure 2: Power of test as a function of sample size for two node features (a) gender
and (b) total nominations of friends. The black, red, blue and yellow lines are the power
under scenarios 1-4 respectively.
1. The inverse probability weighted estimator
µ̂IP W =

2. The sample average µ̂ =

1
n

n
n
1 X y(Xi )
1 X y(Xi )d¯
.
=
n i=1 N πXi
n i=1 deg(Xi )

Pn

i=1

y(Xi ).

3. The bias adjusted estimator

µ̂IP W
µ̂BA =
µ̂

if {Xi }1≤i≤n ∈ W,
if {Xi }1≤i≤n ∈
/ W,

introduced in Section 3.
Figure 3 shows that for gender, the true bias of the sample average is small. As
such, the MSE of µ̂ (solid) is always smaller than that of µ̂IP W (dotted) (for sample
sizes n < 500). For total nominations, the bias is much larger. So, when n > 20, the
MSE of µ̂ is larger than the MSE of µ̂IP W . The MSE of the bias adjusted estimator
µ̂BA (longdash) lies between that of µ̂ and µ̂IP W . In particular, when µ̂IP W drastically
outperforms µ̂ (i.e. on the right of panel b), the null hypothesis is typically rejected and
µ̂BA performs similarly to µ̂IP W .

7

Discussion

A recent review of the RDS literature counted over 460 studies which used RDS
[White et al., 2015]. Many of these studies seek to estimate the prevalence of HIV
15

10.0
2.0
0.1

0.5

Mean square error

0.100
0.010
0.001

Mean square error

1

2

5 10

50

200

1

Sample size n

2

5 10

50

200

Sample size n

(a) Gender

(b) Total nominations

Figure 3: The mean square error of µ̂ (solid), µ̂IP W (dotted), and µ̂BA (longdash) for
gender and total nominations.
or other infectious diseases; for these studies, a point estimate of the prevalence is
insufficient. These studies have used confidence intervals constructed from bootstrap
procedures and from estimates of the standard error [Handcock et al., 2016]. These
standard error intervals implicitly rely on a central limit theorem and this paper provides a partial justification for such techniques, so long as m ≤ 1/λ22 . This paper
makes a first step at studying the distributional properties of two simple estimators in
this regime.
Figure 1 suggests that if m is larger than 1/λ22 , then the simple estimators (µ̂ and
µ̂V H ) are no longer normally distributed. Interestingly, under the simulation setting
where the estimators are no longer normally distributed, the Q-Q plots are flatter than
the line x = y. This indicates that a confidence interval constructed from the standard
errors would be conservative; a nominally 90% confidence interval would cover µ more
than 90% of the time. As such, a properly constructed interval should be narrower than
the interval constructed from the standard error. If one pursues this path, then care
must be taken in estimating the standard errors. For example, a bootstrap procedure
proposed in [Salganik, 2006] has become very popular. However, for reasons beyond
the inequality in Proposition 5, this bootstrap procedure drastically underestimates the
actual standard errors [Goel and Salganik, 2010, Rohe, 2015].
There are many reasons to suspect the construction of the sampling weights in RDS
studies. At the most basic level, the justification for
{selection probability for node i} ∝ deg(i)
comes from a Markov model which has several problematic pieces (e.g. replacement
sampling, uniform selection of friends, referral process has reached equilibrium, and
all network relationships are reciprocated). While these assumptions are all troublesome, they are merely sufficient conditions. It is conceivable that deg(i) is still an

16

adequate approximation of the selection probabilities (up to scaling) even when the assumptions do not hold. Perhaps the most difficult problem is that deg(Xi ) is estimated
via self-report. Taken together, there are many reasons to doubt the sampling weights.
In a related context, the third section of the paper discusses the bias-variance tradeoff
between µ̂ and µ̂IP W . The results in Proposition 2 and Theorem 3 presume that the
sampling weights are known exactly. However, given that the presumed model has
several deficiencies, the stationary distribution of the (presumed) Markov process does
not necessarily reveal the actual sampling probabilities. As such, µ̂IP W (and by extension µ̂V H ) are constructed with incorrect and noisy measurements of the sampling
probabilities. This will likely make the estimator biased and more variable. Because of
this, in practice we should be less inclined towards the weighted estimator (i.e. µ̂V H )
than the proposed estimator µ̂BA suggests. Section 4 and the data analysis with the
AddHealth network suggest that the sample average is perhaps less biased than was
previously considered. While there are certainly situations where bias corrected estimators should be used, it also seems sensible to first estimate the bias; when the bias is
large, this is a relatively easy task.
The theorems in this paper do not apply to general trees, only to m-trees. If T is a
Galton-Watson tree with E(η(σ)) < λ−2
2 , then the simulations support the following
conjecture:
1 X
√
(y(Xσ ) − Eπ (y)) → N (0, σ 2 ),
n
σ∈T

2

where σ could be computed from the results in [Rohe, 2015]. To prove this result
requires a more careful study of the structure of {Xσ }σ∈T . We leave this problem to
future investigation.

Acknowledgement
This material is based upon work supported in part by the U. S. Army Research Office under grant number W911NF1510423 and the National Science Foundation under
grant number DMS-1309998. We thank Zoe Russek for helpful comments.

A

Proof of Theorem 1

In the appendix we give a proof of the theorems and propositions in the paper. First we
give an outline of the proof of our main theorem. Consider the martingale
X
E(Yh |Fh−1 ) − Yh ,
h

where {Fh } is a filtration to be defined later. Using the Markov property and the
estimation of var(Yh ), we show that the martingale difference sequence satisfies the
condition of the martingale central limit theorem. In this section, P will be a reversible
transition matrix with eigenvalues 1 P
= λ1 ≥ |λ2 | ≥ ... ≥ |λN | and corresponding
eigenfunctions f1 , ..., fN satisfying k fi (k)fj (k)πk = δij for any i, j. We refer
17

to [Levin et al., 2009] for the existence of such eigen-decomposition. Unless stated
otherwise, expectations are calculated with respect to the tree indexed random walk on
the graph.
We begin with some lemma.
Lemma 1. (Lemma 12.2 in [Levin et al., 2009]) Let P be a reversible Markov transition matrix on the nodes in G with respect to the stationary distribution π. The
eigenvectors of P , denoted as f1 , . . . , fN , are real valued functions of the nodes i ∈ G
and orthonormal with respect to the inner product
X
hfa , fb iπ =
fa (i)fb (i)πi .
(5)
i∈G

If λ is an eigenvalue of P , then |λ| ≤ 1. The eigenfunction f1 corresponding to
the eigenvalue 1 is taken to be taken to be the constant vector 1. If X(0), . . . , X(t)
represent t steps of a Markov chain with transition matrix P , then the probability of a
transition from i ∈ G to j ∈ G in t steps can be written as
P (X(t) = j|X(0) = i) = Pijt = πj + πj

N
X

λt` f` (i)f` (j).

`=2

Lemma 2. For any nodes σ, τ in T,
cov(y(Xσ ), y(Xτ )) =

N
X

d(σ,τ )

λ`

< y, f` >2π ,

l=2

where < y, f` >π =

PN

i=1

y(i)f` (i)πi .

Proof. From the reversibility of the Markov chain and Lemma 1, we have
d(σ,τ )

P (Xσ = j|Xτ = i) = Pij

= πj + πj

N
X

d(σ,τ )

λ`

f` (i)f` (j).

l=2

Therefore,
cov(y(Xσ ), y(Xτ )) =

X

=

X

X
y(i)y(j)πi P (Xσ = j|Xτ = i) − (
πi yi )2

i,j

i

y(i)y(j)πi πj

i,j

=

N
X

N
X

d(σ,τ )

λ`

l=2
d(σ,τ )

λ`

< y, f` >2π ,

l=2

and the lemma is proved.
The next lemma gives the expression of var(Yh ).

18

f` (i)f` (j)

(6)

Lemma 3 (Variance of Yh ). Suppose that |λ2 | > 0.

if
 O(1)
var(Yh ) =
O(h)
if

O((mλ22 )h ) if

Then as h → ∞,
m < λ−2
2
.
m = λ−2
2
m > λ−2
2

Proof. For k = 0, 1, ..., h, denote by shk the number of ordered pairs (σ, τ ) such that
|σ| = |τ | = h and d(σ, τ ) = 2k. Then sh0 = mh , and
shk = mh−k m(m − 1)(mk−1 )2 = mh+k − mh+k−1
for k ≥ 1. By Lemma 2,
var(Yh ) =

h
N
N
h
X
X
X
1 X
2k
2
2
(s
λ
<
y,
f
>
)
=
(<
y,
f
>
(1
+
mk−1 (m − 1)λ2k
hk
`
`
`
π
π
` ))
mh
k=0

l=2

l=2

k=0

N
h
X
X
≤
(< y, f` >2π
(mλ22 )k ).
l=2

k=0

Thus

if m < λ−2
 O(1)
2
.
var(Yh ) =
O(h)
if m = λ−2
2

O((mλ22 )h ) if m > λ−2
2

Corollary 2. For any function y on the state space,

√ 1 Yh
mh

→ 0 in probability.

Proof. It follows from Lemma 3 that
var( √

1
mh

Yh ) = O(λ2h
2 ) → 0.

The next lemma is a convergence argument which we will use in the proof of Theorem 1.
Lemma 4 (Slutsky’s lemma). If Xh → X in distribution and Yh → 0 in probability,
then Xh + Yh → X in distribution.
The following theorem from [Durrett, 2010] is essential to the proof of our main
theorem.
Theorem 4 (Martingale central limit theorem). Suppose that {Zh }h≥1 is adapted to
Ph
the filtration {Fh }h≥1 and that E(Zh+1 |Fh ) = 0 for all h ≥ 1. Let Sh = i=1 Zi
Ph
and Vh = i=1 E(Zi2 |Fi−1 ). If
(1) Vh /h → σ 2 > 0 in probability and
Ph
(2) h−1 i=1 E(Zi2 1{|Zi |>√h} ) → 0 for every  > 0,
√
then Sh / h → N (0, σ 2 ) in distribution.
19

Now we are ready to prove our main theorems.
Proof of Theorem 1. Define Yh in the same way as Theorem
1. Without loss of gen√
erality, √
suppose that Eπ (y) = 0. Since m < λ−2
,
mP
− I is invertible. Let
2
y 0 = ( mP − I)−1 y. Then y 0 is also a function on the state space. We will first
argue on the new node feature y 0 and then convert back to y. Define
Yh0 = √
Let zhk =
and

P

σ:|σ|=h

1
mh

X

y 0 (Xτ ).

τ ∈T:|τ |=h

1{Xσ =k} for h ≥ 0 and k = 1, ..., N . Define zh = {zh1 , ..., zhN },
Fh = σ(Xτ : |τ | ≤ h)

for h ≥ 1. It is obvious that {Yh }h≥1 is adapted to the filtration {Fh }h≥1 . Let
Zh = E(Yh0 |Fh−1 ) − Yh0 .
Then {Zh , Fh }h≥1 is a martingale difference sequence. We will verify that {Zh , Fh }h≥1
satisfies (1) and (2) in Theorem 4.
We have
√
√
z T ( mP )y 0
zT y0
√
− √h .
Zh = mzh−1 P y 0 / mh − Yh0 = h−1
mh−1
mh
For any σ ∈ T, denote by p(σ) the parent node of σ. Zh can also be expressed as
Zh =

X 1p(σ) P y 0 − y 0 (Xσ )
X
y 0 (Xσ )
y 0 (Xσ )
√
|Fh−1 ) − √
=
=
Wσ ,
E( √
mh
mh
mh
σ:|σ|=h
σ:|σ|=h
σ:|σ|=h
X

where 1p(σ) is the 1 × N vector with 1p(σ),i = 1 if Xp(σ) = i and 0 otherwise.
We have
E(Wσ |Fh−1 ) = 0
and

V arpi (y 0 )
mh
for
where pi is the ith row of the transition matrix P and V arpi (y 0 ) =
P i = Xp(σ) , P
2
j pij (y(j) −
j pij y(j)) . From the definition of tree indexed Markov process, if
|σ| = |τ | = h, then Wσ , Wτ are independent given {Xσ : |σ| = h − 1}. Using
E(Wσ |Fh−1 ) = 0, we have
E(Wσ2 |Fh−1 ) =

E(Zh2 |Fh−1 ) =

X

E(Wσ2 |Fh−1 ) =

σ:|σ|=h

20

N
X
zhi
V arpi (y 0 ).
h
m
i=1

zhi
2h
2
From Corollary 2, var( m
h ) = O(λ2 ) → 0 for every i. Thus var(E(Zh |Fh−1 )) =
P
h
2
O(λ2h
2 ) → 0(h → ∞) and
i=1 var(E(Zi |Fi−1 )) converges. It follows from the
definition of Vh and the Cauchy-Schwarz inequality that
h

lim var(Vh /h) = lim var(

h→∞

h→∞

h

1X
1X
E(Zi2 |Fi−1 )) ≤ lim
var(E(Zi2 |Fi−1 )) = 0.
h→∞ h
h i=1
i=1

Therefore
Vh /h → σ 2
in probability, where
2

σ =

E(E(Zh2 |Fh−1 ))

=

N
X

0

πi V arpi (y ) =

i=1

N
X

N
N
X
X
0
2
πi (
pij y (j) − (
pij y 0 (j))2 )

i=1

j=1

j=1

0

0

= varπ (y ) − varπ (P y ),
and condition (1) in Theorem 4 is satisfied.
Similarly, we have
X
E(Zh4 |Fh−1 ) =
E(Wσ4 |Fh−1 ) +

X

E(Wσ2 Wτ2 |Fh−1 )

σ,τ :|σ|=|τ |=h

σ:|σ|=h

mh (mh − 1)
C0
≤ C,
≤ h + C1
m
m2h
where C0 , C1 , C are constants. Thus E(Zh4 ) ≤ C for any h, and
h−1

h
X
i=1

E(Xi2 1{|Xi |>√h} ) ≤ h−1

h
X

E(Xi2

i=1

h
Xi2
1 X
C
)
=
E(Xi4 ) ≤
→ 0.
2 h
2 h2 i=1
h

Condition (2) is also satisfied. From Theorem 4, we have
√
h
h
T
( mP )y 0
1 X zi−1
zT y0
1 X
√
√
Zi = √
(
− √i ) → N (0, σ 2 )
h i=1
h i=1
mi−1
mi
in distribution. If m < λ−2
2 , then from Lemma 3,
Lemma 4 and the definition of y 0 ,

T √
zh
( mP )y 0
√
mh−1

→ 0 in probability. From

√
h
h
1 X ziT ( mP − I)y 0
1 X
√
√
=√
Yi → N (0, σ 2 )
h i=1
h i=1
mi
√
in distribution,
where σ 2 = varπ (y 0 ) − varπ (P y 0 ) = varπ (( mP − I)−1 y) −
√
varπ (P ( mP − I)−1 y). The proof is now complete.

21

B

Proof of Theorem 2 and Corollary 1

We provide a proof of the central limit theorem using the moment method. It involves
a careful study of all the moments of µ̂h . The following proposition is essential to our
proof.
Proposition 6. (Moment continuity theorem) Let Xh be a sequence of uniformly subgaussian real random variables, and let X be another subgaussian random variable.
Then the following statements are equivalent:
(1)EXhk → EX k for all k
(2)Xh → X in distribution
Following this proposition, we can break down our proof to two parts. We will first
prove that all the moments of µ̂h converge to the moments of some normal distribution.
Then we will verify that µ̂h is a uniformly subgaussian sequence.

B.1

Proof of moments convergence

Let Xr be the root of the 2-tree. Define
γk,h (i) = E[µ̂kh |Xr = i],
and
γk,h = E[µ̂kh ].

√
Let ρ = 2|λ2 | < 1. We will prove that there exist γk , k ≥ 1 such that |γk,h (i)−γk | =
O(ρh ) for all k, i, h and that γk = E(ξ k ) for ξ ∼ N (0, γ2 ). Our key observation is that
the left and right subtree can be seen as i.i.d copies of the whole tree given the left and
right child of the seed, which makes it possible to build a relationship between γk,h (i)
and γk,h−1 (i). Only condition (c3) is needed throughout the proof.
We need the following Lemma.
Lemma 5. Let {ah } be a sequence satisfying
ah = ch (cah−1 + C + dh ),
where |1 − ch | = O(ρh ), |dh | = O(ρh ), C is a constant and c < ρ < 1. Then
|ah − C/(1 − c)| = O(ρh ).
Proof. Without loss of generality, suppose that ch 6= 0 and C = 0. Since |1 − ch | =
Qh
O(ρh ) and |dh | = O(ρh ), there exists M such that k=1 |ck | ≤ M and |dh | ≤ M ρh
for all h. Therefore,
|ah | =

h
X
k=0

ch−k dk

h
Y
i=h−k

and the lemma is proved.

22

|ci | ≤

M2 h
ρ ,
ρ−c

We use an induction on k. First, we will prove that γ1 = 0. In fact, from Lemma 1,
E(y(Xσ )|Xr = i) =

N
X

|σ|

y(j)Pij = O(|λ2 ||σ| ).

j=1

Therefore

h
1 X k
2 O(|λ2 |k ) = O(ρh )
E(µ̂h |Xr = i) = √ h
2 k=1

for all i, and |γ1,h (i) − γ1 | = O(ρh ) for γ1 = 0.
Now we move from k − 1 to k. Without loss of generality, suppose that γ2,h (i) > 1
for all h, i (or we can multiply y with a large constant). It follows that γ2k,h (i) ≥
(γ2,h (i))k > 1 for all k. We can decompose γk,h (i) into
X
1
γk,h (i) = E[µ̂kh |Xr = i] − E[( √ h
Xσ )k |Xr = i]
2 σ∈T,0<|σ|≤h
X
1
Xσ )k |Xr = i]
+ E[( √ h
2 σ∈T,0<|σ|≤h

(7)

:= I1 + I2
If k is even, then from Holder’s inequality we know that
i(s)
|I1 | =|E[µ̂kh |Xr = i] − E[(µ̂h − √ )k |Xr = i]|
2h


k
X
√ −h
k
|Xr = i]
=
(−i 2 )m E[µ̂k−m
h
m
m=1
k  
X
√ −h
k
(|i| 2 )m E[µ̂kh |Xr = i]
≤
m
m=1
√ −h
≤[(1 + M 2 )k − 1]E[µ̂kh |Xr = i].

(8)

Likewise, If k is odd, we have
i
|I1 | =E[µ̂kh |Xr = i] − E[(µ̂h − √ )k |Xr = i]
2h
√ −h k
≤[(1 + M 2 ) − 1]E[µ̂k−1
h |Xr = i]

(9)

Since k is fixed, E[µ̂k−1
h |Xr = i] is bounded from our assumption on γk−1,h (i), and
√ −h
√ −h
(1 + M 2 )k − 1 = O( 2 ) = O(ρh ). Hence,
√ −h
h
|I1 | ≤ [(1 + M 2 )k − 1]E[µ̂k−1
h |Xr = i] = O(ρ ).

23

Let Xlc and Xrc be the left and right child of the root and Tl and Tr the left and
right subtree, we have
i
E[(µ̂h − √ )k |Xr = i]
2h
X
X
1
1
1
Xσ + √ h−1
Xσ ))k |Xr = i)
=E(( √ ( √ h−1
2 2
2
σ∈Tr ,|σ|≤h
σ∈Tl ,|σ|≤h
X
X
X
1
1
1
=
piu piv E(( √ ( √ h−1
Xσ + √ h−1
Xσ ))k |Xr = i, Xlc = u, Xrc = v)
2
2
2
u,v
σ∈Tr ,|σ|≤h
σ∈Tl ,|σ|≤h
X
X
X
1
1
1
Xσ + √ h−1
Xσ ))k |Xlc = u, Xrc = v)
=
piu piv E(( √ ( √ h−1
2
2
2
u,v
σ∈Tr ,|σ|≤h
σ∈Tl ,|σ|≤h
X
1
=
piu piv √ k (γk,h−1 (u) + γk,h−1 (v)) + S1
2
u,v
X
1
= √ k−2
piu γk,h−1 (u) + S1 ,
2
u
(10)
where
k−1
X 1 kX
S1 =
piu piv γm,h−1 (u)γk−m,h−1 (v).
√ k m
2
u,v
m=1
If k = 2, Equation 8 and 10 reduce to
X
γk,h (i) =
piu γk,h−1 (u) + δh (i),
u

where
δh (i) =

X
y(i)2
2y(i)
− √ h γ1,h (i) + (
piu γ1,h (u))2 = O(ρh ).
h
2
2
u

Write νh = {γh2 }0 and δh = {δh }0 . For n ≥ 2,
νh = P νh−1 + δh .
Ph
Thus by setting δ1 = 0 we have νh = P h ν1 + k=1 P k δh−k ,and it is not hard
to verify
P∞
that all the components of νh (i.e., every γk,h (i)) converge to γ2 = π t ν1 + h=1 π t δh
with rate ρh .
Now suppose that k > 2. There are a fixed number of terms in S1 . Since |γl,h (i) −
γl | = O(ρh ) for all i ∈ S and l < k − 1, we have
k−1
X

 
k
1
h
|S1 −
√ k m γm γk−m | = O(ρ ).
2
m=1
Thus,
I2 = √

1

X

k−2

2

u

k−1
X

 
1
k
h
piu γk,h−1 (u) +
√ k m γm γk−m + O(ρ ).
2
m=1
24

(11)

Combining Equation 8, 9 and 11 we arrive at the final equation for γk,h (i):
k−1
X

 
k
1
h
√ k m γm γk−m +O(ρ )),
k−2
2
2
u
m=1
(12)
√ −h k
h
where ck,h = 1 if k is odd and ck,h = 1 + O(ρ ) ∈ [2 − (1 + M 2 ) , (1 +
√ −h
Qh
M 2 )k ] if k is even. Since √ 1k−2 < ρ and i=1 ck,h converges, we conclude
2
from Lemma 5 that
|γk,h (i) − γk | = O(ρh ),

Pk−1
k
where γk = m=1 √1 k m
γm γk−m /(1 − √ 1k−2 ).
γk,h (i) = ck,h I2 = ck,h ( √

1

X

piu γk,h−1 (u)+

2

2

We have proved that |γk,h (i) − γk | = O(ρh ) for all k, i, h. Let h tend to infinity in
Equation 11, we have
k−1
X

 
1
k
γk = √ k−2 γk +
√ k m γm γk−m .
2
2
m=1
1

(13)

Now suppose that ξi is a sequence of i.i.d N (0, γ2 ) variables. Let
γ̃k = lim E((
h→∞

ξ1 + · · · + ξh k
√
) ),
h

then {γ̃k }, k ∈ N also follows Equation 13. Since γ1 = γ̃1 = 0 and γ2 = γ̃2 we have
γk = γ̃k for every k, and the argument is proved.

B.2

Proof of uniform subgaussianity

To prove that µ̂h are uniformly subgaussian for all h, we need to show that there exists
some θ such that
γ2`,h (i) ≤ θ2` γ2`
for all `.
c1 be a large constant to be defined later and ch+1 = (1 + M (2λ02 )−h )(1 +
√ Let
0 h
0
( 2λ2 ) )c√
h , where λ2 = max{|λ2 |, 2/3} and M = ||y||∞ . Let s`,h = ||γ`,h (i)||∞ .
0
Since 0 < 2λ2 < 1 < 2λ02 , ch+1 > ch and θ = limh→∞ ch exists. Thus it suffices
to prove that
s2`,h ≤ c2`
(14)
h γ2`
and

√
s2`−1,h ≤ c2`−1
( 2λ02 )h γ2`
h

(15)

for all ` and h.
Again we use an induction on `. Since γ1,h (i) = O(|λ2 |h ), we can choose c1 large
enough such that the inequalities in Equation 14 and 15 hold for all (h, `) with h = 1
or ` = 1. Suppose that 14 and 15 are verified for all ` ≤ k. We will prove that they are
also true for ` = k + 1.

25

From condition (c1) and (c2), we know that
X
||
piu γ2k+1,h (u)||∞ ≤ λ2 ||γ2k+1,h (i)||∞ = λ2 s2k+1,h .
u

From our assumption of induction we have
√

√

(| 2λ02 |)h
c2k+1
h−1 √ 2k+2 (



2k + 1
s2k+1,h ≤ [(1 + M 2 )
γ2k+2
−
+
0
2
 



k 
X
2k + 1
2k + 1
2k + 1
+
(
+
)γ2k+2−2m γ2m +
γ2k+2 )
2m − 1
2m
2k + 1
m=1
√ −h
2k+1
= [(1 + M 2 )2k+1 − 1]c2k
h−1 γ2k + ch−1
√





k 
X
(| 2λ02 |)h 2k + 2
2k + 2
2k + 2
γ2k+2 +
γ2k+2−2m γ2m +
γ2k+2 )
√ 2k+2 (
0
2m
2k + 2
2
m=1
√ −h
√ 0 h
2k+1
= [(1 + M 2 )2k+1 − 1]c2k
h−1 γ2k + ch−1 (| 2λ2 |) γ2k+2
√ −h
√
√ 0 h
≤ ([(1 + M 2 )2k+1 − 1](| 2λ02 |)−h + 1)c2k+1
h−1 (| 2λ2 |) γ2k+2
√ 0 h
≤ (1 + M (2λ02 )−h )2k+1 c2k+1
h−1 (| 2λ2 |) γ2k+2
√
≤ c2k+1
(| 2λ02 |)h γ2k+2 ,
h
(16)
and Equation 15 is true for 2k + 1.
Now we move from 2k + 1 to 2k + 2. Recall that
−h 2k+1

2k+2
X



2k + 2 X
piu piv γm,h−1 (u)γ2k+2−m,h−1 (v))
√ 2k+2
m
2
u,v
1

γ2k+2,h (i) = c2k+2,h (

m=0
−h 2k+2

≤(1 + M 2

)

2k+2
X

(

√

m=0

1]c2k
h−1 γ2k



1
2

2k+2


2k + 2
sm,h−1 s2k+2−m,h−1 ).
m

Thus,
−h 2k+2

s2k+2,h ≤ (1 + M 2

)

2k+2
X

(

m=0

Let
I1 =

k+1
X

√

m=0

and
I2 =

k
X
m=0

√

2k+2

2

1
2k+2

2



1



√

1
2

2k+2




2k + 2
sm,h−1 s2k+2−m,h−1 ).
m


2k + 2
s2m,h−1 s2k+2−2m,h−1
2m


2k + 2
s2m+1,h−1 s2k+1−2m,h−1 .
2m + 1

Then s2k+2,h ≤ (1 + M 2−h )2k+2 (I1 + I2 ).
26

We have
k+1
X


2k + 2 2m
2k+2−2m
ch−1 γ2m ch−1
γ2k+2−2m
√ 2k+2
2m
2
m=0


k+1
X
2k + 2
1
γ2m γ2k+2−2m
= c2k+2
√ 2k+2
h−1
2m
2
m=0

I1 ≤

1



(17)

= c2k+2
h−1 γ2k+2 ,
where the last equality follows from Equation 13. On the other hand
k
X


√
2k + 2 2m+1 √ 0 h−1
γ2m+2 c2k+1−2m
( 2λ02 )h−1 γ2k+2−2m
√ 2k+2 2m + 1 ch−1 ( 2λ2 )
h−1
2
m=0


k
√ 0 2(h−1) X
1
2k + 2
2k+2
= ch−1 ( 2λ2 )
√ 2k+2 2m + 1 γ2m+2 γ2k+2−2m .
2
m=0
(18)
It can be directly verified that for all m,


 
2k + 2
k
γ2m+2 γ2k+2−2m ≤ 2(k + 1)
γ2k+2 .
2m + 1
m
I2 ≤

1



Thus,
 
k
√ 0 h−1 X
1
k
I2 ≤ c2k+2
(
2λ
)
2(k
+
1)
γ2k+2
√ 2k+2
2
h−1
m
2
m=0
√ 0 2(h−1)
(k + 1)γ2k+2
= c2k+2
h−1 ( 2λ2 )
√ 0 h
2k+2
≤ ch−1 ( 2λ2 ) 2kγ2k+2 .

(19)

Combining Equation 17 and 19, we have
√ 0 h
I1 + I2 ≤ c2k+2
h−1 γ2k+2 (1 + 2k ∗ ( 2λ2 ) )
√ 0 h 2k+2
≤ c2k+2
.
h−1 γ2k+2 (1 + ( 2λ2 ) )

(20)

Therefore,
sh,2k+2 ≤ c2k+2,h (I1 + I2 )

√
−h 2k+2
≤ c2k+2
)
(1 + ( 2λ02 )h )2k+2 γ2k+2
h−1 (1 + M 2
≤

(21)

c2k+2
γ2k+2 ,
h

and the theorem is proved.

B.3

Proof of Corollary 1

Proof. By Theorem 2 and Slutsky’s lemma, it suffices to prove that dˆ¯ → d¯ in proba1
bility. Let D = max1≤i≤N deg(Xi ). For any σ ∈ T, Eπ deg(X
= d1¯. Thus E 1ˆ¯ = d1¯,
σ)
d

27

ˆ¯ d¯ ≤ D, we have
and it follows from Theorem 2 that dˆ¯ → d¯ in probability. Since d,
¯ − d|
¯ > ) ≤ P (| 1 − 1¯| > /D2 ) → 0 for all  > 0, and the corollary is proved.
P (|dˆ
ˆ
¯
d
d

C

Proof of Proposition 2 and Theorem 3

Proof of Proposition 2. Since y(1), ..., y(N ) are uncorrelated, we have

PN
µ2 i=1 πi2 if i 6= j
E(y(Xi )y(Xj )) =
.
µ2
if i = j
Therefore,
N
N
X
X
µ2
µ2
n−1
− 1+
(µ2
πi2 + µ21 (1 −
πi2 ) − µ21 ).
n
n
n
i=1
i=1

var(µ̂) =
Similarly, we have

var(µ̂IP W ) =
=
Let

y(X1 )
n−1
y(X1 ) y(X2 )
1
var(
)+
cov(
)
n
N πX1
n
N πX1 N πX2
N
µ2 X 1
µ2
n − 1 µ2
1
− 1+
(
+ µ21 (1 − ) − µ21 ).
2
n i=1 N πi
n
n
N
N

V D = var(µ̂IP W ) − var(µ̂)
=

N
N
X
1
µ2 X 1
n−1
(
var(y)(
−
πi2 ).
−
1)
+
n i=1 N 2 πi
n
N
i=1

Since max1≤i≤N N πi = C1 , N 2 πi πj ≤ C12 for all i, j. Therefore
√
√
N
N N
N N
X
πj
πi 2
1
1 XX
1 X X (πi − πj )2
(
−
1
=
−
)
=
√
√
N 2 πi
2 i=1 j=1 N πi
N πj
2 i=1 j=1 N 2 πi πj
i=1
N N
1 XX
(πi − πj )2 .
≥
2C12 i=1 j=1

(22)

(23)

We have
N N
N N
N X
N
N
X
X
1 XX
1 XX 2
(πi − πj )2 =
(πi + πj2 ) −
πi πj = N
πi2 − 1
2 i=1 j=1
2 i=1 j=1
(24)
i=1 j=1
i=1

= N 2 var(π).
From Eq.22, 23 and 24,
VD ≥(

µ2 N 2
− N var(y))var(π),
C12 n

and the proposition is proved.
28

Proof of Theorem 3. For the (T, P )-walk on G, the variance of µ̂IP W and µ̂ is given
by
var(µ̂IP W ) =

N
N
d(σ,τ )
X
µ2 X 1
µ21
1 X
pii
var(y)
−
+
n i=1 N 2 πi
n
n2
N 2 πi
i=1
σ6=τ

and
N
X
µ2
µ21
1 X
d(σ,τ )
var(µ̂) =
−
+ 2
var(y)
πi pii
.
n
n
n
i=1
σ6=τ

Let
V D = var(µ̂IP W ) − var(µ̂)
=

N
N
d(σ,τ )
XX
1
pii
µ2 X 1
d(σ,τ )
(
−
1)
+
var(y)
(
− πi pii
).
2π
n i=1 N 2 πi
n2
N
i
i=1

(25)

σ6=τ

Since C1 N ≤ di ≤ C2 N , we have πi ≤
Thus,

C2
N C1

(t)

for all i and pij ≤

1
C1 N

for all i, j, t.

N
N
d(σ,τ )
1 XX
C2
1 X X pii
d(σ,τ )
d(σ,τ )
(
−
π
p
)
≥
πi pii
≥
.
i ii
2π
2
n2
N
n
N
C12
i
i=1
i=1
σ6=τ

σ6=τ

From equation 23 and 24, we know that
N
X
i=1

1
N 2 C12
var(π).
−
1
≥
N 2 πi
C22

Therefore,
V D ≥ µ2 (

N 2 C12
C2
var(π) −
).
nC22
N C12

We conclude the proof by taking C =

D

C23
.
C3 C14

Proofs of Propositions 4 and 5

Before the proofs, some more notation is necessary. Let |λ1 | ≥ |λ2 | ≥ · · · ≥ |λN |
denote the eigenvalues of P . Denote the f1 , . . . , fN : V → R as the corresponding the
eigenfunctions of P .
The following is a proof of Proposition 5.
Proof. Proposition 1 in [Khabbazian et al., 2016] says that
Cov(ỹ(Xp(τ ) ), ỹ(Xτ )) =

N
X
hỹ, f` i2π λ` .
`=2

29

Thus, from the lemma and the definition of R, R is a convex combination of the eigenvalues
N
X
hỹ, f` i2π
R=
λ` .
varπ (y)
`=2

Lemma 6.

N
X
hỹ, f` i2π
= 1.
varπ (y)
`=2

The proof of Lemma 6 is given on page 342 of [Levin et al., 2009] and is repeated
below for completeness.
From Theorem 1 in [Rohe, 2015],
σµ̃2

N
X
=
hỹ, f` i2π G(λ` ).

(26)

`=2

Applying Jensen’s inequality yields
σµ̃2 = varπ (ỹ)

N
X
hỹ, f` i2π
G(λ` ) ≥ varπ (ỹ)G(R).
varπ (ỹ)

(27)

`=2

The proof of Proposition 4 is similar.
Proof. In the case when ỹ(i) = µ + σf (i), this implies that ỹ = µf1 + σfj . By
hỹ,f i2

π
the orthonormality of the eigenvectors, varπ`(ỹ)
= 1{j = `} for ` > 1. As such, the
inequality in equation (27) holds with equality.

The following is a proof of Lemma 6
Proof.
varπ (ỹ) = Eπ (ỹ 2 ) − (Eπ ỹ)2 =

N
N
X
X
hỹ, fj i2π − (Eπ ỹ)2 =
hỹ, fj i2π .
j=1

j=2

Proposition 5 presumes that G is convex. Figure 4 plots G for twenty different Galton-Watson trees with offspring distribution p(0) = .1, p(1) = .1, p(2) =
.3, p(3) = .5. This offspring distribution has expected value 2.2. The construction
of each tree was stopped when it reached 5000 nodes; if it failed to reach 5000 nodes,
then the process was started over. In these simulations and in others not shown, G is
often convex. When it is not convex, its second derivative is positive away from −1.
This simulation was selected because it shows that even when the trees are sampled
from the same distribution, even when there is nothing strange about the offspring distribution (e.g. all moments are finite), even when it is a very big tree, even under all of
these nice conditions some of the trees have a convex G and some of the trees have a
non-convex G. Similar results hold when the trees have 500 nodes; the only thing that
changes is that the red regions extend slightly further away from −1.
30

Some G functions are not convex.
However, they are all convex away from −1.
1e+00

●

1e−01
1e−02
1e−03
1e−04

●

−1.0

−0.5

0.0

0.5

1.0

Figure 4: Each line corresponds to a function G for a random Galton-Watson tree.
Some of the curves are not convex; the regions of these functions which have a negative
second derivative are highlighted in red. While some of the black lines appear to have
a negative second derivative, they do not; this illusion is due to the log transformation
on the vertical axis.

References
[WHO, 2013] (2013).
Introduction
To
HIV/AIDS
And
Sexually
Transmitted
Infection
Surveillance
Module
4:
Introduction
to
Respondent-drive Sampling.
World Health Organization & UNAIDS.
http://applications.emro.who.int/dsaf/EMRPUB 2013 EN 1539.pdf.
[Benjamini and Peres, 1994] Benjamini, I. and Peres, Y. (1994). Markov chains indexed by trees. The Annals of Probability, pages 219–243.
[Biemer and Christ, 2007] Biemer, P. P. and Christ, S. (2007). Weighting survey data.
International Handbook of Survey Methodology.
[Bollen et al., 2016] Bollen, K. A., Biemer, P. P., Karr, A. F., Tueller, S., and Berzofsky, M. E. (2016). Are survey weights needed? a review of diagnostic tests in
regression analysis. The Annual Review of Statistics and Its Application, 3(1):375–
392.
[Clauset et al., 2009] Clauset, A., Shalizi, C. R., and Newman, M. E. (2009). Powerlaw distributions in empirical data. SIAM review, 51(4):661–703.

31

[Durrett, 2010] Durrett, R. (2010). Probability: theory and examples. Edition 4.1.
Cambridge university press.
[Gile, 2011] Gile, K. J. (2011). Improved inference for respondent-driven sampling
data with application to hiv prevalence estimation. Journal of the American Statistical Association, 106(493).
[Goel and Salganik, 2009] Goel, S. and Salganik, M. J. (2009). Respondent-driven
sampling as markov chain monte carlo. Statistics in medicine, 28(17):2202–2229.
[Goel and Salganik, 2010] Goel, S. and Salganik, M. J. (2010). Assessing respondentdriven sampling. Proceedings of the National Academy of Sciences, 107(15):6743–
6747.
[Handcock et al., 2016] Handcock, M. S., Fellows, I. E., and Gile, K. J. (2016). RDS:
Respondent-Driven Sampling. Los Angeles, CA. R package version 0.7-5.
[Heckathorn, 1997] Heckathorn, D. D. (1997). Respondent-driven sampling: a new
approach to the study of hidden populations. Social problems, pages 174–199.
[Holland et al., 1983] Holland, P., Laskey, K., and Leinhardt, S. (1983). Stochastic
blockmodels: First steps. Social Networks, 5(2):109–137.
[Jones et al., 2004] Jones, G. L. et al. (2004). On the markov chain central limit theorem. Probability surveys, 1(299-320):5–1.
[Khabbazian et al., 2016] Khabbazian, M., Hanlon, B., Russek, Z., and Rohe, K.
(2016). Novel sampling design for respondent-driven sampling. arXiv preprint
arXiv:1606.00387.
[Levin et al., 2009] Levin, D. A., Peres, Y., and Wilmer, E. L. (2009). Markov chains
and mixing times. American Mathematical Soc.
[Lorrain and White, 1971] Lorrain, F. and White, H. C. (1971). Structural equivalence
of individuals in social networks. The Journal of mathematical sociology, 1(1):49–
80.
[Lu et al., 2012] Lu, X., Bengtsson, L., Britton, T., Camitz, M., Kim, B. J., Thorson,
A., and Liljeros, F. (2012). The sensitivity of respondent-driven sampling. Journal
of the Royal Statistical Society: Series A (Statistics in Society), 175(1):191–216.
[Pfeffermann, 1996] Pfeffermann, D. (1996). The use of sampling weights for survey
data analysis. Statistical Methods in Medical Research, 5(3):239–61.
[Rohe, 2015] Rohe, K. (2015). Network driven sampling; a critical threshold for design effects. arXiv preprint arXiv:1505.05461.
[Rohe et al., 2011] Rohe, K., Yu, B., and Chatterjee, S. (2011). Spectral clustering and
the high dimensional stochastic blockmodel. The Annals of Statistics, 39(4):1878–
1915.

32

[Salganik, 2006] Salganik, M. J. (2006). Variance estimation, design effects, and sample size calculations for respondent-driven sampling. Journal of Urban Health,
83(1):98–112.
[Salganik and Heckathorn, 2004] Salganik, M. J. and Heckathorn, D. D. (2004). Sampling and estimation in hidden populations using respondent-driven sampling. Sociological methodology, 34(1):193–240.
[Strogatz, 2001] Strogatz, S. H. (2001).
410(6825):268–276.

Exploring complex networks.

Nature,

[Valliant et al., 2013] Valliant, R., Dever, J. A., and Kreuter, F. (2013). Practical Tools
for Designing and Weighting Survey Samples. Springer New York.
[Volz and Heckathorn, 2008] Volz, E. and Heckathorn, D. D. (2008). Probability
based estimation theory for respondent driven sampling. Journal of Official Statistics, 24(1):79.
[White et al., 2015] White, R. G., Hakim, A. J., Salganik, M. J., Spiller, M. W.,
Johnston, L. G., Kerr, L., Kendall, C., Drake, A., Wilson, D., Orroth, K., et al.
(2015). Strengthening the reporting of observational studies in epidemiology for
respondent-driven sampling studies:?strobe-rds? statement. Journal of clinical epidemiology, 68(12):1463–1471.

33

