Behaviour quantification of public health policy adoption - the case of non-pharmaceutical measures during COVID-19

1.   Introduction

The rapid worldwide spread of the novel coronavirus SARS-CoV-2 and its related disease COVID-19 quickly shifted the priorities of public health experts and policymakers in early 2020. Despite strict initial protective measures applied in Canada, such as the closure of borders and non-essential services ^[1], new cases continued to sweep across provinces. For a novel pathogen with unidentified transmission dynamics, a population's behaviour is a significant factor in predicting the spread of infection, which can inform health policies, such as resource allocation and regional guidelines. However, many existing mathematical models do not adequately represent the relation between individual decision making and disease prevalence that occurs in a rapidly evolving public health crisis like a pandemic ^[2]. Personal values and perceived transmission risk to individuals played a large role in the evolution of the COVID-19 pandemic in Canada, a better understanding of which could help decision makers to model improved disease mitigation strategies in the future. In this paper, we investigate the adoption of non-pharmaceutical control measures in the province of Ontario, as well as in the US states of Georgia and Pennsylvania, and use a game theoretic framework to infer how the adoption of non-pharmaceutical interventions influenced the outcomes across the regions, as well as at what levels of efficacy.

2.   Background

On January 25$^{th}$, 2020, the first positive case of COVID-19 in Canada was reported in Toronto, Ontario ^[3]. In the absence of pharmaceutical prophylactic options, the national number of reported cases rose to nearly 10,000 by the end of March ^[4], at which point the Public Health Agency of Canada began recommending, and later mandating, the implementation of non-pharmaceutical interventions (NPIs) to reduce and delay further contact transmission events, including social/physical distancing, shelter-in-place or "lockdown" periods, service closures, intensified sanitation procedures, and face masks ^[1,5]. However, even recent Canadian pandemic planning guidance notes that the impact and limitations of such NPIs in an outbreak situation were still not well studied leading up to 2020 ^[6]. While the success of a public health intervention is known by epidemiologists to depend upon adherence by the target population, the COVID-19 response presents an opportunity to further investigate the spatiotemporal factors affecting individuals' behaviour and decision making regarding NPIs during a pandemic.

Mathematical modeling has come a long way in terms of utility for public health planning and policy, thanks in part to the extension and refinement of seminal works like the Kermack-McKendrick susceptible-infected-recovered (SIR) compartmental model as well as the sophistication of modern data and computing abilities. For example, in ^[7], the authors implemented additional quarantine and isolation compartments and highly detailed population contact rates (developed by ^[8]) to simulate early COVID-19 transmission in Ontario. However, even the most detailed of contact patterns represent population averages that occur under "normal" conditions and cannot be assumed static in the context of epidemic modeling ^[9]. Moreover, deviations from behavioural norms are rarely uniform across entire populations in extreme circumstances ^[10]. Growing evidence from recent public health crises (e.g., 1980s HIV/AIDS, 2003 SARS, 2009 H1N1) demonstrates that people's behaviour evolves with the epidemic progression as they respond to new information and form perceptions of risk (see for example: ^{[11,12,13,14,15,16]}). The sum of these factors may result in underestimated pandemic effects and motivates further refinement of human behaviour modeling in epidemic simulation models.

Earlier works approached behaviour and transmission as having a cause and effect relationship, whereby one is modelled operationally and the other as a basis to interpret the results ^[2]. For example, ^[17] was the first to vary the contact rate $\beta$ over time according to disease prevalence, so that a rise in infection risk would result in a decrease in population mobility. These ideas lead to "behaviour-incidence" models that incorporated both behaviour and transmission submodels as complementary inputs to better reflect real-world dynamics ^[2,18]. In ^[19,20] the authors tackle the issue of determining $\beta$ by using a hypothesis-free machine-learning algorithm to estimate the transmission rate from data on non-pharmaceutical policies, and in turn forecast the confirmed cases using a mechanistic disease model. Dynamic adjustments to contact rates remain a common technique for modeling NPI-related behavioural changes across many infectious diseases, though there are differences involving the degree of heterogeneity in population characteristics and the number of strategies considered; the COVID-19 modeling challenge was no exception. One 2020 study used an economic-based optimization approach coupled with a compartmental model where the proportion of individuals voluntarily adopting NPIs changed depending upon their age-based risk of mortality and the time-varying disease prevalence ^[21]. Another considered differences in NPI compliance simultaneously by allowing for interaction of individuals with behaviour-varying likelihoods of infectiousness.

