key: cord-0667324-x51cdraa authors: Morina, David; Fern'andez-Fontelo, Amanda; Cabana, Alejandra; Puig, Pedro title: New statistical model for misreported data with application to current public health challenges date: 2020-03-20 journal: nan DOI: nan sha: 898dc3e60242acbc183d3d976ef549fae0fe2b64 doc_id: 667324 cord_uid: x51cdraa The main goal of this work is to present a new model able to deal with potentially misreported continuous time series. The proposed model is able to handle the autocorrelation structure in continuous time series data, which might be partially or totally underreported or overreported. Its performance is illustrated through a comprehensive simulation study considering several autocorrelation structures and two real data applications on human papillomavirus incidence in Girona (Catalunya, Spain) and COVID-19 incidence in the Chinese region of Heilongjiang. There has been a growing interest in the past years to deal with data that is only partially registered or underreported in the time series literature. This phenomenon is very common in many fields, and has been previously explored by different approaches in epidemiology, social and biomedical research among many other contexts [6, 3, 23, 1, 27] . The sources and underlying mechanisms that cause the underreporting might differ depending on the particular data. Some authors consider a situation where the registry is updated with time and therefore the underreporting issue is mitigated [13] . That leads to temporary underreporting while this work is focused on permanent underreporting, where the registered data are never updated in order to become more accurate. From the methodological point of view, several alternatives have been explored, from Markov chain Monte-Carlo based methods [27] to recent discrete time series approaches [9, 10] . Several attempts to estimate the degree of underreporting in different contexts have been done [11] , although there is a lack of models incorporating continuous time series structures and handling underreporting. One of the fields where the interest in addressing the underreporting issues is higher is the epidemiology of infectious diseases. In the last few years, many approaches to deal with underreported data have been suggested with a growing level of sophistication from the usage of multiplication factors [25] to several Markov-based models [4, 20] or even spatio-temporal modelling [26] . Even a new R [22] package able to fitting endemic-epidemic models based on approximative maximum likelihood to underreported count data has been recently published [17] . This work presents two examples where such phenomenon appears. Human papillomavirus (HPV) is one of the most prevalent sexually transmitted infections. It is so common that nearly all sexually active people have it at some point in their lives, according to the information provided by the United States' Centers for Disease Control and Prevention (CDC) [8] . Generally, the infection disappears on its own without inducing any health problem, but in some cases it can produce an abnormal growth of cells on the surface of the cervix that could potentially lead to cervical cancer. HPV infection is also related to other cancers (vulva, vagina, penis, anus, . . . ) and other diseases like genital warts (GW). The fact that most cases of HPV infection are asymptomatic causes that public health registries might be potentially underestimating its incidence. The underreporting phenomenon in HPV data from the discrete time series point of view has been recently studied [9] . There is an enormous global concern around 2019-novel coronavirus (SARS-CoV-2) infection in the last few months, leading the World Health Organization (WHO) to declare public health emergency [24] . As the symptoms of this infection can be easily confused with those of similar diseases like Middle East Respiratory Syndrome Coronavirus (MERS-CoV) or Severe Acute Respiratory Syndrome Coronavirus (SARS-CoV), its incidence has been notably underreported, especially at the beginning of the outbreak in Wuhan (Hubei province, China) by December 2019. Consider an unobservable process with an AutoRegressive Moving Average (ARM A(p, r)) structure defined by where ǫt is a Gaussian white noise process with ǫt ∼ N (µǫ, σ 2 ǫ ). The ARMA processes belong to the family of so called linear processes. Their importance relies on the fact that any stationary nondeterministic process can be written as a sum of a linear process and a deterministic component [7] . These models are very well known, have been used in many applications since their introduction in the early 1950's and are general and flexible enough to be useful in a wide range of different contexts. Most used statistical software packages include functions that allow straightforward fitting of this family of models, so it seems a natural choice in the present work. In our setting, this process Xt cannot be directly observed, and all we can see is a part of it, expressed as The interpretation of the parameters in Eq. (2) is straightforward: q is the overall intensity of misreporting (if 0 < q < 1 the observed process Yt would be underreported while if q > 1 the observed process Yt would be overreported). The parameter ω can be interpreted as the overall frequency of misreporting (proportion of misreported observations). Consider that the unobserved process Xt follows an ARM A(p, r) model as defined in Eq. (1). As can be seen in Appendix 1 (Supplementary Material), the observed process has mean E(Yt) = The autocorrelation function of the observed process can be written in terms of the features of the hidden process Xt as where ρX is the autocorrelation function of the unobserved process Xt. A situation of particular interest is the case ω = 1, meaning that all the observations might be underreported and that a simpler model for Yt excluding the parameter ω might be suitable In this case, however, the observed process Yt would be a non-identifiable ARM A(p, r) model as the parameter q cannot be estimated on the basis of the methodology described in the following section. The likelihood function of the observed process Yt is not easily computable but the parameters of the model can be estimated by means of an iterative algorithm based on its marginal distribution, using the R packages mixtools [5] and forecast [16, 15] . The main steps are described in detail below: (i) Following Eq. (2), the observed process Yt can be written as Yt = (1−Zt)·Xt +q ·Zt ·Xt, where Zt is an indicator of the underreported observations, following a Bernoulli distribution with probability of success ω (Zt ∼ Bern(ω)). The marginal distribution of Yt is a mixture of two normal random variables N µ, σ 2 and N q · µ, q 2 · σ 2 . This fact can be used to obtain initial estimates for q and ω. Using the EM algorithm (specifically on the E-step), the posterior probabilities (conditional on the data and the obtained estimates) can be computed. This can be done using, for instance, the R package mixtools. (ii) Using the indicatorẐt obtained in the previous step, the series is divided in two: One including the underreported observations (treating the non-underreported values as missing data) and another with the non underreported observations (treating the underreported values as missing data). An ARM A model is fitted to each of these two series and a newq is obtained by dividing the fitted means. (iii) A mixture of two normals is fitted to the observed series Yt with mean and standard deviation fixed to the corresponding values obtained from the previous step, and a new ω is estimated. (iv) Steps (ii) and (iii) are repeated until the quadratic distance between (v) Once the parameter estimates are stable according to the previous criterion, the underlying process Xt is reconstructed asXt = (1 − Zt)·Yt + 1 q ·Ẑt ·Yt, and an ARM A model is fitted to the reconstructed process to obtainαj , j = 1, . . . , p,θ k , k = 1, . . . , r andσǫ 2 . To account for potential trends or seasonal behaviour, covariates can be included in the described estimation process. Additionally, a parametric bootstrap procedure with 500 replicates is used to estimate standard errors and build confidence intervals based on the percentiles of the distribution of the estimates. The R codes used to estimate the parameters and build the confidence intervals are available from the authors upon request. The results of the proposed methodology over a comprehensive simulation study and an application on two real data sets are shown in this Section. A thorough simulation study has been conducted to ensure that the model behaves as expected, including AR(p), M A(r) and ARM A(p, r) for 1 ≤ p, r ≤ 3 structures for the hidden process Xt with values for the parameters α, θ, q and ω ranging from 0.1 to 0.9 for each parameter (some combinations of parameters have been omitted for p > 1 or r > 1 to ensure stationarity). For ARM A(p, r) structures with p > 1 or r > 1 the parameters covered the same range (0.1 to 0.9) but with a difference of 0.2 instead of 0.1 for computational feasibility. Only average absolute bias, interval coverage and 95% confidence interval corresponding to p = r = 1 are shown in Table 1 , as higher order models behave in a very similar manner (see Supplementary Material for details) . These values are averaged over all combinations of parameters. Additionally, standard AR(1), M A(1) and ARM A(1, 1) models were fitted to the same simulated series without accounting for their underreporting structure. For each autocorrelation structure and parameters combination, a random sample of size n = 1000 has been generated using the function arima.sim from R package forecast [16, 15] . Different sample sizes (n = 50, 100, 500) have also been considered to study the impact of sample size on accuracy and the results are reported in the Supplementary Material. The performance of the proposed methodology is summarised in Tables S1 to S4 for n = 50, 100, 500 and 1000 respectively. Average absolute bias is similar regardless of the sample size, while average interval lengths (AIL) are higher and interval coverages are poorer (around 75% for n = 50) for lower sample sizes as could be expected. Several bootstrap sizes (b = 20, 50, 100, 500) were also considered and the difference between them were negligible, so only results corresponding to b = 500 bootstrap replicates are reported. It is clear from Table 1 that ignoring the underreported nature of data (labeled as Standard models in the table) leads to highly biased estimates with extremely low coverage rates, even with larger average interval lengths. This is especially relevant when the intensity or frequency of underreported observations is high. The series of weekly cases of HPV infection in Girona in the period 2010-2014 was previously analysed as a discrete IN AR(1) hidden Markov process [9] . In a similar way, we aim to analyse the corresponding series of incidence, and an AR process of order 1 seems to be adequate (see Figure 1) . Additionally, the AR(1) structure has the lowest AIC when compared to similar alternative models like AR (2) In particular, in this case we can write log(ρY (k)) = log(c)+k ·log(α), so a statistically significant intercept of this linear regression model (estimating the parameters by ordinary least squares method) could be interpreted as an evidence of underreporting, as in this case (p − value = 0.0016). It is clear from Figure 1 that the estimated regression line does not cross the origin, so the behavior of the observed process is consistent with an underlying underreported AR(1) process. By means of the estimation method described in Section Estimation, it can be seen that the estimated model for the hidden process is Xt = 0.109 · Xt−1 + ǫt, being the observed process Yt, Yt = Xt with probability 0.240 0.242 · Xt with probability 0.760 (5) The estimated parameters are reported in Table 2 . These results are highly consistent with those previously reported in the literature for the number of HPV cases obtained through a discrete time series approach [9] and can be interpreted in a very straightforward way. Moreover, this new methodology can be used to model the incidence of the disease instead of the number of cases, accounting for potential changes in the underlying population. The estimated intensity of underreporting isq = 0.242, with 95% confidence interval (0.176, 0.308). The registered and estimated evolution of HPV incidence within the study period (2010-2014) can be seen in Figure 2 . These results indicate that only 34% of the HPV incidence in the considered period of time was actually recorded. Taking into account that public health cervical cancer prevention strategies are often designed on the basis of simulation models which are calibrated to registered HPV data [21] , it is clear that providing decision makers with accurate data on HPV incidence is key to ensure optimal allocation of scarce public health funds. The betacoronavirus SARS-CoV-2 has been identified as the causative agent of an unprecedented world-wide outbreak of pneumonia starting in December 2019 in the city of Wuhan (China) [24] , named as COVID-19. Considering that many cases run without developing symptoms beyond those of MERS-CoV, SARS-CoV or pneumonia due to other causes, it is reasonable to assume that the incidence of this disease has been underregistered, especially at the beginning of the outbreak [28] . This work focuses on the COVID-19 incidence registered in Heilongjiang province (north-eastern China) in the period (2020/01/22-2020/02/26), and it can be seen in Figure 3 that the registered data (black color) reflect only a fraction of the actual incidence (red color). A disease with a similar behavior (MERS-CoV) has been modeled in a previous work as an ARM A(3, 1) [2] , so we evaluated the performance of this model and similar ones. Probably due to the shortness of the available data this autoregressive structure was not observed and in our case the best performing model was an M A(1) (AIC of -151.04 against -148.49 for the ARM A(3, 1)), consistently with the residuals profile shown in Figure 4 , obtained from fitting an M A(1) model to the most likely process Xt reconstructed following step (v) in Section Estimation. By means of the estimation method described in Section Estimation, it can be seen that the estimated model for the hidden process is Xt = 0.481 · ǫt−1 + ǫt, being the observed process Yt, Yt = Xt with probability 0.483 0.194 · Xt with probability 0.517 (6) The estimated parameters are reported in Table 3 . In biomedical and epidemiological research, the usage of disease registries in order to analyse the impact and incidence of health issues is very common. However, the accuracy and data quality of such registries is in many cases at least doubtful. This is the case, for instance, for rare diseases [18] or health issues that clear asymptomatically in most cases like HPV infection. In the case of HPV incidence in Girona in the period 2010-2014, the registered weekly average is 0.17 cases per 100,000 individuals, while the reconstructed process has a weekly average of 0.50 cases per 100,000 individuals, showing that only 34% of the real incidence is recorded by the public health system. It must be considered that HPV infection is related to subsequent complications such as cervical cancer in some cases and that public health cervical cancer prevention strategies are often designed on the basis of simulation models which are calibrated to registered HPV data [21] and therefore the optimal allocation of scarce public health resources cannot be ensured if the under-reporting issue is not accounted for. This result is very consistent with that of [9] , where the authors claim that only 38% of the HPV cases were registered in the same area and period of time. The Heilongjiang region COVID-19 data reveal that in average about 66% of the actual incidence in the period 2020/01/22-2020/02/26 was reported. The unavailable data estimated by the proposed methodology are crucial to provide public health decision-makers with reliable information, which can also be used to improve the accuracy of dynamic models aimed to estimate the spread of the disease [28] . In China and almost globally afterwards, different non-pharmaceutical interventions were undertaken in order to minimise the impact of the disease on the general population and especially over the health systems, which were put to the limit of their capacity by the pandemic. In this context, one of the main challenges in predicting the evolution of the disease or evaluating the impact of these strategies is to use data as accurate as possible, taking into account that many COVID-19 cases are asymptomatic or with mild symptoms and a generalized shortage of testing kits [14] , and therefore knowing that the registered number of affected individuals might be severely underestimated. The concerns around accuracy of registered data have recently led to the publication of recommendations to improve data collection to ensure accuracy of registries (see for instance [19, 12] ). Nonetheless, these recommendations are very recent and may be difficult for the public health services of many countries to fully implement them, due to operational or structural issues. The proposed methodology is able to deal with underreported (or overreported) data in a very natural and straightforward way, estimating its intensity and frequency on a continuous time series, and allowing to reconstruct the most likely unobserved process. It is also flexible enough to handle covariates straightforwardly, and therefore it is simple to introduce trends or seasonality if necessary, so it can be useful in many contexts, where these issues might arise. The simulation study shows that the proposed methodology behaves as expected and that the parameters used in the simulations, under different autocorrelation structures, are properly recovered, regardless of the intensity and frequency of the underreporting issues. Work-related skin diseases in Norway may be underreported: data from Temporal dynamics of Middle East respiratory syndrome coronavirus in the Arabian Peninsula Reporting of foodborne illness by U.S. consumers and healthcare professionals On the estimation of the reproduction number based on misreported epidemic data mixtools : An R Package for Analyzing Finite Mixture Models Estimating the under-reporting of norovirus illness in Germany utilizing enhanced awareness of diarrhoea during a large outbreak of Shiga toxinproducing E. coli O104: H4 in 2011 -a time series analysis Time Series: Theory and Methods. chapter Stationary CDC grand rounds: Reducing the burden of HPVassociated cancer and disease Under-reported data analysis with INAR-hidden Markov chains Untangling serially dependent underreported count data for gender-based violence Measuring underreporting and underascertainment in infectious disease datasets: a comparison of methods Towards a Core Set of Indicators for Data Quality of Registries Bayesian nowcasting during the STEC O104:H4 outbreak in Germany Rapid asymptomatic transmission of COVID-19 during the incubation period demonstrating strong infectivity in a cluster of youngsters aged 16-23 years outside Wuhan and characteristics of young patients with COVID-19: A prospective contact-tracing study and F Yasmeen. forecast: Forecasting functions for time series and linear models Automatic time series forecasting: The forecast package for R hhh4u: Fit an endemic-epidemic model to underreported data Data quality in rare diseases registries Hanns Lochmüller, and Domenica Taruscio. Recommendations for improving the quality of rare disease registries The parameter identification problem for SIR epidemic models: identifying unreported cases Impact of model calibration on cost-effectiveness analysis of cervical cancer prevention R: A Language and Environment for Statistical Computing How much work-related injury and illness is missed by the current national surveillance system? Ahmed Al-Jabir, Christos Iosifidis, and Riaz Agha. World Health Organization declares Global Emergency: A review of the 2019 Novel Coronavirus (COVID-19) Model selection and parameter estimation for dynamic epidemic models via iterated filtering: application to rotavirus in Germany A Hierarchical Framework for Correcting Under-Reporting in Count Data Markov Chain Monte Carlo analysis of underreported count data with an application to worker absenteeism Estimating the Unreported Number of Novel Coronavirus (2019-nCoV) Cases in China in the First Half of January 2020: A Data-Driven Modelling Analysis of the Early Outbreak This work was supported by grant COV20/00115 from Instituto de Salud Carlos III (Spanish Ministry of Health). This work was partially supported by grant RTI2018-096072-B-I00 from the Spanish Ministry of Science and Innovation.