Quantifying the under-reporting of genital warts cases

Genital warts are a common and highly contagious sexually transmitted disease. They have a large economic burden and affect several aspects of quality of life. Incidence data underestimate the real occurrence of genital warts because this infection is often under-reported, mostly due to their specific characteristics such as the asymptomatic course. Genital warts cases for the analysis were obtained from the catalan public health system database (SIDIAP) for the period 2009-2016, covering 74\% of the Catalan population. People under 15 and over 94 years old were excluded from the analysis as the incidence of genital warts in this population is negligible. This work introduces a time series model based on a mixture of two distributions, capable of detecting the presence of under-reporting in the data. In order to identify potential differences in the magnitude of the under-reporting issue depending on sex and age, these covariates were included in the model. This work shows that only about 80\% in average of genital warts incidence in Catalunya in the period 2009-2016 was registered, although the frequency of under-reporting has been decreasing over the study period. It can also be seen that the under-reported issue has a deeper impact on women over 30 years old. The registered incidence in the Catalan public health system is underestimating the real burden in almost 10,000 cases in Catalunya, around 23\% of the registered cases. The total annual cost in Catalunya is underestimated in at least about 10 million Euros respect the 54 million Euros annually devoted to genital warts in Catalunya, representing 0.4\% of the total budget of the public health system.

from the analysis as the incidence of genital warts in this population is negligible. This work introduces a time series model based on a mixture of two distributions, capable of detecting the presence of under-reporting in the data. In order to identify potential differences in the magnitude of the under-reporting issue depending on sex and age, these covariates were included in the model. This work shows that only about 80% in average of genital warts incidence in Catalunya in the period 2009-2016 was registered, although the frequency of under-reporting has been decreasing over the study period. It can also be seen that the under-reported issue has a deeper impact on women over 30 years old. The registered incidence in the Catalan public health system is underestimating the real burden in almost 10,000 cases in Catalunya, around 23% of the registered cases. The total annual cost in Catalunya is underestimated in at least about 10 million Euros respect the 54 million Euros annually devoted to genital warts in Catalunya, representing 0.4% of the total budget of the public health system.

Introduction
Health information systems are essential to ensure the safety and quality of health care and improve adherence to clinical practice guidelines, but they are also a very powerful tool concerning resources management and control, decision making, and effective and efficient planning of prevention and control interventions [1,2]. However, the incompleteness and inaccuracy of the information is common in this type of registries and can lead to problems at a clinical level, but also at a population level such as the underestimation of some diseases. In Catalunya (Spain), the Information System for Research in Primary Care (SIDIAP) was launched in 2010 with the integration of data from the clinical work station of primary care (ECAP) of the Catalan Health Institute (ICS), which started in 1998, and other complementary sources [3]. The ICS is the main provider of health services in Catalunya and manages 283 out of 370 Primary Care Teams with a catchment of 5,564,292 people, approximately 74% of the Catalan population 1 . Nevertheless, it is reasonable to assume that the incidence of genital warts (GW) will be very similar among the Catalan population not covered by ICS. In the particular case of sexually transmitted diseases, it is even more important to have reliable information due to their remarkable morbidity, and therefore, the importance of controlling trends over time and priority setting (see [4] for a comprehensive discussion focused on developing countries). GW are a common and highly contagious sexually transmitted disease in Catalunya (in 2016 the incidence was about 107 cases per 100,000 women and 139 cases per 100,000 men [5]) caused by a subset of HPV types, with the most common being genotypes 6 and 11. They are usually benign, or non-cancerous, skin growths that develop on the genital area. However, they have an important negative impact on the health service and the individual, in addition to have a large economic burden and affect several aspects of quality of life [6,7,8]. It is well known that incidence data underestimate, to some degree, the real occurrence of genital warts because this infection is often underreported, mostly due to their specific characteristics such as the asymptomatic course of the disease [9]. Further, the SIDIAP database only includes data from the public healthcare sector and around 28% of the general population in Catalunya have a double health insurance coverage, public and private [10], so this fact can also explain why GW incidence rates are underestimated, although this source of under-reporting cannot be detected by the proposed model as we only have data from the public health system. There has been a growing interest in the past recent years to deal with data that are only partially registered or under-reported in the biomedical literature [11,12,13,14,15,16]. Most of these previous works deal with discrete-valued time series, whereas this paper is focused on the incidence of a disease, which should be treated as a continuous-valued time series. Therefore, the aim of this work is to quantify the under-reporting of genital warts cases in Catalunya and the reconstruction of the actual incidence in the period 2009-2016 on the basis of the mixture model described in Section 2.

