Study info

The WCIS was jointly run by the UK Health Security Agency and the Office for National Statistics. A cohort of randomly sampled households was recruited and asked to complete questionnaires monthly from 13 November 2023 to 7 March 2024. Participants were cohorted into weekly waves, with staggered start dates for each survey window. There were 111,752 participants, with 386,343 submissions completed during the study. Not all participants responded in every wave of the survey (Supplementary Table 3). The cohort consisted of adults aged three and over in England, living in residential settings. Participants were representative of all regions, though skewed towards older and wealthier individuals (Supplementary Table 4).

Along with demographic information, participants were asked questions about their symptomatic status, such as ‘Have you had any of the following symptoms in the last 7 days?’ and ‘Which of these symptoms are new or worse than usual for you?’ The count and percentage of the total symptomatic statuses reported by respondents are given in Supplementary Table 5. Multiple modelling steps are required to produce a representative symptom incidence from this data.

To compare symptom incidence, we also use the symptom search dataset (SSD) from Google Research [11]. The data represent a daily national-level relative search volume for over 400 symptoms, originally set up to monitor SARS-CoV-2. The SSD represents a signal of individuals in aggregate searching for information about symptoms, relative to total search volume.

In this work, we explore the following symptoms relating to winter respiratory illnesses: ‘cough’, ‘fever’, ‘loss of smell’, ‘shortness of breath’, and ‘sore throat’.

Google incidence

Throughout, we assume that a Google search for a symptom is a proxy for symptom onset as an individual seeks to learn more about their new symptom. Whilst true incidence itself can vary smoothly over time at a population level, observations of symptom incidence are influenced by behaviour changes and reporting effects. Namely, healthcare-seeking behaviour varies by day of the week and can be impacted by holiday periods. We therefore fit a generalized additive model (GAM) to adjust for these known biases in the reporting of syndromic surveillance Google SSD.

The relative search volume (RV) at time t, V(t), is continuous and bounded by (0,100). The rescaled search volume, denoted S(t) = V(t)/100 ∈ (0,1), is assumed to follow a Beta distribution

$$ S(t)\sim \mathrm{Beta}\left(\unicode{x03BC} \left(\mathrm{t}\right)\unicode{x03D5}, \left[1-\unicode{x03BC} \left(\mathrm{t}\right)\right]\unicode{x03D5} \right), $$

$$ \mathrm{logit}\left(\unicode{x03BC} \left(\mathrm{t}\right)\right)={\beta}_0+\mathrm{GP}(t)+{f}_{\mathrm{DoW}}\left(\mathrm{DoW}(t)\right)+{\beta}_hh(t), $$

where,

• • $ \mu (t) $ is the $ \mathrm{Beta} $ mean on the logit scale; $ \unicode{x03D5} $ is a single precision parameter estimated by restricted maximum likelihood (REML).

• • $ \mathrm{GP}(t) $ is a Gaussian-process spline that imposes a squared-exponential correlation structure decaying with the squared time difference between days.

• • $ {f}_{\mathrm{DoW}} $ is a cyclic cubic spline of day-of-week.

• • $ h(t)=1 $ on UK public holidays, 0 otherwise.

The terms $ {f}_{\mathrm{DoW}} $ and $ {\beta}_hh(t) $ are marginalized out, providing a smooth estimate of RV over time. Prediction samples are generated from the model coefficients using a Metropolis-Hastings approach, giving $ J=1000 $ samples of the symptom incidence, $ {V}_j(t) $ , from which the median and 90% prediction intervals are generated. The model is fit using mgcv [Reference Wood29] and prediction samples extracted with gratia [Reference Simpson30].

WCIS incidence method

To compare Google Trends with a reliable baseline, we require the daily incidence of new symptoms from WCIS. Note that, whilst we could straightforwardly estimate the incidence of symptoms within the study participants from their symptom onset dates, it would not be possible to directly extend this to a symptom incidence rate across the whole population. This would require knowledge of participation rates, which cannot be defined with respect to symptom onset dates. In contrast, we can more easily estimate prevalence rates because we have data on both individuals who do have the symptom and those who do not at each time point. Furthermore, the symptom onset date provided in the study is a single onset date for all symptoms, which forces the assumption of concurrent onsets of each symptom.

Therefore, incidence is not observed directly but can be reconstructed from the prediction draws of prevalence and symptom duration.

Let $ {\pi}_{a,j}(t) $ be the $ j $ -th prevalence prediction draw for age group $ a $ on day $ t $ , and let $ {\overline{d}}_{a,j} $ be the corresponding draw of the mean symptom duration. If we approximate every episode as lasting exactly $ {\overline{d}}_{a,j} $ days, then prevalence and incidence are linked by

$$ {\pi}_{a,j}(t)=\sum_{\tau =0}^{{\overline{d}}_{a,j}-1}{I}_{a,j}\left(t-\tau \right). $$

Therefore,

