Calculating incidence of Influenza-like and COVID-like symptoms from Flutracking participatory survey data

We define three different Symptom Groupings:

• 1.

  CLI1+: Responses indicate any one or more of the above COVID-like symptoms. CLI1+ meets the Aotearoa New Zealand Ministry of Health advice to seek a COVID-19 test[19].

• 2.

  CLI2+: Responses indicate two or more of the above COVID-like symptoms. CLI2+ allows us to be slightly more discerning, given that many non-infectious illnesses such as allergies or asthma can cause a new onset of the COVID-like symptoms in the survey.

• 3.

  ILI: Responses indicate at least both cough and fever, which are symptoms of influenza-like illness. This Symptom Grouping is used for public Flutracking reporting [3].

We wish to estimate ‘incidence’, not ‘prevalence’ from the Flutracking data. Specifically, we wish to identify any new onset of symptoms that meet the threshold for an ‘incident’ according to a chosen Symptom Grouping. The method we use follows the same method employed by the Aotearoa New Zealand Ministry of Health and Flutracking Australia for estimating the incidence of ILI (cough and fever), but extends to the Symptom Groupings CLI1+ and CLI2+ above, in addition to the Symptom Grouping for ILI.

Flutracking data contains each respondent’s reported symptoms for a given survey week. For a selected Symptom Grouping, we first check whether an individual’s symptoms in the given week meet the criteria for that Symptom Grouping. If so, we mark it as an Incident. For example, if we wanted to investigate the number of CLI1+ Incidents, we would record all weeks where an individual recorded one or more of the six COVID-like symptoms as an Incident.

Figure 2 provides a hypothetical example of a respondent who reported different symptoms over the course of 8 weeks.

To distinguish new onset of illness from the continuation of a previously reported Incident, each Incident is assigned a unique ID, with this same Incident ID given to all consecutive weeks that meet the Symptom Grouping. If there is no response for one week between two weeks in which the Symptom Grouping was met, it is assumed that the interstitial week is a continuation of the previous week’s Incident and it is allocated the same Incident ID. For any gaps in responses between Incidents longer than one week we allocate the second Incident a new Incident ID. This procedure is consistent with the existing methodology used by Flutracking Australia and Flutracking New Zealand.

The incidence calculation is applied independently for each Symptom Grouping. This means that if a respondent reported one symptom in their first week of illness and two or more symptoms in the second week, the second week would be recorded as an ‘Incident’ of CLI2+, as though it were the onset of a new illness. This means that if a participant met the criteria for CLI1+ for three consecutive weeks, but met the criteria for CLI2+ only in the first and last weeks, this would be recorded as two Incidents of CLI2+, as shown in Figure 2. This is capturing useful information, as the advice is that someone with new or worsening symptoms should seek a COVID-19 test.

While our approach allows for flexible Symptom Grouping criteria to be implemented within the methodology, future work may seek to apply more sophisticated, unsupervised learning approaches when defining Incidents. While not in the scope of the current research, previous work has demonstrated that unsupervised learning approaches provide a useful technique for modelling heterogeneous symptom experiences among respondents [20].

Flutracking participants are more likely to respond to the survey when they are experiencing symptoms than when they are well[6]. This symptom reporting bias can lead to overestimating symptom incidence rates if it is not accounted for. Consequently, it is important to define a subset of consistent responses to use for incidence rate calculations. These are defined as being the responses, in a given week, from participants who are deemed to have responded consistently to the survey, as assessed at that week.

Previous studies have determined consistent respondents by only considering responses from those participants who have completed more than some minimum fraction of all surveys (e.g., [13, 8]). We improve on this approach by determining a subset of consistent responses for each week of the analysis, based on the response history of each survey participant for a defined window prior to the week under consideration. This allows a respondent’s ‘consistent’ status to change from week to week. We implement this by looking back over a specified number of weeks (the ‘window size’) and requiring the participant to have responded in all those weeks, with an allowance for a certain number of ‘missing weeks’ within the window. This makes it possible to adjust for the known reporting bias at any given point in time, while still maximising the sample size at the corresponding time. Our approach also means that analysis can be performed on a weekly basis throughout the year, rather than only at the end of a year (or season) of data collection.

