Moriña et al. BMC Medical Research Methodology (2021) 21:277 https://doi.org/10.1186/s12874-021-01427-2 RESEARCH Open Access Analysis of zero inflated dichotomous variables from a Bayesian perspective: application to occupational health David Moriña1,2*, Pedro Puig2,3 and Albert Navarro4,5 Abstract Background: Zero-inflated models are generally aimed to addressing the problem that arises from having two differ- ent sources that generate the zero values observed in a distribution. In practice, this is due to the fact that the popula- tion studied actually consists of two subpopulations: one in which the value zero is by default (structural zero) and the other is circumstantial (sample zero). Methods: This work proposes a new methodology to fit zero inflated Bernoulli data from a Bayesian approach, able to distinguish between two potential sources of zeros (structural and non-structural). Results: The proposed methodology performance has been evaluated through a comprehensive simulation study, and it has been compiled as an R package freely available to the community. Its usage is illustrated by means of a real example from the field of occupational health as the phenomenon of sickness presenteeism, in which it is reason- able to think that some individuals will never be at risk of suffering it because they have not been sick in the period of study (structural zeros). Without separating structural and non-structural zeros one would be studying jointly the general health status and the presenteeism itself, and therefore obtaining potentially biased estimates as the phe- nomenon is being implicitly underestimated by diluting it into the general health status. Conclusions: The proposed methodology is able to distinguish two different sources of zeros (structural and non- structural) from dichotomous data with or without covariates in a Bayesian framework, and has been made available to any interested researcher in the form of the bayesZIB R package (https://​cran.r-​proje​ct.​org/​packa​ge=​bayes​ZIB). Keywords: Presenteeism, Bayesian methods, Zero-inflation, Simulation study, Bernoulli mixture models Background and the other is circumstantial (sample zero). An exam- In general, zero-inflated models are aimed to address- ple could be the study of sickness presenteeism (SP), ing the problem that arises from having two different i.e. attending work while sick [1]. If it is not previously sources that generate the zero values observed in a dis- restricted, the population is made up, among others, of tribution. In practice, this is due to the fact that the popu- workers who are zero because they have never been sick lation studied actually consists of two subpopulations: (structural zeros) and workers who, having been sick, did one in which the value zero is by default (structural zero) not attend their work place (sample zeros). Note that the difference is important: roughly the first zero informs us exclusively about the status of health, the second about *Correspondence: dmorina@ub.edu 1 Department of Econometrics, Statistics and Applied Econometrics, the exercise of the right to take a sick leave. Riskcenter‑IREA, Universitat de Barcelona (UB), Avinguda Diagonal 690, The most commonly used zero-inflated models are 08034 Barcelona, Spain those that are related to count- ing variables, where it is Full list of author information is available at the end of the article © The Author(s) 2021. 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://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/. The Creative Commons Public Domain Dedication waiver (http://​creat​iveco​ mmons.​org/​publi​cdoma​in/​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:277 Page 2 of 9 assumed that the zero value has a dichotomous source SP, and the results obtained are compared with those that determines whether or not the subject is at risk of of adjusted logistic regressions on the total population suffering the event of interest and another source, only or only in those individuals at risk. In the literature, the for the individuals at risk, that corresponds to the num- SP registry is carried out in a self-reported way, asking ber of episodes (counts) that have been experienced by about the episodes in the last year and later, recorded each individual at risk. In this context, the most com- in a dichotomized way (no SP: 0 episodes; yes SP: 1 or mon available models would be the well known Zero- more episodes). The justification for this dichotomiza- Inflated Poisson (ZIP) and Negative Binomial (ZINB). tion is fundamentally based on two aspects: one, the A good introduction to the mathematical properties of possible memory bias; second, the excessive influence these models can be found in [2], and they have been of workers who report a very high number of episodes. used in many fields such as quality control ([3]), epi- demiology ([4]) or medicine ([5]) among many others. Methods Some guidelines on how to proceed when dealing with Let Y be the variable that indicates occurrence of the count outcomes potentially overdispersed or zero- phenomenon under study. The proposed model has a inflated have been published recently ([6, 7]), based on probability function defined by classical procedures like Vuong’s test ([8]) to check for P (Y = 0) = (1 − 𝜔) + 𝜔 ⋅ (1 − p) overdis- persion ([9]) and zero-inflation ([10]), although these guides cannot be applied to the case studied here (2) due to the dichotomous nature of the outcome. In gen- P (Y = 1) = 𝜔 ⋅ p, eral, zero-inflated models can be expressed as where ω is the probability of exposure and p is the probability of occurrence of the phenomenon of interest    = 0) P(Y  =g + 1 − g · f (0) (1) P Y = j = 1 − g · f (j), j > 0 among exposed individuals, as shown in Fig. 1. Accord- ing to this scheme, the proportion of structural zeros will where g is the structural zero probability and f (0) be 1 − ω and the proportion of non-structural zeros will is the zero probability of an appropriate distribution be ω · (1 − p). (Poisson, negative binomial or Bernoulli as in our case). To overcome the impossibility of models without covar- In practice, zero-inflated models with both dichoto- iates based on the frequen- tist approach to differentiate mous sources (a mixture of two Bernoulli random vari- between structural and non-structural zeros, in this work ables, one with probability of success ω and the other a model within the Bayesian framework is proposed. In with probability of success p) have received far less this context, we as- sume that the prior distribution of attention. This is due, in large part, to the fact that the the parameter of the first Bernoulli ω is uniform between resulting distribution is once again a Bernoulli with 0 and 0.5 while the prior distribution of the probability probability. of success in the second Bernoulli p is uniform between of success ω · p, so that the proportion of structural zeros (1 − ω) and sample zeros (1 − p) are indistin- guishable from the point of view of frequentist statis- tics. However, from the Bayesian perspective and using known reasonable information about these propor- tions, it is possible to distinguish the two sources of zeros and estimate ω and p. Some authors have recently suggested, in other areas such as the classification or identification of images, the usage of Bernoulli-mixture models, based on numerical algorithms such as Expectation-Maximiza- tion (EM) to estimate the parameters [11, 12], given the complexity of the likelihood functions involved. In these cases, however, the inclusion of covariates or adjustment variables is virtually impossible. Also in other areas there are some recent developments in a similar line, such as [13]. In this article we illustrate the use of Zero Inflated Bernoulli (ZIB) models by means of a real dataset on Fig. 1 Model schema. Only variable Y is observed Moriña et al. BMC Medical Research Methodology (2021) 21:277 Page 3 of 9 0.5 and 1. In this way, the proposed model will be able where m is the frequency of occurrence of the phenome- to distinguish the two sources of zeros. Obviously, these non of interest and n is the total number of observations. hypotheses can be modified based on prior knowledge From here, the joint posterior could be obtained as of the parameters that govern the phenomenon under investigation by making simple changes to the posterior f (p, ω) ∼ (p · ̟ )m · (1 − p · ω)n−m · U[0,1/2] (ω) · U[1/2,1] (p) (5) distributions defined in Eq. (6) and recalculating the mar- ginals shown in Eq. (5). This distinction is not necessary From here the posterior marginal distributions of the if covariates are included in the model since the covari- two parameters can be obtained as ates allow the origin of the zeros to be distinguished. To ensure that the estimates are kept within the appropriate 1 f (ω) ∼ ωm 1/2 pm · (1 − p · ω)n−m dp ∼ parameter space, the logit link, commonly used in logistic 1  ω m n−m ω · ω/2 t · (1 − t) dt ∼ regression, has been used. ω F (ω,m+1,n−m+1)−F ( 2 ,m+1,n−m+1) ω    1/2 (6) logit(ω) = log 1−ω ω = θ0 + θ1 X1 + · · · + θk Xk f (p) ∼ pm o ωm · (1 − p · ω)n−m dω ∼ (3) 1 p/2 m n−m  p · o t · (1 − t) dt ∼   p logit(p) = log 1−p = β0 + β1 Z1 + . . . βm Zm, p F ( 2 ,m+1,n−m+1) p where X1,. .., Xk are the covariates that have a hypotetical impact over the zero inflated part and Z1,. .., Zm are the where F is the beta distribution function with parameters covariates that might have an influence over the non zero m + 1 and n − m + 1, implemented in the R function inflated part. The parameters θi, i = 0,. .., k and βj, j = 0,. .., pbeta. m are assumed to follow a normal distribution with mean The methods used to analyse the real data example in 0 and variance σ2 and σθ2 and σβ2 respectively, modeled as the following sections are in accordance with relevant hyperparameters. guidelines and regulations, in particular with the Inter- The models proposed to analyse the data described in national Labour Organization criteria, also used to define the following section and in the simulation study have the target population in the European Working Condi- been written in the programming language Stan, within tions Survey [17] or the EU Labour Force Survey [18]. the R environment [14] and are freely available from the Participation in the considered study was voluntary and authors as a package called bayesZIB [15]. To the best confidential, and informed consent was obtained from of our knowledge, this is the only package available in all subjects in order to be included. The data were anal- R able to fit zero-inflated Bernoulli regression models. ysed anonymously and all procedures were approved by The use of the package is very similar to other packages the Ethics Committee on Animal and Human Experi- that implement zero-inflated models, such as pscl [16], mentation of the Autonomous University of Barcelona to facilitate the interpretation of the results, while more (CEEAH/3445). advanced users could easily adapt the code to their spe- cific requirements. If necessary, appropriate priors for the Results parameters ω and p can be defined in the function bayes- This section presents the results of the analyses using the ZIB using the argument priors (only uniforms with differ- proposed methodology over a real data set and they are ent parameters are implemented so far in the package). compared to the most common alternatives. The perfor- mance of the method is also studied by means of a com- No covariates prehensive simulation study, with and without covariates. In the particular case in which the interest is in estimat- ing the proportion of struc- tural (1-ω) and sample (1-p) Real data zeros without accounting for the effect of any covariate, In the database used to exemplify the use of the pro- the posterior distributions of ω and p can be obtained posed methodology, we have a total of n = 1564 workers. analytically assuming some a priori knowledge of their Among these, it is known that 946 (around 61%) were not distributions. As mentioned before, one could set ω to at risk of being presenteeist because they were not ill on be uniform distributed on [0, 0.5] and p to be uniform any day during the study period. These observations cor- distributed on [0.5, 1]. Be- cause the observations are respond to the concept of structural zeros (1 − ω = 0.61), Bernoulli(p · ω) distributed, the likelihood function can and an estimate of their proportion can be obtained by be written as using zero inflated models, even taking into account the values of the variables used as explanations in the regres- L ∼ (p · ω)m · (1 − p · ω)n−m , (4) sion model. The proportion of presenteeists among those exposed is p = 0.70. Globally, a total of m = 430 workers Moriña et al. BMC Medical Research Methodology (2021) 21:277 Page 4 of 9 Fig. 2 Prior (left column) and posterior (right column) distributions for ω and p experienced the event of interest. In the following sub- interpretation since it is significantly underestimating sections CI is used as an abbreviation of the confidence the proportion of presenters if the subjects who have not interval for frequentist analyses, and CrI is used for cred- been at risk of being present are excluded, since it is ulti- ibility interval when referring to the proposed Bayesian mately an estimate of ωˆ · pˆ, being impossible to identify model. the two parameters. Additionally, the proposed model allows incorporat- Including the whole population ing covariates in both pro- cesses. To illustrate how it Taking all the population into consideration (n = 1564), i. works, we will consider here the self-perceived general e., including individuals at risk and not at risk (those who state of health (categorized as good or bad) and the were not at risk during the study period), we fit a Bayes- feeling of being replaceable, which is an item included ian zero inflated Bernoulli model, where the proportion in the vulnerability dimension of the Employment Pre- of structural zeros 1 − ω is greater than 0.5 (prior uni- cariousness Scale [20], with categories “Always”, “Some- form for ω at [0, 0.5]) and the proportion of sample zeros times” and “Never”. The hypothesis is that the general 1 − p is less than 0.5 (prior uniform for p at [0.5, 1]). This state of health would be related to the risk of be- ing information is extracted from [19]. In this case, without present (zero inflated part of the model) and the feel- using covariates, the model allows estimating the values ing of being replaceable would be related to presenting of ωˆ = 0.37 (95% CrI: 0.27–0.49) and pˆ = 0.74 (95% CrI: the phenomenon once a worker is exposed (non zero 0.55–0.99). Here ωˆ and pˆ indicate the median of the inflated part in the model), so the model is including marginal posterior of ω and p respectively. The a priori one covariate in each part (k = m = 1). The model was and a posteriori marginal distributions of both param- fitted using 5 Markov chains, 5000 iterations in each eters are shown in Fig. 2. The different shapes between chain (half used for warmup and half for inference), a marginal priors and posteriors show that the models target average acceptance probability of 0.999 and a learn from the data. maximum allowed treedepth of 25. Notice that these On the other hand, analysing these data as is traditional technical values might need to be changed depending in the literature, using a logistic regression model without on the analysed data. As can be seen in the next sec- taking into account that there are subjects who have not tion, the results of the model shown in Table 1 largely been at risk, the proportion obtained from presentists co- incide with the standard logistic analyses reported is 0.27 (95% CI: 0.25–0.30), a value with a controversial in Table 2, particularly in the direction and impact of Moriña et al. BMC Medical Research Methodology (2021) 21:277 Page 5 of 9 Table 1 Bayesian analysis on whole population. CrI stands for effect magnitudes, the following values for each param- credible interval eter were considered: Covariate Coefficient. (95% CrI) • β0 = 0.5, 1, 2 Struct Intercept −0.40 (− 0.83, 0.22) • β1 = 2, 3, 4 Bad 1.43 (0.83, 2.58) • β2 = 3 Non-struct. Intercept 1.03 (0.16, 2.32) • θ0 = − 0.5, − 1, − 2 Sometimes 0.03 (−0.73, 1.37) • θ1 = − 2, − 3, − 4 Never −0.68 (− 1.38, − 0.15) • θ2 = − 3 It is important to notice that this is an extreme situa- tion, in which we consider that all the mass probability of Table 2 Logistic regression on all population and only on exposed individuals. CI stands for confidence interval the parameter distributions is concentrated in one point, the “true” value of the parameter. Population Covariate Coefficient. (95% CI) For each random sample, the posterior marginal distri- Whole population Intercept 1.07 (− 1.19, − 0.95) butions of the parameters have been summarised by their Bad 1.05 (0.68, 1.42) median and percentiles 2.5 and 97.5%. Only exposed Intercept 1.16 (0.79, 1.53) Tables 3 and 4 show, for each combination of param- Sometimes −0.02 (− 0.60, 0.56) eters, the average estimates and upper and lower limits of Never − 0.54 (− 0.97, − 0.11) the 95% credibility intervals. As no relevant differ- ences were observed regarding sample sizes, Tables 3 and 4 shows only the results corresponding to n = 1500. The results corresponding to n = 500 are available as Sup- associations. All R codes used in this paper are available plementary Material. It can be seen that in all cases the as Supplementary Material. original parameters used to generate the simulations can Similarly to other regression models, the effect of never be properly recovered by the fitted models. having the feeling of being replaceable over the odds of The R code used for the simulation is available as Sup- attending work when sick for someone who is at risk plementary Material. An additional simulation was con- compared to workers who always have that sensation can ducted to evaluate the performance of the proposed be quantified by e− 1.04 = 0.35. methodology when there are no covariates involved, the details and results of this simulation can also be find in Excluding healthy population the Supplementary Material (Appendix A, Table S2). If the information about which subjects are actually exposed to the phenomenon of interest is available (an Discussion ideal but unusual situation in practice), unexposed sub- The proposed methodology is able to distinguish two jects could be excluded and a logistic regression model different sources of zeros (structural and non-struc- could be adjusted to the exposed individuals. Using the tural) from dichotomous data in a Bayesian framework same explanatory variables described in the previous sec- by assuming priors with different parameters on pro- tion, the corresponding coefficients and their 95% confi- portion of structural and non- structural zeros. Fur- dence intervals are described in Table 2. thermore, since it is freely available as an R package, it is easily usable for any researcher who needs to adjust this Simulation study type of data and easily modifiable for more advanced In order to check the performance of the proposed meth- users who need to adapt the model to their context, for odology, 100 random samples were generated for each example with different choices of the prior distributions considered sample size (n = 500, 1500), and combina- of ω and p. tion of parameters. The zero inflated part was build upon The approach used to analyse the SP is an important the logistic regression model logit(P (X = 1)) = θ0 + θ1 · topic. Some studies include all working population to x1 + θ2 · x2, where x1 and x2 are two independent covari- estimate SP, whilst other exclude “healthy” workers. As ates, each following a standard normal distribution. The result, different conclusions in terms of prevalence and non zero inflated part was build upon the logistic regres- associated factors are obtained [21]. SP is an outcome sion model logit(P (Y = 1 | X = 1)) = β0 + β1 · x3 + β2 · resulting from mixing two phenomena, i.e. health status x4, where x3 and x4 are two independent covariates, each and exercise of rights. Health status plays a role regard- with a standard normal distribution. To cover different ing the fact of being exposed; and, among the exposed, Moriña et al. BMC Medical Research Methodology (2021) 21:277 Page 6 of 9 Table 3 Simulation study results including covariates (I) β0 β1 θ0 θ1 β0 (95% CrI) β1 (95% CrI) β2 (95% CrI) θ0 (95% CrI) θ1 (95% CrI) θ2 (95% CrI) −2 0.5 (0.2, 0.9) 2 (1.6, 2.5) 3 (2.5, 3.7) −0.5 (− 0.8, − 0.2) − 2 (− 2.5, − 1.7) − 3.1 (− 3.6, − 2.6) −0.5 −3 0.5 (0.2, 0.9) 2 (1.6, 2.5) 3 (2.5, 3.7) −0.5 (− 0.8, − 0.2) − 3 (− 3.6, − 2.5) − 3 (− 3.6, − 2.5) −4 0.5 (0.2, 0.8) 2 (1.6, 2.5) 3 (2.5, 3.6) − 0.5 (− 0.8, − 0.2) − 3.9 (− 4.7, − 3.3) − 2.9 (− 3.5, − 2.4) −2 0.5 (0.2, 0.9) 2 (1.6, 2.5) 3 (2.4, 3.7) − 1 (− 1.3, − 0.7) − 2 (− 2.5, − 1.7) −3.1 (− 3.6, − 2.6) 2 −1 −3 0.5 (0.2, 0.9) 2 (1.6, 2.5) 3 (2.4, 3.7) −1 (− 1.3, − 0.7) −3 (− 3.6, − 2.5) −3 (− 3.6, − 2.5) −4 0.5 (0.2, 0.9) 2 (1.6, 2.5) 3.1 (2.5, 3.7) −1 (− 1.3, − 0.7) −3.8 (− 4.6, − 3.2) − 2.9 (− 3.5, − 2.4) −2 0.5 (0.1, 1) 2 (1.5, 2.6) 3 (2.3, 3.8) − 2 (− 2.4, − 1.7) −2 (− 2.5, − 1.6) −3.1 (− 3.6, − 2.6) −2 −3 0.5 (0.1, 0.9) 2 (1.5, 2.6) 3 (2.4, 3.8) −2 (− 2.3, − 1.6) − 3 (− 3.6, − 2.4) − 3 (− 3.6, − 2.5) −4 0.5 (0.2, 0.9) 2 (1.5, 2.5) 3 (2.4, 3.7) −1.9 (− 2.3, − 1.6) − 3.9 (− 4.7, − 3.3) − 2.9 (− 3.5, − 2.4) −2 0.5 (0.2, 0.9) 3 (2.4, 3.7) 3 (2.4, 3.7) −0.5 (− 0.8, − 0.2) − 2 (− 2.4, − 1.7) − 3 (− 3.6, − 2.6) −0.5 −3 0.5 (0.2, 0.9) 3 (2.5, 3.7) 3 (2.4, 3.6) −0.5 (− 0.8, − 0.2) − 3 (− 3.6, − 2.5) − 3 (− 3.6, − 2.5) −4 0.5 (0.2, 0.8) 3 (2.5, 3.6) 3 (2.4, 3.6) − 0.5 (− 0.8, − 0.2) − 3.9 (− 4.7, − 3.3) − 2.9 (− 3.5, − 2.4) −2 0.5 (0.1, 0.9) 2.9 (2.3, 3.7) 2.9 (2.3, 3.7) −1 (− 1.2, − 0.7) − 2 (− 2.4, − 1.6) − 3 (− 3.5, − 2.5) 0.5 3 −1 −3 0.5 (0.2, 0.9) 2.9 (2.4, 3.6) 3 (2.4, 3.6) −1 (− 1.3, − 0.7) − 3 (− 3.6, − 2.5) − 3 (− 3.6, − 2.5) −4 0.5 (0.2, 0.8) 3 (2.4, 3.6) 3 (2.4, 3.6) −1 (− 1.3, − 0.7) −3.9 (− 4.7, − 3.3) − 2.9 (− 3.5, − 2.4) −2 0.4 (0, 0.9) 2.9 (2.2, 3.7) 2.9 (2.2, 3.7) −2 (− 2.4, − 1.7) −2 (− 2.4, − 1.6) − 3 (− 3.6, − 2.5) −2 −3 0.5 (0.1, 1) 3 (2.3, 3.8) 3 (2.3, 3.8) −2 (− 2.4, − 1.6) − 3 (− 3.6, − 2.5) − 3 (− 3.6, − 2.5) −4 0.5 (0.1, 0.9) 3 (2.4, 3.7) 3 (2.4, 3.7) −2 (− 2.4, − 1.6) − 3.9 (− 4.7, − 3.3) −2.9 (− 3.5, − 2.4) −2 0.4 (0.1, 0.8) 3.9 (3.2, 4.8) 2.9 (2.4, 3.6) −0.5 (− 0.7, − 0.2) − 2 (− 2.4, − 1.6) − 3 (− 3.5, − 2.5) −0.5 −3 0.5 (0.1, 0.8) 3.9 (3.2, 4.7) 2.9 (2.3, 3.6) −0.5 (− 0.8, − 0.2) − 3 (− 3.5, − 2.5) − 3 (− 3.5, − 2.5) −4 0.5 (0.2, 0.9) 4 (3.2, 4.8) 3 (2.4, 3.6) − 0.5 (− 0.8, − 0.2) − 3.9 (− 4.6, − 3.3) −2.9 (− 3.5, − 2.4) −2 0.5 (0.1, 0.9) 3.8 (3.1, 4.8) 2.9 (2.3, 3.6) −1 (− 1.3, − 0.7) − 2 (− 2.4, − 1.6) − 3 (− 3.5, − 2.5) 4 −1 −3 0.5 (0.1, 0.9) 3.9 (3.1, 4.8) 2.9 (2.3, 3.6) −1 (− 1.3, − 0.7) −3 (− 3.5, − 2.5) −3 (− 3.5, − 2.5) −4 0.5 (0.1, 0.9) 3.9 (3.2, 4.8) 2.9 (2.4, 3.6) − 1 (− 1.3, − 0.7) − 4 (− 4.7, − 3.3) − 3 (− 3.6, − 2.5) −2 0.4 (0, 0.9) 3.8 (2.9, 4.9) 2.9 (2.2, 3.7) −2 (− 2.4, − 1.7) −2 (− 2.5, − 1.7) −3 (− 3.6, − 2.5) −2 −3 0.5 (0.1, 0.9) 3.9 (3.1, 5) 2.9 (2.3, 3.7) −2 (− 2.4, − 1.6) − 3 (− 3.6, − 2.5) − 3 (− 3.6, − 2.5) −4 0.5 (0.1, 0.9) 3.9 (3.1, 4.8) 2.9 (2.3, 3.7) −2 (− 2.4, − 1.6) − 3.9 (− 4.7, − 3.3) −2.9 (− 3.5, − 2.4) −2 1 (0.7, 1.5) 2 (1.6, 2.5) 3 (2.5, 3.7) −0.5 (− 0.8, − 0.3) −2 (− 2.4, − 1.7) − 3 (− 3.6, − 2.6) −0.5 −3 1 (0.7, 1.4) 2 (1.6, 2.4) 3 (2.5, 3.6) −0.5 (− 0.7, − 0.2) − 3 (− 3.6, − 2.5) − 3 (− 3.5, − 2.5) −4 1 (0.7, 1.4) 2 (1.6, 2.5) 3 (2.5, 3.6) − 0.5 (− 0.8, − 0.2) − 3.9 (− 4.6, − 3.3) −2.9 (− 3.5, − 2.4) −2 1 (0.6, 1.4) 2 (1.5, 2.5) 3 (2.4, 3.7) −1 (− 1.2, − 0.7) −2 (− 2.4, − 1.7) − 3 (− 3.6, − 2.6) 2 −1 −3 1 (0.6, 1.4) 2 (1.6, 2.5) 3 (2.4, 3.6) −1 (− 1.3, − 0.7) − 3 (− 3.6, − 2.5) − 3 (− 3.6, − 2.5) −4 1 (0.7, 1.4) 2 (1.6, 2.5) 3 (2.5, 3.7) −1 (− 1.3, − 0.7) − 3.9 (− 4.6, − 3.3) − 2.9 (− 3.5, − 2.4) −2 1 (0.5, 1.6) 2 (1.4, 2.6) 3 (2.3, 3.8) −2 (− 2.3, − 1.7) −2 (− 2.4, − 1.6) − 3 (− 3.6, − 2.5) −2 −3 1 (0.6, 1.5) 2 (1.5, 2.6) 3 (2.4, 3.8) −2 (− 2.4, − 1.7) − 3 (− 3.5, − 2.5) − 3 (− 3.6, − 2.5) −4 1 (0.6, 1.5) 2 (1.5, 2.5) 3 (2.4, 3.7) −1.9 (− 2.3, − 1.6) − 3.9 (− 4.6, − 3.3) − 2.9 (− 3.5, − 2.4) −2 1 (0.6, 1.4) 2.9 (2.4, 3.6) 2.9 (2.3, 3.6) −0.5 (− 0.7, − 0.2) −2 (− 2.4, − 1.7) −3.1 (− 3.6, − 2.6) −0.5 −3 1 (0.6, 1.4) 3 (2.4, 3.6) 3 (2.4, 3.6) −0.5 (− 0.8, − 0.2) − 3 (− 3.5, − 2.5) −3 (− 3.5, − 2.5) −4 1 (0.7, 1.4) 3 (2.5, 3.6) 3 (2.4, 3.6) − 0.5 (− 0.8, − 0.2) − 4 (− 4.7, − 3.4) − 3 (− 3.6, − 2.5) −2 0.9 (0.5, 1.4) 2.9 (2.3, 3.6) 3 (2.4, 3.7) −1 (− 1.2, − 0.7) − 2 (− 2.4, − 1.7) − 3 (− 3.6, − 2.6) 1 3 −1 −3 1 (0.6, 1.4) 3 (2.4, 3.7) 3 (2.4, 3.7) −1 (− 1.3, − 0.7) − 3 (− 3.6, − 2.5) − 3 (− 3.5, − 2.5) −4 1 (0.6, 1.4) 3 (2.4, 3.6) 3 (2.4, 3.6) −1 (− 1.3, − 0.7) −3.9 (− 4.6, − 3.3) − 2.9 (− 3.5, − 2.4) −2 1 (0.5, 1.6) 3 (2.3, 3.9) 2.9 (2.3, 3.8) −2 (− 2.4, − 1.7) −2 (− 2.4, − 1.7) − 3 (− 3.6, − 2.6) −2 −3 0.9 (0.5, 1.5) 2.9 (2.3, 3.7) 2.9 (2.3, 3.7) −2 (− 2.3, − 1.6) − 3 (− 3.6, − 2.5) − 3 (− 3.6, − 2.5) −4 0.9 (0.5, 1.4) 2.9 (2.3, 3.7) 2.9 (2.3, 3.7) −1.9 (− 2.3, − 1.6) −3.9 (− 4.7, − 3.3) − 2.9 (− 3.5, − 2.4) −2 0.9 (0.6, 1.4) 3.9 (3.1, 4.7) 2.9 (2.3, 3.6) −0.5 (− 0.7, − 0.3) − 2 (− 2.4, − 1.7) −3 (− 3.5, − 2.6) −0.5 −3 1 (0.6, 1.4) 3.9 (3.2, 4.8) 2.9 (2.4, 3.6) −0.5 (− 0.8, − 0.2) − 3 (− 3.6, − 2.5) −3 (− 3.6, − 2.6) −4 0.9 (0.6, 1.3) 3.9 (3.2, 4.7) 2.9 (2.4, 3.6) − 0.5 (− 0.7, − 0.2) − 4 (− 4.7, − 3.4) − 3 (− 3.5, − 2.5) −2 0.9 (0.5, 1.4) 3.9 (3.1, 4.8) 2.8 (2.2, 3.6) −1 (− 1.3, − 0.7) − 2 (− 2.4, − 1.7) −3 (− 3.6, − 2.6) Moriña et al. BMC Medical Research Methodology (2021) 21:277 Page 7 of 9 Table 3 (continued) β0 β1 θ0 θ1 β0 (95% CrI) β1 (95% CrI) β2 (95% CrI) θ0 (95% CrI) θ1 (95% CrI) θ2 (95% CrI) 4 −1 −3 1 (0.6, 1.4) 3.9 (3.2, 4.8) 2.9 (2.4, 3.7) −1 (− 1.2, − 0.7) − 3 (− 3.6, − 2.5) − 3 (− 3.6, − 2.5) −4 1 (0.6, 1.4) 3.9 (3.2, 4.7) 2.9 (2.3, 3.6) −0.9 (− 1.2, − 0.7) −4 (− 4.7, − 3.3) −3 (− 3.5, − 2.5) −2 0.9 (0.4, 1.5) 3.8 (2.9, 4.9) 2.9 (2.2, 3.7) −2 (− 2.3, − 1.7) −2 (− 2.4, − 1.7) − 3 (− 3.5, − 2.5) −2 −3 1 (0.5, 1.5) 3.9 (3, 4.9) 2.9 (2.2, 3.7) −2 (− 2.4, − 1.7) −3 (− 3.5, − 2.5) −3 (− 3.6, − 2.5) −4 0.9 (0.5, 1.4) 3.9 (3.1, 4.9) 2.9 (2.3, 3.7) −1.9 (− 2.3, − 1.6) − 3.9 (− 4.6, − 3.3) − 3 (− 3.5, − 2.5) Table 4 Simulation study results including covariates (II) β0 β1 θ0 θ1 β0 (95% CrI) β1 (95% CrI) β2 (95% CrI) θ0 (95% CrI) θ1 (95% CrI) θ2 (95% CrI) −2 1.9 (1.5, 2.5) 2 (1.5, 2.5) 3 (2.4, 3.7) −0.5 (− 0.7, − 0.2) −2 (− 2.4, − 1.7) −3.1 (− 3.6, − 2.6) −0.5 −3 2 (1.6, 2.6) 2 (1.6, 2.5) 3 (2.4, 3.6) −0.5 (− 0.8, − 0.3) −3 (− 3.6, − 2.6) −3 (− 3.5, − 2.6) −4 2 (1.6, 2.5) 2 (1.6, 2.5) 3 (2.4, 3.6) − 0.5 (− 0.8, − 0.2) −3.9 (− 4.6, − 3.4) −3 (− 3.5, − 2.5) −2 1.9 (1.4, 2.5) 1.9 (1.5, 2.5) 2.9 (2.3, 3.6) −1 (− 1.2, − 0.8) − 2 (− 2.4, − 1.7) − 3.1 (− 3.6, − 2.6) 2 −1 −3 2 (1.5, 2.6) 1.9 (1.5, 2.5) 2.9 (2.4, 3.6) −1 (− 1.2, − 0.7) − 3 (− 3.5, − 2.6) −3 (− 3.5, − 2.5) −4 2 (1.6, 2.6) 2 (1.6, 2.5) 3 (2.4, 3.7) −1 (− 1.3, − 0.7) − 4 (− 4.7, − 3.4) − 3 (− 3.5, − 2.5) −2 1.9 (1.3, 2.7) 2 (1.5, 2.7) 2.9 (2.2, 3.8) −2 (− 2.3, − 1.7) −2 (− 2.4, − 1.7) −3 (− 3.5, − 2.6) −2 −3 2 (1.4, 2.7) 2 (1.5, 2.6) 2.9 (2.3, 3.8) −2 (− 2.3, − 1.7) −3 (− 3.5, − 2.5) −3 (− 3.5, − 2.5) −4 1.9 (1.4, 2.6) 1.9 (1.5, 2.5) 2.9 (2.3, 3.7) −2 (− 2.3, − 1.7) −4 (− 4.7, − 3.4) − 3 (− 3.5, − 2.5) −2 2 (1.5, 2.6) 3 (2.4, 3.7) 3 (2.4, 3.7) −0.5 (− 0.7, − 0.2) −2 (− 2.4, − 1.7) −3 (− 3.5, − 2.6) −0.5 −3 2 (1.5, 2.5) 3 (2.4, 3.6) 2.9 (2.4, 3.6) −0.5 (− 0.7, − 0.3) −3 (− 3.5, − 2.6) −3 (− 3.5, − 2.5) −4 1.9 (1.5, 2.5) 3 (2.4, 3.6) 2.9 (2.4, 3.6) −0.5 (− 0.7, − 0.2) −3.9 (− 4.6, − 3.4) − 2.9 (− 3.5, − 2.5) −2 1.9 (1.4, 2.6) 2.9 (2.3, 3.7) 2.9 (2.3, 3.7) −1 (− 1.3, − 0.8) − 2 (− 2.4, − 1.7) −3 (− 3.5, − 2.6) 2 3 −1 −3 1.9 (1.4, 2.5) 2.9 (2.3, 3.6) 2.9 (2.3, 3.6) −1 (− 1.3, − 0.7) −3 (− 3.5, − 2.5) −3 (− 3.5, − 2.6) −4 1.9 (1.5, 2.5) 2.9 (2.3, 3.6) 2.9 (2.4, 3.6) −1 (− 1.3, − 0.7) −4 (− 4.7, − 3.4) − 3 (− 3.6, − 2.5) −2 1.9 (1.3, 2.6) 2.9 (2.2, 3.7) 2.9 (2.2, 3.8) −2 (− 2.3, − 1.7) −2 (− 2.4, − 1.7) −3 (− 3.5, − 2.6) −2 −3 1.9 (1.3, 2.6) 2.9 (2.3, 3.7) 2.9 (2.2, 3.7) −2 (− 2.3, − 1.7) −3 (− 3.5, − 2.6) −3 (− 3.5, − 2.5) −4 2 (1.4, 2.6) 2.9 (2.3, 3.7) 2.9 (2.3, 3.7) −1.9 (− 2.3, − 1.6) − 3.9 (− 4.5, − 3.3) −2.9 (− 3.4, − 2.5) −2 1.9 (1.4, 2.5) 3.8 (3.1, 4.7) 2.9 (2.3, 3.6) −0.5 (− 0.7, − 0.3) −2 (− 2.4, − 1.7) −3.1 (− 3.5, − 2.6) −0.5 −3 1.9 (1.4, 2.5) 3.8 (3.1, 4.7) 2.9 (2.3, 3.6) −0.5 (− 0.7, − 0.2) −3 (− 3.6, − 2.6) −3 (− 3.6, − 2.6) −4 1.9 (1.4, 2.4) 3.8 (3.1, 4.7) 2.9 (2.3, 3.6) −0.5 (− 0.7, − 0.2) −4 (− 4.7, − 3.4) −3 (− 3.5, − 2.5) −2 1.9 (1.4, 2.6) 3.8 (3, 4.8) 2.9 (2.3, 3.7) −1 (− 1.2, − 0.8) − 2 (− 2.4, − 1.7) − 3 (− 3.5, − 2.6) 4 −1 −3 1.9 (1.4, 2.5) 3.9 (3.1, 4.8) 2.9 (2.3, 3.6) −1 (− 1.3, − 0.7) − 3 (− 3.5, − 2.6) −3 (− 3.5, − 2.6) −4 2 (1.5, 2.5) 3.9 (3.2, 4.8) 2.9 (2.3, 3.6) − 1 (− 1.2, − 0.7) −3.9 (− 4.6, − 3.4) −3 (− 3.5, − 2.5) −2 1.8 (1.2, 2.6) 3.8 (2.9, 4.9) 2.8 (2.1, 3.7) −2 (− 2.3, − 1.7) −2 (− 2.4, − 1.7) −3 (− 3.5, − 2.6) −2 −3 1.8 (1.3, 2.5) 3.7 (2.9, 4.7) 2.8 (2.1, 3.6) −2 (− 2.3, − 1.7) − 3.1 (− 3.6, − 2.6) −3 (− 3.6, − 2.6) −4 1.9 (1.4, 2.6) 3.8 (3, 4.8) 2.9 (2.2, 3.7) −1.9 (− 2.3, − 1.6) − 3.9 (− 4.6, − 3.3) −2.9 (− 3.5, − 2.5) the lack of the exercise of the right to take a sick leave reasonable estimates for the parameters involved in determines SP. Using the proposed ZIB approach one both the zero inflated and non zero inflated processes. could describe, in a single analysis, both phenomena: As expected, the credibility intervals length diminishes first, which factors are associated to the exposure to with sample size while their coverage grows. presenteeism (to be “sick”, factors related to health sta- tus), and after that, which factors increase the probabil- Conclusions ity of being presenteeist among the exposed workers. The proposed method is a reliable alternative for the The simulation study shows that, even with relatively analysis of zero inflated di- chotomous outcomes, as small sample sizes the model is capable of producing shown by the simulation study, and can be very useful Moriña et al. BMC Medical Research Methodology (2021) 21:277 Page 8 of 9 in situations when there are two potential and indistin- Barcelona (UAB), 08193 Cerdanyola del Vallès, Spain. 3 Departament de Matemàtiques, Universitat Autònoma de Barcelona (UAB), Edifici C, Campus de guishable sources of zeros. If there are covariates to be Bellaterra, 08193 Cerdanyola del Vall̀es, Spain. 4 Universitat Autònoma de Barce- included in the model, the method is able to use them in lona (UAB), Cerdanyola del Vallès, Spain. 5 Unitat de Bioestadística, Facultat de order to identify the subpopulation at risk, and the Bayes- Medicina, Universitat Autònoma de Barcelona, Cerdanyola del Vallès, Spain. ian strategy assures that the two sources of zeros may be Received: 7 June 2021 Accepted: 7 October 2021 detected even when there are no covariates by utilizing different priors for the probability of success of each Ber- noulli variable. The proposed model has been compiled in the form of the bayesZIB R package [15], so it is pub- References licly available for any researcher facing this issue. 1. Johns G. Presenteeism in the workplace: a review and research agenda. J Organ Behav. 2010;31(4):519–42. https://​doi.​org/​10.​1002/​job.​630. 2. Cameron AC, Trivedi PK. Regression analysis of count data. New York: Cam- Abbreviations bridge University Press; 1998. https://​doi.​org/​10.​1017/​CBO97​80511​814365. CI: Confidence interval; CrI: Credible interval; EM: Expectation-Maximization; 3. Lambert D. Zero-inflated poisson regression, with an application to SP: Sickness presenteeism; ZIB: Zero inflated Bernoulli; ZINB: Zero inflated defects in manufacturing. Technometrics. 1992;34(1):1. https://​doi.​org/​10.​ negative binomial; ZIP: Zero inflated Poisson. 2307/​12695​47. 4. Abiodun G, Makinde O, Adeola A, Njabo K, Witbooi P, Djidjou-Demasse R, et al. A dynamical and zero-inflated negative binomial regression Supplementary Information modelling of malaria incidence in Limpopo Province, South Africa. Int The online version contains supplementary material available at https://​doi.​ J Environ Res Public Health. 2019;16(11). https://​doi.​org/​10.​3390/​IJERP​ org/​10.​1186/​s12874-​021-​01427-2. H1611​2000. 5. Zhao H, Pan Y, Wang C, Guo Y, Yao N, Wang H, et al. The effects of metal exposures on Charlson comorbidity index using zero-inflated negative Additional file 1. Additional tables and code. Additional tables obtained binomial regression model: NHANES 2011-2016. Biol Trace Elem Res. from the simulation study and R code to reproduce the analyses. 2021;199(6):2104–11. https://​doi.​org/​10.​1007/​S12011-​020-​02331-4. 6. Perumean-Chaney SE, Morgan C, McDowall D, Aban I. Zero-inflated and overdispersed: what’s one to do? J Stat Comput. 2013;83(9):1671–83. Acknowledgements https://​doi.​org/​10.​1080/​00949​655.​2012.​668550 http://​dx.​doi.​org.​sire.​ub.​ Not applicable. edu/​10.​1080/​00949​655.​2012.​668550. 7. Paulo Favero L, de Freitas Souza R, Belfiore P, Luiz Corrêa H, Haddad MF, Authors’ contributions Paulo L, et al. Count data regression analysis: concepts, overdispersion DM, PP and AN developed the statistical model. DM and PP derived the detection, zero-inflation identification, and applications with R. Pract described properties, and DM implemented the model in R software and Assess Res Eval. 2021;26. https://​doi.​org/​10.​7275/​44nn-​cj68. conducted the analyses. AN provided the context on occupational health and 8. Vuong QH. Likelihood ratio tests for model selection and non-nested sickness presenteeism. All authors have read and approved the manuscript. hypotheses. Econometrica. 1989;57(2):307. https://​doi.​org/​10.​2307/​19125​57. 9. Cameron AC, Trivedi PK. Regression-based tests for overdispersion in the Funding Poisson model. J Econ. 1990;46(3):347–64. https://​doi.​org/​10.​1016/​0304-​ This work was partially supported by grant RTI2018–096072-B-I00 from the 4076(90)​90014-K. Spanish Ministry of Science and Innovation and by the Spanish State Research 10. Desmarais BA, Harden JJ. Testing for zero inflation in count models: bias Agency, through the Severo Ochoa and María de Maeztu Program for Centers correction for the Vuong test. Stata J. 2013;13(4):810–35. https://​doi.​org/​ and Units of Excellence in R&D (CEX2020–001084-M). The funding sources had 10.​1177/​15368​67X13​01300​408 https://​doi-​org.​sire.​ub.​edu/​10.​1177/​15368​ no role in the data collection, analysis or interpretation of the results. 67X13​01300​408. 11. Barbu A, Wu T, Wu YN. Learning mixtures of bernoulli templates by two- Availability of data and materials round EM with performance guarantee. Electron J Stat. 2015;8:3004–30. The dataset analysed during the current study is available in the GitHub https://​doi.​org/​10.​1214/​14-​EJS981.​1305.​0319v6. repository, https://​github.​com/​dmori​nya/​Bayes​ZIB. R codes used to analyse 12. Diop A, Diop A, Dupuy JF. Simulation-based inference in a zero- the real data and to generate the data used in the simulation study are avail- inflated Bernoulli regression model. Commun Stats: Simul Comput. able as supplementary material. 2016;45(10):3597–614. https://​doi.​org/​10.​1080/​03610​918.​2014.​950743. 13. Sun Z, Rosen O, Sampson AR. Multivariate Bernoulli mixture models with application to postmortem tissue studies in schizophrenia. Biometrics. Declarations 2007;63(3):901–9. https://​doi.​org/​10.​1111/j.​1541-​0420.​2007.​00762.x. 14. R Core Team: R: A Language and environment for statistical computing. R Ethics approval and consent to participate Foundation for statistical computing, Vienna, Austria (2021). R Foundation All procedures were approved by the Ethics Committee on Animal and for statistical computing. https://​www.R-​proje​ct.​org/. Human Experimentation of the Autonomous University of Barcelona 15. Moriña Soler D, Puig P, Navarro A. bayesZIB: Bayesian zero-inflated Ber- (CEEAH/3445). Participation in the study was voluntary and informed consent noulli regression model. In: R package version 0.0.2; 2021. was obtained from all subjects. 16. Zeileis A, Kleiber C, Jackman S. Regression models for count data in R. J Stat Softw. 2008;27(8)1–25. Consent for publication 17. Eurofound. EWCS 2015 – Methodology [Internet]. https://​www.​eurof​ Not applicable. ound.​europa.​eu/​surve​ys/​europ​ean-​worki​ng-​condi​tions-​surve​ys/​ sixth-​europ​ean-​worki​ng-​condi​tions-​survey-​2015/​ewcs-​2015-​metho​dol- Competing interests ogy. Accessed 10 Sep 2021. The authors declare no competing interests. 18. Eurostat. EU labour force survey - methodology [Internet]. http://​ec.​ europa.​eu/​euros​tat/​stati​stics-​expla​ined/​index.​php/​EU_​labour_​force_​ Author details surve​y_-_​metho​dology. Accessed 10 Sep 2021. 1 Department of Econometrics, Statistics and Applied Econometrics, Risk- 19. Navarro A, Salas-Nicás S, Moncada S, Llorens C, Molinero-Ruiz E. center‑IREA, Universitat de Barcelona (UB), Avinguda Diagonal 690, 08034 Bar- Prevalence, associated factors and reasons for sickness presenteeism: celona, Spain. 2 Centre de Recerca Matemàtica, Universitat Autònoma de a cross-sectional nationally representative study of salaried workers in Moriña et al. BMC Medical Research Methodology (2021) 21:277 Page 9 of 9 Spain, 2016. BMJ Open. 2018;8(7):021212. https://​doi.​org/​10.​1136/​bmjop​ en-​2017-​021212. 20. Vives A, Amable M, Ferrer M, Moncada S, Llorens C, Muntaner C, et al. Employment precariousness and poor mental health: evidence from Spain on a new social determinant of health. J Environ Public Health. 2013;2013:978656. https://​doi.​org/​10.​1155/​2013/​978656. 21. Navarro A, Salas‐Nicás S, Llorens C, Moncada S, Molinero‐Ruíz E, Moriña D. Sickness presenteeism: are we sure about what we are studying? A research based on a literature review and an empirical illustration. Am J Ind Med. 2019;62(7). https://​doi.​org/​10.​1002/​ajim.​22982. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in pub- lished maps and institutional affiliations. Ready to submit your research ? Choose BMC and benefit from: • fast, convenient online submission • thorough peer review by experienced researchers in your field • rapid publication on acceptance • support for research data, including large and complex data types • gold Open Access which fosters wider collaboration and increased citations • maximum visibility for your research: over 100M website views per year At BMC, research is always in progress. Learn more biomedcentral.com/submissions