Mining Google Trends data for nowcasting and forecasting colorectal cancer (CRC) prevalence

Data

CRC prevalence data were extracted from the publicly available database of the Our World in Data (OWID) platform via the owid() function embedded in the dedicated R software package (https://ourworldindata.org/). At its turn, OWID sources CRC prevalence data from the Global Burden of Disease (GBD) database of the Institute for Health Metrics and Evaluation (IHME) (https://www.healthdata.org/gbd), which is thus the original cancer statistics source. The CRC prevalence data covers the 2004–2019 period (i.e., T = 16) for a total of 202 individual countries and country panels. First, heterogeneity across world countries is assessed through “plotmeans” style plots with the R’s “gplots” package of Gregory, Warnes & Lodewijk (2016). This method carries the non-trivial advantage of plotting group means along with confidence intervals (CI), which, together with descriptive statistics reported in Table 3, offers an insightful view of the extent of heterogeneity present in the data. Of note, the t distribution is used by default to compute confidence intervals within the plotmeans() function. However, as Rovetta (2023) explains, confidence (or credible) intervals concern the mean distribution and generally are not the best indicators for the dataset variability.

Moreover, the heterogeneity in worldwide CRC prevalence is further explored by employing both income-based and geography-based panels, as per the Word Bank (WB) classification.

This approach is particularly helpful for informing global equitable health policies. Consequently, four income-based panels and their corresponding CRC prevalence rates are extracted. Subsequently, seven geography-based panels are also delineated, and the relevant observations are sourced. Table 2 reflects the income and geographical panels extracted from OWID, along with the abbreviations employed for the latter category.

The Google Trends platform (https://www.google.com/trends/) reports the relative search volume index (RSI) of a keyword in a region, over a specified time period, with an RSI of 100 reflecting the point of peak popularity (Silva et al., 2019). As such, the RSI is estimated as per Eq. (1) and subsequently scaled from 0 to 100 (Tudor & Sova, 2022).

(1) $$RSI={\displaystyle \frac{Totalsearchvolumefor“keyword”submittedinthespecifiedgeographyoverthespecifiedperiod}{Totalsearchessubmittedinthespecifiedgeographyoverthespecifiedperiod}}$$

Google Trends data is sourced in this study by making use of R’s “gtrendsR” package (Massicotte & Eddelbuettel, 2022), calling the gtrends() function, specifying alternatively “colorectal cancer” and “colonoscopy” as the search terms of interest, the default geography of “all” corresponding to global queries, and specifying the “health” category (i.e. category no. 45) to assure reliable results. Moreover, the string specifying the time span of the query is set to “all,” which implies that monthly RSI data starting on January 2004 (i.e., the beginning of Google Trends) and spanning to December 2022 is extracted. Consequently, a time series of 229 monthly RSI observations is available and will be employed for modeling and forecasting purposes in this study. Tudor (2022a) offers further relevant details on the methodology behind the RSI index. The “gtrendsR” package is not only a great tool for extracting GT data, but allows the specification of multiple relevant elements, is capable of performing an interest-by-region analysis and is useful for making searches reproducible (however, as explained later, sample instability is an issue that must be accounted for). It should be mentioned that the ISO language code only influences the data returned by related topics, as per the package description, and thus does not impact the sourced data.

Method

This study develops a three-step framework to forecast monthly CRC prevalence evolution based on Google Trends relative search index.

First, a monthly model of CRC prevalence is estimated based on RSI for the world panel (as per World Bank), such that:

(2) $$CRC=f\left(RSI\right)$$where the CRC prevalence in year i is modelled as a linear function f of the same-year RSI, plus some white noise.

To perform the estimation, the monthly RSI index is first converted to annual averages to match the frequency of the CRC prevalence series imported from the OWID database. Moreover, the parameters are estimated via the ordinary least squares (OLS) estimator. Of note, the main research interest does not lie on fine-tuning the CRC-RSI relationship, but rather on the second step of the framework, i.e., accurate forecasting of RSI. Additionally, particular care is given to the choice of the estimation period, juggling the need to include enough observations and the fact that the Google platform has taken off as a main search engine only in the most recent decade. As per Britannica (https://www.britannica.com/topic/Google-Inc), by 2004, Google was processing daily 200 million searches, increasing to over 3 billion by the end of 2011. As of March 2023, Google’s market share of the global search market across all devices is nearly 84%, making it the dominant search engine worldwide (Statista, 2023), which translates into 9 billion daily searches (EarthWeb, 2023). In light of these figures and aiming to increase the relevancy of current findings, the estimations performed within the first step of the framework (i.e., to assess the link between distinct RSI pertaining to alternative relevant keywords and official CRC prevalence statistics), have been limited to the most recent decade of available data. However, this in turn reduces the number of observations in the data sample employed in the first step within the proposed framework.

Second, the RSI time series is modeled in-sample, tested out-of-sample, and ultimately forecasted and validated against official statistics and/or results of other studies, where available. To this end, multiple tools are employed.

As such, the nonlinearity of monthly RSI data is examined by two alternative tests: the Tsay test for quadratic nonlinearity (Tsay, 1986), as well as the likelihood ratio test for threshold autoregression (LR) proposed by Chan (1991). In both cases, the null hypothesis is that the time series includes an AR process (Munim, 2022). Next, for modeling and forecasting purposes, the RSI index is split into a training set (made up of 80% of observations) used for in-sample fitting of three alternative models (i.e., ETS, ARIMA, and FNNAR) and a testing set (containing the remaining 20% of the observations) on which the forecasting ability of the best-fit model in-sample in each category (i.e., best ETS, best ARIMA, best FNNAR) is comparatively assessed. For this purpose, several forecasting accuracy statistics (Hyndman & Athanasopoulos, 2018) are calculated (shown below). Moreover, the Diebold and Mariano (DM) test for superior forecasting ability (Diebold & Mariano, 1995) is estimated to confirm the superior forecasting ability of the best forecasting model relative to the second ranked candidate.

Root mean squared error:

(3) $$RMSE=\sqrt{mean(e_t^2})$$

Mean absolute error:

(4) $$MAE=mean\left(\left|e_t\right|\right)$$

Mean absolute percentage error:

(5) $$MAPE=mean\left(\left|p_t\right|\right)$$where $p_t={\displaystyle \frac{100e_t}{y_t}}$

Mean absolute scaled error:

(6) $$MASE=mean\left(\left|q_j\right|\right),$$where $q_j={\displaystyle \frac{e_t}{\frac{1}{N-1}{\sum }_{i=2}^N\left|y_t-y_{t-1}\right|}}$ in the case of non-seasonal series and $q_j={\displaystyle \frac{e_t}{\frac{1}{N-m}{\sum }_{i=m+1}^N\left|y_t-y_{t-m}\right|}}$ when the time series is seasonal.

In all cases, the forecast error is defined as:

(7) $$e_{T+h}=Y_{T+h}-{\hat{Y}}_{T+h|T}$$

To sum up, the data splitting rule for the training-testing RSI modeling and forecasting is given in the script snippet in Box 1.

Simple statistical methods have dominated the field of time series forecasting for many decades, both in academia and industry, due to their relatively accurate, fast-to-compute, and interpretable forecasts, with neural networks emerging as especially accurate in univariate time series forecasting (Semenoglou, Spiliotis & Assimakopoulos, 2023). However, due to their capability of revealing the influence of many parameters, multivariate time series models have uncontested advantages, but their goodness of fit is intrinsically related to the quality and availability of data related to these parameters (Atchadé & Sokadjo, 2022). Thus, given that the more complex the model, the more the need for data, and the lesser the availability and quality of cancer prevalence statistics, especially in less-developed economies and during the COVID-19 pandemic, we rely on the best-identified univariate model in the current research. Our choice is further backed by the findings of Ziel & Weron (2018) that conclude that, for the task of electricity price forecasting, the multivariate modeling approach does not uniformly outperform the univariate one across 12 distinct electricity spot price datasets and is even often outperformed by the latter. Furthermore, we rely on the findings of Jaidka et al. (2021) that study 208 Designated Market Areas (DMAs) throughout the United States and find that information-seeking behavior from Google Trends data is 19% more accurate than sociodemographic and regional controls at predicting the prevalence of noncommunicable diseases. Consequently, as per Ziel, Steinert & Husmann (2015), we argue that our approach can surpass the complexity of cancer prevalence forecasting by offering several advantages relative to multivariate settings, such as being capable to model the search volume index and to provide accurate forecasts without the need for any data manipulation nor for any additional information, and, not in the least, by employing efficient and rapid state of the art estimation techniques.

As mentioned earlier, artificial neural networks (ANNs) have emerged as particularly useful for modeling and forecasting complex systems while accounting for non-linearities. The FNNs, are the most extensively utilized neural network (NN) models (Hsieh, 2004). A neural network is made up of several small computational units or nodes that are stacked in layers and run in parallel, the, connections among nodes, and an activation function (Chen & Billings, 1992; Allende, Moraga & Salas, 2002). Figure 3 provides the basic architecture of an ANN, which is made up of an input layer X, linked to at least one hidden H, which is, in turn, linked to the output layer Y (Chen & Billings, 1992).

The feedforward network is a type of neural network in which the input feeds forward through the network layers to the output. In the case of a feedforward neural network autoregression (FNNAR) model, the input layer consists of the past observations of a series (Munim, Shakil & Alon, 2019). Consequently, the vector of predicted values Y is found as follows:

(8) $$\mathrm{Y}=\mathrm{f}(\mathrm{H})$$where

(9) $$\mathrm{H}=\{\mathrm{w}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{x}[(\mathrm{p}\ast \mathrm{k})]\ast \mathrm{X}+\mathrm{b}\mathrm{i}\mathrm{a}\mathrm{s}\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r}$$and f is the activation function, with p standing for the nonseasonal inputs for the linear AR process, and k for the number of nodes in the hidden layer.

FNNAR models are fitted through automated algorithms in the “forecast” package in R software. In particular, the nnetar() function makes 25 repetitions and finds parameters through AIC minimization. By default, the number of nodes in the hidden layer is calculated as k = (p + P + 1)/2.

Additionally, to assure the robustness of findings, ETS and ARIMA models are automatically fitted in the R environment.

The exponential smoothing state space model (ETS) is based on the work of Holt (1957), Winters (1960), and Brown (1959). Predictions given by ETS are weighted averages of prior data, with greater weights assigned to more recent observations and weights decreasing exponentially (Hyndman & Athanasopoulos, 2018; Silva et al., 2019). ETS estimates comprise three components: the trend (T), seasonal (S), and error (E) components. In addition, the trend includes both a level term (l) and a growth term (b). According to Yang et al. (2015), the trend and seasonal components might be none (N), additive (A), additive damped (Ad), multiplicative (M), or multiplicative damped (Md). Hence, an ETS model is represented by a three-character string (Z,Z,Z) (Perone, 2021), with the first Z representing the error assumption of the state-space model and the second and third Zs representing the trend type and the season.

The use of the automated ETS technique via the ets() function inside the forecast package in R (Hyndman & Khandakar, 2008; Hyndman et al., 2022) assesses 30 ETS equations during the modeling procedure i.e., all 30 formulae are provided by Hyndman & Khandakar (2008).

The approach optimizes the smoothing parameters and the starting state variable and employs a penalized likelihood, i.e., the corrected Akaike’s Information Criterion (AICc), to choose the best model on the training set, which is then used to generate point predictions.

Autoregressive integrated moving average (ARIMA) models were created by Box & Jenkins (1970) and are among the most widely used parametric time series analysis and forecasting methods (Silva et al., 2019).

In equation form, according to Hyndman & Athanasopoulos (2018), an ARIMA (p,q,d) (P,Q,D) seasonal model is given as:

(10) $$\begin{array}{rl}(1-{\phi }_{1}B-...-{\phi }_p{B}^p)& (1-{\mathrm{\Phi }}_1{B}^s-...-{\mathrm{\Phi }}_P{B}^{sP})(1-B{)}^d(1-B^s{)}^D{Y}_t\\ & =(1-{\theta }_{1}B-...-{\theta }_q{B}^q)(1-{\mathrm{\Theta }}_1{B}^s-...-{\mathrm{\Theta }}_P{B}^{sQ}){\epsilon }_t\end{array}$$where s is the seasonal period, the lowercase and the capital letters represent nonseasonal and seasonal parameters, whereas ${\epsilon }_t$ is a zero-mean random variable with the standard deviation $\sigma$. The “auto.arima” function from the “forecast” package in R software conducts ARIMA modeling using an automated and optimized method that is capable of efficiently traversing the space of models in order to identify the ideal model. The function determines if seasonal differencing is required for the training data, calculates unit root tests, and selects model parameters using a step-by-step AICc reduction procedure (Tudor, 2022a, 2022b).

The best forecasting model on the out-of-sample set is used to issue point forecast for the RSI index setting a forecasting horizon h of 24 months, thus reaching the end of 2024. Of note, h is chosen by complying with two important criteria in time series forecasting. First, ideally the forecasting horizon must be shorter than the testing window (Hyndman & Athanasopoulos, 2018). Second, longer-term forecasts have been shown to fail do provide superior information relative to the mean of the explained variable (Breitung & Knüppel, 2021). Consequently, setting the length of h equal to 24 assures that both criteria are met.

Third, the function f estimated at step one is applied to the RSI series forecasted at step two, under the no-change hypothesis, in order to produce forecasts for the CRC prevalence rate up to the end of 2024, as follows:

(11) $$\widehat{CRC}=\widehat{\text{f}}\left(\widehat{\text{RSI}}\right)$$

To sum-up, the proposed framework embeds the sequential main steps depicted in Fig. 4.

Worldwide trends and heterogeneity in CRC prevalence

Figure 5 reveals that the global panel presents a rapidly increasing CRC prevalence rate over the last decades, which is mostly driven by the upper-middle income panel, confirming that the rising number of CRC cases is creating an increasing global public health concern (Xi & Xu, 2021).

The graphical depiction further shows that the prevalence rate of CRC in high-income countries (i.e., as per the World Bank classification) has consistently been about 10 times that of low- and middle-income economies. Moreover, it should be noted that the growth rate of CRC prevalence is significantly steeper in upper-middle income countries than in both rich and poor countries, suggesting increased exposure to CRC risk factors, albeit at a high basal level. Nevertheless, it is increasing rapidly in less developed countries due to increased exposure to CRC risk factors.

Furthermore, there is also high geographical heterogeneity in CRC prevalence at the world level (Fig. 6). Thus, data indicates that the East-Asia& Pacific region is registering both a high basal CRC prevalence level, as well as a high growth rate, confirming that immediate measures are paramount to offer some relief against this burden.

Additionally, the heterogeneity amongst individuals is also apparent in Fig. 7, which denotes the mean CRC prevalence level over the sample period for each of the 202 countries included in the analysis. The plot is issued through the plotmeans() function in the “gplots” package, which uses the t distribution to compute confidence intervals, with a 0.95 confidence level for error bars.

The high country heterogeneity also emerges from the descriptive statistics for the entire sample, which are reported in Table 3. Estimations show a very high standard deviation (657,663.06) and huge range (11,457,627) for the levels of CRC prevalence for world countries during 2004–2019, whereas its mean level was 86,722.88.

Relevancy of alternative RSIs for CRC prevalence

Figure 8, which plots the trends of the CRC prevalence and of the RSI “colonoscopy” in panel a, and the scatterplot between the two variables in panel b reveals a disaggregation between the two series during the first years in the sample, when the RSI series presents a decreasing trend, whereas the CRC prevalence series shows an uninterrupted increasing trend throughout the analysis period. This in turn leads next to the implementation of a joinpoint (or change-point) regression analysis (Gillis & Edwards, 2019) to confirm and assess the changing trends that emerge from the visual inspection. Of note, this kind of analysis has been routinely applied to detect changing trends in cancer mortality and/or incidence series (see, among others, Qiu et al., 2009; Crispo et al., 2013; Sarakarn et al., 2017; Wilson, Bhatnagar & Townsend, 2017). Results of the joint-point regression analysis, graphically reflected in Fig. 9, confirm that the RSI series presents one join point in 2014, leading to two periods with distinct trends, as follows: a negative trend until 2014, and an increasing trend during 2014–2019 (Fig. 9). On the other hand, the joinpoint regression analysis also found one joinpoint in CRC prevalence series around 2014, but in this case the joinpoint delineates two positive trends. All estimations have been performed R’s software “segmented” package. Hence, given the disaggregation between the series during the first segment, to comply with the linearity assumption, we subset the secondary segment to implement the linear models and estimate the statistics of interest.

Estimation results for the two linear models are summarized in Table 4. Of note, as the RSI indexes are scaled from zero to 100 by construction, to make the data comparable and obtain more consistent estimates, the CRC prevalence rates are divided by 100,000 prior to estimating the model. First, it should be noted that effect sizes are the paramount outcome of empirical statistical research (Lakens, 2013). As explained by Lakens (2013), which in turn cites Rosenthal (1994), effect sizes are usually classified into two families: the d family (i.e., standardized mean differences) and the r family (i.e., measures of strength of association), whereas the d family effect sizes are conceptually determined by the difference between observations divided by the standard deviation of these observations, and the r family effect sizes characterize the proportion of variance that can be attributed to group membership. Moreover, as Lee (2016) explains, p-values can only inform of the statistical significance, but not the magnitude of the effect, whereas CIs mitigate, but do not solve, this issue by providing a range of possibility. Furthermore, it should be mentioned that the Shapiro-Wilk test (Shapiro & Wilk, 1965) estimated for the three series could not reject the null hypothesis that all samples stem from a normal distribution.

In light of these considerations, we estimate and report the Pearson correlation coefficient and the value of Cohen's d statistics (Cohen, 1988) together with its correction for sample bias (i.e., known as Hedges’ g) between the CRC prevalence series and each RSI series, including the 95% CIs for all estimates. Both Cohen’s d and Hedges’ g statistics are estimated by calling the cohen.d() function within the Efficient Effect Size Computation or “effsize” package (Torchiano, 2020) in R software. Results agree that the RSI for the keyword “colonoscopy” is the best indicator for cancer prevalence, showing a 95% CI for the Pearson correlation coefficient of (0.9324; 0.9992), and a goodness-of-fit (R-squared) of 98.56%, indicating that it can explain most of the variability in CRC prevalence at the world level during the most recent years of available data. In turn, the RSI for the keyword “colorectal cancer” can also be a reliable indicator for cancer prevalence, although not a substitute, contributing to explain close to 73% and showing lower correlation with the CRC prevalence. Moreover, although the effect size magnitude is large for both relationships according to the thresholds provided in Cohen (1992), it is nonetheless significantly smaller for the RSI corresponding to the keyword “coloscopy” than for the RSI sourced for the keyword “colorectal cancer”. Consequently, estimations, albeit based on a small sample, agree with previous findings, indicating that online information-seeking for the keyword “colonoscopy” is the best proxy for the prevalence of CRC at the world level.

Forecasting CRC prevalence from Google search interest

The results of the nonlinearity tests performed on the RSI time series are reported in Table 5. Results agree that the web-query index is nonlinear, as both tests strongly reject the null hypothesis.

Consequently, accounting for the proven existence of data nonlinearity, this study proceeds to estimate a univariate neural network forecasting models that can properly handle this characteristic, i.e., FNNAR. For comparative purposes and to further increase the reliability of results, the widely used ETS and ARIMA models are also estimated, and their respective forecasting capability assessed. Results enclosed in Table 6 reveal that the ANN model, NNAR, is acknowledged as best performing in terms of out-of-sample forecasting accuracy by all scale- and scale-free metrics. It should also be mentioned that the estimated DM test further confirmed the superiority of the neural network model when compared to the second ranked forecasting model, ARIMA.

To assure the robustness of current findings, we first mitigate the sample instability that characterizes GT-sourced RSI series by following the technique proposed by Medeiros & Pires (2021), which solves the sampling instability of RSIs by sourcing and averaging across multiple samples with the same specifications. The aforementioned study shows that this approach improves both model selection and forecast accuracy. Here, we source and average across five different RSI samples constructed by indicating the same keyword, category, and geography. Results, reported in Table 7, again indicate the autoregressive neural network model as best performing in terms of out-of-sample forecasting accuracy and also reveal that employing the consolidated sample does mitigate forecast errors over the testing set for both FNNAR and ETS, offering some support for the findings of Medeiros & Pires (2021).

Moreover, we also perform a secondary robustness check by repartitioning the original RSI sample through the implementation of a 174-55 data splitting rule, as per Box 2.

Table 8 contains the accuracy measures estimated over the newly defined test-set that were issued by the three competing models, now trained over the first 174 observations in the sample. The nonlinear FNNAR has again been most able to capture variations in data and thus to issue superior out-of-sample forecasts for the RSI.

In light of the consistent superior performance over alternative test sets, the autoregressive NN model is subsequently reestimated over the entire series spanning January 2004 to December 2022 to further issue the expected RSI for CRC over the forecasting horizon of 24 months corresponding to the January 2023-December 2024 interval. Projections (visually represented in Fig. 10) show a continuation of population web-search interest in CRC, with an overall growth rate of 16.42 percent relative to the year 2019 (i.e., the point corresponding to the most recent observation for registered CRC prevalence rates at the world level), which in absolute terms would indicate a total burden of disease of 12,361,225 cases, or an additional number of 903,598 CRC “survivors” by the end of 2024. Of note, current projections reinforce the estimations of the Global Cancer Observatory (https://gco.iarc.fr/), reported by Xi & Xu (2021) and Morgan et al. (2023), among others, which expect a surge in colorectal cancer incidence and consequently report an expected count of 3.2 million new CRC cases by 2040, or an increase of 63% relative to 2020 levels. Furthermore, current projections also agree with the estimations of the American Cancer Society (2022a) cited by the National Cancer Institute (2022), which expect the number of US all-cancers survivors to increase by 24.4%, thus reaching 22.5 million, by 2032.

As a final step, we attempt to validate the capability of the proposed framework by employing it for the task of nowcasting and forecasting the level of worldwide CRC prevalence and compare its projections with the last year of available statistics at the world level (i.e., 2019) provided by the Global Burden of Disease (GBD) database of the Institute for Health Metrics and Evaluation (IHME) (https://www.healthdata.org/gbd).

Figure 11 attests that the linear model constructed within the step 1 of the framework manages to accurately nowcast the level of global CRC prevalence for 2019, indicating a global burden of 11,452,157 compared to 11,457,627 reported by the GBD database of IHME (i.e., 5,470 residual cases at the world level for year 2019).

For forecasting purposes, we subset the first 192 observations in the sample (spanning January 2004–December 2019) and employ an alternative data partitioning strategy (i.e., 180–12). Thus, the FNNAR is trained over the first 180 observations (training window spanning from December 2018) and used to make predictions for the relative search index for “coloscopy” through December 2019. Figure 12 shows the graphical representation of the FNNAR predictions over January to December of 2019, indicating that the model has been able to detect trends fairly well, although it has slightly underestimated the relative search interest over the testing set, issuing an average forecast for RSI of 71.5 over 2019 relative to the actual average RSI level of 78.33. Then, projections issued by FNNAR for the end of 2019 are fed into the linear model constructed during the first step of the framework to finally estimate the global CRC prevalence level corresponding to the year 2019, which is then compared to statistics published by IHME. Results show that the proposed framework has been able to capture the global CRC prevalence level for 2019, indicating an interval of (10, 114, 105–11, 812, 163) for the total number of cases globally (i.e., 10,963,134-point forecast), whereas the GBD database of IHME reports a global burden for colorectal cancer of 11,457,627 for 2019. However, to obtain a clearer picture of the framework's ability to project the global burden of the disease, additional validation is required; this will be accomplished with the release of new official statistics corresponding to more recent years.

Limitations

Research that employs Google Trends data, as is the case with the current study, is based on the central premise that people turn to Google to find subject-related information when they need it, and, due to its visibility, this information demand can subsequently be used as a reliable predictor in a variety of settings (Bleher & Dimpfl, 2022). Moreover, compared to official statistics, GT data are available in real-time and cover a wide sample of countries, regions, and even some large cities, offering valuable insights especially when other data are lacking or are only available with important time lags (Eichenauer et al., 2022), as is the case with cancer prevalence data that makes the subject of the current research. Consequently, Google Trends data have proven particularly helpful in shedding light on population health behaviors (Arora, McKee & Stuckler, 2019) and hold the potential to improve the capabilities of public health surveillance systems (Althouse et al., 2015; Tkachenko et al., 2017). However, there are also some non-trivial limitations to the GT data that should be properly acknowledged (Tudor, 2022a) and, where possible, overcome.

First, Google Trends RSI is aggregated from a sample of search queries (Narita & Yin, 2018), thus suffering from sample bias (Medeiros & Pires, 2021) which, together with sample instability that arises from its construction (Eichenauer et al., 2022), can have serious consequences for the reliability of research findings that use GT data. In fact, this makes any study employing GT data inherently irreproducible (Rovetta, 2021). To mitigate, albeit not overcome, these issues, we construct a consistent sample by additionally employing a robust sampling procedure that involves sourcing multiple samples at various intervals and averaging across to construct a “consolidated” RSI. Moreover, it should also be mentioned that the aforementioned weaknesses related to GT data are prominent for less searched topics or keywords (Medeiros & Pires, 2021), which was not the case for the main keyword of interest in this study (i.e., “colonoscopy”) over the analyzed period. Moreover, previous research has shown that “worldwide” RSIs, as is the case in the current study, are the most reliable GT samples (Rovetta, 2021).

Second, Google Trends data is susceptible to being influenced by external events, such as news events or media coverage (Sato et al., 2021; Satpathy, Kumar & Prasad, 2023). Thus, given that the period under consideration in the current study contains the COVID-19 pandemic, it should be noted that the media coverage may have influenced the search volume in pandemic-related terms during the analyzed period (Rovetta & Castaldo, 2021), and, given the construction of the RSI as explained in the method section, this in turn may have direct consequences for the relative interest in unrelated terms, such as the one used in this research (i.e., “colonoscopy”). Nonetheless, it should also be mentioned that even allowing for the possibility that media coverage may affect population web searches, GT data still offers a way to measure web interest in a particular subject more effectively than any other approaches previously employed (Rovetta, 2021). Furthermore, as the COVID-19 pandemic significantly affected CRC screening rates, with direct consequences for the accuracy of CRC prevalence data and their usefulness for cancer research, we align with previous studies and argue that digital traces are the best available tools for assessing disease prevalence.