Game theoretic models ^[22] have become an attractive behavioural pairing to disease transmission models, providing a framework to view a population as being made of a group of individuals ("players") who make strategic decisions to try and maximize their own health benefits during a public health threat, while also being influenced by the decisions of others. The point where no player can further improve their outcome by choosing a different strategy is called a Nash equilibrium - and is the most common way to define the solution to a non-cooperative game involving two or more players ^[23]. Game theory has seen prominent application in voluntary vaccine uptake models ^{[24,25,26,27,28]}, and to a lesser extent NPIs ^[16].

3.   The game behavioural model and quantification of perceived risks/discomforts

In this work, we design a population game to serve as a behavioural model that is paired with an existing compartmental transmission model for COVID-19 used in ^[7,29,30]. The game model helps to investigate spatio-temporal population adoption of NPI mitigation measures across the province of Ontario, Canada in 2020 (where Ontario is used as our case study), whereas the pairing between the game and the transmission model is used to infer efficacy of NPIs over time in Ontario. It is important to differentiate here a game model from traditional decision models where individuals simply optimize their well-being, or minimize their cost, without allowing for influence from the decisions of others. In a game context, when population NPI adoption is high, those who choose not to follow NPI's may still see benefit from reduced transmission odds without having to incur the "costs" of avoidance measures.

