Introduction

Within epidemiology, the importance of heterogeneity, whether that host, population, statistical, or environmental, has long been recognized1,2,3,4,5. For example, when designing targeted interventions, it is crucial to understand and account for differences that may exist within populations6,7,8. These differences can present in a variety of forms: heterogeneity in susceptibility, transmission, response to guidance, and treatment effects etc.; all of which affect the dynamics of an infectious disease1,2,6,9,10,11,12,13,14. While heterogeneity may exist on a continuous spectrum, it can be difficult to incorporate into analysis and interpretation, so individuals are often placed in discrete groups according to a characteristic that aims to represent the true differences15,16,17,18,19. When examining optimal influenza vaccination policy in the United Kingdom, Baguelin et al.20 classified individuals within one of seven age groups. Explicitly accounting for, and grouping, individuals by whether they inject drugs can help target interventions to reduce human immunodeficiency virus (HIV) and Hepatitis C Virus incidence21. Similarly, epidemiological models have demonstrated the potential for HIV pre-exposure prophylaxis to reduce racial disparities in HIV incidence22. Therefore, heterogeneity can be used to inform more complete theories of change, increasing intervention effectiveness23.

When discretizing a population for the purposes of inclusion within a mechanistic model, three properties need to be defined: (1) the number of groups, (2) the size of the groups, and (3) the differences between the groups. Typically, as seen in the examples above, demographic data is used e.g., age, sex, race, ethnicity, socio-economic status, etc., often in conjunction with the contact patterns and rates7,9,15,17,20,22,24. There are several reasons for this: the data is widely available, and therefore can be applied almost universally; it is easily understandable; and there are clear demarcations of the groups, addressing properties 1) and 2). However, epidemiological models often aim to assess the effects of heterogeneity with respect to infection, e.g., “how does an individual’s risk tolerance affect their risk of infection for influenza?”. When addressing questions such as these, demographic data does not necessarily provide a direct link between the discretization method and the heterogeneous nature of the exposure and outcome, particularly if behavioral mechanisms are a potential driver. Instead, it relies on assumptions and proxy measures e.g., an individual’s age approximates their contact rates, which in turn approximates their risk of transmission. This paper demonstrates an alternative approach to discretizing populations for use within mechanistic models, highlighting the benefits of an interdisciplinary approach to characterize heterogeneity in a manner more closely related to the risk of infection.

In early 2020, shortly after the World Health Organization (WHO) declared the SARS-CoV-2 outbreak a public health emergency of international concern25, universities across the United States began to close their campuses and accommodations, shifting to remote instruction26,27. By Fall 2020, academic institutions transitioned to a hybrid working environment (in-person and online), requiring students to return to campuses28,29,30. In a prior paper31 we documented the results of a large prospective serosurvey conducted in State College, home to The Pennsylvania State University (PSU) University Park (UP) campus. We examined the effect of 35,000 returning students (representing a nearly 20% increase in the county population32) on the community infection rates, testing serum for the presence of anti-Spike Receptor Binding Domain (S/RBD) IgG, indicating prior exposure33. Despite widespread concern that campus re-openings would lead to substantial increases in surrounding community infections28,34,35, very little sustained transmission was observed between the two geographically coincident populations31.

Given the high infection rate observed among the student body (30.4% seroprevalence), coupled with the substantial heterogeneity in infection rates between the two populations, we hypothesized that there may be further variation in exposure within the student body, resulting from behavioral heterogeneity. Despite extensive messaging campaigns conducted by the University36, it is unlikely that all students equally adhered to public health guidance regarding SARS-CoV-2 transmission prevention. We use students’ responses to the behavioral survey to determine and classify individuals based on their intention to adhere to public health measures (PHMs). We then show that these latent classes are correlated with SARS-CoV-2 seroprevalence. Finally, we parameterize a mechanistic model of disease transmission within and between these groups, and explore the impact of public health guidance campaigns, such as those conducted at PSU36. We show that interventions designed to increase student compliance with PHMs would likely reduce overall transmission, but the relatively high initial compliance limits the scope for improvement via PHM adherence alone.

