Essoham Ali
email: ali.essoham@ugb.edu.sn

Zero-inflated Poisson regression model for a new class of flexible link functions: A case study on healthcare utilization

Keywords: Flexibility, Excess of zeros, Generalized extreme value distribution, Health-care demand, Simulations

Many disciplines produce count data that contain many zeros. Zero inflation count models such as ZIP and ZIB have been widely used to model count data, in particular, to model the latent structure in a ZIP regression model that allows a non-linear functional relationship between covariates and the expected count outcome. A critical problem in modeling the count response data is the appropriate choice of links functions. Commonly used link functions such as logit link have fixed skewness but lack in flexibility to allow the data to determine the degree of the skewness. To overcome this limitation, we propose a flexible ZIP regression model that combines a generalized extreme value link function with the other link functions. The maximum likelihood estimator is used in the estimation of the model. Maximum likelihood estimation is effective in this model in a series of scenarios. Through simulated data sets and analysis of the healthcare utilization, we show that the proposed link function is quite flexible and performs better against link misspecification than standard link functions.

Introduction

Statistical modelling is an important step of data analysis in many fields of scientific research or decision-making purpose. To carry out this approach, one needs to specify a probability distribution that accounts as accurate as possible the variability observed in data. Given the plethora of discrete or continuous distributions available (e.g., [START_REF] Johnson | Univariate Discrete Distributions[END_REF], guidelines are needed to identify not randomly the one or two-parameter family of distributions suited for modelling data on hand. The practice in this procedure is to consider in addition certain phenomenon such as: over-/under-dispersion or zero inflation/deflation for count data ( (e.g., [START_REF] Bonat | Extended Poisson-Tweedie: properties and regression models for count data[END_REF]) and, over-/under-variation or zero mass for continuous data ( (e.g., [START_REF] Abid | Geometric Tweedie regression models for continuous and semicontinuous data with variation phenomenon[END_REF]).

The Poisson dispersion phenomenon is well-known and very widely used in practice; see, (e.g., Kokonendji, C.C. ) for a review of count (or discrete integer-valued) models. Various models have been developed to address zero-inflation, such as zero-inflated (ZI) models which mix a degenerate distribution at zero with a standard count model. Zero-inflated Poisson (ZIP) regression model was proposed by [START_REF] Lambert | Zero-inflated Poisson regression, with an application to defects in manufacturing[END_REF] and further developed by [START_REF] Dietz | On estimation of the Poisson parameter in zero-modifed Poisson models[END_REF], [START_REF] Lim | Zero-inflated Poisson regression mixture model[END_REF] and [START_REF] Monod | Random effects modeling and the zero-inflated Poisson distribution[END_REF], among many others. Recent variants of ZIP regression include random-effects ZIP models [START_REF] Hall | Zero-inflated Poisson and binomial regression with random effects: a case study[END_REF]; [START_REF] Min | Random effect models for repeated measures of zero-inflated count data[END_REF]) and semiparametric ZIP models [START_REF] Lam | Semiparametric analysis of zero-inflated count data[END_REF]. A Zero-inflated negative binomial (ZINB) regression model was proposed by [START_REF] Ridout | A score test for testing a zero-inflated Poisson regression model against zero-inflated negative binomial alternatives[END_REF], see also [START_REF] Moghimbeigi | Multilevel zeroinflated negative binomial regression modeling for over-dispersed count data with extra zeros[END_REF].

Thus, [START_REF] Hall | Zero-inflated Poisson and binomial regression with random effects: a case study[END_REF] introduced the zero-inflated binomial (ZIB) model, see also Hall and Berenhaut (2002), [START_REF] Diop | Maximum likelihood estimation in the logistic regression model with a cure fraction[END_REF], and [START_REF] Diallo | Asymptotic properties of the maximum likelihood estimator in zero-inflated binomial regression[END_REF]. [START_REF] Deng | Score tests for both extra zeros and extra ones in binomial mixed regression models[END_REF] proposed a zero-one inflated binomial regression model for such data. In [START_REF] Nguyen | Asymptotic results in censored zero-inflated Poisson regression[END_REF], authors proposed a Zero-inflated Poisson regression with right-censored data. The usual way to model the response variable is to use a Generalized Linear Model (GLM), where we model the latent probability of "success" by a linear function of covariates through a link function [START_REF] Mccullagh | Generalized linear models[END_REF]. The logit, probit and Student t link functions are three of the common links used in GLM. However, the link functions mentioned above are "symmetric" links in the sense that they assume that the latent probability of a given response variable approaches 0 with the same rate as it approaches 1. Equivalently, the probability density function that corresponds to the inverse cumulative distribution function of the link function is symmetric. However, this may not be a reasonable assumption in many cases. One commonly adopted asymmetric link function is the complementary loglog (cloglog) link function. However, the cloglog link has a fixed negative skewness. As a result, it lacks both the flexibility to let the data display how much skewness should be incorporated and also the ability to allow positive skewness. In short, count data might often be better modeled with flexible link functions that allow for both positive and negative skewness and that allow the data to determine the amount of skewness required.

Many research works have been conducted which introduce flexibility into the link functions. [START_REF] Aranda-Ordaz | On two families of transformations to additivity for binary response data[END_REF] proposed two separate one-parameter models for additional flexibility in the logistic model. [START_REF] Guerrero | Use of the Box-Cox transformation with binary response models[END_REF] used Box-Cox transformation on the odds ratio to form a more flexible class of model. [START_REF] Jones | Reply to Comments on "Families of distributions arising from distributions "of order statistics[END_REF] proposed a family of flexible distributions based on the distribution of order statistics. [START_REF] Stukel | Generalized logistic models[END_REF] proposed a twoparameter class of generalized logistic models. Stukel's model approximates many standard symmetric and asymmetric link functions quite well, but in a Bayesian framework, it may result in improper posteriors when the usual improper uniform prior is used in regressions [START_REF] Chen | A new skewed link model for dichotomous quantal response data[END_REF] . Recently, [START_REF] Wang | Generalized extreme value regression for binary response data:An application to B2B electronic payments system adoption[END_REF] proposed the generalized extreme value link function giving more flexible skewness controlled by the shape parameter. But the standard logistic and probit links are not among the special cases of this family.

A critical problem in modeling the count response data is the appropriate choice of links functions. To overcome this limitation, we propose a flexible Zero-Inflated Poisson regression model that combines a generalized extreme value link function with the other link functions. In the extreme value theory, the GEV distribution is used to model the tail of a distribution [START_REF] Coles | An Introduction to Statitical Modelling of Extreme Values[END_REF]. Currently, the logistic regression model, with its convenient interpretation and implementation, has been routinely employed to estimate and predict.As in this work, we focus on the Poisson parameters we have chosen to vary several link functions in order to see the flexibility of the GEV distribution with respect to the others. In the GLM, [START_REF] Agresti | Categorical Data Analysis[END_REF], log-log and complementary log-log link functions are used since they are asymmetric functions. In particular, the log-log link function is the quantile function of the Gumbel random variable. The inverse function of the complementary log-log is equal to one minus the cumulative distribution function of the Gumbel random variable. Consequences of link misspecification have been studied by numerous authors in the literature. In particular, for independent binary observations, [START_REF] Czado | The effect of link misspecification on binary regression inference[END_REF] demonstrated that falsely assuming a logistic link leads to a substantial increase in the bias and the mean squared error of the parameter estimates as well as the predicted probabilities, both asymptotically and in finite samples. Moreover, these undesirable effects have greater magnitude when the misspecification involves skewness than when it involves kurtosis (or tail weight). [START_REF] Wu | On the relationship between links for binary response data[END_REF] showed also that under certain conditions there exist linear relationships between the regression coefficients though the choice of links is important for goodness of fit. To build an appropriate and extremely flexible model for the count data and to overcome the constraint for the skewed generalized link models, we propose the cloglog, probit and generalized extreme value (GEV) distribution as a link function. In this paper we then suggest a new class of link functions to model count data, and apply it to healthcare utilization data. This paper is organized as follows. In Section 2, we recall the defnition of ZIP regression model, we describe the maximum likelihood estimation under different link functions and we introduce some useful notations. In Section 3, we report the results of our simulations. An application to a health-care utilization dataset is described in Section 4. Some concluding remarks are given in Section 5.

Notations and likelihood calculation

Let us first specify the notation used throughout this paper. Let Z i ∼ π i δ 0 + (1π i )P(λ i ) denote the count of interest and X = (1, X 2 , • • • , X p ) be a p-vector of covariates ( denotes the transpose operator and let J i = 1 {Z i =0} ). π i is the probability of success for the ith observation. We assume that the conditional distribution of Z given X is given by a Poisson regression model with parameter λ i = e β X i , where β ∈ R p is a vector of unknown parameters. We associate π i and W i through a cumulative distribution function F as follows:

π i = F (γ W i ) (2.1)
where F is a cumulative distribution function and F -1 determines the link function.

W i = (1, W i2 , • • • , W iq
) be a q-vector of covariates and γ ∈ R q is a vector of unknown parameter.

Susceptibility probability function with different links functions 2.1.1. Zero-inflated ZIP regression model

The ZIP model assumes that the response variable Z i (where the lower indice i indicates the individual) is such that

Z i ∼ 0 with probability π i , P(λ i ) with probability 1 -π i , (2.2)
where P(λ i ) denotes Poisson distribution with parameter λ i > 0. Obviously, the ZIP model reduces to a standard Poisson distribution if π i = 0. In ZIP regression, the mixing probability π i and parameter λ i are usually modeled by logistic and log-linear models respectively, that is:

F -1 (π i ) = logit(π i ) = γ W i (2.3) and log(λ i ) = β X i , (2.4)
Suppose that we observe a sample of n independent copies

(Z i , X i , W i ), i = 1, . . . , n of (Z, X, W). For i = 1, . . . , n, the log-likelihood of θ = (β , γ ) in the latent class ZIP model (2.6)-(2.3)-(2.4) is : n (θ) = n i=1 J i log e γ W i + e -exp(β X i ) + (1 -J i ) Z i β X i -log (Z i !) -log 1 + e γ W i .
The maximum likelihood estimator of (β, γ) is obtained by maximizing this function. The ML estimator is consistent and asymptotically normally distributed (see [START_REF] Czado | Consistency and asymptotic normality of the maximum likelihood estimator in a zero-inflated generalized Poisson regression[END_REF]).

ZIP-GEV regression model

A key component of the model given in (2.3)-(2.4) is the specification of the link function. The commonly used logit link is specified as π i = F (γ W i ), where F is a cumulative distribution function (cdf) and F -1 determines the link function. The symmetry in the normal distribution leads to the symmetry in the logit link. [START_REF] Wang | Generalized extreme value regression for binary response data:An application to B2B electronic payments system adoption[END_REF] showed that the symmetric link has an inferior performance when the data structure requires a skewed response probability function. They proposed a link function based on the GEV distribution. The distribution function of GEV (µ, σ, ξ) is given by:

G(x|µ, σ, ξ) =      exp -1 + ξ (x-µ) σ -1/ξ + , ξ = 0, exp -exp( (x-µ) σ ) , ξ = 0, (2.5)
where µ ∈ R, σ ∈ R + and ξ ∈ R are, respectively, the location, scale and shape parameters, and x+ = max(0, x). The shape of this distribution function is very flexible with the tail behavior controlled by the shape parameter ξ. When ξ = 0, it is the Gumbel distribution and decays exponentially. When ξ < 0, it reduces to the negative Weibull distribution with a finite short upper endpoint. When ξ > 0, it becomes the Fréchet distribution with a heavy tail behavior. The GEV link is the inverse of F which is assumed as

π i = F (W i |ξ) = 1 -GEV(-γ W i ; ξ) =    1 -exp -(1 -ξγ W i ) -1/ξ + , ξ = 0, 1 -exp -exp( (x-µ) σ ) , ξ = 0, (2.6)
where GEV(x; ξ) represents the cumulative probability at x for the GEV distribution with parameters φ = (µ = 0, σ = 1, ξ). Note µ and σ are set to fixed constants for model identifiability. [START_REF] Wang | Generalized extreme value regression for binary response data:An application to B2B electronic payments system adoption[END_REF] showed that the GEV link model specified in ( 6) is negatively skewed for ξ < log 2 -1 and positively skewed for ξ > log 2 -1. The link function is approximately symmetric at ξ = log 2 -1. The cloglog link, specified as

F -1 (π i ) = -log(-log(π i )) = γ W i , is a special case of the GEV link with ξ = 0.
The GEV regression model proposed by [START_REF] Calabrese | Modelling Small and Medium Enterprise Loan Defaults as Rare Events: The Generalized Extreme Value Regression Model[END_REF] is defined by a link function that corresponds to the inverse cumulative function of the GEV distribution, that can be called GEV regression model or "gevit, in analogy with the "logit". The ZIP regression model under the GEV link is then given by gevit

(π i ) = [-log(π i )] -ξ -1 ξ = γ W i = γ 1 + q j=1 γ j W ij (2.7) log(λ i ) = β X i = β 1 + p k=1 β k X ik (2.8)
where ξ ∈ R is the shape parameter for GEV distribution. According to (2.6)-(2.7)-(2.8), the log-likelihood of θ = (β , γ )

GEV n

(θ) = n i=1 J i log exp -(1 + ξγ W i ) -1 ξ + (1 -exp -(1 + ξγ W i ) -1 ξ )e -exp (β X i ) +(1 -J i ) Z i β X i -e β X i + log 1 -exp -(1 + ξγ W i ) -1 ξ -log (Z i !) .
The MLE θn = ( β n , γ n ) of θ is obtained by solving the score equation

∂ GEV n (θ) ∂θ = 0, (2.9)
which can be achieved by nonlinear optimization

ZIP-cloglog regression model

The asymmetric cloglog link is specified as

F -1 (π i ) = -log(-log(π i )) = γ W i .
(2.10)

Assume that we observe n independent vectors

(Z 1 , X 1 , W 1 ), • • • , (Z n , X n , W n ) from the model (2.4)-(2.
10), all defined on the probability space (Ω, C, P). The log-likelihood of θ = (β , γ ) based on these observations is

cloglog n (θ) = n i=1 J i log e -exp (-γ W i ) + (1 -e -exp (-γ W i ) )e -exp (β X i ) + n i=1 (1 -J i ) Z i β X i -e β X i + log 1 -e -exp (-γ W i ) -log (Z i !) , = n i=1 i (θ).
The maximum likelihood estimator θn

= ( β n , γ n ) of θ is solution of the k-dimensional score equation ˙ i (θ) = ∂ cloglog n (θ) ∂θ = 0 (2.11)
where k = p + q.

ZIP-probit regression

The zero-inflated Poisson model using the probit link function can be defined in the same way as the ZIP model, where the probability of zero inflation is modelled by the probit model. When risk factors are available, the mixing probability π i is usually modeled by a probit model :

F -1 (π i ) = Φ(γ W i ),
(2.12)

where Φ is the distribution function of N (0, 1). According to (2.4)-(2.12) the log-

likelihood of θ = (β , γ ) probit n (θ) = n i=1 J i log Φ(γ W i ) + (1 -Φ(γ W i ))e -exp (β X i ) +(1 -J i ) Z i β X i -e β X i + log 1 -Φ(γ W i ) -log (Z i !) . The MLE θn = ( β n , γ n ) of θ is the solution of the k-dimensional score equation ∂ probit n (θ) ∂θ = 0 (2.13)
Solving this (non-linear) equation is relatively straightforward using standard mathematical softwares Remark 1. A rigorous assessment of the asymptotic properties of θn is presented in the censored ZIP model [START_REF] Nguyen | Asymptotic results in censored zero-inflated Poisson regression[END_REF]. In this paper, such properties can be expected in the ZIP model regardless of the link function used to model the probability of susceptibility. However, leaving aside the distribution theory, we propose to study these properties by means of simulations.

A simulation study

In this section, we compare, via simulations, the performance of four links functions (2.4)-(2.7)-(2.10)-(2.12) used to model the probability of zero-inflated. We generate 2 covariates for our simulation study

X i = (1, X i2 , . . . , X ip ) and W i = (1, W i2 , • • • , W iq ) , where X i1 = W i1 = 1 and the X i2 , • • • , X i6 , W i4 , W i5 are independently drawn from nor- mal N (0, 1),binomial B(1, 0.3), normal N (1, 1.5), exponential E(1), uniform U(2, 5), normal N (-1, 1
) and binomial B(1, 0.5) distributions respectively. Linear predictors are allowed to share common terms by letting

W i2 = X i2 et W i3 = X i3 .
The regression parameter β is chosen as β = (0.7, 0.1, 0.4, 0.85, -0.5, 0) for all simulations. With the same value of β, we carry out our studies under four scenarios based on four true models as follows

Simulation scenario

Scenario 1: The following ZIP regression model is used to simulate data :

log(λ i (β)) = β 1 X i1 + β 2 X i2 + β 3 X i3 + β 4 X i4 + β 5 X i5 + β 6 X i6 , logit(π i ) = γ 1 W i1 + γ 2 W i2 + γ 3 W i3 + γ 4 W i4 , +γ 5 W i5 ,
The regression parameter γ is chosen as γ = (-0.9, -0.65, -0.2, 0.65, 0) . In this setting, the average proportion of zero-inflated data is 0.20.

Scenario 2: The count data are generated from the cloglog link model with

F -1 (π i ) = -log(-log(π i )) = γ W i .
We consider two values for γ, namely : γ = (0.5, -0.60, -0.2, 0.75, 0) and γ = (0.25, -0.9, 0.60, -0.45, 0) . With these values, the average proportion c of zeroinflated data in the simulated data sets is 0.20 and 0.60 respectively.

Scenario 3: We simulate the data according to the ZIP-GEV model (2.7)-(2.8) defined by: log

(λ i (β)) = β 1 X i1 + β 2 X i2 + β 3 X i3 + β 4 X i4 + β 5 X i5 + β 6 X i6 , gevit(ω i ) = γ 1 W i1 + γ 2 W i2 + γ 3 W i3 + γ 4 W i4 , +γ 5 W i5 ,
The advantage of the GEV link model we are talking about here is that it integrates a wide range of asymmetries with the shape parameter ξ. But in our simulations we choose ξ = 0.5 belonging to Frechet's domain. The regression parameter γ is chosen as:

• case 1: γ = (-0.95, 0.5, -0.4, -0.65, 0)

• case 2: γ = (0.8, 0.2, 0.4, -0.8, 0) . Using these values, in case 1 (respectively case 2), the average percentage of zero-infation in the simulated data sets is 0.25 (respectively 0.75).

Scenario 4: In a second set of simulation scenarios, the data sets simulated from the ZIP-probit model with F -1 (π i ) = Φ -1 (γ W i ).Two values for γ, namely : γ = (-0.5, 0, 0.65, 0.8, 0) and γ = (-0.1, 0.85, 0.9, -0.2, 0) .The parameter vector γ ∈ R 5 is chosen to yield various average proportions of zero-inflation within each sample, namely : 0.20 and 0.60.

We consider the following sample sizes: n = 500; 2000. For each combination of the simulation design parameters (sample size, proportions of zero-inflation),We simulate N = 1000 replications for each combination [sample size × proportion of zero-inflation] of the design parameters. Simulations are conducted using the statistical software R R Core Team (2018). We use the package maxLik [START_REF] Henningsen | maxLik: A package for maximum likelihood estimation in R[END_REF] to solve the score equation (2.9)-(2.11)-(2.13) via a Newton-Raphson algorithm.

Results

For each configuration sample size× zero-inflation proportion of the simulation design parameters, we calculate the average bias, standard deviation, average standard error and root mean square error of the estimate over the N simulated samples. We also obtain the empirical coverage probability and average length of 95%-level Wald confidence intervals for the β j . The results are described in Table 1 (first scenario), Table 2 andTable 3 (second scenario), Table 4 andTable 5 (third scenario), Table 6 and Table7 (fourth scenario) for n = 500 and n = 2000 respectively. Through simulations, we also assess the normal approximation by plotting estimated densities obtained from the N normalized estimates ( βj,n -β j )/standard error( βj,n ) j, k = 1, . . . , 6, and by comparing with the density of the standard normal distribution. Standard errors are obtained as the square roots of the diagonal elements of the estimated variance matrix for ours models. Figures 1, 2 and 3 provide results for ZIP-GEV model (n = 500 ,25% of zeroinflation), ZIP-cloglog (n = 500 and 20% of zero-inflation), and ZIP-probit model (n = 500 ,30% of zero-inflation) . Plots for the other scenarios are similar and thus are not given. From these results, it appears, as expected, that the bias, variability and length of confidence intervals of all estimates decrease as the sample size increases. For fixed n, we observe that performances of the βj,n s remain stable when the proportion of zero-inflation varies from small to moderate values.

These observations illustrate the general fact that accurate estimation in a zeroinflated regression model requires a balance between susceptible and non susceptible subpopulations (that is, a sufficient amount of zero and non zero observations should be available to accurately estimate the zero-inflation probabilities and count submodel). Also, empirical coverage probabilities are close to the nominal level, which indicates that the normal approximation of the distribution of the MLE is appropriate, even when the sample size is moderate. This is confirmed by Figures 1, 2 and 3.

A comparison of the four models

In this section, we compare, through simulations, the performance of four models . We obtain the MLE in the four models, for the four scenarios. In the first scenario, our real model is the ZIP where the probability of zero inflation is modeled by the logit link. The other models ( ZIP-probit, ZIP-cloglog and ZIP-GEV ) misspecifies the susceptibility probability π i . In the second scenario, the true model is the ZIP-cloglog . In the third scenario, the true model is the ZIP-cloglog. In the fourth scenario, the true model is the ZIP-probit. In all four cases, the γ estimates are assumed to be biased in the misspecified model. This is confirmed by the simulation results. However, in all four scenarios, the interest is generally on the β , which relates the covariates to the λ i intensity of the account response. For this reason, we provide results only for β. Moreover, since the proposed models adopt the same specification for π i , a comparison of the β's estimates of the four models is fair. The results are described in Table 1 (Scenario 1), Table 2 andTable 3 (Scenario 2), Table 4 andTable 5 (Scenario 3), Table 6 andTable7 (Scenario 4).

It appears that in all four models, the estimate of β is quite robust to a misspecification of the probability of susceptibility. That is, when the logit model is used to generate the data (scenario 1), the β estimates in the ZIP-Probit and ZIP-GEV models are of good quality. Referring to scenario 2 and scenario 3 described above, we validate the β estimates in the ZIP-probit and ZIP-GEV models.

Conversely, when the ZIP-probit and ZIP-GEV models are used to simulate the data (Scenario 3 and Scenario 4), the β estimates in these models perform equally well and better than the others proposed. We also observe that the estimates obtained from the ZIP-probit and ZIP-GEV models behave almost systematically better than the estimates based on the other models, even when the ZIP-logit or ZIP-cloglog is used to simulate the data.

Applications with real-life data

Data description and competing models

The data are obtained from the National Medical Expenditure Survey (NMES) which was conducted in 1987 and 1988 to provide a comprehensive picture of how Americans use and pay for health services. The NMES is based upon a representative, national probability sample of the civilian, non-institutionalized population and individuals admitted to longterm care facilities during 1987. Under the household survey of the NMES, more than 38000 individuals in 15000 households across the United States were interviewed quarterly about their health insurance coverage, the services they used, and the cost and source of payments of those services. These data were verified by cross-checking information provided by survey respondents with providers of health-care services. In addition to health-care data, NMES provides information on health status, employment, sociodemographic characteristics, and economic status.

In this paper we consider a subsample of individuals ages 66 and over (a total of 4406 observations) all of whom are covered by Medicare, a public insurance programme that offers substantial protection against health care costs. Residents of the United States are eligible for Medicare coverage at age 65. Some individuals start receiving Medicare benefits a few months into their 65th year primarily because they fail to apply for coverage at the appropriate time. Virtually all individuals who are 66 or older are covered by Medicare.

In addition, most individuals make a choice of supplemental private insurance coverage shortly before or in their 65th year because the price of such insurance rises sharply with age and coverage becomes more restrictive. The response variable is the number of visits to a physician in an office setting (denoted by ofp in what follows). Available covariates include: i) socio-economic variables: gender (1 for female, 0 for male), age (in years, divided by 10), marital status, educational level (number of years of education), income, ii) various measures of health condition: number of chronic conditions (cancer, arthritis, gallbladder problems • • • ) and a variable indicating self-perceived health level (poor, average, excellent) and iii) a binary variable indicating whether individual is covered by medicaid or not (medicaid is a US health insurance for individuals with limited income and resources, we code it as 1 if the individual is covered and 0 otherwise). Self-perceived health is recoded as two dummy variables denoted by "health1" (1 if health is perceived as poor, 0 otherwise) and "health2" (1 if health is perceived as excellent, 0 otherwise).

