Stochastic variable selection strategies for zero-inflated models

Abstract

When count data exhibit excess zero, that is more zero counts than a simpler parametric distribution can model, the zero-inflated Poisson (ZIP) or zero-inflated negative binomial (ZINB) models are often used. Variable selection for these models is even more challenging than for other regression situations because the availability of p covariates implies 4^p possible models. We adapt to zero-inflated models an approach for variable selection that avoids the screening of all possible models. This approach is based on a stochastic search through the space of all possible models, which generates a chain of interesting models. As an additional novelty, we propose three ways of extracting information from this rich chain and we compare them in two simulation studies, where we also contrast our approach with regularization (penalized) techniques available in the literature. The analysis of a typical dataset that has motivated our research is also presented, before concluding with some recommendations.

Keywords

excess zeros ZI model Hurdle model variable selection stochastic search

1 Introduction

Excess zeros in count data is a situation where we observe more zeros than expected under the assumption that the response variable follows a count distribution such as the Poisson, negative binomial or even Poisson log-normal. Count data with excess zeros are often encountered in health or biomedical sciences, but also in other fields. A non-exhaustive list of typical examples include modelling of dental caries (Mwalili et al. 2008; Preisser et al. 2012), health care utilization (Deb and Holmes 2000), sexual risk behaviour (Yu et al. 2013), highway safety (Lord et al. 2007) and abundance in marine species conservation (Wenger and Freeman 2008). In Section 4, we analyze a popular dataset from Riphahn et al. (2003) on the demand of health care, measured as the number of visits to the doctor in the last three months before the survey, which represents the type of data that have motivated our work. The proportion of zero responses is $41.2 %$ .

For independent data, the zero-inflated regression model (Lambert 1992) and zero-altered or hurdle model (Mullahy 1986) are both commonly used. The zero-inflated (ZI) model is a mixture model which add a point mass distribution on top of a count distribution, whereas the hurdle model describes the zero counts and the positive counts separately. In both models, covariates are used parametrically to describe the two components of the model. Nonparametric extensions to the hurdle model have been proposed by Barry and Welsh (2002), whereas parametric versions with random effects are a current research topic (e.g., Cantoni et al. 2017). For a review of these (and related) models and their properties, we refer to Ridout et al. (1998) and Zuur et al. (2009).

Variable selection is important in applications because it allows to explain data in the simplest way for better interpretability (e.g., efficient identification of risk factors). It is also a mean for cost management, if the model is next used for prediction. From the statistical point of view, parsimonious models have to be preferred for better prediction accuracy and better interpretation. In addition, collinearity-related problems are mitigated when there are fewer variables involved. We focus on the more challenging ZI model even though our methodology can be applied to the hurdle model as well: for $p$ available covariates, there are as many as $4^{p}$ possible different ZI models, whereas this number is $2 \times 2^{p}$ for the hurdle model, due to the orthogonality of its two parts, which can be fitted separately.

The ZI model being a regression model fitted by maximum likelihood, all of the classical techniques for variable selection based on the likelihood are a priori available. These include several kinds of information criteria, stepwise and subset selection. We refer to Burnham and Anderson (2002) as a general reference on model selection. The mentioned selection methods either imply exploring all possible models, which is unfeasible in presence of a large number of covariates, or require selecting a prior set of promising covariates (or models), which make them prone to fail in identifying the correct model which has generated the data. This is one of the reasons why stochastic random search procedures are used (e.g., George and McCulloch 1993; Qian and Field 2002; Cantoni et al. 2007; Jochmann 2013). They have the advantage to avoid an exhaustive search. Markov Chain Monte Carlo (MCMC) techniques are used to screen potential models and ensure efficiency in selecting the set of candidate models. Alternatives that do not require fitting of all possible models are regularization techniques which simultaneously fit the model and perform variable selection; for example, Buu et al. (2011), Wang et al. (2014a), Tang et al. (2014) and Wang et al. (2015) for ZIP and/or ZINB models.

We extend to ZI models the approach by Qian and Field (2002), which combines the use of an MCMC procedure with a criterion for model assessment. We will use an information criterion as in Qian and Field (2002), but also an estimated prediction error based on cross-validation (CV). Within this framework where the algorithm produces an MCMC chain, we pay particular attention to how we can take advantage of the rich output produced and suggest to look at three new concepts: (a) the consensus model, which selects the variables whose marginal frequencies of appearance are above a certain threshold, (b) the minimum criterion, which retains the model with minimum value of the criterion used for model assessment and (c) the model vote, which selects the model which has been most visited.

We contrast our proposal to recent regularization techniques. The simulation show that our proposal is effective in identifying the variables generating the data and either outperforms or is comparable to its competitors. The benefits of our technique are also evident in our real data example, where the number of possible models is $2^{54}$ .

The manuscript has the following structure. In Section 2, we introduce the ZI model and our proposal for stochastic variable selection. Section 3 presents two simulations settings, where two versions of our approach are contrasted with existing alternatives. A real data analysis is presented in Section 4 and the article ends with a concluding section.

2 Methodology

The ZI model assumes that for a random variable $Y_{i}$ , we either observe a zero with probability $π_{i}$ or we observe a realization of a count random variable with mean $μ_{i}$ , with probability $1 - π_{i}$ . More formally, for independent individuals $i = 1, \dots, n$ , we have

