Quantifying the under-reporting of uncorrelated longitudal data: the genital warts example

Abstract

Background

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.

Methods

Genital warts cases for the analysis were obtained from the Catalan public health system database (SIDIAP) for the period 2009-2016. 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.

Results

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 this issue has a deeper impact on women over 30 years old.

Conclusions

Although this study shows that the quality of the registered data has improved over the considered period of time, the Catalan public health system is underestimating genital warts real burden in almost 10,000 cases, around 23% of the registered cases. The total annual cost is underestimated in about 10 million Euros respect the 54 million Euros annually devoted to genital warts in Catalunya, representing 0.4% of the total budget.

Background

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 (http://ics.gencat.cat/es/lics/). 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–8]. A higher risk of CIN2+ lesions in women following a GW diagnose has been reported in a comprehensive recent study, even more than four years after the GW diagnose [9]. It is well known that incidence data underestimate, to some degree, the real occurrence of genital warts because this infection is often under-reported, mostly due to their specific characteristics such as the asymptomatic course of the disease [10]. This issue might be even more severe in specific vulnerable populations as imprisoned women [11]. 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, so this fact can also explain why GW incidence rates are underestimated [12], 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 [13–18]. 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 the next Section.

Methods

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 σ². In our setting, this process cannot be directly observed, and all we can see is a part of it, expressed as

$$ Y_{t}=\left\{ \begin{array}{ll} X_{t} & \text{ with probability}\ 1-\omega_{t}\\ q \cdot X_{t} & \text{ with probability}\ \omega_{t}\\ \end{array} \right. $$

(1)

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=Y_1t with probability (1−ω_t) and Y_t=Y_2t with probability ω_t, where Y_1t coincides with the unobserved process X_t and Y_2t is a normal random variable with mean q·μ and variance q²·σ². The parameter ω_t is modeled as logit(ω_t)=α₀+α₁·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, both taking values 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 μ_1,t=β₀+β₁·t+β₂·a+β₃·s+β₄·a∗s (where a is the age, s is the sex and a∗s is the interaction between age and sex). The average of the second component Y_2t can be recovered as μ_2,t=q·(β₀+β₁·t+β₂·a+β₃·s+β₄·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: $f(t)=\beta _{5} \cdot sin\left (\frac {2 \cdot \pi \cdot t}{3}\right)+\beta _{6} \cdot cos\left (\frac {2 \cdot \pi \cdot t}{3}\right)$ on the terms μ_1,t and μ_2,t. Other similar models were considered and the best fitting one according to the validation process described in the next Section was chosen. In particular, as coefficients β₁ and β₆ are not significant, models without linear trend and with only one periodicity term were considered but the resulting validations were not satisfactory. Therefore, the final expressions were $\mu _{1,t}=\beta _{0}+\beta _{1} \cdot t+\beta _{2} \cdot a+\beta _{3} \cdot s+\beta _{4} \cdot a*s+\beta _{5} \cdot sin\left (\frac {2 \cdot \pi \cdot t}{3}\right)+\beta _{6} \cdot cos\left (\frac {2 \cdot \pi \cdot t}{3}\right)$ and $\mu _{2,t}=q \cdot \left (\beta _{0}+\beta _{1} \cdot t+\beta _{2} \cdot a+\beta _{3} \cdot s+\beta _{4} \cdot a*s+\beta _{5} \cdot sin\left (\frac {2 \cdot \pi \cdot t}{3}\right)+ \beta _{6} \cdot cos\left (\frac {2 \cdot \pi \cdot t}{3}\right)\right)$. The estimates and their associated standard errors were obtained by maximizing the log-likelihood function described in Eq. (2) and from its Hessian matrix respectively, using the nlm procedure in R [19].

$$ \begin{aligned} l(Y,\theta) = \sum_{t=1}^{n} \log \left((1- \omega_{t}) \frac{1}{\sqrt{2 \pi} \sigma} e^{\frac{(y_{t}-\mu_{1, t})^{2}}{2 \sigma^{2}}} \right.\\ \left.+ \omega_{t} \frac{1}{\sqrt{2 \pi} q \sigma} e^{\frac{(y_{t}- \mu_{2, t})^{2}}{2 q^{2} \sigma^{2}}} \right), \end{aligned} $$

(2)

where Y=y₁,…,y_n is the observed series, $\theta =(\alpha _{0}, \alpha _{1}, \gamma, \beta _{0}, \ldots, \beta _{6}, \sigma), \omega _{t} = \frac {e^{\alpha _{0}+\alpha _{1} t}}{1+e^{\alpha _{0}+\alpha _{1} t}}, q = \frac {e^{\gamma }}{1+e^{\gamma }}$ and μ_1,t and μ_2,t are as defined before.

In order to get proper initial values for the maximization routine, an Expectation-Maximization (EM) algorithm for mixtures of linear regressions was used, through the R package mixtools [20]. The estimates provided by the EM algorithm could have been used directly, but although this methodology is widely used when dealing with mixtures of distributions, it is unable to produce standard errors directly [21], and this is an important drawback in our context and in many other situations. If the main focus was not on quantifying the under-reporting issue, an alternative approach to analyze these data might be a hierarchical generalized linear model with random effects [22], implemented in the R package HGLMM [23]. By means of this methodology the most likely unobserved real GW incidence process is reconstructed based on the classification (underreported or not underreported) given by the posterior probabilities for the observations, provided by the output of the mixtools procedure, and on the estimates of the parameters. All the R code used to fit the models and to obtain the reported results and figures is available as Supplementary material.

Validation

The goodness of fit of the proposed mixture approach can be assessed by means of the Akaike’s Information Criterion (AIC) compared to a single normal model. In this case, this measure favors the proposed model (AIC: 1717.9) in front of the single normal model (AIC: 1826.1). 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 as

$$ \begin{aligned} \hat{r_{t}}&=Y_{t}-\left(\hat{\omega_{t}} \cdot \hat{q} \cdot \left(\hat{\beta_{0}}+\hat{\beta_{1}} \cdot t+\hat{\beta_{2}} \cdot a+\hat{\beta_{3}} \cdot s +\hat{\beta_{4}} \cdot a*s+\hat{\beta_{5}} \cdot sin \left(\frac{2 \cdot \pi \cdot t}{3}\right)\right.\right.\\ &\quad \left.+\hat{\beta_{6}} \cdot cos \left(\frac{2 \cdot \pi \cdot t}{3}\right)\right) + (1-\hat{\omega_{t}}) \cdot \left(\hat{\beta_{0}} \right. +\hat{\beta_{1}} \cdot t+\hat{\beta_{2}} \cdot a+\hat{\beta_{3}} \cdot s+\hat{\beta_{4}} \cdot a*s\\ &\quad\left.\left.+ \hat{\beta_{5}} \cdot sin\left(\frac{2 \cdot \pi \cdot t}{3}\right)+\hat{\beta_{6}} \cdot cos\left(\frac{2 \cdot \pi \cdot t}{3}\right)\right)\right) \end{aligned} $$

(3)

Fig. 1

Auto-correlations and partial auto-correlations of the model residuals

where Y_t is the total observed GW incidence at time t, and the letters with a hat ($\hat {}$) represent the estimated parameters.

If we were dealing with counts (number of cases) instead of incidence, the underlying distribution might be a Poisson, although the monthly number of GW cases is large enough to be approximated by a normal distribution. Additionally, the assumption that the underlying distributions of the two processes are Gaussian seems reasonable considering the qqplot of the residuals shown in Fig. 2.

Fig. 2

Residuals qqplot

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 1 Registered and estimated GW monthly average incidence (number of cases x 100,000 individuals) in the period 2009-2016

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.

Table 2 Parameter estimates

Figure 3 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 β₁, which is not significantly different from zero. Increasing trends are mainly explained by the significant coefficient α₁, which leads to a decreasing frequency of under-reporting ω_t.

Fig. 3

Registered (solid black line) and estimated underlying series (dashed red line) for each of the considered sub-populations

The under-reporting frequency is about 95% in 2009 (ω₁) and around 21% in 2016 (ω₉₆). This is measured by parameter α₁ in model (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.

Table 3 Registered, estimated and projected number of GW cases in Catalunya

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.

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) [12], 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.

Conclusions

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, 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 [24, 25]. This work presents a methodology that 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, similarly to [13, 18].