We fit the following four models : i) a ZIP regression model where λ i and π i are specified as in (2.3)-(2.4); ii) the ZIP model with cloglog link (denoted by ZIP-cloglog thereafter)where π i is as in (2.10) , iii) the ZIP model with with probit link (denoted by ZIP-probit), where π i is as in (2.12) and iv) the ZIP model with GEV link (denoted by ZIP-GEV),where π i is as in (2.7). Selection of regressors for inclusion in π i requires some care. Indeed, it was previously observed in various other zero-inflated models that including all available regressors in both count and zero-inflation probabilities can yield lack of identification of model parameters. See for example [START_REF] Diop | Maximum likelihood estimation in the logistic regression model with a cure fraction[END_REF], who suggest to solve this issue by letting at least one of the covariates included in the count model to be excluded from the zero-inflation model (or the converse). Such condition is not required in the ZIP model. Using the Wald testing, we identify five significant predictors : age, gender, educational level, number of chronic conditions and medicaid status, that are included in π i .

Results for the fours competing models (ZIP, ZIP-cloglog , ZIP-probit and ZIP-GEV ) are displayed in Table 8. We report estimate, standard error and significance level of Wald test for each parameter. For purpose of comparison, we also report AIC and BIC values for the four models. ZIP-GEV appears as the best model in terms of both AIC . A closer look at the results from the widely used logit link regression model in the healthcare utilization research and our GEV regression model reveals some difference in the estimation of the covariates effects. Gender, educational level and medicaid status are identified by ZIP-GEV as the most influencing factors for being a permanent non-user, with medicaid recipients being more likely to be permanent non-users. The four models identify the same subset of influent factors for healthcare utilization, with similar parameter estimates.