P (Y_{i} = y_{i}) = \{\begin{matrix} π_{i} + {1 - π_{i}} P_{F_{μ_{i}, θ}} (Y_{i} = 0) & y_{i} = 0 \\ {1 - π_{i}} P_{F_{μ_{i}, θ}} (Y_{i} = y_{i}) & y_{i} = 1, 2, \dots, \end{matrix}

(2.1)

where

F_{μ_{i}, θ}

is a count distribution with expectation

μ_{i}

and possibly additional (nuisance) parameters

θ

. Supplementary Section 1 gives the explicit expressions of model given in equation (2.1) when

F_{μ_{i}, θ}

is the Poisson or the negative binomial distribution.

Let $x_{i}$ and $z_{i}$ be two sets of covariates measured for subject $i$ (including an intercept term), with possibly some regressors in common. The parameters $π_{i}$ and $μ_{i}$ can be modelled in terms of $x_{i}$ and $z_{i}$ by different link functions. The logit link is a common choice for $π_{i} = π (x_{i})$ and the logarithmic link is the most popular choice for $μ_{i} = μ (z_{i})$ :

log \{\frac{π (x_{i})}{1 - π (x_{i})}\} = x_{i}^{t} β and log {μ (z_{i})} = z_{i}^{t} γ,

(2.2)

where

β

and

γ

are parameters, with

dim (β) = p + 1

and

dim (γ) = q + 1

We define $Y = (Y_{1}, \dots, Y_{n})^{t}$ , $X$ and $Z$ of dimension $n \times (p + 1)$ and $n \times (q + 1)$ , the matrices such that row $i$ contains $x_{i}^{t}$ and $z_{i}^{t}$ , respectively. We also define a model by a vector $υ = (υ_{1}, \dots, υ_{p + q})$ , where the first $p$ components represent the covariates in $x_{i}$ and the following $q$ components represent the covariates in $z_{i}$ . We code $υ_{i} = 1$ if the corresponding variable is included in the model, and $0$ otherwise.

The parameter estimates $\hat{β}$ and $\hat{γ}$ of the ZI model given in equation (2.1) with equation (2.2) are obtained by maximizing the log-likelihood $ℓ (β, γ; Y, X, Z) = \sum_{i = 1}^{n} log P (Y_{i} = y_{i})$ . These estimates inherit the properties of maximum likelihood estimators, in particular, the asymptotic normality with asymptotic variance equal to the inverse of the Fisher information matrix.

Stochastic search procedures consist of moving randomly through the model space, containing all possible models, and draw a subset of promising models. The advantage of such procedures is that they avoid selecting a prior set of candidate models based on potential misleading assumptions. As Qian and Field (2002) put it, the key idea is ‘to convert the model selection procedure into a problem of random sample generation from a finite population’. To address the main concern that the procedure moves well through the model space and guarantees to select the optimal model or a set of models very close to it, MCMC techniques such as the Gibbs sampler or the Metropolis—Hastings algorithm are generally implemented. The idea is to generate from the finite population of all possible models a ‘sample of candidate models which contain the optimal model $α^{0}$ with the highest possible probability’ Qian and Field (2002). We define the probability distribution of model $υ$ on the set $V$ of all possible models as $P (υ) = B exp {- S (υ | Y, X, Z)}$ , where $B = {[\sum_{v \in V} exp {- S (υ | Y, X, Z)}]}^{- 1}$ is a normalizing constant and $S (υ | Y, X, Z)$ is a utility measure, see Section 2.2. The term $B$ is usually non-explicitly computable, but one can nevertheless generate a sample of ‘candidate’ models from $P (υ)$ . We expect $P (υ)$ to become the stationary distribution of the chain, after a certain number of initial steps called the burn-in period, ensuring an efficient move through the model space. A consistent sample of promising models to be considered for the ultimate variable selection is obtained. We refer to Roberts and Rosenthal (1998) as a general reference about MCMC techniques.

2.1 The Metropolis–Hastings algorithm

We will use the Metropolis–Hastings (MH) algorithm (Metropolis et al. 1953; Hastings 1970), which can be summarized with the three following steps:

Choose an initial model $υ^{(0)}$ .

At iteration $s$ , given $υ^{(s - 1)}$ , generate a candidate $\tilde{υ}$ from a set $Ψ_{υ^{(s - 1)}}$ of neighbouring models of $υ^{(s - 1)}$ , using an operating transition kernel $q (υ υ^{(s - 1)})$ for $υ$ . Then draw $u$ from $U \sim$ Uniform(0,1) and compare it with the threshold

\begin{matrix} r (υ^{(s - 1)}, \tilde{υ}) = min [{\frac{P (\tilde{υ})}{P (υ^{(s - 1)})}}^{k} \times \frac{q (υ^{(s - 1)} \tilde{υ})}{q (\tilde{υ} υ^{(s - 1)})} 1] = \\ min [exp \{k (S (υ^{(s - 1)} | Y, X, Z) - S (\tilde{υ} | Y, X, Z))\} \frac{q (υ^{(s - 1)} \tilde{υ})}{q (\tilde{υ} υ^{(s - 1)})} 1], \end{matrix}

(2.3)

where

k

is a scaling parameter varying according to the form of

P (υ)

. Finally, set

υ^{(s)} = \tilde{υ}

u \leq r (υ^{(s - 1)}, \tilde{υ})

and

υ^{(s)} = υ^{(s - 1)}

otherwise.

Repeat the procedure for $s = 1, \dots, S$ , returning the model sequence ${υ^{(0)}, \dots, υ^{(S)}}$ .

In Step 2, a candidate $\tilde{υ}$ is always retained (accepted) if $r (υ^{(s - 1)}, \tilde{υ}) = 1$ . If $r (υ^{(s - 1)}, \tilde{υ}) < 1$ , a candidate $\tilde{υ}$ is accepted with probability $P (U < r) = r$ .

We recommend to choose the initial model $υ^{(0)}$ based on an educated guess (e.g., using a preliminary selection method), so that a more interesting subspace of $V$ is visited.

The set of neighbouring models $Ψ_{υ^{(s - 1)}}$ is defined as a collection of models which differ by a fixed number of covariates with respect to model $υ^{(s - 1)}$ . The amount of differing variables impacts the jumps of the MH algorithm. If $Ψ_{υ^{(s - 1)}}$ equals the set of all possible models, the process is unstable and can hardly converge to the stationary distribution of the chain. With a fixed amount of differing variables, the process is stable. The larger the amount of differing variables, the further the model visited can be from the previous one. Our experience is that a small amount of differing variables (e.g., one variable) is preferable since as soon as we will visit a model close to the optimal one, we will stay close to it and visit rapidly a set of interesting models.

The operating transition kernel $q (υ υ^{(s - 1)})$ is the probability of moving from a model to another. We define it as a function of the $p$ -values of a $z$ -test on the coefficients of each covariate in a full (i.e., including all the available covariates) ZI model. Details of the technicalities are given in Supplementary Section 2, where we exemplify the computation of the kernel. In situations where it is not appropriate to use the $p$ -values from a full model (e.g., large number of covariates and multicollinearity), informative alternatives include the use of the $p$ -value of a marginal (univariate) model, measures of effect size (Cohen's measures, $η^{2}$ , etc.) or VIF (for models for which these quantities are available). In the worst case of absence of relevant information, equal probabilities for each variable can be used.

The scaling parameter $k$ controls the jumps of the MH algorithm through the model space. It depends on the transition kernel $q$ and the utility criterion $S (υ | Y, X, Z)$ , and has to be calibrated empirically so that the probability of accepting a candidate model (also called the acceptance rate) reaches a satisfactory level. Recommendations of Gelman et al. (1996) define a good acceptance rate as approximately $0.5$ for the case of one-dimensional parameters and $0.25$ for the case of multivariate parameters.

The length of the Markov chain $S$ represents the number of submodels to be visited. If $S$ should be moderate to save computational time compared to an exhaustive search, it still has to be sufficiently large to attain stationarity. Qian and Field (2002) use a length of 10 000 and infer that an initial burn-in period of 1 000 might be enough to obtain a stationary process. Cantoni et al. (2007) use a chain of size 5 000. Depending on the use of the resulting chain, stationarity may not be required (see Section 2.3).

The consistency properties of this stochastic search approach have been derived in Qian and Field (2002) in a general setting. They proved that, under some ergodic conditions, the probability that the ‘true’ model appears at least once in the chain tends to one with the length of the chain. With high probability, the ‘true’ model will appear more frequently and early in the chain, and therefore the model with the highest marginal frequency and the one with the smallest utility measure are consistent and converge to the true model almost surely.

2.2 Alternatives for the utility measure

The utility measure $S (υ | Y, X, Z)$ can take several forms. One approach is to use an information criterion like the Akaike's information criterion (AIC; Akaike 1973) $AIC (υ) = - 2 ℓ ({\hat{β}}_{υ}, {\hat{γ}}_{υ}; Y, X, Z) + 2 (p_{υ} + q_{υ} + 2)$ , for a ZI model $υ$ of dimension $p_{υ} + q_{υ}$ , where ${\hat{β}}_{υ}$ and ${\hat{γ}}_{υ}$ are the maximum likelihood estimates for model $υ$ . When implementing $AIC (υ)$ in the MH algorithm, the threshold $r (υ^{(s - 1)}, \tilde{υ})$ in (2.3) becomes

r (υ^{(s - 1)}, \tilde{υ}) = min (exp [k {AIC (υ^{(s - 1)}) - AIC (\tilde{υ})}] \times \frac{q (υ^{(s - 1)} \tilde{υ})}{q (\tilde{υ} υ^{(s - 1)})} 1) .

The Mallow's $C_{p}$ (Mallows 1973), the Bayesian information criterion (BIC; Schwarz 1978), the stochastic complexity (SSC; Rissanen 1996) or the focused information criterion (FIC; Claeskens and Hjort 2008) could also be implemented.

A second approach is to use an estimate of prediction error computed with CV, as in Cantoni et al. (2007). This alternative can also be used in situations where the likelihood is not available. Besides an estimate of the prediction error, one can also obtain an estimate of its standard error and use this information to select a subset of good models in the spirit of the one-standard-error rule (Friedman et al. 2001).

Following Shao (1993), we define $M$ random splits which divide the data into a construction sample $C$ of size $n_{C}$ and a validation sample $T$ of size $n_{T} = n - n_{C}$ . Step 2 of the MH algorithm of Section 2.1 is modified as follows. For a candidate model $\tilde{υ}$ :

a. For each split $m = 1, \dots, M$ , estimate the parameters of the model $\tilde{υ}$ with the construction sample, predict $\hat{y_{i}}$ ( $i \in T 2 pt$ ) for the validation sample and compute the prediction error for this split with the loss function $(n_{T})^{- 1} \sum_{i \in T 2 pt} (y_{i} - \hat{y_{i}})^{2}$ .

b. Compute the average prediction error $PE (\tilde{υ})$ and its empirical standard error $σ {PE (\tilde{υ})}$ across splits.

c. Draw $u$ from $U \sim$ Uniform(0,1) and compare it with

r (υ^{(s - 1)}, \tilde{υ}) = min (exp [k {PE (υ^{s - 1}) - PE (\tilde{υ})}] \times \frac{q (υ^{(s - 1)} \tilde{υ})}{q (\tilde{υ} υ^{(s - 1)})}, 1),

Set,

υ^{(s)} = \tilde{υ}

u \leq r (υ^{(s - 1)}, \tilde{υ})

and

υ^{(s)} = υ^{(s - 1)}

otherwise.

We recommend to use $M = 50$ splits of the data which should be sufficiently large to ensure a valid estimate of the prediction error. The choice of $n_{C}$ depends on the model. Following the recommendation of Shao (1993) for linear models, one could use $n_{C} = n^{3 / 4}$ . For excess zeros models where the proportion of zeros is very large, the choice of $n_{C}$ has to be larger (e.g., $n_{C} = n / 2$ ) to guarantee that both zero and non-zero responses are present in the construction sample.

We can reparametrize $k = cc / σ {PE (υ^{(0)})}$ , yielding a a combination of the jumping tuning constant $cc$ and the standard error of the prediction error computed for the initial model $υ^{(0)}$ . The constant $cc$ can be interpreted as the probability of accepting a new model whose average prediction error is within $σ {PE (υ^{(0)})}$ of the previous one.

The choice of the utility measure can also be driven by the computational burden implied. With respect to the use of an information criterion, the CV approach involves, for each candidate model $\tilde{v}$ , the refitting and prediction of $M$ splits of the data. As a counterpart, it allows to obtain an estimator of $σ (PE)$ , to use, for example, for a one-sigma-rule. Roughly speaking, CV is $M$ times more computationally intensive than the AIC.

2.3 How to summarize results

We describe three ways of summarizing the output of the MCMC chain obtained with our stochastic search. This rich output is a by-product of the procedure and is an important feature from which we can take advantage. The performance of these alternatives will be compared in our simulations.

The consensus model: The model minimizing the utility measure, even if a good one, is not necessarily the true one generating the data. If we are interested in a single best model, we propose to define a consensus model identified based on the marginal frequency of appearance of each variable within the MCMC chain: it includes the variables that appears in a proportion equal or above a certain threshold, named the inclusion level. The higher the inclusion level, the smaller the resulting consensus model. In theory, the computation of the marginal frequencies has to be done after the removal of the initial burn-in period in order not to bias the results. Our practical experience shows that results do not differ much if the entire chain is used instead. This consensus model is close in spirit to the median probability model of Barbieri and Berger (2004).

The minimum criterion: The minimum criterion technique consists of defining the optimal model as the one which optimizes the utility criterion across the chain. The resulting best model is the optimal one among the visited models in the chain, but it is possibly not the model with optimal criterion among all possible models. As the aim of the minimum criterion technique is to find the model with smallest criterion, the removal of a burn-in period is not required and we can use the information of the entire chain.

Model votes: The model votes technique amounts to select the model which has been visited the most in the Markov chain. We consider each visit as a vote and we retain the most ‘popular’ model, that is, the model receiving the largest number of votes. The count of the votes has to be done on the stationary part of the process, after the removal of the burn-in period.

3 Simulation study

We assess the performance of our proposed approach and contrast it to competing alternatives introduced in Supplementary Section 3.1. In Study 1, we consider a moderated-sized design ( $p = q = 10$ , therefore $2^{20}$ possible models) for the ZIP model, whereas in Study 2, we consider a larger dimension setting ( $p = q = 20$ , therefore $2^{40}$ possible models) for the ZINB model. Our approach is implemented in R (R Development Core Team 2012), with the parameter estimates of the ZI model obtained from the zeroinfl function of package pscl (Jackman 2011). Dr A. Buu kindly provided her R code (Buu et al. 2011) and we used the R package mpath (Wang et al. 2014b) to implement the approach of Wang et al.(2014a)) and Wang et al.(2015). In both penalized approaches, we use the default parameter setting.

3.1 Study 1

We consider the same set of potential covariates for both parts of a ZIP model, so that $X = Z$ and $p = q = 10$ , but with two scenarios: one with independent covariates and one with a dependence structure. For the independent dataset, we generate independently $p$ variables for $X$ with $x_{1}$ - $x_{5}$ $\to$ N(0,1), $x_{6}$ and $x_{7}$ $\to$ Bernoulli(0.5) and $x_{8}, x_{9}$ and $x_{10}$ $\to$ Uniform(0,1). For the correlated dataset, we generate independently $x_{1}$ – $x_{8}$ according to the same distributions used for the independent dataset and introduce a correlation between $x_{9}$ and $x_{10}$ using a ‘compound symmetry’ design, according to the following procedure. We generate three i.i.d. variables $W_{1}, W_{2}$ and $U$ from a Uniform(0,1) and set for $t = 0.9$ that $x_{9} = (W_{1} + tU) / (1 + t)$ and $x_{10} = (W_{2} + tU) / (1 + t)$ . The correlation between $x_{9}$ and $x_{10}$ is approximately $Corr (x_{9}, x_{10}) = t^{2} / (1 + t^{2}) = 0.45$ . The nature of the variables has been chosen to represent the variety of situations one can encounter in practice.

For $i = 1, \dots, n$ ( $n = 300$ and $n = 600$ ), our true generating model is defined by

\begin{matrix} logit {π (x_{i})} & = & 0.2 + 0.2 x_{i 3} - 1 x_{i 7} - 1 x_{i 10}, \\ log {λ (z_{i})} & = & 1 + 0.5 z_{i 1} + 0.2 z_{i 6} - 1 z_{i 10} . \end{matrix}

(3.1)

The level of zero-inflation varies between

30 %

and

32 %

across the four situations (

n = 300

and

n = 600

, independent and correlated dataset). For each setting, 300 replications are considered.

The coefficients in model (3.1) are defined such that the $p$ -values obtained from a $z$ -test when we fit the full model (i.e., including all the available variables) give an interesting range of significance levels and henceforth an interesting range of signal strength. In Supplementary Table 1, we summarize for each covariate included in (3.1) the median of $p$ -values across the simulated samples. Note that the effects are in general stronger for $n = 600$ . By including weak and strong effects in the model generating the data, we expect to evaluate the sensitivity of each algorithm in identifying the true explanatory variables with respect to their signals.

Our preliminary calibration has found that $k = 0.25$ gives an average acceptance rate around $30 %$ when AIC is used as utility measure in our simulation setting. For the CV method, we found that an appropriate calibration is $cc = 20$ for an average acceptance rate just below $30 %$ . We set $M = 50$ with the same splits used throughout the MH algorithm, and $n_{C} = n / 2$ . The Markov chains considered are of length $S = 10 000$ for each dataset, with a burn-in period of $5 000$ . As initial model $υ^{(0)}$ , we chose to start with the full model including all the covariates in $X$ and $Z$ .

We discuss in detail the results for $n = 300$ in Sections 3.1.1–3.1.5. The results for $n = 600$ are discussed in Supplementary Section 3.4. A general summary of the simulation results is provided in Section 3.1.6.

Table 1:

Summary of the marginal frequency of appearance computed over the 10 000 visited models of the Markov chain, $n = 300$ . The variables included in the true data-generating model are highlighted in bold

	Binary part
	Independent data						Correlated data
	CV prediction error			AIC			CV prediction error			AIC
	1st Q	Median	3rd Q	1st Q	Median	3rd Q	1st Q	Median	3rd Q	1st Q	Median	3rd Q
$x_{1}$	0.004	0.078	0.537	0.097	0.152	0.275	0.009	0.112	0.515	0.102	0.153	0.268
x ₂	0.005	0.048	0.312	0.097	0.145	0.268	0.012	0.075	0.402	0.104	0.163	0.325
$x_{3}$	0.022	0.392	0.935	0.134	0.307	0.746	0.010	0.160	0.837	0.156	0.319	0.668
$x_{4}$	0.004	0.040	0.349	0.096	0.152	0.250	0.015	0.079	0.426	0.103	0.161	0.328
$x_{5}$	0.008	0.096	0.457	0.106	0.162	0.304	0.010	0.061	0.418	0.101	0.162	0.294
$x_{6}$	0.027	0.151	0.526	0.099	0.167	0.318	0.037	0.175	0.612	0.099	0.176	0.365
$x_{7}$	0.940	1.000	1.000	0.908	0.992	0.999	0.888	0.999	1.000	0.840	0.989	0.998
$x_{8}$	0.003	0.026	0.226	0.098	0.151	0.283	0.006	0.031	0.344	0.104	0.156	0.262
$x_{9}$	0.006	0.031	0.224	0.098	0.138	0.239	0.016	0.108	0.441	0.107	0.164	0.282
$x_{10}$	0.113	0.772	0.991	0.250	0.642	0.907	0.077	0.479	0.888	0.144	0.335	0.755
	Count part
	Independent data						Correlated data
	CV prediction error			AIC			CV prediction error			AIC
	1st Q	Median	3rd Q	1st Q	Median	3rd Q	1st Q	Median	3rd Q	1st Q	Median	3rd Q
$z_{1}$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000
$z_{2}$	0.004	0.051	0.399	0.103	0.156	0.275	0.005	0.034	0.310	0.100	0.159	0.279
$z_{3}$	0.012	0.095	0.488	0.097	0.152	0.283	0.004	0.033	0.256	0.094	0.146	0.316
$z_{4}$	0.004	0.039	0.324	0.097	0.149	0.304	0.005	0.046	0.341	0.103	0.170	0.288
$z_{5}$	0.005	0.064	0.385	0.099	0.155	0.294	0.005	0.031	0.200	0.097	0.154	0.260
$z_{6}$	0.111	0.878	0.998	0.264	0.674	0.965	0.114	0.856	0.998	0.282	0.696	0.975
$z_{7}$	0.035	0.212	0.540	0.102	0.152	0.281	0.043	0.179	0.613	0.106	0.188	0.347
$z_{8}$	0.003	0.025	0.256	0.091	0.146	0.239	0.004	0.036	0.227	0.097	0.155	0.271
$z_{9}$	0.007	0.057	0.271	0.097	0.137	0.281	0.006	0.067	0.377	0.106	0.162	0.350
$z_{10}$	1.000	1.000	1.000	1.000	1.000	1.000	0.993	1.000	1.000	0.971	0.998	1.000

3.1.1 Marginal frequency of appearance

The marginal frequency of appearance of each variable gives a sense of whether our procedure visit the most interesting models. Table 1 summarizes the frequencies for $n = 300$ . These are computed for each variable on the full Markov chain of each simulation run (i.e., without removing a burn-in period for the computations). Overall, by looking at the medians, we see that the true variables generating the data (in bold) appear in a higher proportion relatively to the superfluous ones, which gives some prior indications about the performances of the algorithms in identifying the true variables. The variables with weaker effect (i.e., $x_{3}, x_{10}$ and $z_{6}$ ) appear relatively less than the variables with stronger effect (i.e., $x_{7}, z_{1}$ and $z_{10}$ ). The variable $z_{1}$ appears in all cases regardless to the setting, which is not surprising as its signal is very strong (see Supplementary Table 1).

If we compare CV and AIC for the same type of data (either independent or correlated), we observe that CV provides larger median marginal frequencies of appearance for the true variables and smaller ones for the variables not included in the data-generating model, but that the distribution of the marginal frequencies is wider for CV than for AIC (as expressed by the first and third quartile).

Focusing on the correlated setting, we observe that the correlation between the pairs of variables $(x_{9}, x_{10})$ and $(z_{9}, z_{10})$ does not affect the marginal frequency of appearance of $z_{10}$ , but affects the marginal frequency of appearance of $x_{10}$ , most probably because the information is sparser in a binary regression.

Figure 1:

Proportion of simulated samples (over the last 5 000 visited models of each MCMC chain) including the variable in the consensus model. In dark grey are the variables included in the true data-generating model

3.1.2 The consensus model

We compare three different inclusion levels for the consensus model, namely, $30 %$ , $50 %$ and $70 %$ . In Figure 1 (independent dataset) and Supplementary Figure 1 (correlated dataset) the proportion of inclusion of each variable in the consensus model are presented. The analysis of these figures shows that, for the three inclusion levels, the true variables (in dark grey) are included in the consensus model in a large proportion of simulated samples relatively to the other variables. The variables with strong effect (i.e., $x_{7}, z_{1}$ and $z_{10}$ ) are included in the consensus model in a large proportion of simulated samples, even when the inclusion level is $70 %$ . This indicates that our algorithm (both with CV and AIC) performs well in identifying these kind of covariates regardless of the threshold defined for the consensus model. We see that the proportion of simulated samples in which we identified the variables with weak effect (i.e., $x_{3}, x_{10}$ and $z_{6}$ ) tends to decrease as the inclusion level increases, showing that the level of the consensus model impacts the identification of these lower signal covariates. This drop is more important when AIC is used (a consequence of the larger spread of the marginal frequencies, as observed in Table 1), but it is joint with a decrease of the proportion of superfluous variables retained in the model. In Figure 1, we see that the true variables are included in the $30 %$ consensus model for the majority of the simulated samples within each setting. A low inclusion level increases the probability of identifying the entire set of true variables (even the ones with weak effect), even though the probability of inclusion of other (superfluous) variables also increases, leading to a less parsimonious model. Figure 1 shows that a high inclusion level does not prevent a high probability of keeping the variables with a strong effect whilst the probability of keeping superfluous covariates is lower. Thus, an inclusion level equal of $50 %$ seems to be a good compromise between accuracy and simplicity of the model. Focusing on the $50 %$ consensus model, we see that the proportions of inclusion of the true data-generating variables are quite close when we compare CV and AIC. For the binary part of the model, we see that the two versions exhibit better performance for independent data than for correlated data. For the count model, the two methods seem to perform as good for both kinds of data, indicating some robustness in recognizing the true data-generating variables. Globally speaking, we notice that the proportion of inclusion of the superfluous variables are constantly higher with CV than with AIC, showing that the CV procedure is less selective in excluding the superfluous covariates.

3.1.3 Regularizations techniques

We also did the same plots for the penalized approach. We consider three versions of the Buu et al. (2011) approach (BJLT) for each of the LASSO and SCAD penalty, see options 1), 2) and 3) in Supplementary Section 3.1. For the approach of Wang et al.(2014a) and (Wang et al. (2015) (WMW), we consider the three penalties: LASSO, SCAD and MCP. The results are given in Figure 2 (independent data), and in Supplementary Figure 2 (correlated data). For the independent data setting, focusing on BJLT, the different behaviour of the technique depending on the choice of the penalty (either LASSO or SCAD) is apparent at first sight. The LASSO tends to include lots of superfluous variables, in particular with versions 2) and 3), where the tuning parameters depend on the standard error of the coefficients. At the contrary, SCAD is very selective: it captures the stronger signals, but often misses the weakest. There is less difference between the three versions of SCAD than there is between the three versions of LASSO. We also notice that LASSO, and less so SCAD, seem affected by the nature of the covariates: see, for example, the behaviour for $x_{8}$ – $x_{10}$ , which are uniformly distributed, compared to $x_{1}$ – $x_{5}$ , which are normally distributed. If we move on and look at the results for WMW, we observe that the behaviour is globally much more satisfactory than the BJLT approach. For all three penalties, WMW is good in identifying the strong signals. WMW–MCP is less effective in discarding the irrelevant covariates. Altogether, the WMW results with either LASSO or SCAD are comparable with those obtained with our $50 %$ or $70 %$ consensus model.

