Moriña et al. BMC Medical Research Methodology (2021) 21:6 https://doi.org/10.1186/s12874-020-01188-4 RESEARCH ARTICLE Open Access Quantifying the under-reporting of uncorrelated longitudal data: the genital warts example David Moriña1,2* , Amanda Fernández-Fontelo3 , Alejandra Cabaña4 , Pedro Puig4 , Laura Monfil5 , Maria Brotons5 and Mireia Diaz5 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. Keywords: Genital warts, Estimation, HPV, Under-reporting, Time series Background erful tool concerning resources management and con- Health information systems are essential to ensure the trol, decision making, and effective and efficient planning safety and quality of health care and improve adherence to of prevention and control interventions [1, 2]. However, clinical practice guidelines, but they are also a very pow- the incompleteness and inaccuracy of the information is common in this type of registries and can lead to prob- *Correspondence: dmorina@ub.edu lems at a clinical level, but also at a population level such 1 Department of Econometrics, Statistics and Applied Economics, as the underestimation of some diseases. In Catalunya Riskcenter-IREA, Universitat de Barcelona (UB), Avinguda Diagonal, 690, 08034 Barcelona, Spain (Spain), the Information System for Research in Primary 2 Centre de Recerca Matemàtica, Universitat Autònoma de Barcelona (UAB), Care (SIDIAP) was launched in 2010 with the integra- Edifici C, Campus de Bellaterra, 08193 Cerdanyola del Vallès, Spain tion of data from the clinical work station of primary Full list of author information is available at the end of the article care (ECAP) of the Catalan Health Institute (ICS), which © The Author(s). 2020 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Moriña et al. BMC Medical Research Methodology (2021) 21:6 Page 2 of 7 started in 1998, and other complementary sources [3]. The assigned to an ICS primary care center (74% of the Catalan ICS is the main provider of health services in Catalunya population). Monthly GW incident cases for the analysis and manages 283 out of 370 Primary Care Teams with a were obtained from the SIDIAP database for the period catchment of 5,564,292 people, approximately 74% of the 2009-2016. Episodes of GW were classified as incident if Catalan population (http://ics.gencat.cat/es/lics/). Never- they were preceded by at least 12-month period without theless, it is reasonable to assume that the incidence of any episode. People under 15 and over 94 years old were genital warts (GW) will be very similar among the Catalan excluded from the analysis as the incidence of GW in this population not covered by ICS. In the particular case population is negligible (averages of 0.24 cases and 0.22 x of sexually transmitted diseases, it is even more impor- 100,000 individuals over the period of study respectively). tant to have reliable information due to their remarkable morbidity, and therefore, the importance of controlling Model trends over time and priority setting (see [4] for a com- Consider Xt the series of real GW incidence, where t = prehensive discussion focused on developing countries). 1, 2, . . . is the time, following a normal distribution with GW are a common and highly contagious sexually trans- mean μ and variance σ 2 . In our setting, this process can- mitted disease in Catalunya (in 2016 the incidence was not be directly observed, and all we can see is a part of it, about 107 cases per 100,000 women and 139 cases per expressed as 100,000 men[5]) caused by a subset of HPV types, with the  Xt with probability 1 − ωt most common being genotypes 6 and 11. They are usu- Yt = (1) q · Xt with probability ωt ally benign, or non-cancerous, skin growths that develop on the genital area. However, they have an important The series Yt represents the registered values correspond- negative impact on the health service and the individ- ing to GW incidence in the part of Catalunya covered ual, in addition to have a large economic burden and by ICS. According to Eq. (1), the registered observations affect several aspects of quality of life [6–8]. A higher series Yt is a mixture of two normally distributed ran- risk of CIN2+ lesions in women following a GW diag- dom variables Yt = Y1t with probability (1 − ωt ) and nose has been reported in a comprehensive recent study, Yt = Y2t with probability ωt , where Y1t coincides with even more than four years after the GW diagnose [9]. It the unobserved process Xt and Y2t is a normal random is well known that incidence data underestimate, to some variable with mean q · μ and variance q2 · σ 2 . The degree, the real occurrence of genital warts because this parameter ωt is modeled as logit(ωt ) = α0 + α1 · t and infection is often under-reported, mostly due to their spe- can be interpreted as the frequency of under-reporting cific characteristics such as the asymptomatic course of at a time t, while q can be interpreted as the intensity the disease [10]. This issue might be even more severe of such under-reporting, both taking values between 0 in specific vulnerable populations as imprisoned women and 1. When q = 0 the observed incidence is Yt = 0 [11]. Further, the SIDIAP database only includes data and when q = 1 there is no under-reporting. A value of from the public healthcare sector and around 28% of ωt equal to 0 indicates that the observed value at time the general population in Catalunya have a double health t is not under-reported, and a value of ωt equal to 1 insurance coverage, public and private, so this fact can means that under-reporting is for sure happening. In also explain why GW incidence rates are underestimated order to detect potential differences in GW incidence [12], although this source of under-reporting cannot be depending on sex (men and women) and age (16-29 and detected by the proposed model as we only have data from 30-94), these covariates were included in the model, so the public health system. There has been a growing inter- the mean of the observed process Y1t was modeled as est in the past recent years to deal with data that are only μ1,t = β0 + β1 · t + β2 · a + β3 · s + β4 · a ∗ s (where a is the partially registered or under-reported in the biomedical age, s is the sex and a ∗ s is the interaction between age literature [13–18]. Most of these previous works deal with and sex). The average of the second component Y2t can be discrete-valued time series, whereas this paper is focused recovered as μ2,t = q · (β0 + β1 · t + β2 · a + β3 · s + β4 · a ∗ s). on the incidence of a disease, which should be treated as a After fitting the previous model and performing residuals continuous-valued time series. Therefore, the aim of this examination, a seasonal behavior with period 3 months work is to quantify the under-reporting of genital warts was observed. Hence the model was updated by includ- cases in Catalunya and the reconstruction of the actual ing the following trigonometric function   to reflect this  2·π·t incidence in the period 2009-2016 on the basis of the periodic behavior: f (t) = β5 · sin 2·π·t3 + β 6 · cos 3 mixture model described in the next Section. on the terms μ1,t and μ2,t . Other similar models were considered and the best fitting one according to the Methods validation process described in the next Section was Population and incidence estimation chosen. In particular, as coefficients β1 and β6 are not The study population included all residents in Catalunya significant, models without linear trend and with only Moriña et al. BMC Medical Research Methodology (2021) 21:6 Page 3 of 7 one periodicity term were considered but the resulting an alternative approach to analyze these data might be a validations were not satisfactory. Therefore, the final hierarchical generalized linear model with random effects expressions were μ1,t = β0 + β1 · t + β2 · a + β3 · s + [22], implemented in the R package HGLMM [23]. By β4 · a ∗ s + β5 · sin 2·π·t + β6 · cos 2·π·t and μ2,t = q · means of this methodology the most likely unobserved  3 3   β0 + β1 · t + β2 · a + β3 · s + β4 · a ∗ s + β5 · sin 2·π·t + real GW incidence process is reconstructed based on the  2·π·t  3 β6 · cos 3 . The estimates and their associated classification (underreported or not underreported) given standard errors were obtained by maximizing the log- by the posterior probabilities for the observations, pro- likelihood function described in Eq. (2) and from its vided by the output of the mixtools procedure, and on the Hessian matrix respectively, using the nlm procedure in R estimates of the parameters. All the R code used to fit the [19]. models and to obtain the reported results and figures is ⎛ available as Supplementary material.  n (yt −μ1,t )2 1 l(Y , θ ) = log ⎝(1 − ωt ) √ e 2σ 2 Validation 2π σ t=1 The goodness of fit of the proposed mixture approach ⎞ (2) (yt −μ2,t )2 can be assessed by means of the Akaike’s Information 1 2 2 ⎠, +ωt √ e 2q σ Criterion (AIC) compared to a single normal model. 2π qσ In this case, this measure favors the proposed model (AIC: 1717.9) in front of the single normal model (AIC: where Y = y1 , . . . , yn is the observed series, θ = eα0 +α1 t eγ 1826.1). The model has been validated by analyzing its (α0 , α1 , γ , β0 , . . . , β6 , σ ), ωt = 1+eα0 +α1 t , q = 1+eγ and μ1,t residuals. Figure 1 shows that they behave like white and μ2,t are as defined before. noise as expected and that there are no significant auto- In order to get proper initial values for the maximization correlations that should be accounted for. The residuals rt routine, an Expectation-Maximization (EM) algorithm for have been estimated as 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 2·π ·t rˆt = Yt − ω̂t · q̂ · βˆ0 + βˆ1 · t + βˆ2 · a + βˆ3 · s + βˆ4 · a ∗ s + βˆ5 · sin 3 this methodology is widely used when dealing with mix- 2·π ·t tures of distributions, it is unable to produce standard +βˆ6 · cos + (1 − ω̂t ) · βˆ0 + βˆ1 · t + βˆ2 · a + βˆ3 · s + βˆ4 · a ∗ s 3 errors directly [21], and this is an important drawback 2·π ·t 2·π ·t +βˆ5 · sin + βˆ6 · cos in our context and in many other situations. If the main 3 3 focus was not on quantifying the under-reporting issue, (3) ACF PACF 0.2 0.2 0.1 0.1 0.0 0.0 −0.1 −0.1 −0.2 −0.2 2.5 5.0 7.5 10.0 2.5 5.0 7.5 10.0 Lag Lag Fig. 1 Auto-correlations and partial auto-correlations of the model residuals Moriña et al. BMC Medical Research Methodology (2021) 21:6 Page 4 of 7 1 Sample quantiles 0 −1 −2 −2 −1 0 1 2 Teoretical quantiles Fig. 2 Residuals qqplot where Yt is the total observed GW incidence at time monthly average of 5.9 cases per 100,000 men for 7.1 t, and the letters with a hat (ˆ) represent the estimated cases per 100,000 men on the reconstructed series, 21.8% parameters. higher. Regarding males under 30 years old, the recon- If we were dealing with counts (number of cases) instead structed series is 13.3% higher (monthly averages of 18.4 of incidence, the underlying distribution might be a Pois- and 20.8 cases per 100,000 men for the registered and son, although the monthly number of GW cases is large reconstructed processes respectively). For women under enough to be approximated by a normal distribution. 30 years old, the monthly average registered incidence of Additionally, the assumption that the underlying distribu- GW in Catalunya is 19.0 per 100,000 women, while the tions of the two processes are Gaussian seems reasonable reconstructed hidden process has an average of 23.0 cases considering the qqplot of the residuals shown in Fig. 2. per 100,000 women, about 21.0% larger. This information is summarized in Table 1 and described in more detail in Results the Supplementary material (Table S1). Our analysis estimates that, globally, only around 80% Table 2 shows the estimated effect of the age and sex of actual GW incidence was registered in the SIDIAP over the under-reporting issue. In particular, it can be seen database in the period 2009-2016. For women over 30 that the GW incidence is higher among younger popu- years old, the monthly average registered incidence is 3.9 lations and men. It can also be noticed that a significant cases per 100,000 women, while the estimated monthly interaction between sex and age group is found, which can incidence is 4.9 cases per 100,000 women, 24.9% higher. be interpreted as a distinguishable impact of sex on GW On males over 30 years old, the registered series has a incidence depending on the age group. Table 1 Registered and estimated GW monthly average incidence (number of cases x 100,000 individuals) in the period 2009-2016 Sex Age Incidence (registered) Incidence (estimated) Difference (%) 15-29 19.0 23.0 21.0% Females 30-94 3.9 4.9 24.9% Average 6.8 8.4 23.2% 15-29 18.4 20.8 13.3% Males 30-94 5.9 7.1 21.8% Average 8.3 9.8 18.3% Global 7.6 9.1 19.9% Moriña et al. BMC Medical Research Methodology (2021) 21:6 Page 5 of 7 Table 2 Parameter estimates The under-reporting frequency is about 95% in 2009 Covariate Parameter Estimate (95% CI) (ω1 ) and around 21% in 2016 (ω96 ). This is measured by α0 2.99 (1.77; 4.20) parameter α1 in model (1), and should not be confused to t α1 -4.31 (-6.53; -2.09) overall under-reporting of the data, as its intensity (mea- sured by parameter q in the model) also plays a crucial β0 13.76 (7.11; 20.40) role. For instance, all observations in a certain period of t β1 0.36 (-12.75; 13.46) time could be slightly under-reported (ω = 1, q near age β2 -13.53 (-14.13; -12.92) to 1), resulting in small differences between registered sex β3 -1.60 (-2.24; -0.95) and estimated values or just a few observations might be age ∗ sex β4 3.25 (2.44; 4.06) under-reported (ω near to zero) but with a high intensity β5 4.16 (0.44; 7.88) (q near to zero), potentially resulting on large differences between registered and estimated values. Table 3 shows β6 0.52 (-5.59; 6.64) the total number of GW cases registered in the SIDIAP in q 0.75 (0.72; 0.77) 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 Figure 3 shows the registered (solid black line) and area outside ICS coverage is the same. reconstructed unobserved (dashed red line) processes for each of the considered sub-populations. Although this Discussion figure shows increasing trends for all series, they are The results of this work show that in relative terms, not well explained by coefficient β1 , which is not signif- the under-reporting issue has a deeper impact on people icantly different from zero. Increasing trends are mainly over 30 years old (where GW incidence is lower), espe- explained by the significant coefficient α1 , which leads to cially among women. Nonetheless, the relative difference a decreasing frequency of under-reporting ωt . between registered and estimated annual averages range Women 15 29 years old Women 30 94 years old 8 30 7 Incidence x 100,000 Incidence x 100,000 6 25 5 20 4 3 15 2 10 1 2010 2012 2014 2016 2010 2012 2014 2016 Time Time Men 15 29 years old Men 30 94 years old 12 30 10 25 Incidence x 100,000 Incidence x 100,000 8 20 6 15 4 10 2 5 2010 2012 2014 2016 2010 2012 2014 2016 Time Time Fig. 3 Registered (solid black line) and estimated underlying series (dashed red line) for each of the considered sub-populations Moriña et al. BMC Medical Research Methodology (2021) 21:6 Page 6 of 7 Table 3 Registered, estimated and projected number of GW cases in Catalunya Sex Age SIDIAP (registered) SIDIAP (estimated) Catalunya (registered projection) Catalunya (estimated projection) 15-29 8,051 9,769 10,280 12,460 Females 30-94 7,625 9,520 9,062 11,337 Total 15,676 19,289 19,342 23,797 15-29 7,967 9,097 10,166 23,797 Males 30-94 10,774 13,842 12,914 16,598 Total 18,741 22,939 23,080 28,182 Global 34,417 42,228 42,422 51,979 between 13.3% and 24.9%. It is also remarkable that the although about 2.8 million Euros would correspond to pri- quality of SIDIAP register regarding GW in Catalunya has vate insurances. It is, therefore, clear that knowing the true been significantly improving during the study period, as burden of GW at the general population level is impor- the frequency of under-reported observations has been tant for health policy makers, especially after the intro- decreasing over time. Facing under-reported information duction of prophylactic vaccines against HPV in many from public health registers is very common in many countries, as it plays a crucial role in developing and eval- situations, especially regarding potentially asymptomatic uating prevention strategies [24, 25]. This work presents diseases like GW. The proposed methodology considers a methodology that opens a wide field for future research the potential under-reporting in continuous time series lines. In particular, if temporal correlations are found in data in a very flexible way, estimating its frequency and the data, an appropriate model should take this structure intensity, and it is general enough to be appropriate in a into account, similarly to [13, 18]. wide range of real situations in the public health context. Additionally, the most likely non-observed process can be Supplementary Information reconstructed on the basis of estimated posterior proba- The online version contains supplementary material available at https://doi.org/10.1186/s12874-020-01188-4. bilities. Moreover, the GW data show that these models can deal with time-dependent under-reporting parame- ters, seasonal behavior, trends and also incorporate the Additional file 1: Additional tables and R code to reproduce the analyses. effect of other factors by including covariates. Abbreviations One of the potential limitations of this study is that AIC: Akaike’s information criterion; CI: Confidence interval; ECAP: e-Centre the database used included data from the public health- d’atenció primària (e-Primary attention center); EM: Expectation-maximization; care setting and not from the private sector. In Catalunya, ICS: Institut català de la salut (Catalan health institute); GW: Genital warts; SIDIAP: Sistema d’informació per al desenvolupament de la investigació en it is estimated that 33% of women and 25% of men atenció primària (Information system for the development of primary aged 15 to 44 years have a double health insurance attention research) coverage (i.e. the public health insurance and a pri- vate insurance plan) [12], so the rates estimated in our Acknowledgements We acknowledge the SIDIAP, with special thanks to Maria Aragon for her help study are likely still underestimating the real incidence in data collection. We thank CERCA Programme / Generalitat de Catalunya for of GW. One of its strengths is that the same method- institutional support. ology (possibly with minor model modifications) could Authors’ contributions be used to analyze the frequency and intensity of poten- DM, AFF, AC and PP developed the statistical model. DM implemented it in R tial under-reporting issues for any condition or set- software. LM, MB and MD helped providing an epidemiological context for the ting in the absence of temporal dependence among the methodology in the framework of genital warts and enriched the discussion focusing on the relevance on public health practice. All authors have read and observations. approved the manuscript. Conclusions Funding David Moriña acknowledges financial support from the Spanish Ministry of The GW incidence registered in SIDIAP is underestimat- Economy and Competitiveness, through the María de Maeztu Programme for ing the real burden in almost 10,000 cases in Catalunya, Units of Excellence in R&D (MDM-2014-0445) and Fundación Santander around 23% of the registered cases. The annual per person Universidades. This work has partial funding promoted by the Department of Health of the Generalitat de Catalunya for the execution of the project cost of GW was around 1000 Euros [8], so the poten- “Monitorización y evaluación del impacto de la introducción de nuevas tial total annual cost is underestimated in at least about estrategias preventivas del cáncer de cuello de útero en Catalunya” (reference 10 million Euros respect the 54 million Euros devoted to 0599S/7613/2010). This work was partially funded by the Instituto de Salud Carlos III-ISCIII (Spanish Government) through the projects PIE16/00049, GW in Catalunya annually, representing 0.4% of the total PI16/01254, PI16/01056, PI19/01118, Río Hortega CM15/00061 (Co-funded by budget of the Catalan Government intended for health, FEDER funds / European Regional Development Fund. ERDF, a way to build Moriña et al. BMC Medical Research Methodology (2021) 21:6 Page 7 of 7 Europe), Agència de Gestió d’Ajuts Universitaris i de Recerca (2017SGR1718) 10. Hsueh P-R. Human papillomavirus, genital warts, and vaccines. J Microbiol and by grant RTI2018-096072-B-I00 from the Spanish Ministry of Science, Immunol Infect Wei mian yu gan ran za zhi. 2009;42(2):101–6. Innovation and Universities. The funding sources had no role in the data 11. Escobar N, Plugge E. Prevalence of human papillomavirus infection, collection, analysis or interpretation of the results. cervical intraepithelial neoplasia and cervical cancer in imprisoned women worldwide: A systematic review and meta-analysis. BMJ Publ Availability of data and materials Group. 2020. https://doi.org/10.1136/jech-2019-212557. R codes used in the described analyses are available as Supplementary 12. Enquesta de salut de Catalunya (ESCA 2017) [Internet]. Departament de material. Salut. 2017 [cited 2019 Mar 31]. https://salutweb.gencat.cat/web/. content/_departament/estadistiques-sanitaries/enquestes/Enquesta-de- Ethics approval and consent to participate salut-de-Catalunya/Resultats-de-lenquesta-de-salut-de-Catalunya/ This study was approved by the Clinical Research Ethics Committee and the documents/els40prals_2017_web.xlsx. Institutional Review Board of the University Institute for Primary Care Research 13. Fernández-Fontelo A, Cabaña A, Puig P, Moriña D. Under-reported data (IDIAP) Jordi Gol (P15/106). analysis with INAR-hidden Markov chains. Stat Med. 2016;35(26):4875–90. https://doi.org/10.1002/sim.7026. Consent for publication 14. Bernard H, Werber D, Höhle M. Estimating the under-reporting of Not applicable. norovirus illness in Germany utilizing enhanced awareness of diarrhoea during a large outbreak of Shiga toxin-producing E. coli O104:H4 in Competing interests 2011–a time series analysis. BMC Infect Dis. 2014;14:116. https://doi.org/ The Cancer Epidemiology Research Programme, where LM, MB and MD are 10.1186/1471-2334-14-116. affiliated to, has received institutional sponsorship for grants from Merck and 15. Alfonso JH, Løvseth EK, Samant Y, Holm J.-Ø.. Work-related skin diseases GlaxoSmithKline. in Norway may be underreported: data from 2000 to 2013. Contact Dermatitis. 2015;72(6):409–12. https://doi.org/10.1111/cod.12355. Author details 16. Rosenman KD, Kalush A, Reilly MJ, Gardiner JC, Reeves M, Luo Z. How 1 Department of Econometrics, Statistics and Applied Economics, much work-related injury and illness is missed by the current national Riskcenter-IREA, Universitat de Barcelona (UB), Avinguda Diagonal, 690, 08034 surveillance system? J Occup Environ Med. 2006;48(4):357–65. https://doi. Barcelona, Spain. 2 Centre de Recerca Matemàtica, Universitat Autònoma de org/10.1097/01.jom.0000205864.81970.63. Barcelona (UAB), Edifici C, Campus de Bellaterra, 08193 Cerdanyola del Vallès, 17. Arendt S, Rajagopal L, Strohbehn C, Stokes N, Meyer J, Mandernach S. Spain. 3 Chair of Statistics, School of Business and Economics, Reporting of foodborne illness by U.S. consumers and healthcare Humboldt-Universität zu Berlin, Berlin, Germany. 4 Barcelona Graduate School professionals,. Int J Env Res Public Health. 2013;10(8):3684–714. https:// of Mathematics (BGSMath), Departament de Matemàtiques, Universitat doi.org/10.3390/ijerph10083684. Autònoma de Barcelona (UAB), Edifici C, Campus de Bellaterra, 08193 18. Fernández-Fontelo A, Cabaña A, Joe H, Puig P, Moriña D. Untangling Cerdanyola del Vallès, Spain. 5 Unit of Infections and Cancer - Information and serially dependent underreported count data for gender-based violence. Interventions (UNIC - I&I), Cancer Epidemiology Research Program (CERP), Stat Med. 2019;38(22):4404–22. https://doi.org/10.1002/sim.8306. Catalan Institute of Oncology (ICO)-IDIBELL, L’Hospitalet de Llobregat, Spain. 19. R Core Team. R: A Language and Environment for Statistical Computing. 2019. https://www.r-project.org/. 20. Benaglia T, Chauveau D, Hunter DR, Young D. mixtools : An R Package Received: 2 July 2020 Accepted: 4 December 2020 for Analyzing Finite Mixture Models. J Stat Softw. 2009;32(6):1–29. https:// doi.org/10.18637/jss.v032.i06. 21. Jamshidian M, Jennrich RI. Standard errors for EM estimation. J R Stat Soc References Ser B (Stat Methodol). 2000;62(2):257–70. https://doi.org/10.1111/1467- 1. Groseclose SL, Buckeridge DL. Public Health Surveillance Systems: Recent 9868.00230. Advances in Their Use and Evaluation. Ann Rev Public Health. 2017;38(1): 22. Lee Y, Nelder JA. Double Hierarchical Generalized Linear Models. J R Stat 57–79. https://doi.org/10.1146/annurev-publhealth-031816-044348. Soc Ser B (Methodol). 1996;58(4):619–78. 2. Ford MA, Spicer CM. Monitoring HIV Care in the United States: Indicators 23. Molas M, Lesaffre E. Hierarchical generalized linear models: The R and Data Systems; 2012. http://www.nap.edu/catalog.php?record_id= package HGLMMM. J Stat Softw. 2011;39(13):1–20. https://doi.org/10. 13225. Accessed 12 Apr 2019. 18637/jss.v039.i13. 3. SIDIAP, Information System for Research in Primary Care [Internet]. [cited 24. Kjær S, Tran T, Sparen P, Tryggvadottir L, Munk C, Dasbach E, Liaw K, 2019 Mar 29]. https://www.sidiap.org/. Nygård J, Nygård M. The Burden of Genital Warts: A Study of Nearly 70,000 4. McCormack D, Koons K. Sexually Transmitted Infections. W.B. Saunders. Women from the General Female Population in the 4 Nordic Countries. J 2019. https://doi.org/10.1016/j.emc.2019.07.009. Infect Dis. 2008;196(10):1447–54. https://doi.org/10.1086/522863. 5. Brotons M, Monfil L, Roura E, Duarte-Salles T, Casabona J, Urbiztondo L, 25. Kostaras D, Karampli E, Athanasakis K. Vaccination against HPV virus: a Cabezas C, Bosch FX, de Sanjosé S, Bruni L. Impact of a single-cohort systematic review of economic evaluation studies for developed HPV vaccination strategy with quadrivalent vaccine in northeast Spain: countries: Taylor and Francis Ltd; 2019. https://doi.org/10.1080/14737167. Population-based analysis of genital warts in men and women. In: 2019.1555039. EUROGIN; 2018. https://www.eurogin.com/content/dam/Informa/ eurogin/previous/Abstracts-Eurogin-2018.pdf. Publisher’s Note 6. Woodhall SC, Jit M, Soldan K, Kinghorn G, Gilson R, Nathan M, Ross JD, Springer Nature remains neutral with regard to jurisdictional claims in Lacey CJN, study group Q. The impact of genital warts: loss of quality of published maps and institutional affiliations. life and cost of treatment in eight sexual health clinics in the UK. Sex Transm Infect. 2011;87(6):458–63. https://doi.org/10.1136/sextrans-2011- 050073. 7. Sénécal M, Brisson M, Maunsell E, Ferenczy A, Franco EL, Ratnam S, Coutlée F, Palefsky JM, Mansi JA. Loss of quality of life associated with genital warts: baseline analyses from a prospective study,. Sex Transm Infect. 2011;87(3):209–15. https://doi.org/10.1136/sti.2009.039982. 8. Castellsagué X, Cohet C, Puig-Tintoré LM, Acebes LO, Salinas J, Martin MS, Breitscheidel L, Rémy V. Epidemiology and cost of treatment of genital warts in Spain. Eur J Public Health. 2009;19(1):106–10. https://doi. org/10.1093/eurpub/ckn127. 9. Blomberg M, Dehlendorff C, Kjaer SK. Risk of CIN2+ following a diagnosis of genital warts: A nationwide cohort study. Sex Transm Infect. 2019;95(8): 614–8. https://doi.org/10.1136/sextrans-2019-054008.