From Table 8, we observe that in the overall population, significant determinants of the decision to consult a non-physician when visiting in an office setting include health status, age, gender, educational level and medicaid status. Patients with poor health will favor office visits to a physician over office visits to a non-physician, which seems a natural observation. Women and people with higher education have higher probability to consult a non-physician, while medicaid recipients are more likely to visit physicians than nonphysicians. The probability of visiting a non-physician when consulting in an office setting decreases with age. This may be due to several factors, such as decreasing mobility associated with ageing (aged patients will tend to limit their consultations to those considered as the most necessary, that is, to physician visits) and worsening of the health condition with ageing (patients whose health declines are likely to favor visits to a physician). Table 8: : Health-care data analysis: estimates,standard errors and significance codes: *** significant at the 0.1% level, ** significant at the 1% level, * significant at the 5% level, . significant at the 10% level..

Concluding Remarks

In this paper, we study the properties of MLE in ZIP regression models when the susceptibility probability function is modeled with different links functions. Our simulations suggest that the MLE works well and that reliable statistical inferences about the parameters of interest in the different models can be based on the normal approximation of the MLE distribution. Maximum likelihood estimation is shown to perform well in this model, under a range of scenarios. Moreover, in our analysis of health-care utilization, the proposed model provides plausible explanations and interpretations and gives useful insight of the decision of using or not available healthcare services. Several issues now require more attention, such as estimation in the bivariate ZIP-GEV regression in various forms. Investigating the estimation of a flexible ZIP regression model that combines a generalized extreme value link function with a Gaussian process is also desirable. All these issues will be tackled in future works.