In Figure 3 we plot trade-off curves showing the median fraction of responses excluded and the median relative change in the weekly incidence estimates as a function of ‘Window size’ and ‘Missing weeks allowed’. This helps us to determine suitable parameter values for the selection criteria of consistent responses. For more information on how these effects vary throughout the year, see Figures 12 and 13 in the Additional Information section.

Increasing the number of prior weeks for which a participant must have previously responded (the response window size) reduces the number of participants included in the ‘consistent response’ subset for any given week — that is, the fraction of responses excluded increases with increasing window size. Allowing respondents to have skipped some of the weeks in the consistent response window means that the responses of some participants are re-included in the consistent response subset.

When we consider the impact of window size and missing weeks allowed on the incidence estimate we find that a more stringent (i.e. larger) window size reduces the rate of symptom incidence calculated for any given week. This is expected from previous work looking at the impact of reporting bias [6]. The biggest jump is seen when requiring any response in the window prior to reporting symptom onset (i.e. a one-week window, or a two-week window with one missing week allowed, or a three-week window with two missing weeks allowed). Imposing this minimal consistent response window results in an immediate decrease of approximately 5% in the incidence estimate for all Symptom Groupings. We attribute this initial decrease to a correction for the known symptom reporting bias — if participants were equally likely to respond in a given week, independent of symptom incidence, then applying this criterion would not affect the estimated incidence rate.

Increasing the consistent response window, beyond the minimum, initially causes a decrease in the number of consistent responses but does not significantly change incidence estimates. However, further increasing the window size beyond about four weeks causes both the number of consistent responses and the estimated incidence rate to decrease. We attribute this to the more stringent consistent response window introducing a different form of bias by increasing the over-representation of some groups of participants in the sample population — specifically older Pākehā/NZ European cohorts who tend to have lower incidence rates — see for example Figure 5.

While our method allows for the user to choose any consistent respondent window size and any number of missed weeks within that, we have chosen in our subsequent analysis to use a window size of four weeks with up to one missing week. This defines a consistent response subset where small changes in the window size do not result in significant changes in the resulting incidence rate and where the fraction of responses included in the sample population for a given week is not too severely reduced, in part due to the fact that the Flutracking survey allows respondents to retrospectively enter results for up to four weeks into the past. We provide two examples of the selected values for window size and allowed missing responses using hypothetical response series illustrated in Figure 4.

There are several demographic biases in the Flutracking cohort with registration rates differing by age, ethnicity and location of respondents. If any of these self-selection biases also align with a corresponding variation in incidence rates along the same demographic factor, then it will contribute to an under- or over-estimate of incidence for the total population.

The strongest factor affecting incidence estimates is age (see Figure 5), with younger ages, particularly 0–4 years, tending to have higher incidence across CLI1+, CLI2+, and ILI Symptom Groupings. Therefore, we adjust our incidence estimates to account for the relative number of respondents in the survey population according to age each week.

Although there is a relationship between ethnicity and incidence of several respiratory illnesses in Aotearoa New Zealand [21, 22], it can be unhelpful to treat ethnicity as a predictive factor for these illnesses. Research has found that factors including socioeconomic deprivation, household crowding, and housing quality, are all associated with increased incidence of respiratory infections [22, 23, 24], and that these are highly correlated with ethnicity [23]. These same underlying factors are also likely to contribute to a participant selection bias within the Flutracking survey. Under-representation of lower socio-economic groups has also been noted within similar participatory surveys for Influenza. [9]