$$ \hskip0.24em {I}_{a,j}\left(t-\frac{{\overline{d}}_{a,j}}{2}\right)\approx \frac{\pi_{a,j}(t)}{{\overline{d}}_{a,j}}. $$

This is calculated by dividing each prevalence draw by its paired mean duration and shifting the resulting incidence curve backwards by $ \mathrm{round}\left({\overline{d}}_{a,j}/2\right) $ . This produces a set of incidence draws $ {I}_{a,j} $ that match the daily resolution of the WCIS data and carry the full uncertainty from both prevalence and duration estimates. Further information on the derivation of incidence can be found in Supplementary Section 2.

WCIS prevalence

We estimated the proportion of participants who were currently symptomatic by combining a multilevel-regression-and-poststratification (MRP) framework [Reference Downes31] with Gaussian-process smoothers previously developed for influenza-like illness [Reference Mellor32]. Sub-groups are defined by age group $ a $ , region $ r $ (nine English regions), and invitation window index, $ w $ , defined to be the time between the response window start date and the survey submission date.

For each calendar day, $ t $ , the survey yields $ {n}_{a,r,w}(t) $ symptomatic responses out of $ {N}_{a,r,w}(t) $ returned questionnaires. We assume the number of symptomatic responses follows a binomial distribution with $ {N}_{a,r,w}(t) $ as the number of trials and prevalence $ {\pi}_{a,r,w}(t) $ as the probability of success:

$$ {n}_{a,r,w}(t)\sim \mathrm{Binomial}\left({N}_{a,r,w}(t),{\pi}_{a,r,w}(t)\right), $$

$$ \mathrm{logit}\left({\pi}_{a,r,w}(t)\right)={\beta}_0+{f}_1(t)+{f}_2^a(t)+{b}_r+{\gamma}_w, $$

where,

• • $ {f}_1(t) $ is the national GP smoother

• • $ {f}_2^a(t) $ is an age-specific GP smoother

• • $ {b}_r $ is a Gaussian random intercept for region

• • $ {\gamma}_w $ is a categorical fixed effect for time from response window start to participant submission date.

Prediction sampling and marginalization

We draw $ J=\mathrm{1,000} $ prediction samples from the model using gratia, which randomizes all smoothing parameters via a Metropolis–Hastings step. To obtain window-marginal prevalences, we weight each $ \pi (t) $ by the frequency of its window:

$$ {\pi}_{a,r,j}(t)={\sum}_w{\delta}_w{\pi}_{a,r,j}(t), $$

where,

$$ {\delta}_w=\frac{\mathrm{number}\ \mathrm{of}\ \mathrm{responses}\ \mathrm{in}\ \mathrm{window}\ \mathrm{in}\mathrm{dex}\;w}{\mathrm{total}\ \mathrm{responses}}. $$

Poststratification

Census population totals $ {p}_{a,r} $ are treated as fixed. For every draw $ j $

$$ {\pi}_{\mathrm{nat},j}=\frac{\sum_a{\sum}_r{\pi}_{a,r,j}(t)\;{p}_{a,r}}{\sum_a{\sum}_r{p}_{a,r}}, $$

$$ {\pi}_{a,j}(t)=\frac{\sum_r{\pi}_{a,r,j}(t){p}_{a,r}}{\sum_r{p}_{a,r}}. $$

The resulting national time-series $ {\pi}_{\mathrm{nat},j}(t) $ provide 1,000 prediction draws of daily national and age-specific prevalence, used for downstream incidence calculation.

WCIS symptom duration

We modelled the distribution of symptom durations reported by WCIS participants using the R package epidist (v0.3.0 [Reference Howes, Park and Abbott33]) via Bayesian MCMC, allowing for interval censoring and right truncation [Reference Ward34, Reference Park35]. Let $ {D}_i $ be the true duration (in days) of episode $ i $ reported by a participant in age group $ a $ . We assume a log-normal duration model:

$$ \log \left({D}_i\right)\sim N\left({\mu}_a,{\sigma}_{,a}^2\right), $$

with age group specific parameters

$$ {\mu}_a={\beta}_{\mu, 0}+{\beta}_{\mu, \mathrm{a}}, $$

$$ \log \left({\sigma}_{\mathrm{a}}\right)={\beta}_{\sigma, 0}+{\beta}_{\sigma, a}, $$

where $ {\beta}_{\mu, 0} $ and $ {\beta}_{\sigma, 0} $ are global intercepts, and $ {\beta}_{\mu, \mathrm{a}},{\beta}_{\sigma, a} $ are fixed effects for each of the seven age bands. We select priors for the intercept terms with mean and standard deviation parameterization of a normal distribution,

$$ {\beta}_{\mu, 0}\sim \mathrm{Normal}\left(\log (4),0.5\right), $$

$$ {\beta}_{\sigma, 0}\sim \mathrm{Normal}\left(\log (0.7),0.5\right), $$

and uninformative priors for the age-specific fixed effects,