Figure 2:

Proportion of simulated samples including the variable in the final model. In dark grey are the variables included in the true data-generating model

Comparing Figure 2 with Supplementary Figure 2 we see that, similarly as for the consensus model, the performance of the penalized approaches in recognizing the signal associated with $x_{10}$ is affected by the correlation between $x_{9}$ and $x_{10}$ .

3.1.4 The minimum criterion

In Table 2 are displayed the proportion (across the entire MCMC chain) of inclusion of each variable in the model with minimum criterion. We notice that the minimum criterion technique (either with CV or AIC) seems globally less effective than the $50 %$ consensus model counterpart in identifying the true variables (compare with results in Figure 1 and Supplementary Figure 1). It is particularly the case concerning the variables with a weak effect ( $x_{3}, x_{10}$ and $z_{6}$ ). On the other hand, the minimum criterion technique performs better than the consensus model in excluding the superfluous variables given their low proportion of inclusion, in particular with AIC. The CV version perform slightly better than AIC in identifying the variables with weak effects, but on the other hand, it is less selective in excluding the uninformative variables. When we compare the independent and the correlated datasets, we see that the performance in capturing the signal of $x_{10}$ and $z_{10}$ is lower in the correlated dataset.

Table 2:

Proportion of simulated samples for which each variable is included in the optimal model that minimizes the criterion across the entire MCMC chain ( $n = 300$ )