The low response rate in Aotearoa New Zealand Flutracking data for ethnic groups other than Pākeha/NZ European, combined with the under-representation of lower socio-economic groups, mean that it is unlikely that the sampled population of non-Pākeha ethnicities in Flutracking is truly representative of those groups in the wider population. Therefore, weighting based on ethnicity alone may actually amplify selection biases rather than mitigating them. For this reason, although it is possible using our method, we do not calculate or re-weight symptom incidence by ethnicity. We note that the current re-weighting methodology is easily able to be applied to ethnicity or socio-economic status, should Flutracking enrolment be sufficiently broad that responses are representative of these groups. In order to avoid this limitation in future studies, we point to previous research on best practice for participatory surveillance data collection [25, 26].

The reference population used for the age weightings presented in this method is the 2021 Estimated Residential Population, calculated by Statistics NZ with 5-year age groups and spatial units of District Health Board (DHB) of residence [27]. However, any other source of reference population data can be used, provided that it can be mapped on to the age (1-year age bands) and location (Aotearoa New Zealand postcode) groupings used in the Flutracking data. In general, it is better to use higher levels of aggregation (e.g. 5-year age groups) and DHB. This is in order to avoid the low response numbers per analysis unit, and consequent lower statistical power and increased uncertainty, that a more granular aggregation will lead to.

In order to weight survey responses according to participant age, we use the rake function provided in R’s survey package [1] to assign a weighting coefficient to each response. The process of raking matches the marginal distribution of ages from the survey respondents with those of a reference population.²²2Since this is specified in 5-year age groups, sets of ages must also consist of combinations of 5-year age groups. The raking can be applied at any level of aggregation appropriate to the reported results, with the caveat that fine-grained aggregation can lead to small cell sizes which may limit the statistical power of results in some cases.

When analysing responses for sub-national spatial units (e.g. DHBs), we calculate the age weighting factor for each age group with respect to the age distributions for each of the corresponding sub-national units. However, for analysis at a national level, we calculate only national level age weightings, rather than combining sub-national weighting factors. This avoids the pit-fall of national level weighting factors inheriting lower statistical power and increased uncertainty from small cell sizes for some location-age combinations.

Constructing confidence intervals for weighted survey data is non-trivial. We construct estimates of mean incidence with 95% confidence intervals for each survey week for a given set of demographic factors and Symptom Groupings using the survey package in R [1], incorporating the age weighting previously computed in Section 1.3.4. The 95% confidence intervals for each survey week can be turned into 95% pointwise confidence bands across survey weeks.

We first create a ‘survey design’ object using our set of consistent responses, their incidence classification and weighting, and other relevant factors such as their survey week, age group, and location. This survey design object is then passed to the svyby function, which independently calls the svyciprop function for each subset in a set of factors and outputs the estimated mean incidence and 95% confidence intervals for each subset. svyciprop has several methods of producing confidence intervals; we use its logistic regression method ("logit"). An example set of factors could be Symptom Grouping and survey week, as in Figure 6. One independent subset used in these calculations is the Symptom Grouping ‘CLI1+’ and the Survey Week ending on Sunday the 31st of May 2020.

The process used by svyby and svyciprop is equivalent to fitting a weighted logistic regression model, where every subset in a set of factors is fitted independently such that the fitted value of any subset in the set of factors is equal to the weighted mean incidence of that subset. Thus the same model could have been fitted using the svyglm function from the survey package and then fitting interaction effects between all factors. In order to do this, a quasi-binomial model family needs to be used in svyglm [1] because the use of survey weights may make the number of ‘Incidents’ non-integer. However, as we were mainly concerned about the fitted means and confidence intervals for each subset rather than fitted values of regressors and their significance, we did not do this in general. Additionally, our process makes it very simple to extract confidence intervals for our estimates, which, while not too difficult (e.g. using a non-parametric bootstrap [28]), is not nearly as straightforward to do with the fitted model returned by svyglm.

The only exception to this, is that in order to consider the significance of differences in incidence for a given Symptom Grouping, as seen in Figure 9, we did need to use svyglm to test whether the difference in the estimated mean incidence between two locations on a given week was statistically significant (i.e. if the p-value for the Auckland Metro DHBs coefficient was less than 0.05).