Population and incidence estimation
The study population included all residents in Catalunya assigned to an ICS primary care center (74% of the Catalan population). Monthly GW incident cases for the analysis were obtained from the SIDIAP database for the period 2009-2016. Episodes of GW were classified as incident if they were preceded by at least 12-month period without any episode. People under 15 and over 94 years old were excluded from the analysis as the incidence of GW in this population is negligible (averages of 0.24 cases and 0.22 x 100,000 individuals over the period of study respectively).

Model
Consider X t the series of real GW incidence, where t = 1, 2, . . . is the time, following a normal distribution with mean µ and variance σ 2 . In our setting, this process cannot be directly observed, and all we can see is a part of it, expressed as The series Y t represents the registered values corresponding to GW incidence in the part of Catalunya covered by ICS. According to Eq. 1, the registered observations series Y t is a mixture of two normally distributed random variables Y t = (1 − ω t )·Y 1t + w t ·Y 2t , where Y 1t coincides with the unobserved process and Y 2t is a normal random variable with mean µ and variance σ 2 . The parameter ω t is modeled as log(ω t ) = α 0 + α 1 ·t and can be interpreted as the frequency of under-reporting at a time t, while q can be interpreted as the intensity of such under-reporting, taking a value between 0 and 1. When q = 0 the observed incidence is Y t = 0 and when q = 1 there is no under-reporting. A value of ω t equal to 0 indicates that the observed value at time t is not under-reported, and a value of ω t equal to 1 means that under-reporting is for sure happening. In order to detect potential differences in GW incidence depending on sex (men and women) and age (16-29 and 30-94), these covariates were included in the model, so the mean of the observed process Y 1t was modeled as µ 1t = β 0 + β 1 ·t + β 2 · a + β 3 · s + β 4 · a * s (where a is the age, s is the sex and a * s is the interaction between age and sex). Similarly, the average of the second component can be recovered as µ 2t = q · (β 0 + β 1 · t + β 2 · a + β 3 · s + β 4 · a * s). After fitting the previous model and performing residuals examination, a seasonal behavior with period 3 months was observed. Hence the model was updated by including the following trigonometric function to reflect this periodic behavior: on the terms µ 1t and µ 2t .
Other similar models were considered and the best fitting one, according to the validation process described in Section 2.3, was chosen. In particular, as coefficients β 1 and β 6 are not significant, models without linear trend and with only one periodicity term were considered but the resulting validations were not satisfactory. The estimates and their associated standard errors were obtained by maximizing the likelihood function using the nlm procedure in R [17]. The R package mixtools[18] was used to obtain proper initial values for the maximization algorithm. All the data and code used are available as supplementary material. If the main focus is not on quantifying the under-reporting issue, an alternative approach to analyze these data might be a hierarchical generalized linear model with random effects [19], implemented in the R package HGLMM [20]. By means of this methodology the most likely unobserved real GW incidence process is reconstructed using the components of the estimated mixture, provided by the output of the mixtools procedure.

Validation
The model has been validated by analyzing its residuals. Figure 1 shows that they behave like white noise as expected and that there are no significant auto-correlations that should be accounted for. The residuals r t have been estimated aŝ r t = Y t − ω t ·q · β 0 +β 1 · t +β 2 · a +β 3 · s +β 4 · a * s + (1 −ω t ) · β 0 +β 1 · t +β 2 · a +β 3 · s +β 4 · a * s where Y t is the total observed GW incidence at time t, and the letters with hat indicate the estimated parameters.