Independent data						Correlated data
Binary part			Count part			Binary part			Count part
	CV	AIC		CV	AIC		CV	AIC		CV	AIC
$x_{1}$	0.200	0.033	$z_{1}$	1.000	1.000	$x_{1}$	0.203	0.053	$z_{1}$	1.000	1.000
$x_{2}$	0.150	0.043	$z_{2}$	0.163	0.053	$x_{2}$	0.143	0.057	$z_{2}$	0.113	0.040
$x_{3}$	0.317	0.197	$z_{3}$	0.163	0.043	$x_{3}$	0.290	0.193	$z_{3}$	0.137	0.053
$x_{4}$	0.140	0.027	$z_{4}$	0.127	0.040	$x_{4}$	0.127	0.053	$z_{4}$	0.113	0.037
$x_{5}$	0.193	0.030	$z_{5}$	0.167	0.043	$x_{5}$	0.163	0.043	$z_{5}$	0.087	0.037
$x_{6}$	0.187	0.047	$z_{6}$	0.520	0.373	$x_{6}$	0.203	0.070	$z_{6}$	0.510	0.410
$x_{7}$	0.760	0.783	$z_{7}$	0.197	0.043	$x_{7}$	0.750	0.783	$z_{7}$	0.213	0.037
$x_{8}$	0.143	0.050	$z_{8}$	0.110	0.037	$x_{8}$	0.127	0.037	$z_{8}$	0.103	0.043
$x_{9}$	0.120	0.020	$z_{9}$	0.137	0.037	$x_{9}$	0.177	0.043	$z_{9}$	0.157	0.060
$x_{10}$	0.483	0.313	$z_{10}$	0.993	0.997	$x_{10}$	0.327	0.210	$z_{10}$	0.780	0.887

