Marginal effects and incremental effects in two-part models for endogenous healthcare utilization in health services research

Abstract

In health services research, endogenous healthcare utilization refers to the notion that the choice of utilizing health services is endogenous due to its correlation with the intensity of utilization outcomes, such as the number of emergency room visits. Greene in (Empir Econ 36:133–173, 2009) extended four conventional two-part models for zero-abundant count utilization to the two-part models that account for endogenous utilization. However, statistical inference on the (average) marginal and incremental effects in these models has not been carefully studied. The present article provides the estimation formulations for (average) marginal and incremental effects of the four two-part models: zero-inflated Poisson and negative binomial models, and hurdle Poisson and negative binomial models with correlated errors that characterize endogenous healthcare utilization. The variance estimation derived from the delta method is provided to facilitate the statistical inference of these effects. We then perform simulation studies to numerically justify our methodology. An empirical study is presented to investigate the average effects of household income and health insurance status on healthcare utilization with the German Scocioeconomic Panel data. The four models give consistent results regarding interpretations for moral hazards and adverse selection in the study.

Appendices

Appendix 1: Density functions and log-likelihood functions

1.1 Zero-inflated Poisson and zero-inflated negative binomial models with endogenous utilization

In the framework of zero-inflated models with endogenous utilization, for the ${\hbox {ZIP}} _{\mathrm{eu}}$ model, the marginal probability function is given by