Table 3 :

 3

					Sample size n = 500				Sample size n = 2000
			β1,n	β2,n	β3,n	β4,n	β5,n	β6,n	β1,n	β2,n	β3,n	β4,n	β5,n	β6,n
		bias	-0.0009 -0.0012 -0.0000 -0.0003 0.0004 0.0002 -0.0014 0.0000 -0.0002 0.0000 0.0002 0.0003
		SD	0.0902 0.0193 0.0373 0.0138 0.0283 0.0223 0.0439 0.0092 0.0183 0.0063 0.0139 0.0107
	GEV SE	0.0911 0.0193 0.0384 0.0137 0.0286 0.0219 0.0441 0.0093 0.0186 0.0064 0.0139 0.0106
		RMSE 0.1282 0.0273 0.0535 0.0195 0.0402 0.0312 0.0623 0.0131 0.0261 0.0090 0.0196 0.0151
		CP	0.9440 0.9530 0.9610 0.9440 0.9470 0.9460 0.9450 0.9530 0.9580 0.9540 0.9500 0.9590
		(CI)	0.3560 0.0754 0.1503 0.0532 0.1115 0.0857 0.1729 0.0363 0.0727 0.0249 0.0542 0.0415
		bias	-0.0346 -0.0026 0.0078 0.0102 -0.0086 0.0001 -0.0324 -0.0013 0.0066 0.0093 -0.0083 0.0004
		SD	0.0917 0.0196 0.0382 0.0143 0.0294 0.0227 0.0449 0.0094 0.0187 0.0067 0.0143 0.0109
	logit	SE	0.0929 0.0196 0.0388 0.0141 0.0294 0.0222 0.0449 0.0094 0.0187 0.0065 0.0142 0.0107
		RMSE 0.1350 0.0278 0.0550 0.0226 0.0425 0.0318 0.0713 0.0133 0.0273 0.0132 0.0218 0.0153
		CP	0.9370 0.9470 0.9580 0.8860 0.9420 0.9470 0.9020 0.9470 0.9400 0.7080 0.9170 0.9490
		(CI)	0.3628 0.0763 0.1519 0.0548 0.1149 0.0868 0.1759 0.0367 0.0734 0.0255 0.0558 0.0420
		bias	-0.0004 -0.0012 -0.0000 -0.0004 0.0005 0.0002 -0.0011 -0.0001 -0.0002 0.0000 0.0002 0.0003
		SD	0.0902 0.0193 0.0373 0.0138 0.0283 0.0223 0.0439 0.0092 0.0183 0.0063 0.0139 0.0107
	cloglog SE	0.0911 0.0193 0.0384 0.0137 0.0286 0.0219 0.0442 0.0093 0.0186 0.0064 0.0139 0.0106
		RMSE 0.1282 0.0273 0.0535 0.0195 0.0402 0.0312 0.0623 0.0131 0.0260 0.0090 0.0196 0.0151
		CP	0.9440 0.9550 0.9620 0.9430 0.9460 0.9460 0.9440 0.9520 0.9580 0.9540 0.9500 0.9600
		(CI)	0.3561 0.0754 0.1503 0.0533 0.1116 0.0857 0.1729 0.0363 0.0727 0.0249 0.0543 0.0415
		bias	-0.0003 -0.0014 0.0002 -0.0005 0.0006 0.0002 -0.0006 -0.0003 0.0000 -0.0002 0.0003 0.0003
		SD	0.0902 0.0192 0.0373 0.0138 0.0283 0.0223 0.0439 0.0092 0.0183 0.0064 0.0139 0.0107
	probit SE	0.0912 0.0193 0.0384 0.0137 0.0286 0.0219 0.0442 0.0093 0.0186 0.0064 0.0139 0.0106
		RMSE 0.1282 0.0273 0.0535 0.0195 0.0402 0.0313 0.0622 0.0131 0.0260 0.0090 0.0196 0.0151
		CP	0.9450 0.9560 0.9620 0.9430 0.9440 0.9450 0.9450 0.9530 0.9580 0.9520 0.9500 0.9600
		(CI)	0.3562 0.0754 0.1503 0.0533 0.1117 0.0857 0.1730 0.0363 0.0727 0.0249 0.0543 0.0415