Results
Our analysis estimates that, globally, only around 80% of actual GW incidence was registered in the SIDIAP database in the period 2009-2016. For women over 30 years  old, the monthly average registered incidence is 3.9 cases per 100,000 women, while the estimated monthly incidence is 4.9 cases per 100,000 women, 24.9% higher. On males over 30 years old, the registered series has a monthly average of 5.9 cases per 100,000 men for 7.1 cases per 100,000 men on the reconstructed series, 21.8% higher. Regarding males under 30 years old, the reconstructed series is 13.3% higher (monthly averages of 18.4 and 20.8 cases per 100,000 men for the registered and reconstructed processes respectively). For women under 30 years old, the monthly average registered incidence of GW in Catalunya is 19.0 per 100,000 women, while the reconstructed hidden process has an average of 23.0 cases per 100,000 women, about 21.0% larger. This information is summarized in Table 1 and described in more detail in the supplementary material (Table S1). Table 2 shows the estimated effect of the age and sex over the under-reporting issue. In particular, it can be seen that the GW incidence is higher among younger populations and men. It can also be noticed that a significant interaction between sex and age group is found, which can be interpreted as a distinguishable impact of sex on GW incidence depending on the age group.  Figure 2 shows the registered (solid black line) and reconstructed unobserved (dashed red line) processes for each of the considered sub-populations. Although this figure shows increasing trends for all series, they are not well explained by coefficient β 1 , which is not significantly different from zero. Increasing trends are mainly explained by the significant coefficient α 1 , which leads to a decreasing frequency of under-reporting ω t . The under-reporting frequency is about 95% in 2009 (ω 1 ) and around 21% in 2016 (ω t ). This is measured by parameter α 1 in Eq 1, and should not be confused to overall under-reporting of the data, as its intensity (measured by parameter q in the model) also plays a crucial role. For instance, all observations in a certain period of time could be slightly under-reported (ω = 1, q near to 1), resulting in small differences between registered and estimated values or just a few observations might be under-reported (ω near to zero) but with a high intensity (q near to zero), potentially resulting on large differences between registered and estimated values. Table 3 shows the total number of GW cases registered in the SIDIAP in the period of study, the reconstructed values according to these registered cases and the projection over the whole Catalan population, assuming that the incidence on the area outside ICS coverage is the same.

Discussion
The results of this work show that in relative terms, the under-reporting issue has a deeper impact on people over 30 years old (where GW incidence is lower), especially among women. Nonetheless, the relative difference between registered and estimated annual averages range between 13.3% and 24.9%. It is also remarkable that the quality of SIDIAP register regarding GW in Catalunya has been significantly improving during the study period, as the frequency of under-reported observations has been decreasing over time. Facing under-reported information from public health registers is very common in many situations, especially regarding potentially asymptomatic diseases like GW. The proposed methodology considers the potential under-reporting in continuous time series data in a very flexible way, estimating its frequency and intensity, and it is general enough to be appropriate in a wide range of real situations in the public health context. Additionally, the most likely non-observed process can be reconstructed on the basis of estimated posterior probabilities. Moreover, the GW data show that these models can deal with time-dependent under-reporting parameters, seasonal behavior, trends and also incorporate the effect of other factors by including covariates. The described methodology opens a wide field for future research lines. In particular, if temporal correlations are found in the data, an appropriate model should take this structure into account. One of the potential limitations of this study is that the database used included data from the public healthcare setting and not from the private sector. In Catalunya, it is estimated that 33% of women and 25% of men aged 15 to 44 years have a double health insurance coverage (i.e. the public health insurance and a private insurance plan) [10], so the rates estimated in our study are likely still underestimating the real incidence of GW. One of its strengths is that the same methodology (possibly with minor model modifications) could be used to analyze the frequency and intensity of potential under-reporting issues for any condition or setting in the absence of temporal dependence among the observations. The GW incidence registered in SIDIAP is underestimating the real burden in almost 10,000 cases in Catalunya, around 23% of the registered cases. The annual per person cost of GW was around 1000 Euros [8], so the potential total annual cost is underestimated in at least about 10 million Euros respect the 54 million Euros devoted to GW in Catalunya annually, representing 0.4% of the total budget of the Catalan Government intended for health 2 , although about 2.8 million Euros would correspond to private insurances. It is, therefore, clear that knowing the true burden of GW at the general population level is important for health policy makers, especially after the introduction of prophylactic vaccines against HPV in many countries, as it plays a crucial role in developing and evaluating prevention strategies [21,22].