Methods

Design, setting, and participants

This research was conducted with PSU Institutional Review Board approval and in accordance with the Declaration of Helsinki, and informed consent was obtained for all participants. The student population has been described in detail previously31, but in brief, students were eligible for the student cohort if they were: ≥ 18 years old; fluent in English; capable of providing their own consent; residing in Centre County at the time of recruitment (October 2020) with the intention to stay through April 2021; and officially enrolled as PSU UP students for the Fall 2020 term. Upon enrollment, students completed a behavioral survey in REDCap37 to assess adherence and attitudes towards public health guidance, such as attendance at gatherings, travel patterns, and nonpharmaceutical interventions. Shortly after, they were scheduled for a clinic visit where blood samples were collected. Students were recruited via word-of-mouth and cold-emails.

Outcomes

The primary outcome was the presence of S/RBD IgG antibodies, measured using an indirect isotype-specific (IgG) screening ELISA developed at PSU38. An optical density (absorbance at 450 nm) higher than six standard deviations above the mean of 100 pre-SARS-CoV-2 samples collected in November 2019, determined a threshold value of 0.169 for a positive result. Comparison against virus neutralization assays and RT-PCR returned sensitivities of 98% and 90%, and specificities of 96% and 100%, respectively39. Further details in the Supplement of the previous paper31.

Statistical methods

To identify behavioral risk classes, we fit a range of latent class analysis (LCA) models (two to seven class models) to the student’s behavioral survey responses, using the poLCA package40 in the R programming language, version 4.3.3 (2024-02-29)41. We considered their answers regarding the frequency with which they intended to engage in the following behaviors to be a priori indicators of behavioral risk tolerance: wash hands with soap and water for at least 20s; wear a mask in public; avoid touching their face with unwashed hands; cover cough and sneeze; stay home when ill; seek medical attention when experiencing symptoms and call in advance; stay at least 6 feet (about 2 arms lengths) from other people when outside of their home; and, stay out of crowded places and avoid mass gatherings of more than 25 people. The behavioral survey collected responses on the Likert scale of: Never, Rarely, Sometimes, Most of the time, and Always. For all PHMs, Always and Most of the time accounted for > 80% of responses (with the exception of intention to stay out of crowded places and avoid mass gatherings, where Always and Most of the time accounted for 78.8% of responses). To reduce the parameter space of the LCA and minimize overfitting, the behavioral responses were recoded as Always and Not Always. Measures of SARS-CoV-2 exposure e.g., IgG status, were not included in the LCA model fitting, as they reflect the outcome of interest. We focused on responses regarding intention to follow behaviors because this information can be feasibly collected during a public health campaign for a novel or emerging outbreak; it has also been shown that intentions are well-correlated with actual behaviors for coronavirus disease 2019 (COVID-19) public health guidelines, as well as actions that have short-term benefits42,43. We examined the latent class models using Bayesian Information Criterion, which is a commonly recommended as part of LCA model evaluation44,45, to select the model that represented the best balance between parsimony and maximal likelihood fit.

Using the best-fit LCA model, we performed multivariate logistic regression of modal class assignment against IgG seropositivity to assess the association between the latent classes and infection. This “three-step” approach is recommended over the “one-step” LCA model fit that includes the outcome of interest as a covariate in the LCA model45,46. The following variables were determined a priori to be potential risk factors for exposure31: close proximity (6 feet or less) to an individual who tested positive for SARS-CoV-2; close proximity to an individual showing key COVID-19 symptoms (fever, cough, shortness of breath); lives in University housing; ate in a restaurant in the past 7 days; ate in a dining hall in the past 7 days; only ate in their room/apartment in the past 7 days; travelled in the 3 months prior to returning to campus; and travelled since returning to campus for the Fall term. Variables relating to attending gatherings were not included in the logistic regression due to overlap with intention variables of the initial LCA fit. Missing variables were deemed “Missing At Random” and imputed using the mice package47, as described in the supplement of the previous paper31. To examine the effect of modal class assignment, we computed the GLM results of 100 simulations where class assignment was drawn according to each individual’s probability of class membership, resulting from the LCA fitting process. The percentage of simulations where class membership p-values ≤ 0.05 were computed, and the mice package was used to produce pooled odds ratios and associated 95% confidence intervals.