Table 3:

Proportion of simulated samples for which each variable is included in the optimal model with largest number of votes over the last 5 000 visited models of each MCMC chain ( $n = 300$ )

Independent data						Correlated data
Binary part			Count part			Binary part			Count part
	CV	AIC		CV	AIC		CV	AIC		CV	AIC
$x_{1}$	0.277	0.155	$z_{1}$	1.000	1.000	$x_{1}$	0.267	0.170	$z_{1}$	1.000	1.000
$x_{2}$	0.200	0.178	$z_{2}$	0.227	0.174	$x_{2}$	0.210	0.196	$z_{2}$	0.197	0.186
$x_{3}$	0.423	0.414	$z_{3}$	0.260	0.207	$x_{3}$	0.373	0.405	$z_{3}$	0.183	0.199
$x_{4}$	0.213	0.132	$z_{4}$	0.193	0.194	$x_{4}$	0.223	0.196	$z_{4}$	0.190	0.186
$x_{5}$	0.240	0.214	$z_{5}$	0.203	0.171	$x_{5}$	0.223	0.183	$z_{5}$	0.147	0.170
$x_{6}$	0.270	0.174	$z_{6}$	0.623	0.658	$x_{6}$	0.313	0.232	$z_{6}$	0.627	0.663
$x_{7}$	0.843	0.941	$z_{7}$	0.277	0.164	$x_{7}$	0.840	0.938	$z_{7}$	0.320	0.229
$x_{8}$	0.177	0.174	$z_{8}$	0.170	0.155	$x_{8}$	0.197	0.180	$z_{8}$	0.180	0.154
$x_{9}$	0.173	0.151	$z_{9}$	0.170	0.151	$x_{9}$	0.243	0.183	$z_{9}$	0.217	0.212
$x_{10}$	0.627	0.622	$z_{10}$	0.993	1.000	$x_{10}$	0.490	0.441	$z_{10}$	0.863	0.971

3.1.5 Model votes

The proportions of inclusion of each variable in the final model obtained by the model votes technique are summarized in Table 3. These proportions are computed after burn-in removal. We notice that the proportion of inclusion of the true and superfluous variables for the CV version are very close to the inclusion probabilities of the $50 %$ consensus model presented in Figure 1 and Supplementary Figure 1. When AIC is used with the model votes techniques, the inclusion probabilities are as good as those of CV. Here again, the algorithm is less effective in recognizing the variables with weak signal (e.g., $x_{3}$ , $x_{10}$ and $x_{6}$ ) compared to the strong ones: both versions (CV and AIC) seem to identify better the true variables when they are applied to an independent dataset than a correlated one, in particular for the binary part of the model.