In what follows, we define Ontario's public health (PH) regions as a group of players competing to minimize personal risks of getting infected with COVID-19 as well to minimize personal discomfort with adoption of NPIs. Each PH region is thought of a set of identical players, where players' decisions at any given time are interdependent. The choice of one player to adopt NPIs provides a net positive health payoff for others by reducing COVID-19 transmission risk, but in turn it may incentivize another player not to use NPIs and still benefit from the conferred protection. Since both options are available at each time step across the population, players face uncertainty and make use of available information to assess their personal risks and discomforts. The concepts of free ridership and information signaling are well studied in vaccine game literature (see for instance https://royalsocietypublishing.org/doi/full/10.1098/rsif.2016.0820).

3.1.   Data driven analysis of perceived risk of infection and personal discomfort with measures in ON public health units

COVID-19 data were coordinated across Ontario's 34 Public Health units (PHUs), which are responsible for the delivery of local health services and whose jurisdictions are generally partitioned according to census municipalities ^[31]. We assume that time-varying individual decisions regarding NPIs are affected by both personal discomfort incurred from NPIs adoption (denoted by $\tilde{r}$) and by the personal perceived risk of getting infected when not adopting NPIs (denoted by $\bar{r}$). To quantify $\bar{r}$ and ${\tilde{r}}$, we set up two multivariable linear regressions (MLR) per PHU per month to express each in terms of relevant empirical data.

First, we collected daily reported counts of COVID-19 incidence and deaths by PHU. Such information was broadcast daily by major provincial news bodies and their regional affiliates (e.g., CBC, CTV News, Ontario Newsroom) via television and the internet to provide ongoing coverage of the pandemic and related developments. Similarly, we collected daily COVID-19-related hospitalizations, but at the province level instead of by PHU. The decision was based partly on reporting trends by news agencies, but also on Ontario's disproportionate population density and concentration of hospital facilities. For example, the Northern health region, encompassing 7 of 34 PHUs and 88% of Ontario's landmass, serves just 6% of the province's population and hosts 47 of 141 hospital sites ^[32,33,34]. Additionally, overwhelmed hospitals coordinated patient transfers to accommodate for the sheer volume of severe COVID-19 cases ^[35], making it difficult to track cases by PHU of origin. We assume that the overall daily number of hospitalizations in Ontario better reflects the severity of the COVID-19 situation that would factor in to an individual's personal risk perception than their local hospital numbers.

Second, we examined daily population mobility changes across all PHUs as a proxy for community contact reduction reflecting voluntary and government-mandated risk avoidance measures before and after the onset of COVID-19. We utilize the COVID-19 Community Mobility Reports made available by Google ^[36] that tracked and partitioned cellphone mobility data according to the following five indicators: grocery and pharmacy stores, parks, transit stations, retail and recreation, and workplaces. Each indicator is calculated as frequency of visits relative to activity during the baseline time period, defined as the median of the corresponding day of the week between January 3 and February 6, 2020. More information regarding the Google mobility data and how it was partitioned according to the PHUs can be found in the Appendix B.

Third, to examine the individual effects of government-mandated health policies, we used the Oxford COVID-19 Government Response Tracker (OxCGRT) which was created to collect and standardize information on COVID-19 government policy measures in over 180 countries and provide tools to compare and measure their impacts ^[37]. We focus specifically on the OxCGRT Stringency Index, which measures the strictness of "lockdown" policies that primarily restrict regular activities and mobility. This data was available at the provincial level. More information on the OxCGRT can be found in the Appendix C.

Finally, using the data provided by the Institute for Health Metrics and Evaluation (IHME), we are able to see the provincial level of mask compliance for Ontario ^[38]. Provincial mask compliance was directly fed in as a 7-day rolling average of mask use mean, given as a percentage.

Estimation of perceived risks of infection, as well as personal discomfort with mask use, are in general difficult problems which necessitate further assumptions alongside available data. In our case here, we will use a combination of regional ($_R$) and province-wide ($_{P}$) freely available data on incidence, hospitalizations, Oxford index, Google mobility index, and observed mask use. Further, we will assume that the observed level of mask use in a regional population is reflected in a perceived discomfort index with wearing a mask, while the work mobility index in a region is reflected in a perceived risk of infection, for an average individual in a PHU.

If we assume a linear relationship between observed mask use (as dependent variable), cases and Oxford index on one hand (as independent variables), and another linear relation between the reduced work mobility (as dependent) and cases and hospitalizations (as independent variables) we have:

Here $\bar{r_i}$ is the value of the work-related mobility index in the absence of new COVID-19 cases and hospitalizations: this means that when cases and hospitalizations increase, the perceived risk will decrease, and vice-versa. We then assume that $\bar{r}_i(t)$ is a proxi-measure of the perceived risk of infection in region $i$ over time. Similarly, we assume that the observed mask wearing adoption level in Ontario is a function of the disease prevalence and previous mask use:

Here $\tilde{r_i}$ is the value of the observed mask use at endemic level of cases ^* and absence of economic measures: this means that when cases and hospitalizations increase, the perceived discomfort will decrease ^†, and vice-versa. We then assume that $\tilde{r}_i(t)$ is a proxi-measure of the perceived discomfort of wearing a mask in region $i$ over time.

^*The use of masks in public was well established pre-COVID in certain countries such as Japan, South Korea, Singapore, etc., and less so in others, during seasonal respiratory viruses spread.

^†As masks reduce transmission, individuals are more likely to tolerate the mask at higher levels of disease prevalence.

Both computations for $\bar{r_i}$ and $\tilde{r_i}$ capture very broad cases which are sufficient for our research, however possibilities for further refinement are discussed in our discussions Section 6. The regression analysis was computed using the statistical software R's lm function which is used for fitting linear models. We repeat for each PHU, collecting the perceived risk of infection and personal discomfort terms into the $\bar{r_i}$ matrix and $\tilde{r_i}$ matrix respectively, where $\bar{r_i}(t)$ for $i \in \{1, 2, \dots, 34\}$ represents the corresponding PHU via column, and $t \in \{1, 2, \dots, 10\}$ is the monthly time step starting from March to December of 2020.

Some interesting insights in the evolution of perceptions over time are highlighted in the Appendix D. We saw some interesting differences in the evolution of the perceived risks of COVID-19 between the PHU regions with the lowest population density and those with the highest population density, while the evolution of the personal discomfort with wearing a mask remained similar amongst the PHU regions. We also discuss some correlations between how the population perceived NPI's and some socio-economic and demographic factors.

3.2.   Transmission model

Following ^[30], COVID-19 transmission in Ontario is simulated via an SEIRL (Susceptible-Exposed-Infectious-Recovered-isoLated) model with daily incidence data across Ontario and at the PHU level. This model assumes that the population is homogeneous and well-mixed, that infected individuals are either symptomatic or asymptomatic, and that most individuals who are aware of their infection isolate for a fixed period of time. The system is described with the following ordinary differential equations:

Our model parameters (Table 1) have been taken from literature, with the exception of the isolation rate $1/\kappa$ which we assume to be 1 day.

Table 1.  Parameter values for our model 3.3.

Symbol	Definition	Initial value	Reference
$N_{total}$	Population size	Table 4
$\sigma$	Rate at which exposed become infectious (days$^{-1}$)	1/2.5	[39]
$a$	Proportion of permanently asymptomatic cases	0.5	[40,41,42]
$\epsilon$	Proportion of compliance with isolation	0.78	[43]
$\kappa$	Isolation rate	1	assumed
$\gamma$	Recovery/removal rate	1/7	[44]

---

 | Show Table

DownLoad: CSV

4.   The NPI adoption game

Our model has 34 interacting groups whose decisions affect each other's perceived relative risk of infection. That is to say we have a game taking place between population groups. A multiplayer Nash game (see first ^[45]) consists of a finite number of players denoted by $N > 0$. A generic player $i \in \{1, \dots, N\}$ has a strategy vector $x_i$ selected from a closed convex set $S_i \subset \mathbb{R}^{n_i}$ so that $n_1+\cdots+ n_N = n$. A player has a payoff function $u_i: K \longrightarrow \mathbb{R}$, where $K := S_1 \times \dots \times S_N$ or the cross product of the strategy sets of all players; evidently: $K\subset \mathbb{R}^{n_1+...+n_N = n}$. In general we strive to find a Nash equilibrium strategy vector $\underline{x^*} := (x_{1}^*, x_{2}^*, \dots, x_{N}^*)$, such that when player $i \in \{1, \dots, N\}$ plays $x_{i}^*$, there is no reason to switch strategies and as such, these strategies should be stable equilibrium solutions of our game. Mathematically, we define a Nash vector of strategies as follows:

Definition 1. Assuming each player is rational and wants to minimize their payoff function $U_i : K \longrightarrow \mathbb{R}$. Then aNash equilibrium is a vector $x^* \in K$ which satisfies the inequalities:

Now recall from our previous section that $\bar{r}_i(t)$ is a measure of the perceived risk of COVID-19 infection to an individual in region $i$. Similarly, recall that $\tilde{r}_i(t)$ is a measure for the personal discomfort to an individual when wearing a mask. Let us now denote by $r_i(t)$ the relative perceived risk estimate of an individual residing in one of the 34 PHUs, defined as the ratio of these two factors:

That is to say, whenever the relative risk exceeds $r_i(t) > 1$, the personal discomfort with masks is higher that the perceived risk of infection (for that region, over time, for an average individual), and vice versa.

Note that all residents in a given PHU are provided with the same information and use this information to asses their risks and discomforts. In actuality, there is a population of individual decision-makers within each PHU whose aggregate decisions constitute an overall probability of adopting NPIs, therefore we will call a PHU a "player" in what follows, and we will call the aggregate decisions of individuals residing in each PHU an "average player" in what follows. The population of each PHU represents a fixed proportion of the total Ontario population $\epsilon_i$, where $\epsilon_i \in (0, 1)$ and ${\sum\limits_{i = 1}^{34} \epsilon_i = 1}$. The strategy of player $i\in\{1, \cdots, 34\}$ is the probability that they will adopt NPI's, denoted by $x_i \in [0, 1]$, hence the strategy set for each player is $[0, 1]$ and the strategy set of our game will be $K:= [0, 1]^{34}$. To define the payoffs to each player, we adapt a utility function that was previously used to determine Nash equilibrium solutions for well-studied vaccine-preventable childhood diseases from ^[25,26]:

Scaling the payoffs by dividing with $\bar(r)_i(t)$ and using (4.1), we get:

In this setting, $U_i$ represents the player's disutility (or cost) of disease prevalence as a function of total population NPI compliance given by $\sum_i \epsilon_i x_i$. A greater perceived coverage in the population means a reduced perceived infection risk for susceptible individuals.

In general, $\delta_i(t)$ signifies the perception of player $i$ on the overall adoption of NPI measures in the province, and it is not something we can estimate from the decision model or from data, but we rely on it to solve the game. Since it is a quantity very difficult to measure, we will solve our decision model using a statistically uniform spread of values of $\delta_i(t)\in[0, 1]$^‡, for each $i$ and $t$, thus generating an envelope of NPI adoption curves for each region and for ON, and reporting the quantiles of the curves envelope in the figures below.

^‡The value of 0 signifies that the player perceives they has no information on the level of NPI adoption in ON, whereas a value of 1 signifies the player perceives they have a very accurate level of information on NPI adoption in ON.

Therefore, we first solve the game:

and obtain the Nash equilibria $x^*_i$ for each PHU in Ontario. Each $x^*_i$ represents a given player in region $i$'s probability of respecting NPIs. We use a projected dynamical systems (PDS) approach to compute the optimal strategy and respective equilibrium NPI-adoption estimates (see ^[26], using MATLAB). We are able to determine the expected NPI-adoption level for Ontario, namely:

We present the determined $eNPI_{adopt}$ for ON in Figure 1:

Figure 1.  Weekly $eNPI$ adoption across Ontario from March to December 2020; note that for the first 2 weeks of March 2020, the adoption levels are 0.

DownLoad: Full-Size Img PowerPoint

5.   Efficacy of adopted measures as a fraction of transmission reduction in ON

In this section we investigate the use of non-pharmaceutical interventions on disease transmission. To reflect the contribution of the NPI's adoption on the reduction of transmission of COVID-19, we further use the epidemiological model for ON (Section 3 above), specifically recalling that, although the transmission rate $\beta$ is generally constant for a given disease, in the absence of any interventions, the presence of NPIs has affected the transmission over time. Therefore, it is more realistic to think of $\beta: = \beta(t)$. Moreover, the transmission intrinsically depends on the number of contacts a single indiviual has, and the probability of transmission per contact between infected and susceptible individuals. Since NPIs were meant to decrease contacts while masks were meant to lower the probability of transmission per contact, $\beta(t)$ should be able to capture the changes in transmission induced by NPI's.

Figure 2.  The interplay between the transmission dynamics (reflected in observables such as cases, deaths, hospitalizations and economic measures), the decision-making happening as their result, and the reflection of the decisions in the transmission dynamics 1 week further into the future.

DownLoad: Full-Size Img PowerPoint

Since our estimated $eNPI_{adopt}(t)$ should be meaningful vis-a-vis the evolution of the transmission $\beta(t)$, we modify the SEIRL model of Section 3.2 to incorporate a time dependent transmission rate given by:

where $\beta_0$ is the value of the transmission in March 2020 (near disease-free equilibrium) deduced above and $\tau$ is the delay between the behaviour of individuals and its effects on transmission (in days). Further, we investigated a few different values for the delay in the effects of NPI-adoption ($\tau$) and found that the closest fit to the Ontario incidence data was $\tau = 1\, \, week$ delay. Note that there is no delay in the first couple weeks, as the lockdown would have an immediate effect. The delay is needed to reflect that respecting NPIs today is only measurable in the disease transmission later in the future.

Equation (5.1) is only useful if we are able to derive/estimate the transmission rate $\beta_0$ near the start of the pandemic. Similar to ^[30], here we compute the value of the transmission rate $\beta_0$ as a function of the exponential growth of the SEIRL model near the disease free equilibrium:

To estimate the exponential growth factor for ON we use weekly case incidence data (which we denote by $inc(t)$) for each region, and for the whole province. We assume that $inc(t)$ is given by an exponential curve of the type (see ^[29,30]):

and we compute $\rho\approx 1.239646$ (the first positive slope in Figure 3 below):

Figure 3.  Exponential growth estimate in March 2020 in Ontario; in the plot, week 0 is the week of March 1st 2020.

DownLoad: Full-Size Img PowerPoint

With this value of $\rho$, we use Eq (5.2) to compute $\beta_0\approx 6.890257$ and an initial reproduction number of $R_0 \approx 1.725$ at province level. Regionally, at the level of each PHU, we computed the emergence and initial growth phase, and depicted the spread in the map graph below, showcasing the spatio-temporal aspect (see Figure 4).

Figure 4.  This is a temporal depiction of the exponential growth at the beginning of the pandemic, seen per PHU. We can see that the earlier PHU was Toronto and surrounding area - March 7, 2020; then followed by immediate neighbouring PHU's (i.e., Southwestern Ontario and Ottawa areas) - March 14, 2020; then further spreading into the lower population density PHU's - March 21, 2020 & April 4, 2020.

DownLoad: Full-Size Img PowerPoint

With all these pieces together, we execute a 2-parameter fitting of simulated disease incidence compared to real data incidence in ON, to derive optimal values of the efficacy of NPI adoption, $eff(t)$ in Ontario (between $[0, 1]$), and optimal values of the $fit(t)$ parameter (between $[0, 1]$) which accounts for other variables that may have affected transmission, such as: temperature, holidays, school closures, etc. The resulting $\beta(t)$, efficacy of NPIs and $fit(t)$ parameters from mid-May 2020 to December 2020 are shown in Figure 5.

Figure 5.  Median and confidence interval for effective $\beta(t)$ values (upper panel). Median and confidence interval for efficacy of NPIs ($eff(t)$) with overlaying trendline in black (middle panel). Median and confidence interval for fit parameter $fit(t)$ (lower panel).

DownLoad: Full-Size Img PowerPoint

5.1.   Model application - other regions

To demonstrate the versatility of our mathematical framework (game and epidemiological model), we show below how we applied it to two other similar regions in the United States: Pennsylvania and Georgia. Table 2 includes the characteristics of these two states, and their distinctions from the Ontario case.

Table 2.  Populations of ON, PA, and GA in millions, and the corresponding numbers of public health units.

	Population	Health regions
ON	14 mil.	34
PA	12 mil.	6
GA	10.6 mil.	18

---

 | Show Table

DownLoad: CSV

In a similar manner, we estimate the expected adoption levels of NPI in PA and GA, from March 2020 to December 2020, and we plot them in Figure 6 below.

Figure 6.  Weekly expected NPI adoption across PA (left) and GA (right) from March to December 2020; note that for the first 2 weeks of March 2020, the adoption levels are 0.

DownLoad: Full-Size Img PowerPoint

Again, we fit our SEILR model to actual incidence data using a 1-week delay starting in April, as we did for Ontario.

Figure 7.  The 2-parameter fit of simulated incidence on real incidence data for Georgia and Pennsylvania in U.S.

DownLoad: Full-Size Img PowerPoint

This resulted in the following $\beta(t)$, efficacy of NPIs, and $fit(t)$ parameters from mid-May 2020 to December 2020 in Georgia, shown in Figure 8, respectively.

Figure 8.  Georgia, U.S.: Median and confidence interval for effective $\beta(t)$ values (upper panel). Median and confidence interval for efficacy of NPIs ($eff(t)$) with overlaying trendline in black (middle panel). Median and confidence interval for fit parameter $fit(t)$ (lower panel).

DownLoad: Full-Size Img PowerPoint

Similarly, the results for Pennsylvania are shown in Figure 9, respectively.

Figure 9.  Pennsylvania, U.S.: Median and confidence interval for effective $\beta(t)$ values (upper panel). Median and confidence interval for efficacy of NPIs ($eff(t)$) with overlaying trendline in black (middle panel). Median and confidence interval for fit parameter $fit(t)$ (lower panel).

DownLoad: Full-Size Img PowerPoint

Table 3 shows the average efficacy of NPIs for Ontario, Georgia, and Pennsylvania.

Table 3.  Average efficacy of NPIs implemented in Ontario, Pennsylvania, and Georgia from mid-May 2020 to December 2020.

Region	Average efficacy
ON	55.43%
PA	42.15%
GA	53.50%

---

 | Show Table

DownLoad: CSV

6.   Discussion and conclusions

In this study, we employed a multi-faceted approach to understand the adoption of non-pharmaceutical interventions (NPIs) during the COVID-19 pandemic across Ontario, Canada. Our findings highlight several key insights into the dynamics of public health policy adoption, behavioural responses, and the efficacy of these measures in mitigating disease transmission. We then applied our approach to Pennsylvania and Georgia, United States, to demonstrate the applicability of our model and to gain further insight of NPI adoption in regions with similar populations.

Our study showed how the perceived risk of infection and personal discomfort with NPIs affect decision-making processes among individuals. We found that these factors varied across PHUs, influenced by local epidemiological conditions (such as cases, hospitalizations, and deaths) and the strictness of policy measures imposed by authorities. We then integrated an SEIRL compartmental model to estimate the effectiveness of NPIs in reducing COVID-19 transmission across Ontario. By incorporating time-dependent transmission rates influenced by NPI adoption levels, we demonstrated a tangible (and quantifiable) reduction in disease transmission attributable to these measures. Moreover, we also showed that quantifying the adoption and effectiveness of NPI measures in a pandemic does not give the full picture of dampening the transmissibility - thus the fit parameter we quantified is a good proxi for dampening effects we did or could not measure (such as weather, types of governance, effectiveness of public policy communications, etc.).

When comparing findings in Ontario, Pennsylvania, and Georgia, we found that Pennsylvania exhibited the lowest average adoption of NPIs. This disparity in adoption rates can be attributed to several factors, including differences in governmental mandates, public compliance, and socio-economic influences. The lower adoption in Pennsylvania suggests potential challenges in garnering widespread adherence to public health measures, which may have contributed to its comparatively higher COVID-19 transmission rates. The effectiveness of NPIs in reducing COVID-19 transmission varied across these three regions, reflecting the differences seen in NPI adoption rates.

There were several limitations during our study. First, our analysis relied on publicly reported daily cases, deaths, and hospitalizations as the primary sources of epidemiological data. These data are subject to potential under-reporting or delayed reporting, which may lead to an underestimation of the true incidence and severity of COVID-19 in the studied regions. The reliance on reported figures could introduce biases, particularly if there were variations in testing rates or reporting practices over time. Thus, the actual disease burden and dynamics could differ from our model's estimations, based on these reported numbers alone. Second, our SEIRL transmission model uses an estimated $\epsilon$ = 78% compliance rate for isolation measures in ON, and 70% rate in Georgia and Pennsylvania, as well as an assumed 1-day isolation rate. These are idealized and may not fully capture the variability in human behaviour and adherence to NPIs observed in real-world settings.

In conclusion, this study has presented a robust framework for analyzing the effect of NPI measures on the spread of COVID-19 in province/state areas in Canada and the US, focusing on the heterogeneity of the public health sub-regions. By integrating game theory with epidemiological models, we gain a direct way to quantify/implement decision-making behaviours during an epidemic, such as COVID-19. For future work, we plan to include vaccination data in our approach. This will allow us to apply our modeling framework to 2021 and 2022, as vaccinations became pivotal in controlling the pandemic's progression. Applying this model to a variety of sub-regions within a country (or state/province) may allow further insight into how regional differences affect decision making and thereby the effectiveness of NPIs.

Use of AI tools declaration

The authors declare they have not used artificial intelligence (AI) tools in the creation of this article.

Acknowledgments

The work was supported to by a joint grant from the National Sciences and Engineering Research Council (NSERC) of Canada and Sanofi under an Alliance Ⅱ grant (Cojocaru, M. G. principal investigator).

Conflict of interest

Thommes, E. W. is a Sanofi employee and may hold stock options. Cojocaru, M. G., Loster, R, Smook, S., Humphrey. L., Lyver, D, Mohammadi, Z. declare no conflict of interests.

Appendix

A.   Google mobility data

Google creates these aggregated and anonymized data sets from users who have turned on and agreed to share the information from their location history setting of their Google accounts on their phones ^[36]. As such, these data sets may not be representative of the population, especially in more rural settings with bad data reception or regions where cellphone usage is lower. Additionally, Google has not publicly shared the precise methodologies for calculating these social mobility scores, thus there is a certain degree of ambiguity as to how representative this data is of true population trends. However, through data provided by Statistics Canada, the number of individuals with personal cell phones in Canada as of 2020 was above $80\%$ ^[46]. Additionally, through data provided by Statista, Google Maps is the most downloaded mapping app in the United States, hence it is the most indicative mobility application for Canadian usage ^[47].

Some manipulation of Google mobility data was required to reconcile differences in geography and data recording. Ontario has 49 census divisions and six health regions (North East, North West, East, Central, Toronto, and West) that are further refined into 34 public health units (PHUs). In general, PHUs correspond to census borders, with smaller and less densely populated regions combining to form singular PHUs. In these cases the social mobility indicators were averaged across the census divisions making up each PHU. Additionally, the large census division of Kenora spans the northern PHUs of Northwestern (NWR), Thunder Bay (THB) and Porcupine (PQP) near equally. Given that the largest proportion of Kenora's population falls within NWR, we assign the mobility data for Kenora to this region and combine it with Rainy River division. Table A1 outlines the 51 Google mobility data collection regions and to which of the 34 Ontario public health units they were assigned.

Table A1.  PHU regions and their corresponding Google Mobility Regions.

Public Health Region	Code	Google Mobility Regions
Algoma Public Health	ALG	Algoma District
Brant County Health Unit	BRN	Brant County Brantford
Chatham-Kent Health Unit	CHK	Chatham-Kent
Durham Region Health Department	DUR	Regional Municipality of Durham
Eastern Ontario Health Unit	EOH	Prescott and Russel Stormont, Dundas and Glengarry
Grey Bruce Health Unit	GBC	Grey County Bruce County
Halton Region Health Department	HAL	Regional Municipality of Halton
Hamilton Public Health Services	HAM	Hamilton
Haldimand- Norfolk Health Unit	HDN	Haldimand County Norfolk County
Haliburton, Kawartha, Pine Ridge District Health Unit	HKP	Haliburton County Kawartha Lakes Northumberland County
Hastings and Prince Edward Counties Health Unit	HPE	Hastings County Prince Edward County
Huron Perth District Health Unit	HPH	Huron County Perth County
Kingston, Frontenac and Lennox and Addington Public Health	KFL	Frontenac County Lennox and Addington County
Lambton Public Health	LAM	Lambton County
Leeds, Grenville and Lanark District Health Unit	LGL	Lanark County Leeds and Grenville United Counties
Middlesex-London Health Unit	MID	Middlesex County
Niagara Region Public Health Department	NIA	Regional Municipality of Niagara
North Bay Parry Sound District Health Unit	NPS	Nipissing District Parry Sound District
Northwestern Health Unit	NWR	Rainy River District Kenora District
Ottawa Public Health	OTT	Ottawa
Southwestern Public Health	OXF	Oxford County Elgin County
Peel Public Health	PEL	Regional Municipality of Peel
Porcupine Health Unit	PQP	Cochrane District
Peterborough Public Health	PTC	Peterborough County
Renfrew County and District Health Unit	REN	Renfrew County
Simcoe Muskoka District Health Unit	SMD	Simcoe County Muskoka District Municipality
Sudbury and District Health Unit	SUD	Sudbury District
Thunder Bay District Health Unit	THB	Thunder Bay District
Toronto Public Health	TOR	Toronto
Timiskaming Health Unit	TSK	Timiskaming District
Region of Waterloo Public Health	WAT	Waterloo Regional Municipality
Wellington-Dufferin-Guelph Public Health	WDG	Wellington Dufferin
Windsor-Essex County Health Unit	WEC	Essex County
York Region Public Health Services	YRK	Regional Municipality of York

---

 | Show Table

DownLoad: CSV

B.   Stringency index

The Bank of Canada (BoC) adapted the OxCGRT stringency index to create a stringency index of their own which is a measure of containment policies and public information campaigns ^[48]. More specifically, this index focuses on the economic impact of COVID-19. Similar to the OxCGRT, the BoC's stringency index does not measure the efficacy of a province's response to COVID-19 nor does it provide a direct measure of the impact of government policies on the economy ^[48]. As of the summer 2020, researchers have been collecting publicly available information on government policies, creating daily government response indexes for all 10 provinces. With this, the bank could systematically measure, track and compare government policy responses. The BoC's stringency index follows the methodology of the OxCGRT with a few adjustments to make it more appropriate for the Canadian context and to capture specific differences in policy responses. In the BoC's stringency index, they added 3 policy indicators to the OxCGRT stringency index which included the following:

● $C9$ - Restrictions on inter-provincial travel which splits the original OxCGRTs "internal movement restrictions" ($C7$) to capture only intra-provincial travel. Therefore adding this new $C9$ term gives greater weight to travel restrictions in the stringency index.

● $C10$ - Enforcement mechanisms for individuals: OxCGRT does not include a measure identifying differences in provinces that have similar containment policies but different punishments or enforcement of these policies. This C10 indicator captures penalties for individuals violating public health orders.

● $C11$ - Enforcement mechanisms for businesses violating public health orders.

By adding more targeted dimensions, they increase the level of detail in the variation across provinces as well as over time. Additional geographical dimensions were included in the BoC index including types of outdoor events, outdoor versus indoor events and exceptions for certain travellers. Additionally, they refined the categories of coding values for several policy indicators based on a more detailed examination of policies in Canada. While they added categories to their indicators, the sub-index calculation decreases as the denominator in this equation increases

● $a_j$ is the ordinal value for the most restrictive measures in place within the province.

● $b_j$ is the maximum coding value.

● $n_j$ is the total number of targeting dimensions (e.g., by geography).

● $f_{ij}$ is the flag for targeting. dimension, $i$, which is equal to 1 if the policy is general or 0 if targeted.

As a result of this increasing $b$ term, often the sub-index is lower when compared with OxCGRT. While the BoC stringency index is better tailored to the Canadian landscape and the resulting economic implications, it does not provide a more refined scale and thus we still only have provincial data. Additionally, these factors which cater to the economic side of the stringency measures do not help us better fit our model to describe the behaviour of individuals in a given PHU. Thus we use the broader scope of the OxCGRT Stringency Index in our determination of the risk parameters.

C.   Evolution of perceptions

There are some interesting insights that we could highlight in the evolution of perceptions over time. For instance, in Figure A1 we depict the evolution of the risk of infection in the PH regions of Algoma, North Bay Parry Sound, Northwestern, Sudbury, Renfrew, Thunder Bay, and Timiskaming, which are some of the lowest density populated regions of ON (upper panel, versus the regions of Durham, Halton, City of Hamilton, Niagara, Ottawa, Peel, Waterloo, York, Toronto, which are some of the highest populated (lower panel). A similar plot is presented in Figure A2, where we see the perceived discomfort with mask wearing among the two subpopulations.

Figure A1.  Change of $\bar{r}$ values over time for eight PHU regions with the lowest population densities (upper panel) and nine PHU regions with the highest population densities (lower panel).

DownLoad: Full-Size Img PowerPoint

Figure A2.  Change of $\tilde{r}$ values over time for eight PHU regions with the lowest population densities (upper panel) and nine PHU regions with the highest population densities (lower panel).

DownLoad: Full-Size Img PowerPoint

Last but not least, we looked at some socio-economic and demographic factors and their possible relation to how the population perceived the NPI's. We found that there are some correlations here, probably due to mobility and incidence, as seen in Figure A3.

Figure A3.  There is somewhat of a positive correlation between median age and the perceived risk of infection, and a somewhat negative correlation between median income and the perceived risk of infection (right panel). There are no correlations of these factors with the personal discomfort (left panel).

DownLoad: Full-Size Img PowerPoint