We parameterized a deterministic compartmental Susceptible-Infected-Recovered (SIR) model using Approximate Bayesian Computation (ABC) against the seroprevalence within each latent class. The recovery rate was set to 8 days (Supplemental Table 1). Diagonal values of the transmission matrix were constrained such that βHH ≤ βMM ≤ βLL (H represents high-adherence to public health guidelines, and M and L represent medium- and low-adherence, respectively), with the following parameters fit: the transmission matrix diagonals, a scaling factor for the off-diagonal values (ϕ), and a scaling factor for the whole transmission matrix (ρ). The off-diagonal values are equal to a within-group value (diagonal) multiplied by a scaling factor (ϕ). This scaling factor can either multiply the within-group beta value of the source group (e.g., βHL = ϕ·βLL; Eq. 1A), or the recipient group (e.g., βLH = ϕ·βLL; Eq. 1B), each with a different interpretation.

(1)

The former assumes that between-group transmission is dominated by the transmissibility of the source individuals, implying that adherence to the PHMs primarily prevents onwards transmission, rather than protecting against infection. The latter assumes that between-group transmission is dominated by the susceptibility of the recipient individuals, implying that adherence to the PHMs primarily prevents infection, rather than protecting against onwards transmission. A range of between-group scaling values (ϕ) were simulated to perform sensitivity analysis for the degree of assortativity. Results are only shown for matrix structure A, but alternative assumptions about between-group mixing can be found in the supplement (Supplemental Figs. 1–4). To examine the effect of an intervention to increase PHM adherence, we redistributed a proportion of low- and medium adherence individuals to the high adherence latent class, i.e., a fully effective intervention is equivalent to a single-group SIR model of high adherent individuals. Model fitting and simulation was conducted using the Julia programming language, version 1.10.5 48.

Results

Demographics

Full details can be found in the prior paper31, but briefly: 1410 returning students were recruited, 725 were enrolled, and 684 students completed clinic visits for serum collection between 26 October and 21 December 2020. Of these, 673 students also completed the behavioral survey between 23 October and 8 December 2020. The median age of the participants was 20 years (IQR: 19–21), 64.5% identified as female and 34.6% as male, and 81.9% identified as white. A large proportion (30.4%) were positive for IgG antibodies, and 93.5% (100) of the 107 students with a prior positive test reported testing positive only after their return to campus.

LCA fitting

Of the 673 participants, most students intended to always mask (81.0%), always cover their coughs/ sneezes (81.9%), and always stay home when ill (78.2%) (Table 1). Two of the least common intentions were social distancing by maintaining a distance of at least 6 feet from others outside of their home, avoiding crowded places and mass gatherings > 25 people (43.4% and 53.1% respectively), and avoiding face-touching with unwashed hands (43.5%).

Table 1 Participants’ intention to always or not always follow 8 public health measures.
Full size table

The four- and the three-class LCA models had the lowest BIC respectively (Table 2). Examining the four-class model, there was minimal difference in the classification of individuals, relative to the three-class model. In the four-class model, the middle class (of the three-class model) was split into two groups with qualitatively similar class-conditional item response probabilities i.e., conditional on class membership, the probability of responding “Always” to a given question, except for hand washing and avoiding face-touching with unwashed hands (Supplemental Table 2).

Table 2 Log likelihood, AIC, and BIC of two to seven class LCA model fits.
Full size table

We fit a logistic regression model to predict binary IgG serostatus that included inferred class membership, in addition to other predictor variables we previously identified in31. The mean and median BIC and AIC indicated similar predictive ability of the three- and four-class LCA models (Table 3). Given these factors, the three-class model was selected for use in simulation for parsimony, requiring fewer assumptions and parameters to fit.