Table 4:

Number of times the true model has been correctly identified (Exact), the number of times one or more data-generating variables are missing from the selected model, but with possibly many extra spurious variables (Underfit) and the number of times a largest model containing the true one (Overfit) has been selected ( $n = 300$ )

	Independent data
	Underfit ( $\leq 3$ )	Underfit ( $- 2$ )	Underfit ( $- 1$ )	Exact	Overfit ( $+ 1$ )	Overfit ( $+ 2$ )	Overfit (>3)
CV $- 30 %$	27	95	121	1	3	4	49
CV $- 50 %$	37	119	101	2	7	10	24
CV $- 70 %$	63	119	97	3	8	5	8
AIC $- 30 %$	10	71	159	0	8	21	31
AIC $- 50 %$	38	114	126	3	10	7	2
AIC $- 70 %$	81	117	91	5	3	3	0
min CV	86	114	84	5	3	3	5
min AIC	130	125	44	1	0	0	0
votes CV	38	109	114	2	7	6	24
votes AIC	24	95	151	3	15	9	7
BJLT–LASSO1	124	103	54	0	1	2	16
BJLT–LASSO2	0	20	194	0	0	0	86
BJLT–LASSO3	1	17	151	0	0	0	131
BJLT–SCAD1	149	120	31	0	0	0	0
BJLT–SCAD2	134	103	51	1	3	1	7
BJLT–SCAD3	90	115	77	1	4	3	10
WMW–LASSO	89	123	71	1	2	9	5
WMW–SCAD	97	124	67	3	6	3	0
WMW–MCP	74	120	90	1	6	4	5

If we look at Supplementary Table 2 which summarizes the median number of votes for the best model across the simulated samples, we notice that the most visited model is generally more popular (i.e., gets relatively more votes) in the CV setting than in the AIC setting. On the other hand, the median number of model visited within the last 5 000 models of the chain is much larger when AIC is used than when CV is used. These observations indicate that the AIC method explore more models than the CV method.

Finally, Table 4 (independent dataset) and Supplementary Table 3 (correlated dataset) give a picture of the correct model identification of each technique considered here. The table gives the number of times the true model has been correctly identified (Exact), the number of times one or more data-generating variables are missing from the selected model, but with possibly many extra spurious variables (Underfit) and the number of times a largest model containing the true one (Overfit) has been selected. As a general comment, we can say that, given the mixture of weak and strong signals in the data-generating model, it is very difficult (for each technique) to identify exactly the true model. Also, due to the presence of weaker signals, the general tendency is to miss some variables rather than to overfit. A notable exception are the BJLT-LASSO2 and BJLT-LASSO3 techniques, in particular for the independent data. This is in line with what we have observed in Figure 2.

For the independent data, the consensus model, either based on CV or AIC, very often selects a model where one or two data-generating variables are missing. The minimum criterion misses the signal more heavily than the consensus model, even more so for AIC, whereas the performance of the approach of model votes is comparable to that of the consensus model. As already mentioned, BJLT–LASSO2 and BJLT–LASSO3 either largely overfit (three extra variables) or miss one variable, whereas BJLT–LASSO1 heavily discard variables that should appear in the selected model (two or more variables missing). The three versions of BJLT–SCAD tend to behave like BJLT–LASSO1 (with some exceptions), which means that they miss two or more true variables, making it less appealing than our consensus model approach, for example. WMW shows a similar behaviour to BJLT–SCAD, irrespective of the penalty used.

For the correlated dataset, the conclusions are similar, but with a stronger tendency to underfit. This is particularly so for BJLT–LASSO2 and BJLT–LASSO3.

3.1.6 Summary

We remark that better performance in identifying the signals, in particular the weak ones, is accompanied by the inclusion of more superfluous variables in the final model. This is quite sensible.

Our simulation exercises show that for the consensus model, both versions which make use of either CV or AIC are good in identifying strong signals. CV is more sensitive than AIC to the inclusion level when it comes to the identification of lower signal variables. As a consequence, CV tends to include more superfluous variable than AIC. Globally speaking, CV and AIC exhibit comparable performances regarding the proportion of inclusion of the true data-generating variables. Both versions are affected by the correlation between covariates, but mainly in the binary part of the model. These same comments apply to the model votes technique, which shows a very similar behaviour to the 50% consensus model. The minimum criterion approach is globally less effective than the 50% consensus model, in particular concerning the weak effects. CV seems to behaves slightly better than AIC, though, in identifying the variables with weak effects. Lower performance of the minimum criterion is observed for the correlated dataset, but only for the correlated variables themselves.

Between the penalized techniques, the three versions of WMW clearly outperform BJLT–LASSO. BLJT–SCAD1 and WMW–SCAD are comparable in terms of signal identification (e.g., Figure 2), but WMW–SCAD tend to miss fewer variables (e.g., Table 4).

3.2 Study 2

In this study, we consider one of the scenarios in Wang et al. (2015), namely their example 1. A ZINB model is considered with $20$ covariates for each regression (2.2), so that $p = q = 20$ , and the number of models to compare is $2^{40}$ . For a sample size $n = 300$ , the predictors are randomly drawn from a $N_{20} (0, Σ)$ distribution, where $Σ$ has elements $ρ^{| i - j |}$ , for $i, j = 1, \dots, 20$ , with $ρ = 0.4$ . The parameters are set to $θ = 2$ , $β = (0.30, - 0.48, 0, 0, 0, 0.4, 0, 0, 0, 0, 0.44, 0, 0.44, 0, 0, 0, 0, 0, 0, 0, 0),$ and $γ = (1.10, 0, 0, 0, - 0.36, 0, 0, 0, 0, 0, 0, 0, 0, - 0.32, 0, 0, 0, 0, 0, 0, 0)$ , which yields a level of zero-inflation of $56 %$ on average. We consider 300 replications.

Based on the conclusions of Study 1, we compare here our consensus model based on AIC (with inclusion levels of $30 %$ , $50 %$ and $70 %$ ) and the WMW approach with the three different penalties. For the MH algorithm with AIC, the tuning parameter $k$ is set to $0.15$ , which produces an acceptance rate approaching $30 %$ .

For consistency with the original reference, the results are summarized by looking at the median of the MSE ratio between the identified model and the ‘full’ ZINB model, the sensitivity (proportion of correctly identified number of non-zero coefficients) and the specificity (proportion of identified number of zero coefficients). In Table 5, the median MSE values show that the accuracy of our AIC approach is better in the binary part of the model, but that WMW performs better in the count part. Increasing the inclusion level of the consensus model seems to have the effect of decreasing MSE ratio and sensitivity, but increasing specificity (note, however, that the standard deviation is larger than the observed differences, in particular for the sensitivity). Our approach has better sensitivity for both the binary and the count part of the model, whereas the performance in terms of specificity is better for the WMW approach. The MCP penalty provides the smaller MSE ratios between the WMW versions and also gives good sensitivity. The LASSO penalty exhibits good performance in terms of specificity. Identifying zero and non-zero coefficients is a competing goal, and for all approaches a trade-off takes place.

Table 5:

Median of MSE ratio between the model identified by each technique and the full ZINB model. Mean and standard deviation (between parentheses) of sensitivity and specificity

	MSE ratio	Sensitivity	Specificity
