Vaccination with BNT162b2 reduces transmission of SARS-CoV-2 to household contacts in Israel

Setting

Vaccination in Israel began on December 20, 2020, mainly using the BioNTech-Pfizer BNT162b2 vaccine, with a few individuals receiving the vaccine earlier. The vaccination campaign first targeted high-risk individuals, including those 60 years of age and older, medical personnel, workers at nursing homes, and individuals with comorbidities. After this phase, which lasted until January 21, age restrictions were lowered. By February 6, 2021, every Israeli citizen above 16 years old was eligible for the vaccine [10]. By the beginning of April 2021, 61% of the population had received at least one dose of the BNT162b2 vaccine [15].

The vaccination roll-out coincided with Israel’s third and largest wave of SARS-CoV-2 registered cases [16]. Consequently, a third national lockdown was issued in Israel starting December 24, 2020, with more severe restrictions (e.g., schools closures) issued starting January 8. These restrictions were progressively lifted starting on February 7, 2021.

Data sources

We used data from Maccabi Healthcare Services (MHS) centralized computerized database, which captures all data on members’ healthcare-related interactions (including demographics, inpatient and outpatient visits, diagnoses, procedures etc). MHS is a nationwide 2.5 million-member state-mandated, not-for-profit sick fund in Israel, representing a quarter of the Israeli population, and is a representative sample of the Israeli population. The individual-level data for cases and household contacts include demographic information (i.e. age, sex), date of any polymerase chain reaction (PCR) tests for SARS-CoV-2 and the result of the test (considering that all such tests of MHS members are recorded centrally), and date of receipt of the first and second doses of the vaccine (if received). Individuals were defined as unvaccinated if they had not received any doses of BNT162b2 and fully vaccinated if at least 10 days had passed since receiving the second dose of the vaccine.

Due to computational constraints, we focused on households with at least one infected individual and two or more members; it was not possible to include all households with no infections. Sensitivity analyses were conducted using a randomly-selected subset of households with no infections to evaluate possible biases of this approach. We restricted our analysis to data from June 15, 2020 to March 24, 2021, since viral testing was not widely available prior to this date.

Data augmentation for inferring transmission rates

For each individual, we observed the date at which the viral test was performed and the outcome of the test, but we do not have data on date of infection. To infer transmission rates, it is necessary to estimate when each individual within a household was infected and the period when each person was infectious. We therefore used a data augmentation approach to impute when a person with a positive PCR test was infected and infectious. This was accomplished using random samples from three different Gamma distributions representing the delay between onset of infectiousness and the date of the PCR test (i.e., τ_report), between the date of infection and the onset of infectiousness (i.e., latent period, τ_latent), and the time from the onset of infectiousness to the end of infectiousness (i.e., infectious period, τ_infections). The values for these distributions were based on prior knowledge derived from observational studies on the latent period, times to seeking a test, and the duration of infectiousness (table S2) [17, 18]. For each individual, a random draw was taken from each of these distributions, and this process was repeated 100 times. For clarity, we refer to these Gamma distributions as the delay distributions.

For each person j in household i with a positive PCR test, let T_ij be the (imputed) time of infection in days since the beginning of the study:, with being the number of days after the start of the study until the PCR test date. We set the beginning of the study as May 29, 2020, since we allow for infections occurring up to 17 days prior to the start of the data (i.e., June 15, 2020). If no infection occurred for person j, T_ij is censored and equal to the total number of days in the study (i.e., t_max). For the purposes of model fitting, we define Y_ijt as a binary variable equal to zero for each day of the study up until the time of infection (t = 1, …, T_ij − 1), equal to one for T_ij (assuming the person is infected), and censored from that point onwards (i.e., that person is relevant only in terms of transmitting to other household members). For people who are never infected, Y_ijt is equal to zero for all days of the study.

Statistical modeling