Table 3 Mean and median AIC and BIC of multiply-imputed logistic regressions for two to seven class LCA models against IgG serostatus.
Full size table

In the three-class model, approximately 15.75% of individuals were members of the group that rarely intended to always follow the PHMs, 38.04% intended to always follow all guidelines, and the remaining 46.21% mostly intended to mask, test, and manage symptoms, but not distance or avoid crowds (Table 4). We have labelled the three classes as “Low”, “High” and “Medium Adherence” groups, respectively, for ease of interpretation. Examining the class-conditional item response probabilities, the Medium Adherence class had a probability of 0.88 of always wearing a mask in public, but a probability of only 0.19 of social distancing when outside of their homes, for example. Calculating the class-specific seroprevalence, the Low Adherence group had the highest infection rates (37.7%, 95% Binomial CI: 28.5–47.7%), the medium adherence the next highest (32.2%, 95% Binomial CI: 27.0-37.7%), and the most adherent group experienced the lowest infection rates (25.4%, 95% Binomial CI: 20.2–31.1%). Incorporating latent class membership into the imputed GLM model described in our previous paper (30) retained the relationship between adherence and infection. Relative to the least adherent group, the Medium Adherence group experienced a non-significant reduction in infection risk (aOR, 95% CI: 0.73, 0.45–1.18), and the most adherent group a significant reduction (aOR, 95% CI: 0.59, 0.36–0.98) (Table 5). When class assignment was determined probabilistically, similar relationships were observed in the pooled results of logistic regression without confounders, relative to the least adherent group: Medium Adherence group (OR, 95% CI: 0.77, 0.47–1.28); and the High Adherence group (OR, 95% CI: 0.60, 0.36–0.99). 66% of the High Adherence group simulations resulted in a p-value ≤ 0.05.

Table 4 Class-conditional item response probabilities shown in the main body of the table for a three class LCA model, with footers indicating the size of the respective classes, and the class-specific seroprevalence.
Full size table
Table 5 Adjusted odds ratio (aOR) for risk factors of infection among the returning PSU UP student cohort.
Full size table

Compartmental model

The ABC distance distributions indicated that low-moderate levels of between-group mixing better fit the data (Fig. 1). After model parameterization (Supplemental Tables 1,4–6), we examined the effect of increasing adherence to public health guidance. Moving all individuals into the High Adherence class resulted in a 77–96% reduction in final size; when low-moderate between-group mixing is simulated, a fully effective intervention results in approximately 96% (95th percentiles: 88–99%) reduction in final seroprevalence, and when between-group mixing is as likely as within-group mixing, a 89% (95th percentiles: 34–99%) reduction is observed (Fig. 2).

Fig. 1
figure 1

Distribution of the distance from the ABC fits, with the minimum and maximum distances illustrated by the whiskers, and the median distance by the point. Between-group mixing of 1.0 equates to between-group mixing as likely as within-group mixing.

Full size image
Fig. 2
figure 2

(A) The reduction in final infection size across a range of intervention effectiveness (1.0 is a fully effective intervention), accounting for a range of assortativity. Between-group mixing of 1.0 equates to between-group mixing as likely as within-group mixing. Each line represents the mean reduction (of 200 simulations) resulting from an intervention, and is associated with a β matrix specific to the degree of between-group mixing illustrated; (B) The relative distribution of group sizes at three levels of intervention effectiveness (0.0, 0.2, 1.0).

Full size image

Discussion

In this interdisciplinary analysis, we collected behavioral data from surveys and integrated it with serosurveillance results. This approach allowed us to use LCA to categorize a population’s transmission potential with measures related to risk tolerance and behavior. The LCA model was fit without inclusion of infection status data, but class membership was correlated with IgG seroprevalence. The classes that were the most adherent to PHMs experienced the lowest infection rates, and the least adherent exhibited the highest seroprevalence. As the logistic model cannot account for indirect effects resulting from between-class interactions, we parameterized a dynamical SIR model to explore the effect of interventions of varying degrees of effectiveness.