Binary part
AIC $- 30 %$	0.418	0.798 (0.190)	0.752 (0.109)
AIC $- 50 %$	0.395	0.721 (0.212)	0.837 (0.107)
AIC $- 70 %$	0.390	0.666 (0.224)	0.900 (0.105)
WMW–LASSO	0.457	0.583 (0.288)	0.957 (0.092)
WMW–SCAD	0.420	0.641 (0.270)	0.935 (0.094)
WMW–MCP	0.368	0.602 (0.246)	0.955 (0.075)
Count part
AIC $- 30 %$	0.548	0.958 (0.150)	0.668 (0.097)
AIC $- 50 %$	0.474	0.942 (0.176)	0.744 (0.095)
AIC $- 70 %$	0.406	0.927 (0.191)	0.800 (0.093)
WMW–LASSO	0.267	0.850 (0.293)	0.921 (0.094)
WMW–SCAD	0.247	0.898 (0.246)	0.874 (0.122)
WMW–MCP	0.215	0.918 (0.214)	0.901 (0.094)

4 Real data analysis

This dataset contains a sub-sample of 1 812 German male individuals in 1994 used in Riphahn et al.(2003) to study the demand for health care. The dataset is publicly available in the R package zic. A description of the variables is given in Supplementary Table 10. The available covariates include demographic information, education and employment information, health status and type of insurance. Supplementary Figure 7 shows a histogram, truncated at 20, of the number of doctor visits in the last three months. We can observe a large peak at zero followed by a certain shape that possibly looks like a negative binomial distribution. The proportion of zeros represents $41.2 %$ of the distribution which is quite important and gives us some indications about an excess zeros situation. We use the same covariates specification (including interactions) as in Wang et al. (2015), which implies $2^{54}$ possible ZINB models.

Table 6:

Variable selection results for the health care demand dataset. The couple $\cdot / \cdot$ indicates whether a variable has been retained (x) or discarded (–) in the binary part and the count part of the model, respectively

	CV			AIC			WMW
	$50 %$	min	votes	$50 %$	min	votes	LASSO	SCAD	MCP
	binary/count
health	x/x	x/x	x/x	x/x	x/x	x/x	x/x	x/x	x/x
age30TRUE	x/–	x/–	x/–	–/–	–/–	–/–	–/–	–/–	–/–
age35TRUE	x/–	–/–	x/–	–/–	–/–	–/–	–/–	–/–	–/–
age40TRUE	–/–	–/–	–/–	–/–	–/–	–/–	–/–	–/–	–/–
age45TRUE	x/–	x/–	x/–	x/x	–/–	–/x	–/–	–/–	–/–
age50TRUE	x/–	x/–	x/–	–/–	–/–	x/–	x/x	x/–	x/–
age55TRUE	–/–	–/–	x/–	–/–	–/–	–/–	–/x	x/x	–/–
age60TRUE	x/–	–/–	–/–	–/–	–/–	–/–	–/–	–/–	–/–
handicap	–/–	x/–	x/–	–/–	–/–	–/–	/x-	–/x	–/–
hdegree	–/–	–/–	–/–	–/–	–/x	–/–	–/–	–/–	–/–
married	x/x	x/x	x/x	–/–	–/x	x/–	–/–	–/–	–/–
schooling	–/–	–/–	–/–	–/–	–/–	–/–	–/–	–/–	–/–
hhincome	x/–	–/–	x/x	–/–	–/–	–/–	–/–	–/–	–/–
children	x/–	x/–	x/–	x/–	x/–	x/–	x/–	x/–	x/–
self	–/x	–/–	x/x	–/x	–/–	–/x	–/–	–/x	–/–
civil	–/x	–/–	–/–	–/x	–/–	–/x	–/–	–/x	–/–
bluec	–/x	–/–	–/x	–/–	–/–	–/–	–/–	–/–	–/–
employed	–/x	–/–	–/x	–/–	–/–	–/–	–/–	–/–	–/–
public	–/–	x/–	x/–	–/–	–/–	–/–	–/x	–/–	–/–
addon	–/x	x/x	–/x	–/–	–/x	–/–	–/–	–/–	–/–
health:age30TRUE	x/–	x/–	x/–	–/–	–/–	x/–	–/–	–/–	–/–
health:age35TRUE	–/–	–/x	x/x	–/–	–/–	–/–	–/–	–/–	–/–
health:age40TRUE	–/–	–/–	–/x	–/–	–/–	–/–	–/–	–/–	–/–
health:age45TRUE	x/x	x/x	x/x	x/x	x/–	x/x	–/–	–/–	–/–
health:age50TRUE	x/–	x/–	x/x	–/–	–/–	–/–	–/–	–/–	–/–
health:age55TRUE	–/x	–/x	x/x	–/x	–/x	–/x	–/–	–/–	–/–
health:age60TRUE	–/x	–/x	–/x	–/–	–/–	–/–	–/–	–/–	–/–

We apply to this dataset our consensus model approach, the minimum criterion technique, the model votes technique (all of them with both CV ( $k = 20$ ) or AIC ( $c = 0.25$ )) and the WMW approach. Table 6 summarizes the results in terms of variable selection. The retained variables by each technique are indicated by a ‘x’ and the discarded variables by a ‘–’.

Few variables have been selected by all the techniques: health and children in the binary part, and for health in the count part. Some variables have been (almost) consistently discarded. It is the case for age40TRUE, age55TRUE, age60TRUE and their interaction with health, schooling, civil, bluec, employed and addon in the binary part of the model, and for almost all the age factors, schooling and children in the count part of the model. For some other variables, there is less agreement between the methods. If we focus only on our proposal, we remark that CV retains more variables than AIC.

5 Conclusions

We propose a technique for variable selection when the number of covariates is large and the screening of all possible models is unfeasible. The technique is simulation based and comes with several options for using the rich computed information. From our two simulation studies, we have learned that our technique is very effective in identifying the correct variables that have generated the data.