$$\begin{aligned} f(y_i|x_i,z_i,\varepsilon _i) = {\left\{ \begin{array}{ll} \varPhi (-A_i) + \varPhi (A_i) e^{-e^{x_i'\beta +\sigma \varepsilon _i}},&{} \mathrm{for~~ } y_i =0,\\ \varPhi (A_i) {e^{-e^{x_i'\beta +\sigma \varepsilon _i}}e^{(x_i'\beta +\sigma \varepsilon _i)y_i}}/{y_i!},&{} \mathrm{for~~ } y_i =1,2,\dots , \end{array}\right. } \end{aligned}$$

where $\frac{\mathrm{var}(y_i|x_i, z_i, \varepsilon _i)}{\mathrm{E}(y_i|x_i, z_i, \varepsilon _i)} =1+\varPhi (-A_i)e^{x_i'\beta +\sigma \varepsilon _i}$, which demonstrates overdispersion of the outcomes. If we denote the parameter by $\theta =(\beta ', \gamma ', \rho ,\sigma )$, and the log-likelihood function is

$$\begin{aligned} \ell (\theta | x_i,z_i)&= {} \sum \limits _{i=1,\ldots ,n,y_i=0} \ln \int _{-\infty }^\infty \left\{ \varPhi (-A_i) + \varPhi (A_i) e^{-e^{x_i'\beta +\sigma \varepsilon _i}}\right\} \phi (\varepsilon _i)d\varepsilon _i\nonumber \\&\quad + \sum \limits _{i=1,\ldots ,n,y_i>0} \ln \int _{-\infty }^\infty \frac{e^{-e^{x_i'\beta +\sigma \varepsilon _i}}e^{(x_i'\beta +\sigma \varepsilon _i)y_i}}{y_i!} \varPhi (A_i) \phi (\varepsilon _i)d\varepsilon _i. \end{aligned}$$

(24)

for the $\hbox {ZINB}_{\mathrm{eu}}$ model, the marginal probability function is as below

$$\begin{aligned} f(y_i|x_i,z_i,\varepsilon _i) = {\left\{ \begin{array}{ll} \displaystyle \varPhi (-A_i) + \varPhi (A_i) \left( \frac{\alpha }{\alpha +e^{x_i'\beta +\sigma \varepsilon _i}}\right) ^\alpha ,&{} \mathrm{for~~ } y_i =0,\\ \displaystyle \frac{\varPhi (A_i) \varGamma (y_i+\alpha )}{\varGamma (\alpha )\varGamma (y_i+1)} \left( \frac{\alpha }{\alpha +e^{x_i'\beta +\sigma \varepsilon _i}}\right) ^\alpha \left( \frac{e^{x_i'\beta +\sigma \varepsilon _i}}{\alpha +e^{x_i'\beta +\sigma \varepsilon _i}}\right) ^{y_i}, &{} \mathrm{for~~ } y_i =1,2,\dots , \end{array}\right. } \end{aligned}$$

where $\frac{\mathrm{var}(y_i|x_i, z_i, \varepsilon _i)}{\mathrm{E}(y_i|x_i, z_i, \varepsilon _i)} =1+\varPhi (-A_i)e^{x_i'\beta +\sigma \varepsilon _i} +\frac{e^{x_i'\beta +\sigma \varepsilon _i} }{\alpha }$. We can see that the overdispersion in the $\hbox {ZINB}_{\mathrm{eu}}$  models is more intense than that in the regular NB models. Let’s denote the parameter by $\theta =(\beta ', \gamma ', \rho ,\sigma ,\alpha )$. Then, the log-likelihood function is

$$\begin{aligned} \ell (\theta | x_i,z_i)&= {} \sum \limits _{i=1,\ldots ,n,y_i=0} \ln \int _{-\infty }^\infty \left\{ \varPhi (-A_i) + \varPhi (A_i) \left( \frac{\alpha }{\alpha +e^{x_i'\beta +\sigma \varepsilon _i}}\right) ^\alpha \right\} \phi (\varepsilon _i)d\varepsilon _i\nonumber \\&\quad + \sum \limits _{i=1,\ldots ,n,y_i>0} \ln \int _{-\infty }^\infty \left( \frac{\alpha }{\alpha +e^{x_i'\beta +\sigma \varepsilon _i}}\right) ^\alpha \left( \frac{e^{x_i'\beta +\sigma \varepsilon _i}}{\alpha +e^{x_i'\beta +\sigma \varepsilon _i}}\right) ^{y_i} \varPhi (A_i) \phi (\varepsilon _i)d\varepsilon _i \nonumber \\&\quad + \sum \limits _{i=1,\ldots ,n,y_i>0} \ln \frac{\varGamma (y_i+\alpha )}{\varGamma (\alpha )\varGamma (y_i+1)}. \end{aligned}$$

(25)

1.2 Hurdle Poisson and hurdle negative binomial models with endogenous utilization

For the $\hbox {HP}_{\mathrm{eu}}$ model, the log-likelihood function for the model is

$$\begin{aligned} \ell (\theta | x_i,z_i)&= {} \sum \limits _{i=1,\ldots ,n,y_i=0} \ln \int _{-\infty }^\infty \varPhi (-A_i) \phi (\varepsilon _i)d\varepsilon _i\nonumber \\&\quad + \sum \limits _{i=1,\ldots ,n,y_i>0} \ln \int _{-\infty }^\infty \frac{e^{-e^{x_i'\beta +\sigma \varepsilon _i}}e^{(x_i'\beta +\sigma \varepsilon _i)y_i}}{\left( 1-e^{-e^{x_i'\beta +\sigma \varepsilon _i}}\right) y_i!} \varPhi (A_i) \phi (\varepsilon _i)d\varepsilon _i. \end{aligned}$$

(26)

For the $\hbox {HNB}_{\mathrm{eu}}$ model, the log-likelihood function is

$$\begin{aligned}&\ell (\theta | x_i,z_i)\nonumber \\&\quad = \sum \limits _{i=1,\ldots ,n,y_i=0} \ln \int _{-\infty }^\infty \varPhi (-A_i) \phi (\varepsilon _i)d\varepsilon _i\nonumber \\&\qquad + \sum \limits _{i=1,\ldots ,n,y_i>0} \ln \int _{-\infty }^\infty \frac{\left( \frac{\alpha }{\alpha +e^{x_i'\beta +\sigma \varepsilon _i}}\right) ^\alpha \left( \frac{e^{x_i'\beta +\sigma \varepsilon _i}}{\alpha +e^{x_i'\beta +\sigma \varepsilon _i}}\right) ^{y_i}}{1- \left( \frac{\alpha }{\alpha +e^{x_i'\beta +\sigma \varepsilon _i}}\right) ^\alpha } \varPhi (A_i) \phi (\varepsilon _i)d\varepsilon _i \nonumber \\&\qquad + \sum \limits _{i=1,\ldots ,n,y_i>0} \ln \frac{\varGamma (y_i+\alpha )}{\varGamma (\alpha )\varGamma (y_i+1)}. \end{aligned}$$

(27)

Appendix 2: Gradients of marginal and incremental effects

The calculation of variance estimates shown in (18), (19), (22), (23) involve the gradients of the marginal and incremental effects. The theoretical formulas of the gradients for our four models are focused on in this section. The corresponding estimates are obtained by substituting the estimate ${\hat{\theta }}$ for the parameter $\theta$.

1.1 Zero-inflated Poisson and zero-inflated negative binomial models with endogenous utilization

Recall that the observed marginal mean formulas for ${\hbox {ZIP}} _{\mathrm{eu}}$ and $\hbox {ZINB}_{\mathrm{eu}}$ models in (5) sare identical: $\mathrm{E}(y_i|x_i, z_i) = \int _{-\infty }^\infty e^{x_i'\beta +\sigma \varepsilon _i}\varPhi (A_i)\phi (\varepsilon _i)\,d\varepsilon _i$, and so are their marginal effects and incremental effect formulas. The difference is that the parameters in ${\hbox {ZIP}} _{\mathrm{eu}}$ models are $\theta =(\beta ',\gamma ',\rho ,\sigma )$, and those in the $\hbox {ZINB}_{\mathrm{eu}}$ models have an additional component $\alpha$: $\theta =(\beta ',\gamma ',\rho ,\sigma ,\alpha )$. The following discussions are the same for both ${\hbox {ZIP}} _{\mathrm{eu}}$ and $\hbox {ZINB}_{\mathrm{eu}}$ models and will not be emphasized again unless indicated otherwise.

For the purpose of notation simplicity, we introduce a pair of smooth functions $p(t)=e^t$ and $q(t)=\varPhi (t)$. Then, ${\dot{p}}(t)=\ddot{p}(t)=e^t$ and ${\dot{q}}(t)= \phi (t)$, $\ddot{q}(t)= -t\phi (t)$, for any $t\in \mathbb {R}$. The similar p, q notations were used in our previous publication (Liu et al. 2018).

Suppose $x_{ij}=z_{ij}$ is a continuous covariate, then we use the short-hand notations $p_i, {\dot{p}}_i,\ddot{p}_i, q_i, {\dot{q}}_i,\ddot{q}_i$ when they are evaluated at $x_i'\beta +\sigma \varepsilon _i$ and $A_i=\left( {z_i'\gamma +\rho \varepsilon _i}\right) /{\sqrt{1-\rho ^2}}$, respectively. Then, the marginal mean of count outcomes $y_i$ is $\mathrm{E}(y_i|x_i, z_i) = \int _{-\infty }^\infty p_iq_i\phi (\varepsilon _i)\,d\varepsilon _i$; the marginal effect with respect to $x_{ij}$ is

$$\begin{aligned} \eta _j(x_i,z_i,\theta )=\int _{-\infty }^\infty \left( \beta _j{\dot{p}}_iq_i+p_i{\dot{q}}_i\frac{\partial A_i}{\partial z_{ij}}\right) \phi (\varepsilon _i)\,d\varepsilon _i, \end{aligned}$$

where $\frac{\partial A_i}{\partial z_{ij}}={\gamma _{j}}/{\sqrt{1-\rho ^2}}$.

If the case is that $x_{ij}=z_{ij}$ is categorical, changing from level $l_1$ to $l_2$, then we use the notations $\{p_{2i}, {\dot{p}}_{2i},\ddot{p}_{2i}\}, \{p_{1i}, {\dot{p}}_{1i},\ddot{p}_{1i}\}, \{q_{2i}, {\dot{q}}_{2i},\ddot{q}_{2i}\}, \{q_{1i}, {\dot{q}}_{1i},\ddot{q}_{1i}\}$ for p, q and their derivatives evaluated at ${\tilde{x}}_{2i}'\beta +\sigma \varepsilon _i$, ${\tilde{x}}_{1i}'\beta +\sigma \varepsilon _i$, $A_{2i}=\left( {{\tilde{z}}_{2i}'\gamma +\rho \varepsilon _i}\right) /{\sqrt{1-\rho ^2}}$ and $A_{1i}=\left( {{\tilde{z}}_{1i}'\gamma +\rho \varepsilon _i}\right) /{\sqrt{1-\rho ^2}}$, respectively, where ${\tilde{x}}_{2i}$, ${\tilde{x}}_{1i}$, ${\tilde{z}}_{2i}$, ${\tilde{z}}_{1i}$ are the covariate vectors $x_i$ and $z_i$ with $x_{ij}=z_{ij}$ equal to $l_2$ and $l_1$, respectively. Thus, the incremental effect of marginal mean with respect to $x_{ij}=z_{ij}$ is

$$\begin{aligned} \pi _j(x_{i(-j)},z_{i(-j)},\theta )=\int _{-\infty }^\infty \left( p_{2i}q_{2i}-p_{1i}q_{1i}\right) \phi (\varepsilon _i)\,d\varepsilon _i. \end{aligned}$$

Then, we can find the derivatives of $\eta _j$ with respect to the parameters $\theta$:

$$\begin{aligned} \begin{array}{llll} &{}&{}\displaystyle \nabla _\beta \,\eta _j(x_i,z_i,\theta ) = \int _{-\infty }^\infty \left\{ \left( \beta _j\ddot{p}_iq_i+{\dot{p}}_i{\dot{q}}_i\frac{\partial A_i}{\partial z_{ij}}\right) x_i + {\dot{p}}_iq_i u_{(j+1)} \right\} \phi (\varepsilon _i)\,d\varepsilon _i,\\ &{}&{}\displaystyle \nabla _\gamma \,\eta _j(x_i,z_i,\theta ) = \int _{-\infty }^\infty \left( \beta _j {\dot{p}}_i{\dot{q}}_i\frac{\partial A_i}{\partial \gamma } + p_i\ddot{q}_i\frac{\partial A_i}{\partial \gamma } \frac{\partial A_i}{\partial z_{ij}} +p_i{\dot{q}}_i \frac{\partial ^2 A_i}{\partial z_{ij}\partial \gamma } \right) \phi (\varepsilon _i)\,d\varepsilon _i,\\ &{}&{}\displaystyle \frac{\partial }{\partial \rho } \,\eta _j(x_i,z_i,\theta ) = \int _{-\infty }^\infty \left( \beta _j {\dot{p}}_i{\dot{q}}_i\frac{\partial A_i}{\partial \rho } + p_i\ddot{q}_i\frac{\partial A_i}{\partial \rho } \frac{\partial A_{ij}}{\partial z_{ij}} +p_i{\dot{q}}_i \frac{\partial ^2 A_i}{\partial z_{ij}\partial \rho } \right) \phi (\varepsilon _i)\,d\varepsilon _i,\\ &{}&{}\displaystyle \frac{\partial }{\partial \sigma } \,\eta _j(x_i,z_i,\theta ) =\int _{-\infty }^\infty \left( \beta _j\ddot{p}_iq_i+{\dot{p}}_i{\dot{q}}_i\frac{\partial A_i}{\partial z_{ij}}\right) \phi (\varepsilon _i)\varepsilon _i\,d\varepsilon _i, \end{array} \end{aligned}$$

(28)

where $u_{(j+1)}$ is the unit vector of length $J_1$ with 1 in the $(j+1)$-st component, $\frac{\partial A_i}{\partial \gamma }={z_{i}}/{\sqrt{1-\rho ^2}}$, $\frac{\partial ^2 A_i}{\partial z_{ij}\partial \gamma }=v_{(j+1)}/\sqrt{1-\rho ^2}$ (length of the unit vector $v_{(j+1)}$ is $J_2$ and its $(j+1)$-st component is 1), $\frac{\partial A_i}{\partial \rho }=(\varepsilon _i+z_i'\gamma \rho )/\sqrt{(1-\rho ^2)^3}$, $\frac{\partial ^2 A_i}{\partial z_{ij}\partial \rho }=\gamma _{j}\rho /\sqrt{(1-\rho ^2)^3}$. The derivatives of $\pi _j$ with respect to $\theta$ are

$$\begin{aligned} \begin{array}{llll} &{}&{}\displaystyle \nabla _\beta \,\pi _j(x_{i(-j)},z_{i(-j)},\theta )=\int _{-\infty }^\infty ( {\dot{p}}_{2i}q_{2i} {\tilde{x}}_{2i}-{\dot{p}}_{1i}q_{1i} {\tilde{x}}_{1i} ) \phi (\varepsilon _i)\,d\varepsilon _i,\\ &{}&{}\displaystyle \nabla _\gamma \,\pi _j(x_{i(-j)},z_{i(-j)},\theta )=\int _{-\infty }^\infty \left( {p}_{2i}{\dot{q}}_{2i}\frac{\partial A_{2i}}{\partial \gamma } - {p}_{1i}{\dot{q}}_{1i}\frac{\partial A_{1i}}{\partial \gamma } \right) \phi (\varepsilon _i)\,d\varepsilon _i,\\ &{}&{}\displaystyle \frac{\partial }{\partial \rho } \,\pi _j(x_{i(-j)},z_{i(-j)},\theta )=\int _{-\infty }^\infty \left( {p}_{2i}{\dot{q}}_{2i}\frac{\partial A_{2i}}{\partial \rho } - {p}_{1i}{\dot{q}}_{1i}\frac{\partial A_{1i}}{\partial \rho } \right) \phi (\varepsilon _i)\,d\varepsilon _i,\\ &{}&{}\displaystyle \frac{\partial }{\partial \sigma } \,\pi _j(x_{i(-j)},z_{i(-j)},\theta )=\int _{-\infty }^\infty \left( {\dot{p}}_{2i}{q}_{2i} - {\dot{p}}_{1i}{q}_{1i} \right) \phi (\varepsilon _i)\varepsilon _i\,d\varepsilon _i, \end{array} \end{aligned}$$

(29)

where $\frac{\partial A_{2i}}{\partial \gamma }={{\tilde{z}}_{2i}}/{\sqrt{1-\rho ^2}}$, $\frac{\partial A_{1i}}{\partial \gamma }={{\tilde{z}}_{1i}}/{\sqrt{1-\rho ^2}}$, $\frac{\partial A_{2i}}{\partial \rho }=(\varepsilon _i+{\tilde{z}}_{2i}'\gamma \rho )/\sqrt{(1-\rho ^2)^3}$, $\frac{\partial A_{1i}}{\partial \rho }=(\varepsilon _i+{\tilde{z}}_{1i}'\gamma \rho )/\sqrt{(1-\rho ^2)^3}$. For $\hbox {ZINB}_{\mathrm{eu}}$ models, $\nabla _\alpha \,\eta _j(x_i,z_i,\theta )=\nabla _\alpha \,\pi _j(x_{i(-j)},z_{i(-j)},\theta ) =0$.

1.2 Hurdle Poisson models with endogenous utilization

For the convenience of computation and notations, we introduce smooth functions $p^{\mathrm{HP}}(t) =\displaystyle \frac{e^{t+e^t}}{e^{e^t}-1}$ and the same $q(t)=\varPhi (t)$ for $\hbox {HP}_{\mathrm{eu}}$ models as used for the ${\hbox {ZIP}} _{\mathrm{eu}}$ and $\hbox {ZINB}_{\mathrm{eu}}$ models. Then

$$\begin{aligned} p^\mathrm{HP}(t) = e^t +\kappa (t),\quad {\dot{p}}^\mathrm{HP}(t) = e^t +{\dot{\kappa }}(t),\quad \ddot{p}^\mathrm{HP}(t) = e^t +\ddot{\kappa }(t), \end{aligned}$$

where $\kappa (t)=\frac{e^t}{e^{e^t}-1}$, ${\dot{\kappa }}(t)=\kappa (t)\{1-e^t-\kappa (t)\}$, $\ddot{\kappa }(t)={\dot{\kappa }}(t)\{1-e^t-2\kappa (t)\} -e^t\kappa (t)$. Then, for $p^\mathrm{HP}$ and q, we use the similar notations $p^\mathrm{HP}_i,q_i,p^\mathrm{HP}_{2i}, q_{2i}, p^\mathrm{HP}_{1i}, q_{1i}$ and their derivatives when evaluated at $x_i'\beta +\sigma \varepsilon _i$ and $A_i=\left( {z_i'\gamma +\rho \varepsilon _i}\right) /{\sqrt{1-\rho ^2}}$ and when $x_{ij}$ takes the values $l_2$ or $l_1$ for the i-th observation, as we have done for the ${\hbox {ZIP}} _{\mathrm{eu}}$ and $\hbox {ZINB}_{\mathrm{eu}}$ models. Therefore, the marginal mean of $\hbox {HP}_{\mathrm{eu}}$ model in (10) can be rewritten as $\mathrm{E}(y_i|x_i,z_i)= \int _{-\infty }^\infty p_i^\mathrm{HP}q_i\,\phi (\varepsilon _i)\,d\varepsilon _i$. Consequently, the marginal effect is

$$\begin{aligned} \eta _j(x_i,z_i,\theta )=\int _{-\infty }^\infty \left( \beta _j{\dot{p}}_i^\mathrm{HP}q_i+p_i^\mathrm{HP}{\dot{q}}_i\frac{\partial A_i}{\partial z_{ij}}\right) \phi (\varepsilon _i)\,d\varepsilon _i, \end{aligned}$$

and the incremental effect with respect to categorical $x_{ij}=z_{ij}$ from $l_1$ to $l_2$ is

$$\begin{aligned} \pi _j(x_{i(-j)},z_{i(-j)},\theta )=\int _{-\infty }^\infty \left( p_{2i}^\mathrm{HP}q_{2i}-p_{1i}^\mathrm{HP}q_{1i}\right) \phi (\varepsilon _i)\,d\varepsilon _i, \end{aligned}$$

where $\theta =(\beta ',\gamma ',\rho , \sigma )$. The remaining lines for the derivatives of $\eta _j$ and $\pi _j$ with respect to $\theta$ have the same forms as (28), (29) in the ${\hbox {ZIP}} _{\mathrm{eu}}$ and $\hbox {ZINB}_{\mathrm{eu}}$ models, with p and its derivatives replaced with $p^{\mathrm{HP}}$ and its derivatives.

1.3 Hurdle negative binomial models with endogenous utilization

For $\hbox {HNB}_{\mathrm{eu}}$ models, our parameter vector is $\theta = (\beta ',\gamma ,\rho ,\sigma ,\alpha )$, and our auxiliary functions are $p^\mathrm{HNB}(t,\alpha ) =\displaystyle \frac{e^t}{1-\delta (t,\alpha )}$, and $q(t)=\varPhi (t)$, where $\delta (t,\alpha ) = \tau ^\alpha (t,\alpha )$, $\tau (t,\alpha ) =\displaystyle \frac{\alpha }{\alpha + e^t}$, and $\alpha >0$.

It is tedious but not hard to get the partial derivatives of $p^\mathrm{HNB}$ with respect to t and $\alpha$. For simplicity of notations, let’s drop the function arguments of $p^\mathrm{HNB}, \delta ,\tau$. Then,

$$\begin{aligned}&{\dot{p}}^\mathrm{HNB}_t(t,\alpha ) = p^\mathrm{HNB}\left( 1-p^\mathrm{HNB}\delta \tau \right) ,\\&\ddot{p}^\mathrm{HNB}_{tt}(t,\alpha ) ={{\dot{p}}^\mathrm{HNB}_t \left( 1-2p^\mathrm{HNB}\delta \tau \right) + (\alpha + 1)\left( p^\mathrm{HNB}\right) ^2 \delta \tau (1-\tau )},\\&\displaystyle {\dot{p}}^\mathrm{HNB}_\alpha (t,\alpha ) = \frac{p^\mathrm{HNB}\delta (\ln \tau + 1-\tau )}{(1-\delta )},\\&\ddot{p}^\mathrm{HNB}_{t\alpha }(t,\alpha ) =\big \{\displaystyle {\dot{p}}^\mathrm{HNB}_\alpha (1-p^\mathrm{HNB}\delta \tau -p^\mathrm{HNB}\tau )\big \}-\frac{(p^\mathrm{HNB})^2\delta e^t}{(\alpha + e^t)^2}, \end{aligned}$$

where we used ${\dot{\tau }}_t(t,\alpha ) = \tau (\tau -1)$, $\ddot{\tau }_{tt}(t,\alpha ) = \tau (\tau -1)(2\tau -1)$, ${\dot{\tau }}_\alpha (t,\alpha ) = {e^t}/{(\alpha + e^t)^2}$, ${\dot{\delta }}_t(t,\alpha ) = \alpha \delta (\tau -1) = -e^t\delta \tau$, $\ddot{\delta }_{tt}(t,\alpha ) = \alpha \delta (\tau -1)\{(\alpha + 1)\tau -\alpha \} = \delta \tau ^2e^t(e^t-1)$, ${\dot{\delta }}_\alpha (t,\alpha ) = \delta (\ln \tau + 1-\tau )$.

At the i-th observation, the function $p^\mathrm{HNB}$ and its partial derivatives evaluated at $(x_i'\beta +\sigma \varepsilon _i,\alpha )$ are denoted by $p_i^\mathrm{HNB}, {\dot{p}}_{ti}^\mathrm{HNB}, \ddot{p}_{tti}^\mathrm{HNB}, {\dot{p}}^\mathrm{HNB}_{\alpha i}, \ddot{p}^\mathrm{HNB}_{t\alpha i}$. When the covariate $x_{ij}$ is fixed at $l_1$ or $l_2$, the values of $p^\mathrm{HNB}, {\dot{p}}_t^\mathrm{HNB}, {\dot{p}}^\mathrm{HNB}_{\alpha }$ are denoted by $p^\mathrm{HNB}_{2i}, {\dot{p}}_{2ti}^\mathrm{HNB}, {\dot{p}}^\mathrm{HNB}_{2\alpha i}$, and $p^\mathrm{HNB}_{1i}, {\dot{p}}_{1ti}^\mathrm{HNB}, {\dot{p}}^\mathrm{HNB}_{1\alpha i}$, respectively. Notations $q_i,{\dot{q}}_i,\ddot{q}_i,q_{2i}, q_{1i},{\dot{q}}_{2i}, {\dot{q}}_{1i}$ for the q function are the same as those in the zero-inflated models and hurdle Poisson models.

Hence, the marginal mean of $\hbox {HNB}_{\mathrm{eu}}$ model in (12) can be rewritten as $\mathrm{E}(y_i|x_i,z_i)= \int _{-\infty }^\infty p_i^\mathrm{HNB}q_i\,\phi (\varepsilon _i)\,d\varepsilon _i$, yielding the marginal effect with respect to continuous $x_{ij}=z_{ij}$, i.e.,

$$\begin{aligned} \eta _j(x_i,z_i,\theta )=\int _{-\infty }^\infty \left( \beta _j{\dot{p}}_{ti}^\mathrm{HNB}q_i+p_i^\mathrm{HNB}{\dot{q}}_i\frac{\partial A_i}{\partial z_{ij}}\right) \phi (\varepsilon _i)\,d\varepsilon _i \end{aligned}$$

and the incremental effect with respect to $x_{ij}=z_{ij}$ from $l_1$ to $l_2$, i.e.,

$$\begin{aligned} \pi _j(x_{i(-j)},z_{i(-j)},\theta )=\int _{-\infty }^\infty \left( p_{2i}^\mathrm{HNB}q_{2i}-p_{1i}^\mathrm{HNB}q_{1i}\right) \phi (\varepsilon _i)\,d\varepsilon _i. \end{aligned}$$

Then, the derivatives of $\eta _j$ with respect to the parameter $\theta$ are

$$\begin{aligned}&\nabla _\beta \,\eta _j(x_i,z_i,\theta ) \\&\quad =\int _{-\infty }^\infty \left\{ \left( \beta _j\ddot{p}_{tti}^\mathrm{HNB}q_i+{\dot{p}}_{ti}^\mathrm{HNB}{\dot{q}}_i\frac{\partial A_i}{\partial z_{ij}}\right) x_i + {\dot{p}}_{ti}^\mathrm{HNB}q_i u_{(j+1)} \right\} \phi (\varepsilon _i)\,d\varepsilon _i,\\&\qquad \nabla _\gamma \,\eta _j(x_i,z_i,\theta ) = \int _{-\infty }^\infty \left( \beta _j {\dot{p}}_{ti}^\mathrm{HNB}{\dot{q}}_i\frac{\partial A_i}{\partial \gamma } + p_i^\mathrm{HNB}\ddot{q}_i\frac{\partial A_i}{\partial \gamma } \frac{\partial A_i}{\partial z_{ij}} +p^\mathrm{HNB}_i{\dot{q}}_i \frac{\partial ^2 A_i}{\partial z_{ij}\partial \gamma } \right) \phi (\varepsilon _i)\,d\varepsilon _i,\\&\qquad \frac{\partial }{\partial \rho } \,\eta _j(x_i,z_i,\theta ) = \int _{-\infty }^\infty \left( \beta _j {\dot{p}}_{ti}^\mathrm{HNB}{\dot{q}}_i\frac{\partial A_i}{\partial \rho } + p_i^\mathrm{HNB}\ddot{q}_i\frac{\partial A_i}{\partial \rho } \frac{\partial A_{ij}}{\partial z_{ij}} +p_i^\mathrm{HNB}{\dot{q}}_i \frac{\partial ^2 A_i}{\partial z_{ij}\partial \rho } \right) \phi (\varepsilon _i)\,d\varepsilon _i,\\&\qquad \frac{\partial }{\partial \sigma } \,\eta _j(x_i,z_i,\theta ) =\int _{-\infty }^\infty \left( \beta _j\ddot{p}_{tti}^\mathrm{HNB}q_i+{\dot{p}}_{ti}^\mathrm{HNB}{\dot{q}}_i\frac{\partial A_i}{\partial z_{ij}}\right) \phi (\varepsilon _i)\varepsilon _i\,d\varepsilon _i,\\&\qquad \frac{\partial }{\partial \alpha } \,\eta _j(x_i,z_i,\theta ) =\int _{-\infty }^\infty \left( \beta _j\ddot{p}_{t\alpha i}^\mathrm{HNB}q_i+{\dot{p}}_{\alpha i}^\mathrm{HNB}{\dot{q}}_i\frac{\partial A_i}{\partial z_{ij}}\right) \phi (\varepsilon _i)\,d\varepsilon _i, \end{aligned}$$

where $u_{(j+1)}$ is the unit vector of length $J_1$. The derivatives of $\pi _j$ with respect to $\theta$ are

$$\begin{aligned}&\nabla _\beta \,\pi _j(x_{i(-j)},z_{i(-j)},\theta )=\int _{-\infty }^\infty ( {\dot{p}}_{2ti}^\mathrm{HNB}q_{2i} {\tilde{x}}_{2i}-{\dot{p}}_{1ti}^\mathrm{HNB}q_{1i} {\tilde{x}}_{1i} ) \phi (\varepsilon _i)\,d\varepsilon _i,\\&\qquad \nabla _\gamma \,\pi _j(x_{i(-j)},z_{i(-j)},\theta )=\int _{-\infty }^\infty \left( {p}_{2i}^\mathrm{HNB}{\dot{q}}_{2i}\frac{\partial A_{2i}}{\partial \gamma } - {p}_{1i}^\mathrm{HNB}{\dot{q}}_{1i}\frac{\partial A_{1i}}{\partial \gamma } \right) \phi (\varepsilon _i)\,d\varepsilon _i,\\&\qquad \frac{\partial }{\partial \rho } \,\pi _j(x_{i(-j)},z_{i(-j)},\theta )=\int _{-\infty }^\infty \left( {p}_{2i}^\mathrm{HNB}{\dot{q}}_{2i}\frac{\partial A_{2i}}{\partial \rho } - {p}_{1i}^\mathrm{HNB}{\dot{q}}_{1i}\frac{\partial A_{1i}}{\partial \rho } \right) \phi (\varepsilon _i)\,d\varepsilon _i,\\&\qquad \frac{\partial }{\partial \sigma } \,\pi _j(x_{i(-j)},z_{i(-j)},\theta )=\int _{-\infty }^\infty \left( {\dot{p}}_{2ti}^\mathrm{HNB}{q}_{2i} - {\dot{p}}_{1ti}^\mathrm{HNB}{q}_{1i} \right) \phi (\varepsilon _i)\varepsilon _i\,d\varepsilon _i,\\&\qquad \frac{\partial }{\partial \alpha } \,\pi _j(x_{i(-j)},z_{i(-j)},\theta )=\int _{-\infty }^\infty \left( {\dot{p}}_{2\alpha i}^\mathrm{HNB}{q}_{2i} - {\dot{p}}_{1\alpha i}^\mathrm{HNB}{q}_{1i} \right) \phi (\varepsilon _i)\,d\varepsilon _i, \end{aligned}$$

where ${\tilde{x}}_{2i},{\tilde{x}}_{1i}$, ${\tilde{z}}_{2i},{\tilde{z}}_{1i}$ are the covariate vectors $x_i$ and $z_i$ with $x_{ij}=z_{ij}$ being equal to $l_2$ and $l_1$, respectively. Formulas of derivatives of $A_i,A_{2i},A_{1i}$ are the same as those for the ${\hbox {ZIP}} _{\mathrm{eu}}$ and $\hbox {ZINB}_{\mathrm{eu}}$ models.

Appendix 3: Box plots of the parameter estimates in the simulation study

The box plots of parameter estimates of the four models with endogenous utilization are given in Fig. 1.

Fig. 1

The box plots of parameter estimates of the four models in the simulation study. The black horizontal line is the true value

Appendix 4: Parameter estimates of the empirical study

This section presents the three tables that summarize the parameter estimates of the two-part models with endogenous utilization in Table (4) and without endogeneity in Table (5) and the estimates of average marginal effect and average incremental effect for the latter (Table 6) for the GSOEP data.

Table 4 Parameter estimates (SEs) given by the four two-part models with endogeneity from fitting to the GSOEP data.

Table 5 Parameter estimates (SEs) given by the four standard two-part models with no endogenous utilization from fitting to the GSOEP data.

Table 6 Estimates (and SEs, p values) of average marginal effect with respect to “hhincome”, average incremental effect with respect to “public” and “addon”, the AIC, BIC of the four standard two-part models with no endogenous utilization from fitting to the GSOEP data