Although a four-class LCA model was a marginally better fit for the data, there were not substantial differences in class assignment relative to the three-class LCA model. The three-class model was selected for use in simulation for parsimony, requiring fewer assumptions and parameters to fit. Upon parametrizing the compartmental model, smaller ABC distance values were observed for low to moderate levels of between-group mixing, implying some degree of assortativity in our population, though the exact nature cannot be determined from our data. Examining the three classes, 38% of individuals already intended to always follow all PHMs. As a result, only 62% of the study population could have their risk reduced with respect to the PHMs surveyed. Further, the infection rates observed in the High Adherence group indicates that even a perfectly effective intervention aimed at increasing adherence to non-pharmaceutical PHMs (i.e., after the intervention, all individuals always followed every measure) would not eliminate transmission in a population, an observation that aligns with prior COVID-19 research49,50,51,52. The extent to which the infection in the High Adherence group is a result of mixing with lower adherence classes cannot be explicitly described, but the sensitivity analysis allows for an exploration of the effect and ABC fits suggest low-moderate levels of between-group mixing occurred. Varying the structure of the transmission matrix yielded very similar quantitative and qualitative results (Supplemental Figs. 1–4).

Examining the impact of increasing adherence to PHMs (modeled as increasing the proportion of the population in the High Adherence class), a fully effective intervention saw between a 77–96% reduction in the final size of the simulation outbreak. We note that the effect at a fully effective intervention is conceptually analogous to the population attributable fraction (PAF) proposed in53; though rather than quantifying the impact of removing one risk group, as in53, we consider the impact of all individuals moving to the low-risk group. Each set of simulations for a given degree of assortativity has a different associated set of parameter values for the transmission matrix. The difference in the magnitude of the achievable reduction at a given level of intervention for the different assortativity levels is attributed to the difference in the corresponding fitted parameters (Supplemental Tables 4–6). With higher levels of between-group mixing, the initial SIR parameterization generally results in lower transmission parameters for the High-High adherence interactions, as more infections in the High Adherence group originate from interactions with Low and Medium Adherence individuals. Increasing adherence, therefore, results in a greater reduction of the overall transmission rate than in simulations with less assortativity.

Limitations and strengths

The student population was recruited using convenience sampling, and therefore may not be representative of the wider population. Those participating may have been more cognizant and willing to follow public health guidelines. Similarly, because of the University’s extensive messaging campaigns and efforts to increase access to non-pharmaceutical measures36, such as lateral flow and polymerase-chain reaction diagnostic tests, the students likely had higher adherence rates than would be observed in other populations. However, these limitations are not inherent to the modeling approach laid out, and efforts to minimize them would likely result in stronger associations and conclusions due to larger differences in the latent behavioral classes and resulting group infection rates.

It is well known that classification methods, like LCA, can lead to the “naming fallacy”44, whereby groups are assigned and then specific causal meaning is given to each cluster, affecting subsequent analyses and interpretation of results. In this paper, this effect is reduced by virtue of the analysis plan being pre-determined, and the relationship with the outcome showing a positive association with the classes in the mechanistically plausible direction (i.e., increasing adherence to PHMs results in reduced infection rates). Our decision to conduct the simulation analysis with the three-class model was, in part, to avoid the potential bias that would arise from naming or assigning an order to the two intermediate risk groups.

Despite these limitations, this work presents a novel application of a multidisciplinary technique, outlining how alternate data sources can guide future model parameterization and be incorporated into traditional epidemiological analysis, particularly within demographically homogeneous populations where there is expected or observed heterogeneity in transmission dynamics. This is particularly important in the design of interventions that aim to target individual behaviors, allowing the categorization of populations into dynamically-relevant risk groups and aiding in the efficient use of resources through targeted actions. Future research should consider including perceived agency and efficacy for PHM adherence.