Show the code
Complex Survey Notes
Preface
Overview
What follows is a coded work through of Thomas Lumley’s “Complex Surveys: A Guide to Analysis in R” (Lumley 2011) and some of the methods underlying design based statistics more broadly. This write-up reflects my understanding of the material with additions made to try and clarify ideas further. These include simulations, derivations or references that I found helpful in working through Lumley’s material. Most of the data sets that I use throughout these notes are from Lumley’s website for the book. Some data sets are not available, or at least I was not able to find them, or the the documentation for the data was not available. I remark on this when relevant in my notes below.
How to use these notes?
I’d imagine there are two general ways one might use these notes:
As a quick reference on any of the code or topics.
If you are working through the book, you may find these notes useful to double check Lumley’s code or your understanding of the principles tested in the exercises.
- Note that Lumley’s code was written at the time of publishing and much has changed, both in the survey package and R more generally. I try to keep this as close to his original code as possible while updating when necessary. I also often show how to produce the same code using the
srvyrpackage, but this is less relevant for the later chapters.- I would not always use the same method if I was writing this kind of code but I try to replicate his to make it easy for any reader with the book to follow along.
- My answers to the exercises are not guaranteed to be correct. I would strongly encourage you to avoid copying them if you are, for example, working through this textbook for a class.
How shouldn’t I use these notes?
These notes are no substitute for buying the book itself. I would encourage you to buy the book if you’re intent on working through these notes.
Drawing Samples in R
It wasn’t very long into reading this book that I found that weighted sampling in R is, unfortunately, not well set-up for complex designs. It would not be possible, for example, to simply use the base R sample or popular tidyverse package dplyr’s function slice_sample() to draw a weighted sample from a population for example with appropriate inclusion probabilities.Further details are in this stats exchange post.
Instead a function from the sampling package would have to be used. I use this function below in any setting where a non-uniform sample with inclusion probabilities is needed.
Its worth further pointing out that the topic of how samples themselves are drawn is a complicated one its own right and that the functions in the sampling package each have pros and cons according to the target estimand or question of interest. Drawing samples with replicate weights – discussed in Chapter 2 – is a similarly complex question which I haven’t yet resolved to my satisfaction.
Chapter 1: The Basics
Design vs. Model
This book focuses on “Design-based” Inference. That is the methods in this book focus on the design from which the data are constructed, rather than the data itself. In a traditional survey setting the data are assumed to be fixed and the probabilities of sampling different entities are used to derive the desired estimate. These inclusion probabilities or their inverse, “sampling weights”, are used to re-balance the data so that they more accurately reflect the target population distribution. Different sampling techniques — clustering, 2-phase, etc. — are used to either decrease the variance of the resulting estimate, the cost associated with the design or both.
Horvitz Thompson Estimation
The Horvitz Thompson Estimator (HTE) is the starting point for non-uniform random estimates. If we observe measure X_i on subject i drawn with probability \pi_i from a population of N total subjects the HTE is formulated as follows:
HTE(X) = \sum_{i=1}^{N} \frac{1}{\pi_i}X_i, \tag{1.1}
which is an unbiased estimator as shown in the Chapter 1 Appendix.
The variance of this estimate is:
V[HTE(X)] = \sum_{i,j} \left ( \frac{X_i X_j}{\pi_{ij}} - \frac{X_i}{\pi_i} \frac{X_j}{\pi_j} \right ), \tag{1.2}
which follows from the Bernoulli covariance using indicator variables R_i=1 if individual i is in the sample, R_i=0 otherwise. A proof is provided in the [Questions From Chapter 1].
Design And Misspecification Effects
(Kish 1965) defined the notion of a design effect as the ratio of a variance of an estimate in a complex sample to the variance of the same estimate in a simple random sample (SRS). The motivation for this entity being that it can guide researchers in terms of how much sample size they may need; If the sample size for a given level of precision is known for a simple random sample, the sample size for a complex design can be obtained by multiplying by the design effect.
While larger sample sizes may be necessary to maintain the same level of variance as a SRS, the more complex may still be more justified because of the lower cost associated. See (Meng 2018) for an example of where design effects are used in a modern statistical setting by comparing competing estimators.
Other preliminary items
From this point Lumley works through an introduction to the datasets used in the book and the idea that we’ll often be taking samples from datasets where we know the “true” population and computing estimates from there. This isn’t always the case and there may some subtlety worth discussing how to interpret results once we get into topics like regression, but for the most part his description makes sense.
One thing I found lacking in this introductory section is the motivation for why we might take non-uniform samples. It isn’t until Chapter 3 that Lumley discusses probability proportional to size (PPS) sampling, but this is very often the reason why a non-uniform sample is used.
If we have some measure that is right skewed in our population of interest and we’d like to estimate the mean, we could take a SRS to estimate the mean but the variance on that item would be lower than if we sampled proportional to the right skew measure itself. I’ll demonstrate with the following quick example, suppose we want to measure the income of a population. Incomes are often right skewed, but we can get a lower variance estimate if we take a weighted sample.
I generate a right skewed population and visualize the distribution.
Here I’ll take a uniform and weighted sample of size 50. Note that the differences in the samples are subtle. They might not look all that different on visual inspection.
Show the code
uniform_sample <- population %>%
slice_sample(n = 50) %>%
transmute(
income = income,
method = "uniform",
pi = 50 / 1E3
)
weighted_sample <- population %>%
mutate(
pi = sampling::inclusionprobabilities(floor(population$income), 50),
in_weighted_sample = sampling::UPbrewer(pi) == 1
) %>%
filter(in_weighted_sample) %>%
transmute(
income = income,
pi = pi,
method = "weighted"
)
rbind(
uniform_sample,
weighted_sample
) %>%
ggplot(aes(x = income, fill = method)) +
geom_histogram() +
ggtitle("Sample Comparisons") +
xlab("Income (USD)") +
theme(legend.title = element_blank(), legend.position = "top")
Finally I’ll estimate the population mean from both samples and include the design effect calculation in the weighted sample estimate.
Show the code
uniform_sample %>%
as_survey_design() %>%
summarize(
mean_income = survey_mean(income)
)# A tibble: 1 × 2
mean_income mean_income_se
<dbl> <dbl>
1 56609. 3918.Show the code
weighted_sample %>%
as_survey_design(probs = pi) %>%
summarize(mean_income = survey_mean(income, deff = TRUE))# A tibble: 1 × 3
mean_income mean_income_se mean_income_deff
<dbl> <dbl> <dbl>
1 50685. 3266. 0.890We see that the weighted estimate standard error is not quite half the uniform estimate. Accordingly the design effect for the weighted sample is less than 1.
Exercises
Doesn’t make sense to reproduce here.
Don’t make sense to reproduce here.
Each visit to the front page of a newspaper’s website has (independently) a 1/1000 chance of resulting in a questionnaire on voting intentions in a forthcoming election. Assuming that everyone who is given the questionnaire responds, why are the results not a probability sample of:
- Voters?
- Readers of the newspaper?
- Readers of the newspaper’s online version?
Lumley lists 4 properties needed for a sample to be considered a probability sample.
- Every individual (unit of analysis) in the population must have a non-zero probability of ending up in the sample (\pi_i>0 \forall i)
- \pi_i must be known for every individual who does end up in the sample.
- Every pair of individuals in the sample must have a non-zero probability of both ending up in the sample (\pi_{i,j} \forall i, j)
- The probability \pi_{i,j} must be known for every pair that does end up in the sample.
- is not guaranteed when considering voters — there are voters who don’t read the paper who have will have \pi_i = 0 — or the broader heading of “readers” of the newspaper - since those who only read the physical paper will have a $_i = 0 $. For “readers of the newspaper’s online version” the sample would only be a probability sample if the time window was further specified, as there could be online readers who do not visit during the survey window, and would thus be assigned a \pi_i=0.
- You are conducting a survey that will estimate the proportion of women who used anti-malarial insecticide-treated bed nets every night during their last pregnancy. With a simple random sample you would need to recruit 50 women in any sub-population where you wanted a standard error of less than 5 percentage points in the estimate. You are using a sampling design that has given design effects of 2-3 for proportions in previous studies in similar areas.
- Will you need a larger or smaller sample size than 50 for a sub-population to get the desired precision?
Larger, a design effect >1 indicates that the variance is larger in the complex design with the same sample size - consequently the sample size will need to be increased to maintain the same level of precision.
- Approximately what sample size will you need to get the desired precision?
100 - 150. Derived from multiplying 50 by 2 and 3.
- Systematic sampling involves taking a list of the population and choosing, for example, every 100th entry in the list.
- Which of the necessary properties of a probability sample does this have?
Items ii-iv from the list enumerated above. The only condition that is not satisfied is that not every item has a nonzero probability of being chosen.
- For systematic sampling with a random start, the procedure would be to choose a random starting point from 1, 2, …, 1000 and then take every 100th entry starting at the random point. Which of the necessary properties of a probability sample does this procedure have?
This satisfies all items from the above list.
- For systematic sampling with multiple random starts we might choose 5 random starting points in 1, 2, ….., 5000 and then take every 500th entry starting from each of the 5 random points. Which of the necessary properties of a probability sample does this procedure have?
Again, this satisfies all items from the above list.
- If the list were shuffled into a random order before a systematic sample was taken, which of the properties would the procedure have.
Again, all of them. The key is adding the known randomness and not excluding any items from selection.
- Treating a systematic sample as if it were a simple random sample often gives good results. Why would this be true?
This would be because the items are not ordered in any particular fashion prior to taking the “systematic sampling”. In this setting a systematic sample is equivalent to a simple random sample.
- Why must all the sampling probabilities be non-zero to get a valid population estimate?
If any of the sampling probabilities are zero, that would introduce bias in shifting the estimate away from the portion of the population that would always be unobserved under repeated sampling.
- Why must all the pairwise probabilities be non-zero to get a valid uncertainty estimate.
This is basically a second order statement equivalent to the previous. If any pair is unable to be observed together that is a form of selection bias that would shift the sample estimate away from the true population value.
- A probability design assumes that people who are sampled will actually be included in the sample, rather than refusing. Look up the response rates for the most recent year of BRFSS and NHANES.
Lumley is highlighting the fact that even though we set up samples thinking that every sample will be observed that is rarely the case. Looking at just the most recent NHANES data I see response rates at ~ 78% for the un-weighted, 6-year household survey.
- In a telephone study using random digit dialing, telephone numbers are sampled with equal probability from a list. When a household is recruited, why is it necessary to ask how many telephones are in the household, and what should be done with this information in computing the weights.
It is necessary to ask how many telephones are in the household to down weight the a priori sampling probability accordingly because every additional telephone line increases the odds that a given house is sampled. For example in a simple population with two houses, where house one has 5 telephones and house two has 2 telephones, and we’re looking to take a n=1 sample, but we don’t know the number of telephones a priori, house one has a \frac{5}{7} probability of being sampled. If that is the house that is chosen its weight needs to go from 2 to \frac{5}{7} to better reflect its sampling probability. In a real sample this would be corrected relative to all the other households number of telephones or perhaps a population average of the number of telephones.
Derive the Horvitz Thompson variance estimator for the total as follows.
- Write R_i = 1 if individual i is in the sample, R_i=0 otherwise. Show that V[R_i] = \pi_i(1-\pi_i) and that Cov[R_i,R_j]=\pi_{ij} - \pi_i\pi_j.
This follows in a straightforward fashion from the assumption that R_i is distributed according to the Bernoulli distribution and R_i \perp R_j. This is an accurate model for sampling with replacement, or sampling from large populations with small sample sizes without replacement, but less true for small sample sizes without replacement.
- Show that the variance of the Horvitz Thompson estimator is:
V[\hat{T}_{HT}] = \sum_{i=1}^{N}\sum_{j=1}^{N} \check{x}_i\check{x}_j(\pi_{ij} - \pi_i \pi_j) We have, \hat{T}_{HT} := \sum_{i=1}^{N} \frac{X_i I(X_i \in {S})}{\pi_i} \\ V[\hat{T}_{HT}] = V\left[\sum_{i=1}^{N}\frac{X_i I(X_i \in {S})}{\pi_i} \right] \\ =\sum_{i=1}^{N}\sum_{j=1}^{N} Cov\left[\frac{X_i I(i \in {S})}{\pi_i}, \frac{X_j I(j \in {S})}{\pi_j}\right]\\ = \sum_{i=1}^{N}\sum_{j=1}^{N} \frac{X_i}{\pi_i}\frac{X_j}{\pi_j}Cov(I(i \in {S}),I(j \in {S})) \\ = \sum_{i=1}^{N}\sum_{j=1}^{N} \frac{X_i}{\pi_i}\frac{X_j}{\pi_j} (\pi_{ij} - \pi_i \pi_j)
which is equivalent to the above, where \check{x_i} = \frac{X_i}{\pi_i}.
- Show that an unbiased estimator of the variance is \hat{V}[\hat{T}_{HT}] = \sum_{i=1}^{N}\sum_{j=1}^{N} \frac{R_i R_j}{\pi_{ij}}\check{x_i}\check{x_j}(\pi_{ij} - \pi_i \pi_j)
To show the expression above is unbiased for \hat{V}[\hat{T}_{HT}] we must show that E\left [\hat{V}[\hat{T}_{HT}] \right] = V[\hat{T}_{HT}] E \left [\hat{V}[\hat{T}_{HT} ] \right] = E \left[\sum_{i=1}^{N}\sum_{j=1}^{N} \frac{R_iR_j}{\pi_{ij}} \check{x}_i \check{x}_j(\pi_{ij} - \pi_i\pi_j) \right] \\ = \sum_{i=1}^{N} \sum_{j=1}^{N} \frac{E[R_iR_j]}{\pi_{ij}} \check{x}_i \check{x}_j (\pi_{ij} - \pi_i \pi_j) \\ = \sum_{i=1}^{N} \sum_{j=1}^{N} \frac{\pi_{ij}}{\pi_{ij}} \check{x}_i\check{x}_j(\pi_{ij} - \pi_i \pi_j) \\ = \sum_{i=1}^{N} \sum_{j=1}^{N} \check{x}_i\check{x}_j(\pi_{ij} - \pi_i \pi_j) \\ \blacksquare
- Show that the previous expression simplifies to equation 1.2
\sum_{i=1}^{N} \sum_{j=1}^{N} \frac{R_iR_jx_ix_j}{\pi_{ij}\pi_i\pi_j}(\pi_{ij} - \pi_i \pi_j) \\ = \sum_{i=1}^{n} \sum_{j=1}^{n} \frac{x_ix_j}{\pi_{ij}\pi_i\pi_j}(\pi_{ij} - \pi_i \pi_j) \\ = \sum_{i=1}^{n} \sum_{j=1}^{n} x_i x_j(\frac{1}{\pi_i \pi_j} - \frac{1}{\pi_{ij}}) \\ = \sum_{i=1}^{n} \sum_{j=1}^{n} \frac{x_ix_j}{\pi_i \pi_j} - \frac{x_i x_j}{\pi_{ij}}\
I’m not sure how the signs switch on the last line to reproduce expression 1.2
- Another popular way to write the Horvitz-Thompson variance estimator is
\hat{V}[\hat{T}_{HT}] = \sum_{i=1}^{n} x_i^{2} \frac{1-\pi_i}{\pi_i^2} + \sum_{i\neq j}x_ix_j\frac{\pi_{ij} - \pi_i \pi_j}{\pi_i\pi_j\pi_{ij}}
Show that this is equivalent to equation 1.2
We need to show that the above is equivalent to \sum_{i,j} \frac{X_iX_j}{\pi_{ij}} - \frac{X_i}{\pi_i}\frac{X_j}{\pi_j}
First we fix i \neq j in the above expression and we find \sum_{i\neq j} \frac{X_iX_j}{\pi_{ij}} - \frac{X_i}{\pi_i}\frac{X_j}{\pi_j} = \sum_{i\neq j} X_iX_j(\frac{1}{\pi_{ij}} - \frac{1}{\pi_i} \frac{1}{\pi_j}) \\ = \sum_{i \neq j} X_iX_j(\frac{\pi_i \pi_j - \pi_{ij}}{\pi_{ij}\pi_i\pi_j})
which is the latter part in the desired expression save for a sign, which again I must be missing somehow or is an error in the book.
Now we take i=j and return to expression 1.2 in which we have,
\sum_{i=j} \frac{X_i X_j}{\pi_{ij}} - \frac{X_i}{\pi_i} \frac{X_j}{\pi_j} \\ = \sum_{i=j} \frac{X_i^2}{\pi_{ii}} - \frac{X_i}{\pi_i} \frac{X_i}{\pi_i} \\ = \sum_{i=j} X_i^2 (\frac{1}{\pi_{ii}} - \frac{1}{\pi_i^2} ) \\ = \sum_{i=j} X_i^2 (\frac{\pi_i^2 - \pi_{ii}}{\pi_i^2 \pi_{ii}}) Clearly we have to formulate \pi_{ii} in terms of \pi_i but isn’t immediately clear to me how to do so. We know that for the two terms to be equal we must have
\frac{\pi_i^2 - \pi_{ii}}{\pi_i^2 \pi_{ii}} = \frac{1 - \pi_i^2}{\pi_i^2} \\ \iff \\ \pi_{ii} = \frac{\pi_i^2}{2- \pi_i^2} Which I suppose we’ll take to be the expression of a co-inclusion probability of an entity sampled with itself (this must assume sampling with replacement) for this expression to be true.
Chapter 1 Appendix
The HTE is an unbiased estimator of the population total - I reproduce the expression from above, but now make explicit the indicator variables that express which observations are included in our sample, S.
HTE := \sum_{i=1}^{N} \frac{X_i I(X_i \in S)}{\pi_i} \\ E[HTE] = E\left [\sum_n \frac{X_i I(X_i \in S)}{\pi_i} \right ] \\ = \sum_n E \left [\frac{X_iI(X_i \in S)}{\pi_i} \right ] \\ = \sum_n \frac{X_iE[I(X_i \in S)]}{\pi_i} \\ = \sum_n \frac{X_i \pi_i}{\pi_i} = \sum_n X_i
Chapter 2: Simple and Stratified Sampling
Starting from Simple Random Samples
When dealing with a sample of size n from a population of size N the HTE of the total value of X_i in the population can be written as
\begin{equation} HTE(X) = \hat{T_X} = \sum_{i=1}^{n} \frac{X_i}{\pi_i}. \end{equation}
For a simple random sample, the variance can be more explicitly written as
\begin{equation} V[\hat{T_X}] = \frac{N-n}{N} \times N^{2} \times \frac{V[X]}{n}, \end{equation}
where \frac{N-n}{N} is the finite population correction factor. This factor is derived from the hypergeometric distribution and explains the reduction in uncertainty that follows from sampling a large portion of the population. Consequently, if the sample is taken with replacement — the same individual or unit has the possibility to be sampled twice — this term is no longer relevant. It should be noted that sampling with replacement is not usually used however, but sometimes this language is used to refer to the fact that the finite correction factor may not be used.
The second term, N^2, rescales the estimate from the mean to the total, while the final term is simply the scaled variance of X.
A point worth deliberating on, that Lumley notes as well, is that while the above equations suggest that a larger sample size is always better that is not always the case in reality. Non-response bias or the cost of surveys can dramatically diminish the quality of the dataset, even if the size is large. I state this is worth deliberating on because it is a matter of increasing importance in the world of “Big Data” - where it can be easy to delude oneself with confidence in their estimates because their sample is large, even when the sample is not well designed. See (Meng 2018) for a larger discussion of this topic.
It follows from the above that the HTE for the population size is defined as \hat{N} = \sum_{i=1}^{n} \frac{1}{\pi_i}. This holds true in the case where, as here \pi_i = \frac{n}{N}, a bit trivial, but also in those where \pi_i may be defined differently.
Confidence Intervals
The sampling distribution for the estimates — typically sample means and sums — across “repeated surveys” is Normal by the Central Limit Theorem, so the typical \bar{x} \pm 1.96 \sqrt{\frac{\sigma^2_X}{n}}, expression is used to calculate a 95% confidence interval. Lumley offers the following example from the California Academic Performance Index (API) dataset to illustrate this idea.
Show the code
data(api)
mn_enroll <- mean(apipop$enroll, na.rm = TRUE)
p1 <- apipop %>%
ggplot(aes(x = enroll)) +
geom_histogram() +
xlab("Student Enrollment") +
geom_vline(xintercept = mn_enroll, linetype = 2, color = "red") +
ggtitle("Distribution of School Enrollment")
p2 <- replicate(n = 1000, {
apipop %>%
sample_n(200) %>%
pull(enroll) %>%
mean(., na.rm = TRUE)
})
mn_sample_mn <- mean(p2)
p2 <- tibble(sample_ix = 1:1000, sample_mean = p2) %>%
ggplot(aes(x = sample_mean)) +
geom_histogram() +
xlab("Student Enrollment Averages") +
geom_vline(
xintercept = mn_sample_mn,
linetype = 2, color = "red"
) +
ggtitle("Distribution of Sample Means")
p1 + p2