Based on the results of Study 1, we recommend the use of either the $50 %$ consensus model or the model votes technique, which perform similarly. These two approaches perform better than the minimum criterion version of our procedure, probably because they better take advantage of entire Markov Chain information. The performance of the regularization approach depends heavily on the implementation and on the choice for smoothing parameter selection. The implementation of Buu et al. (2011) is disappointing (at least in the setting of our Study 1 design, which include also discrete and binary covariates), but the approach of (Wang et al.(2014a) and Wang et al.(2015) gives better results, which are comparable to our $50 %$ consensus model based on AIC, for instance. In summary, in Study 1, our technique either outperforms or behave similarly than the competitors.

In the second simulation study, where we compare more closely our consensus model to the WMW approach in a larger setting, we observe that specificity is improved at the price of sensitivity. The consensus approach has better sensitivity, whereas the WMW approach has better specificity. Whether one should prefer a method that exhibits better specificity or better sensitivity depends on the final goal of the analysis. For instance, sensitivity is more important in studies where one does not want to miss potential effects. On the other hand, specificity is more important to avoid suggesting nonexistent relationships.

This second simulation study confirms the good performance of the $50 %$ consensus model, which we recommend as the better choice between our three approaches.

Footnotes

Acknowledgments

We would like to thank the Editor, the Associate Editor and the anonymous referees for their encouragements and their careful reading of our manuscript. Their comments have led to an improved version of our article.

References

Akaike

(1973) Information theory and an extension of the maximum likelihood principle. In Petrov, BN and Csaki, F (eds), 2nd International Symposium on Information Theory , Vol. 1, pages 267–81. Budapest: Akademiai Kiado.

Barbieri

and Berger

(2004) Optimal predictive model selection. Annals of Statistics , 32, 870–97.

Barry

and Welsh

(2002) Generalized additive modelling and zero inflated count data. Ecological Modelling , 157, 179–88.

Burnham

and Anderson

(2002) Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach . New York, NY: Springer-Verlag Inc.

Buu

Johnson

and Tan

(2011) New variable selection methods for zero-inflated count data with applications to the substance abuse field. Statistics in Medicine , 30, 2326–40.

Cantoni

Field

Mills Flemming

and Ronchetti

(2007) Longitudinal variable selection by cross-validation in the case of many covariates. Statistics in Medicine , 26, 919–30.

Cantoni

Mills Flemming

and Welsh

(2017) A random-effects hurdle model for predicting bycatch of endangered marine species, Annals of Applied Statistics, to appear.

Claeskens

and Hjort

(2008) Model Selection and Model Averaging , Vol. 330. Cambridge: Cambridge University Press.

Deb

and Holmes

(2000) Estimates of use and costs of behavioural health care: A comparison of standard and finite mixture models. Health Economics , 9, 475–89.

10.

Friedman

Hastie

and Tibshirani

(2001) The Elements of Statistical Learning, Vol. 1. New York: Springer.

11.

Gelman

Roberts

and Gilks

(1996) Efficient metropolis jumping rules. Bayesian Statistics , 5, 599–608.

12.

George

and McCulloch

(1993) Variable selection via Gibbs sampling. Journal of the American Statistical Association , 88, 881–89.

13.

Hastings

(1970) Monte Carlo sampling methods using Markov chains and their applications. Biometrika , 57, 97–109.

14.

Jackman

(2011) pscl: Classes and Methods for R Developed in the Political Science Computational Laboratory . Department of Political Science, Stanford University, Stanford, California. URL https://cran.r-project. R package version 1.04.1.

15.

Jochmann

(2013) What belongs where? Variable selection for zero-inflated count models with an application to the demand for health care. Computational Statistics , 28, 1947–64.

16.

Lambert

(1992) Zero-inflated Poisson regression, with an application to defects in manufacturing. Technometrics , 34, 1–14.

17.

Lord

Washington

and Ivan

(2007) Further notes on the application of zero-inflated models in highway safety. Accident Analysis & Prevention , 39, 53–7.

18.

Mallows

(1973) Some comments on C_p. Technometrics , 15, 661–75.

19.

Metropolis

Rosenbluth

Teller

and Teller

(1953) Equation of state calculations by fast computing machines. The Journal of Chemical Physics , 21, 1087.

20.

Mullahy

(1986) Specification and testing of some modified count data models. Journal of Econometrics , 33, 341–65.

21.

Mwalili

Lesaffre

and Declerck

(2008) The zero-inflated negative binomial regression model with correction for misclassification: An example in caries research. Statistical Methods in Medical Research , 17, 123–39.

22.

Preisser

Stamm

Long

and Kincade

(2012) Review and recommendations for zero-inflated count regression modeling of dental caries indices in epidemiological studies. Caries Research , 46, 413–23.

23.

Qian

and Field

(2002) Using MCMC for logistic regression model selection involving large number of candidate models. In Fang, KT, Hickernell, FJ and Niederreiter, J (eds), Monte Carlo and Quasi-Monte Carlo Methods 2000 , pages 460–74. Berlin, Heidelberg: Springer.

24.

R Development Core Team (2012) R: A language and environment for statistical computing . R Foundation for Statistical Computing, Vienna, Austria. URL \mychar http://[www.R-project.org/](last accessed 04 May 2017).

25.

Ridout

Demétrio

CGB

and Hinde

(1998) Models for count data with many zeros. In Proceedings of the 19th International Biometrics Conference, Cape Town , pages 179–90. URL https://www.kent.ac.uk/smsas/personal/msr/webfiles/zip/ibc_fin.pdf

26.

Riphahn

Wambach

Million

(2003) Incentive effects in the demand for health care: A bivariate panel count data estimation. Journal of Applied Econometrics , 18, 387–405.

27.

Rissanen

(1996) Fisher information and stochastic complexity. Information Theory, IEEE Transactions on , 42, 40–7.

28.

Roberts

and Rosenthal

(1998) Markov-chain Monte Carlo: Some practical implications of theoretical results. Canadian Journal of Statistics , 26, 5–20.

29.

Schwarz

(1978) Estimating the dimension of a model. The Annals of Statistics , 6, 461–64.

30.

Shao

(1993) Linear model selection by cross-validation. Journal of the American Statistical Association , 88, 486–94.

31.

Tang

Xiang

and Zhu

(2014) Risk factor selection in rate making: EM adaptive LASSO for zero-inflated Poisson regression models. Risk Analysis , 34, 1112–27.

32.

Wang

C-Y

Zappitelli

Devarajan

and Parikh

(2014a) EM for regularized zero-inflated regression models with applications to postoperative morbidity after cardiac surgery in children. Statistics in Medicine , 33, 5192–5208. ISSN 1097-0258.

33.

Wang Z, with contributions from Achim Zeileis, Jackman S, Ripley B et al. (2014b) mpath: Regularized linear models . URL http://CRAN.R-project.org/package=mpath. R package version 0.1-17. (last accessed 04 May 2017).

34.

Wang

and Wang

(2015) Variable selection for zero-inflated and overdispersed data with application to health care demand in Germany. Biometrical Journal , 57, 867–84.

35.

Wenger

and Freeman

(2008) Estimating species occurrence, abundance, and detection probability using zero-inflated distributions. Ecology , 89, 2953–59.

36.

Chen

Tang

. (2013) Distribution-free models for longitudinal count responses with overdispersion and structural zeros. Statistics in Medicine , 32, 2390–2405.

37.

Zuur

Ieno

Walker

Saveliev

and Smith

(2009) Mixed Effects Models and Extensions in Ecology with R . New York: Springer.

Stochastic variable selection strategies for zero-inflated models

Abstract

Keywords

1 Introduction

2 Methodology

2.3 How to summarize results

3 Simulation study

3.1 Study 1

Summary of the marginal frequency of appearance computed over the 10 000 visited models of the Markov chain, n = 300 . The variables included in the true data-generating model are highlighted in bold

Figure 1:

Proportion of simulated samples (over the last 5 000 visited models of each MCMC chain) including the variable in the consensus model. In dark grey are the variables included in the true data-generating model

3.1.3 Regularizations techniques

Figure 2:

Proportion of simulated samples including the variable in the final model. In dark grey are the variables included in the true data-generating model

Table 2:

Proportion of simulated samples for which each variable is included in the optimal model that minimizes the criterion across the entire MCMC chain ( n = 300 )

Proportion of simulated samples for which each variable is included in the optimal model with largest number of votes over the last 5 000 visited models of each MCMC chain ( n = 300 )

Table 4:

3.2 Study 2

Table 5:

Median of MSE ratio between the model identified by each technique and the full ZINB model. Mean and standard deviation (between parentheses) of sensitivity and specificity

Table 6:

Variable selection results for the health care demand dataset. The couple · / · indicates whether a variable has been retained (x) or discarded (–) in the binary part and the count part of the model, respectively

Footnotes

Acknowledgments

References

Summary of the marginal frequency of appearance computed over the 10 000 visited models of the Markov chain, $n = 300$ . The variables included in the true data-generating model are highlighted in bold

Proportion of simulated samples for which each variable is included in the optimal model that minimizes the criterion across the entire MCMC chain ( $n = 300$ )

Proportion of simulated samples for which each variable is included in the optimal model with largest number of votes over the last 5 000 visited models of each MCMC chain ( $n = 300$ )

Variable selection results for the health care demand dataset. The couple $\cdot / \cdot$ indicates whether a variable has been retained (x) or discarded (–) in the binary part and the count part of the model, respectively