Quantifying population contact patterns in the United States during the COVID-19 pandemic

Survey methodology

We designed and fielded a survey to measure interpersonal interaction in the United States. Following the POLYMOD project (Mossong et al. 2008) and subsequent studies (e.g., Zhang et al. 2020; Jarvis et al. 2020; Dorélien et al. n.d.), survey respondents were asked to report the number of people they had conversational contact with on the day before the interview; in Waves 1 and 2, we also asked about physical contact. Respondents were asked to provide detailed information about up to three of their reported contacts; this detailed information included who those contacts were, how long those contacts lasted, and where they took place. In Wave 0, respondents were asked to report all contacts, and to then report how many of their contacts were not household members. In Waves 1 and 2, respondents were asked to provide a household roster, and then report only contacts outside of the household.

Respondents were recruited using Lucid, an online panel provider. In each wave, we obtained two samples: first, a quota sample that is intended to be representative of the United States; and, second, several smaller quota samples from specific cities: New York, the San Francisco Bay Area, Atlanta, Phoenix, and Boston. In Wave 1, Philadelphia was added.

The project was approved by the UC Berkeley IRB (Protocol 2020-03-13128).

Weighting

City samples

For our city-specific samples, we recruited respondents who lived in the Designated Market Area (DMA) surrounding each city. DMAs are intended to capture media markets, and therefore often include much more than just the urban core of a city.

Our city-specific samples recruited respondents in the DMAs associated with six cities: Atlanta, the San Francisco Bay Area, Boston, New York, Phoenix and, starting in Wave 1, Philadelphia. Figure 6 shows the media markets.

• Download figure
• Open in new tab

Figure 6:

The geographical areas corresponding to the cities in our sample.

Statistical Model

To investigate factors associated with interpersonal contacts, we developed statistical models. We fit separate models to the total reported contacts and to the number of reported non-household contacts. In each case, we model the expected number of contacts using a negative-binomial distribution. The negative binomial distribution is appealing because it allows for overdispersion – that is, it enables us to model count data that exhibit more variance than would be expected under a Poisson distribution. This modeling approach has previously been used to study contact data (e.g., Mossong et al. 2008).

In our model, the log of the expected number of contacts for respondent i is given by where X_i is a vector of covariates that includes age category, gender, household size, survey wave, city, ethnicity (Black / White/ Other), Hispanicity, and whether or not the day being reported about is a weekday. β is a vector of coefficients to be estimated.

Given μ_i, we define λ_i = exp(μ_i) to be the expected number of contacts for respondent i. Then we model the reported number of contacts for respondent i, y_i, as where ϕ ∈ [1, ∞) is a shape parameter that is inversely related to overdispersion; that is, the higher ϕ is estimated to be, the more similar y_i’s distribution is to a Poisson distribution with rate parameter λ_i.

In our data, observations from Wave 0 are censored above 10, because the survey instrument allowed respondents to report up to ‘10 or more’ contacts. Wave 1 and Wave 2 allowed respondents to enter any number of contacts, but we topcoded contacts at 29, following Mossong et al. (2008). Reports that are topcoded or censored in any of the three waves are treated as right-censored in the model. We adopt a Bayesian approach to fitting the model. For all of the regression coefficients β, we assume flat priors. For the intercept and the shape parameter, we assume very weak priors. Specifically, we assume a priori that the intercept α is distributed with mean 0 and a large variance by using pr(α) ~ Student-t(3, 0,10); And we assume a priori that the the shape parameter ϕ, is distributed with mean 1 and a very large variance by using pr(ϕ) ~ Gamma(.01, .01). Figure 2 shows posterior estimates for the mean and 95% credible intervals for key coefficients of interest. Table 1, shows summary statistics for the predictors used in our model.

Accounting for this censoring, in our models the log posterior of the parameters given the data, logpr(α, β, ϕ|y, X), is proportional to where s_nc is the set of responses that are not censored; s_c is the set of responses that are right-censored, with response i ∈ s_c being censored at value c_i; is the PMF of the negative binomial distribution, and F_NB is the cumulative distribution function ; and is the expected number of contacts or non-household contacts for respondent i. (The parameterizations here are the ones used in stan). For each model, we run 4 chains of the sampler; each chain was run for 1000 warmup iterations and then 1000 sampling iterations. All R-hat statistics are 1, suggesting that the chains mixed effectively.

Epidemiological Model

We estimate age-structured contact matrices for each wave of the BICS study using previously described methods (Wallinga, Teunis, and Kretzschmar 2006; Mossong et al. 2008; Klepac, Kissler, and Gog 2018; Jarvis et al. 2020). We group respondents and their contacts into six age bins: 0-18, 18-25, 25-35, 35-45, 45-65, 65+. For each age group, we estimate the average daily number of contacts reported by respondents in that age group with contacts in every age group. In other words, our raw contact matrix, M, has entries m_ij which is the average number of daily contacts between respondents in age group, j, with their reported contacts in age group, i. Adjusting for survey weights, we calculate m_ij as: where w_t,j is the weight for particpant, t, in age group j, and y_i is the weight for their reported contact in age group i. T_j is the total number of respondents in age group, j.

Contacts in the population must be reciprocial but due to differences in reporting in the survey our raw social contact matrix, M, is not. We impose reciprocity by: where c_ij are the entries of the reciprocal contact matrix, C, and N_i and N_j the population size in age classs i and j, respectively. For the youngest age group, for which we have no survey respondents, we assume:

We estimate within age group average number of contacts, c_ii, for the youngest age group using previously developed methods (Klepac, Kissler, and Gog 2018; Jarvis et al. 2020) and data from the United Kingdom POLYMOD study (Mossong et al. 2008). For each wave of the BICS study, we calculate the ratio of the dominant eigenvalue for the contact matrix estimated from the BICS data to the dominant eigenvalue of the contact matrix from the POLYMOD study, with school contacts removed to reflect current school closures, for all age groups that are overlapping between the two studies. The within age group average number of contacts for the [0,18) group in the POLYMOD study is then scaled by this ratio to impute c_[0,18)[0,18) in the BICS contact matrix.

The transmission dynamics of infectious diseases is summarized by the next generation matrix, N, that determines how an infection spreads when a pathogen is first introduced into a fully susceptible population. The basic reproduction number, R₀, is the average number of secondary infections arising from a single infection in a fully susceptible population, and is typically estimated as the spectral radius (dominant eigenvalue), ρ(N) of the next generation matrix, N (Farrington, Kanaan, and Gay 2001). The N matrix is proportional to the population contact matrix, C. The exact relationship between N and C is model- dependent, but for respiratory pathogens such as SARS-CoV-2, N is typically modeled as C scaled by the duration of infectiousness, , and the probability of transmission for a single contact, q. Therefore, the spectral radius of N: where ρ(C) is the dominant eigenvalue of the reicprocal population contact matrix. In other words, R₀ is proportional to the dominant eigenvalue of C.

Since R₀ is proportional to the dominant eigenvalue of C, relative differences in R₀ under different contact patterns is equivalent to the the ratios of the dominant eigenvalues of the different contact matrices. Specifically, if we assume that contact patterns in the population before physical distancing became widespread are equivalent to a baseline contact matrix, and that disease-specific parameters remained unchanged over the course of the survey period, the relative reduction in R₀ during physical distancing, compared to the baseline, is equivalent to the ratios of the dominant eigenvalues of the C matrices from the BICS study, C^BICS, to the dominant eigenvalue of the baseline pre-pandemic contact matrix C^baseline:

Further, if we assume a distribution for R₀ for COVID-19 in the absence of physical distancing, we can estimate the implied R₀ during the study period, by multiplying this ratio with the R₀ value in the absence of physical distancing. We assume that R₀ prior to physical distancing followed a normal distribution with mean 2.5 and standard deviation of 0.54 based on estimates from literature (Anderson et al. 2020; Jarvis et al. 2020). We compare the BICS contact matrices to two baseline business-as-usual scenarios: contact patterns estimated by (Feehan and Cobb 2019), from a probability sample of US Facebook users, and contact patterns from the UK POLYMOD study (Mossong et al. 2008), that has been widely used in many settings. We compute confidence intervals for the estimated R₀ by repeating the age-imputation and relative R₀ estimation on 5,000 bootstrapped samples from the BICS, POLYMOD and 2015 study contact matrices.

Bootstraped contact matrices

Uncertainty estimates for the descriptive results - including Figure 1 and the R₀ analysis summarized in Figure 4 - are based on the bootstrap. To obtain bootstrap resamples, we resampled respondents separately in each stratum (i.e., in each city and in the national sample). For each bootstrap resample, we first drew n_c samples with replacement from among the n_c respondents in city c. For resampled respondents who provided detailed reports about all contacts (i.e., those respondents who had a_i = 1 for all contacts i), we used the reported detailed contacts without resampling. For resampled respondents who did not provide detailed reports about all contacts, but who reported about a subset of contacts (i.e., for each respondent who had a_i > 1 for some contact i), we resampled r_i contacts with replacement from among the respondent’s r_i detailed contacts. This second stage is intended to capture sampling variation due to the respondent choosing which contacts to report about from among her total contacts. Using this approach, we obtained 5,000 bootstrap resamples of our dataset, and these bootstrap resamples are the basis for uncertainty estimates.

Relationships and locations

Here, we present Tables containing the values shown in Figure 1, Panels C and D. These numerical values may be useful as numerical inputs for future modeling work.

Sensitivity analysis for contact definition

In Wave 0 (the pilot study), contact was defined using ‘conversational contact’, as explained in this text:

By in-person conversational contact, we mean a two-way conversation with three or more words in the physical presence of another person.

You might have conversational contact with family members, friends, co-workers, store clerks, bus drivers, and so forth.

(Please do not count people you contacted exclusively by telephone, text, or online. Only consider people you interacted with face-to-face.)

After the pilot study, in Waves 1 and 2, the survey instrument was modified and contact was defined using this text:

Now we would like to ask you some questions about people you had in-person contact with yesterday.

By in-person contact, we mean

EITHER a two-way conversation with three or more words in the physical presence of another person OR physical skin-to-skin contact (for example a handshake, hug, kiss, or contact sports)

You might have in-person contact with family members, friends, co-workers, store clerks, bus drivers, and so forth.

Please do not count people you contacted exclusively by telephone, text, or online. Only consider people you interacted with face-to-face.

In the main text, for Waves 1 and 2 we combine physical and conversational contacts together, since the combination of these two is most relevant for disease transmission. However, 9 percent of non-household detailed contacts in Waves 1 and 2 were reported to be strictly physical (and not conversational). To assess how sensitive our results are to including all contacts in our analysis, we repeated key analyses excluding these contacts reported to be only physical in Waves 1 and 2. Figures 7 and 8 suggest that this decision does not substantively affect our conclusions.

• Download figure
• Open in new tab

Figure 7:

Estimated average number of contacts each person reported to have taken place by contact’s relationship (top panel) and location (bottom panel). Uncertainty estimates are 95% intervals derived from the bootstrap. Contacts are restricted to those who were conversational contacts for Waves 1 and 2.

• Download figure
• Open in new tab

Figure 8:

Estimated R₀ for each survey wave. Contacts are restricted to those who were conversational contacts for Waves 1 and 2.