Using the augmented data, we developed two discrete time-to-event data models of household transmission to estimate vaccine effectiveness against susceptibility to infection and against infectiousness given infection. In both approaches, we model the infection status for person j in household i on study day t (i.e., Y_ijt) using conditionally independent Bernoulli distributions with corresponding probability of infection π_ijt. These probabilities are then defined based on personal demographics, community risk, vaccination status, and characteristics of household transmission, with the approaches differing in how transmission is described. Both models were fit 100 times with different draws from the delay distributions to assess uncertainty in the results due to uncertainty in the unobserved dates of infection and infectiousness.

Household transmission model

For the primary transmission model, we define the probability of infection on a given day as where n_i is the number of members in household i, p_ij0t is the probability that person j in household i is infected by the community on study day t (i.e., community risk of infection), p_ijkt is the probability that person j is infected by household member k (i.e., household risk of infection), and d_ikt is an indicator of whether person k can transmit to j on day t. In other words, , with and equal to the time of onset and end of infectiousness, respectively, for person k, and 1(.) representing the indicator function taking the value of one if the input condition is true and the value of zero otherwise.

The probability that individual j never tested positive for SARS-CoV-2 is given as whereas the probability that individual j is infected on day t^* (before the end of the study) is given as (i.e., the probability of escaping infection up to time (t^* − 1) multiplied by the probability of not escaping infection at time t^*).

Conveniently (with respect to computation), the likelihood function can be written in terms of the introduced binary variables such as with N being the total number of households.

We define the per-person, per-day community risk of infection using the logit link function as where δ₀ is the baseline risk of infection from the community, vax_ijt is a binary variable equal to one if at least 10 days has passed since person j received the second dose of the vaccine, age_1,ij is a binary variable equal to one if the person is between 10 and 60 years old at the start of the study (reference category is the ≤10-year-old age group), age_2,ij is similarly defined for those aged ≥60 years old, and cases_t describes the time-varying risk from the community and is computed as the standardized number of positive PCR tests on day t in the data.

Similarly, the per-person, per-day risk of transmission from an infectious individual k to a susceptible household member j is defined (for k≠j) as where α₀ is the baseline risk of infection from an infected household member and vax_ikt is the vaccination status of household member k. All other terms have been previously described.

Vaccine effectiveness is expressed as a percentage and computed as 100*(1-RR), with RR defined as a risk ratio comparing vaccinated and unvaccinated individuals respectively. Vaccine effectiveness against susceptibility stratified by age group was defined as with a_j being a categorical variable equal to 0 for ages ≤10-year-old, 1 for ages between 10 and 60, and 2 for ages ≥60 years old; vax_ijt and vax_ikt are as previously defined. The explicit equation for age group 0 would then be Including γ₁ in both the numerators and denominators of the VE_S provides the VE with respect to age groups 1 (VE_S,1) and similarly with γ₂ for the VE with respect to age group 2 (VE_S,2).

Vaccine effectiveness against infectiousness given infection is defined as which is based on the vaccination status of infected household member k. Again, for age group 0, this is equivalent to similarly, with γ₁ and γ₂ included for the other age groups.

We also estimated the vaccine effectiveness against transmission as (i.e., the reduction in risk associated with vaccine-derived protection against both infection and infectiousness given infection). This is equivalent to The age-adjusted vaccine effectiveness measures can be derived as , with prop(a_j) being the fraction of people in each age group.

We initially applied the household transmission model to only households with a single occupant to inform the values of the baseline risk (δ₀) and the time-varying risk from the community (δ₁). We found δ₀=-8.25 and δ₁=0.79 and used these values as initial conditions for our analysis on households with at least one infected individual and two or more household members.

Infection-hazard regression model

We compared the primary household transmission model defined by equations (1-3) with an alternative infection-hazard model. The models differ in their definition of the probability of infection, π_ijt. We note that similar to the primary transmission model, the likelihood can be defined using the introduced binary variables, Y_ijt. In this analysis, we use the complementary log-log-link function to connect the probabilities to individual-, household-, and community-level risk factors, such that and where vax_ijt, age_1,ij, and age_2,ij are as previously defined, λ₀ is the intercept parameter, h_t is a smooth function of study time (modeled using splines) that describes the time-varying community risk, m_ijt is a binary variable equal to 1 if at least one other unvaccinated household member is infectious during study time t (not including person j), and z_ijt is a binary variable equal to 1 if at least one other vaccinated household member is infectious (not including person j). Interactions between person j’s vaccination status and the household risk variables are also included.