Simulation results for scenario 2 (data are simulated from the ZIP-cloglog model (2.4)-(2.10), average proportion of zero-inflation = 60%). SD: empirical standard deviation. SE: average standard error. CP: empirical coverage probability of 95%-level confidence intervals. (CI): average length of confidence intervals.

Table

4

: Simulation results for scenario 3 (data are simulated from the ZIP-GEV model (2.7)-(2.8), average proportion of zero-inflation = 25%). SD: empirical standard deviation. SE: average standard error. CP: empirical coverage probability of 95%-level confidence intervals. (CI): average length of confidence intervals.

Table

5

: Simulation results for scenario 3 (data are simulated from the ZIP-GEV model (2.7)-(2.8), average proportion of zero-inflation = 75%). SD: empirical standard deviation. SE: average standard error. CP: empirical coverage probability of 95%-level confidence intervals. (CI): average length of confidence intervals.

Table 6 :

 6 Simulation results for scenario 4 (data are simulated from the ZIP-Probit model (2.4)-(2.12), average proportion of zero-inflation = 20%). SD: empirical standard deviation. SE: average standard error. CP: empirical coverage probability of 95%-level confidence intervals. (CI): average length of confidence intervals.

					Sample size n = 500				Sample size n = 2000
			β1,n	β2,n	β3,n	β4,n	β5,n	β6,n	β1,n	β2,n	β3,n	β4,n	β5,n	β6,n
		bias	-0.0003 -0.0014 0.0019 0.0004 -0.0022 -0.0004 0.0009 0.0001 -0.0012 -0.0001 -0.0013 0.0001
		SD	0.1326 0.0340 0.0730 0.0214 0.0428 0.0329 0.0624 0.0155 0.0313 0.0099 0.0214 0.0156
	probit SE	0.1393 0.0347 0.0723 0.0220 0.0442 0.0336	0.066 0.0162 0.0332 0.0098 0.0208 0.0158
		RMSE 0.1923 0.0486 0.1027 0.0307 0.0615 0.0470 0.0909 0.0224 0.0456 0.0140 0.0299 0.0222
		CP	0.9580 0.9530 0.9550 0.9600 0.9570 0.9520 0.9680 0.9580 0.9620 0.9530 0.9450 0.9480
		(CI)	0.5433 0.1352 0.2804 0.0849 0.1719 0.1310 0.2585 0.0635 0.1297 0.0383 0.0815 0.0619
		bias	-0.2035 -0.0603 -0.0426 0.0649 -0.0576 -0.0018 -0.1848 -0.0489 -0.0374 0.0537 -0.0494 0.0001
		SD	0.1664 0.0455 0.0915 0.0320 0.0524 0.0404 0.0796 0.0205 0.0400 0.0154 0.0264 0.0185
	logit	SE	0.1468 0.0368 0.0769 0.0229 0.0472 0.0350	0.069 0.0170 0.0347 0.0102 0.0221 0.0163
		RMSE 0.3011 0.0840 0.1269 0.0759 0.0910 0.0535 0.2127 0.0556 0.0648 0.0567 0.0602 0.0247
		CP	0.7000 0.6180 0.8750 0.2250 0.7580 0.9120 0.2680 0.2280 0.7770 0.0000 0.4160 0.9160
		(CI)	0.5719 0.1429 0.2969 0.0886 0.1832 0.1365 0.2700 0.0664 0.1355 0.0396 0.0863 0.0639
		bias	-0.0033 0.0002 0.0034 0.0013 -0.0027 -0.0002 -0.0012 0.0018 0.0003 0.0005 -0.0015 0.0001
		SD	0.1334 0.0343 0.0730 0.0216 0.0431 0.0331 0.0625 0.0158 0.0314 0.0100 0.0214 0.0157
	GEV SE	0.1392 0.0348 0.0724 0.0219 0.0442 0.0336	0.066 0.0162 0.0332 0.0098 0.0208 0.0158
		RMSE 0.1928 0.0488 0.1028 0.0308 0.0618 0.0471 0.0909 0.0227 0.0457 0.0140 0.0299 0.0222
		CP	0.9580 0.9470 0.9550 0.9560 0.9560 0.9550 0.9710 0.9540 0.9630 0.9500 0.9460 0.9480
		(CI)	0.5430 0.1353 0.2805 0.0847 0.1718 0.130	0.2584 0.0635 0.1298 0.0383 0.0814 0.0619
		bias	-0.0017 0.0001 0.0034 0.0009 -0.0025 -0.0004 -0.0003 0.0015 0.0001 0.0003 -0.0015 0.0001
		SD	0.1327 0.0341 0.0729 0.0215 0.0429 0.0330 0.0624 0.0156 0.0313 0.0100 0.0214 0.0156
	cloglog SE	0.1393 0.0348 0.0724 0.0219 0.0442 0.0336	0.066 0.0162 0.0332 0.0098 0.0208 0.0158
		RMSE 0.1923 0.0487 0.1028 0.0307 0.0616 0.0470 0.0909 0.0225 0.0456 0.0140 0.0299 0.0222
		CP	0.9580 0.9490 0.9540 0.9580 0.9570 0.9540 0.9680 0.9540 0.9620 0.9480 0.9480 0.9480
		(CI)	0.5431 0.1353 0.2804 0.0848 0.1718 0.1310 0.2584 0.0635 0.1298 0.0383 0.0814 0.0619

Table 7 :

 7 Simulation results for scenario 4 (data are simulated from the ZIP-Probit model (2.4)-(2.12), average proportion of zero-inflation = 60%). SD: empirical standard deviation. SE: average standard error. CP: empirical coverage probability of 95%-level confidence intervals. (CI): average length of confidence intervals.