$$ {\beta}_{\mu, \mathrm{a}},{\beta}_{\sigma, a}\sim \mathrm{Normal}\left(0,3\right). $$

Censoring and truncation

• • Onset times were interval censored to $ \left[{d}_i,{d}_i+1\right) $ days, where $ {d}_i $ is the date that the participant reported symptom onset.

• • End times were right-truncated to lie in $ \left[\mathrm{submissio}{\mathrm{n}}_i+1,\mathrm{submissio}{\mathrm{n}}_i+4\right) $ days, allowing up to 3 days of unreported recovery.

• • Episodes implying durations >60 days (0.4% of records) were excluded to remove potential chronic conditions.

For each posterior draw j and age group a, we computed the mean duration

$$ {\overline{d}}_{a,j}=\mathit{\exp}\left({\mu}_{\log, a,j}+\frac{1}{2}\;{\sigma}_{\log, a,j}^2\right), $$

We summarize the distribution of $ {\overline{d}}_{a,j} $ by its posterior median and 90% credible intervals.

We assumed that symptoms started within a one-day window after the symptom onset date provided by the participant and ended at most 3 days after the survey submission date. We explored the sensitivity to the choice of window length for the symptom end date in Supplementary Figure 5. We found that increasing the window length increased the estimate of the mean duration across age groups, and inclusion of a short window reduced the uncertainty in our estimates.

Furthermore, individuals were surveyed multiple times throughout the study period; however, we assumed these surveys were independent due to being sufficiently far apart (around one month between responses) for individuals to have recovered from an infection event. In addition, some respondents gave symptom onset dates many months in the past, despite the survey focusing on symptoms that have recently developed (or recently become worse than usual); a substantial number gave a symptom onset date of exactly one year ago, the earliest date allowed to be inputted. To avoid these longer-term, chronic symptoms skewing our duration distributions, we excluded individuals with symptom durations longer than 60 days.

WCIS incidence implementation

For each age group, we have $ J=1000 $ prediction draws of

• • Prevalence $ {\pi}_{a,j} $ (t),

• • Mean symptom duration $ {\mu}_{a,j}. $

Because the duration model is fitted independently of the prevalence model, we treat the two sets of draws as independent and pair them by index $ j $ .

For every pair ( $ {\pi}_{a,j}(t),{\overline{d}}_{a,j} $ ), we approximate incidence as

$$ {I}_{a,j}\left(t-\frac{1}{2}\;{\overline{d}}_{a,j}\right)=\frac{\pi_{a,j}(t)}{{\overline{d}}_{a,j}}, $$

and then shift the series backwards $ \mathrm{round}\left({\overline{d}}_{a,j}/2\right) $ days.

The back-shift differs across draws, so we restrict each $ {I}_{a,j}(t) $ to the common window:

$$ {t}_{\mathrm{min}}=\min (t)-\frac{1}{2}\;\max \left({\overline{d}}_{a,j}\right),{t}_{\mathrm{max}}=\max (t)-\frac{1}{2}\;\max \left({\overline{d}}_{a,j}\right). $$

Age-specific incidence is summarized at each $ t $ by the median and 90% prediction interval of the aligned draws, giving $ {\overset{\sim }{I}}_a(t) $ .

National incidence for draw $ j $ is

$$ {I}_j(t)=\sum_a{I}_{a,j}(t), $$

and is summarized to yield $ \overset{\sim }{I}(t) $ with its 90% prediction intervals.

Growth rates

In addition to comparing the Google symptom incidence, $ {V}_j(t), $ against WCIS symptom incidence, $ {I}_j(t) $ , we also compared the growth rates for each sample. Since we have discrete values for the incidence at each time step, $ t $ , we use the following approximation for the WCIS growth rate

$$ {r}_j(t)=\frac{I_j(t)-{I}_j\left(t-1\right)}{I_j\left(t-1\right)}\ast 100, $$

and similarly for Google symptom incidence.

The growth rates across samples were summarized to give median estimates and 90% prediction intervals.

Cross-correlation function (with uncertainty)

Each Google incidence sample, $ {V}_j(t) $ , is paired with a WCIS incidence sample, $ {I}_j(t) $ , for each symptom, with their means defined as $ \overline{V_j} $ and $ \overline{I_j} $ , respectively. The cross-correlation was calculated for each pair of samples, for lags, $ k $ , between −17 and +17 days,

$$ {C}_j(k)=\frac{1}{n}\sum_{t=1}^{n-k}({V}_j(t)-\overline{V_j})({I}_j(t+k)-\overline{I_j}). $$

Here, a positive lag compares WCIS incidence k days in the future against Google Trends. The cross-correlations across samples were then summarized to give a median estimate, $ C(k) $ , and a 90% prediction interval for each lag considered.

Similarly, we calculated cross-correlations for Google and WCIS incidence growth rates, again summarizing to give a median estimate and 90% prediction intervals for each lag considered.