Use of the complementary log-log-link function leads to a hazard ratio interpretation for the exponentiated regression parameters. We define multiple susceptibility vaccine effects using this model output based on comparing different within-individual and within-household scenarios.

First, we compute as the vaccine effectiveness against susceptibility to infection given there are no other infectious household members at that time (i.e., m_ijt = z_ijt = 0) and interpret it as the decrease in the hazard of infection due to vaccination. Similarly, we define and as the vaccine effectiveness against susceptibility to infection given there is at least one unvaccinated (VE_S,u) or vaccinated (VE_S,v) member in the household.

Vaccine effectiveness against infectiousness given the infection of household contacts is defined as and which can be interpreted as the percent reduction in the hazard of infection for individual j when exposed to a vaccinated versus unvaccinated infectious household member k, given individual j is unvaccinated (VE_I,u) or vaccinated (VE_I,v).

For both models, we summarized the vaccine effectiveness estimates by taking the mean over the 100 samples of the delay distributions. We derived the 95% confidence intervals (CI) using the law of total variance. All analyses were carried out in the R statistical software [19].

Sensitivity analysis: including the first vaccine dose

For the primary household transmission model, we included information on the first vaccine dose status by defining community risk as where vax_1,ijt, and vax_2,ijtrepresent mutually exclusive binary variables of vaccination status (i.e., individual j in household i is partially or fully vaccinated at time t, respectively), and all other terms have been previously described.

Similarly, household risk is defined (for k≠j) as: where vax_1,ikt is the first-dose vaccination status of household member k at time t and vax_2,ikt is the second-dose vaccination status of household member k at time t. Vaccine effectiveness estimates are the same as for the primary analysis.

For the infection-hazard model, we now define the variable ψ_ijt from equations 4-5 as where , and represent mutually exclusive binary variables of vaccination level (i.e., partially and fully vaccinated, respectively), m_0,ijt is a binary variable equal to one if at least one other unvaccinated household member is infectious during study time t (not including person j), z_1,ijt is a binary variable equal to one if at least one other partially vaccinated household member is infectious (not including person j), and z_2,ijt is a binary variable equal to one if at least one other fully vaccinated household member is infectious (not including person j). Interactions between person j’s vaccination status and household risk are also included.

The vaccine effectiveness against susceptibility to infection given there are no other infectious household members at that time (i.e., m_0,ijt = z_1,ijt = z_2,ijt = 0) is given as: and is interpreted as the decrease in the hazard of infection due to full vaccination.

Similarly, we define as the vaccine effectiveness against susceptibility to infection given there is at least one unvaccinated (VE_S,u), partially vaccinated (VE_S,pv), or fully vaccinated (VE_S,v) member in the household.

Vaccine effectiveness against infectiousness given the infection of household contacts is defined as which can be interpreted as the effect of being exposed to a fully vaccinated versus unvaccinated infectious household member given individual j is unvaccinated (VE_I,u), partially vaccinated (VE_I,pv), or fully vaccinated (VE_I,v). Results from both models are reported in tables S5 and S6.

Sensitivity analysis: including a subset of households with no infections

We randomly selected 10,000 households with at least two household members and no detected infections and included them along with our original set of households with at least two household members and at least one infection. We compared the vaccine effectiveness results from both models for one iteration of the delay distributions (tables S9 and S10).

Sensitivity analysis: testing the robustness of the results to misclassification of cases

We run a sensitivity analysis to test the robustness of the results to misclassification of individuals who were infected but did not receive a SARS-CoV-2 test. For each individual with a negative PCR test, we randomly selected a new PCR test date. On this new PCR test date, each individual could now have a positive PCR test based on a Bernoulli distribution with probability p = α * prob_test. We define prob_test = (total number of positive PCR cases in the dataset)/(total population in the dataset) = 0.08. With this new dataset, we ran the delay distribution process for one iteration, and we estimated the vaccine effectiveness from both models. We tested two scenarios: (a) α = 0.01 and (b) α = 0.10 (tables S7 and S8).