2.1 Data pre-processing

Prior to carrying out any modeling, two data pre-processing steps were undertaken, the details of which are outlined below.

1. (1) One-hot encoding: for all categorical variables with more than two levels, one-hot encoding was applied. One-hot encoding represents a categorical variable with $l$ categories $c_1, c_2, \ldots, c_l$ using a $l$ dimensional feature vector which is of the form

   (1) \begin{equation} x_j \mapsto \Big (\unicode{x1D7D9}_{\{x_j = c_1\}},\ldots, \unicode{x1D7D9}_{\{x_j = c_{l}\}} \Big )^T \in \mathbb{R}^{l}. \end{equation}

2. (2) Min-max scaler: a min-max scaler which transforms the variable to a $[-1,1]$ scale was implemented for the numerical and binary variables. The transformation was applied using the formula

   (2) \begin{equation} x_j \mapsto x_j^{*} = \frac{2(x_j - m_j)}{M_j - m_j} -1\in [-1, 1] \end{equation}
   where $m_j$ and $M_j$ represent the minimum and maximum values of variable $x_j$ .

A more detailed description of the application of both one-hot encoding and min-max scaler on the admission data set under consideration is given in Jose et al. (Reference Jose, Macdonald, Tzougas and Streftaris2022). Following the data pre-processing steps, a 90:10 random split of the entire data set was created to be used as the learning $\mathcal{D}$ and testing $\mathcal{T}$ data sets. All the models were fitted using the learning data set, and the performance of the models on the testing data set was used to compare their performances.

3.1 Zero-inflated Poisson regression

As mentioned earlier, zero-inflated Poisson models are based on a mixture distribution comprising two components (Lambert, Reference Lambert1992; Cameron & Trivedi, Reference Cameron and Trivedi2013). In practice, a count distribution such as the Poisson or negative binomial distribution is used to represent the count component, and a Bernoulli distribution is assumed for the zero component. Hence a zero-inflated Poisson (ZIP) regression model is of the form

(3) \begin{equation} Y_i\sim \begin{cases} 0& \text{with probability}\; \pi _i\\ Poisson(\lambda _i e_i)&\text{with probability}\; (1-\pi _i)\qquad \ldots \;\text{for}\;\; i = 1, \ldots, n\\ \end{cases} \end{equation}

with $\mu _i=\lambda _ie_i$ being the mean of the Poisson part of the ith record with exposure $e_i$ and rate parameter $\lambda _i$ , while $\pi _i$ represents the probability of only having zero admissions. Zeros arise from both the zero component and the count component, with the two components having probability $\pi _i$ and $(1-\pi _i)$ , respectively. Hence, the probability mass function of the ZIP mixture distribution is

(4) \begin{equation} Pr(Y_i=y_i)= \begin{cases} \pi _i + (1-\pi _i) e^{-\mu _i}, & y_i = 0\\ (1-\pi _i) \frac{e^{-\mu _i}\mu _i^{y_i} }{y_i!} & y_i \gt 0 \\ \end{cases} \end{equation}

and the corresponding mean and variance are given by

(5) \begin{equation} E(Y_i) = (1-\pi _i)\mu _i,\qquad V(Y_i) = (1-\pi _i)\mu _i(1+\pi _i\mu _i). \end{equation}

The ZIP model allows both $\mu _i$ and $\pi _i$ to be modeled using a set of covariates. This would be of the form :

(6) \begin{equation} \log (\mu _i) = o_i + \beta _{0} + \boldsymbol{\beta }_{reg}^{\top }\boldsymbol{x_i} = o_i + \beta _0 + \langle \boldsymbol{\beta }_{reg}, \boldsymbol{x_i}\rangle \end{equation}

and

(7) \begin{equation} logit(\pi _i) = \gamma _{0} + \boldsymbol{\gamma ^\top w_i} \end{equation}

where $o_i=\log (e_i)$ is the offset term, $\beta _{0},\gamma _{0}$ the intercept terms and $\boldsymbol{\{\beta }_{reg}^\top,\boldsymbol{\gamma ^\top \}} = (\beta _{1},\ldots, \beta _{q}, \gamma _{1},\ldots, \gamma _{s})$ being the unknown vector of coefficients to be estimated corresponding to the sets of covariates $\boldsymbol{x}_i=(x_{i,1},\ldots,x_{i,q})^ \top$ and $\boldsymbol{w}_i=(\boldsymbol{w}_{i,1},\ldots,\boldsymbol{w}_{i,s})^ \top$ considered in the two regression functions with dimensions $q\times 1$ and $s\times 1$ , respectively. Regarding the treatment of exposure within the ZIP model, it is possible that the exposure or period at risk could potentially influence the probability of admission, as well as the rates of admission. Consequently, the exposure could be factored into both components, as detailed in Feng (Reference Feng2022). Here, we follow the general practice of treating it as an offset term in the regression function for $\mu$ , as discussed in Lee et al. (Reference Lee, Wang and Yau2001).

The model selection procedure was carried out for both count and zero components of the ZIP regression model. All variables entered the model for the count component, and thus, a full model with all the covariates was utilized for the Poisson regression function. Similarly, for the regression function associated with the zero component, a forward step-wise variable selection process based on the Bayesian information criterion (BIC) was carried out(Vrieze, Reference Vrieze2012). As we are using logistic regression to model the probability parameter, we refer to it as the logistic component from here on. The summary of the variable selection process of the logistic component of the model is given below in Table 3. Under the forward step-wise variable selection process, we start from a model without any terms in the logistic component and consider adding one variable at a time. The outcome of any particular step is the model that yields the lowest BIC value, and the whole process is repeated until there is no further reduction in the BIC value. The model thus identified from the variable selection process had the vector of coefficients $\{\beta _{0}, \boldsymbol{\beta }_{reg},\gamma _{0}, \boldsymbol{\gamma \}} = (\beta _{0}, \beta _{1},\ldots, \beta _{q}, \gamma _{0}, \gamma _{1})$ corresponding to the complete set of covariates as presented in Table 4. The ZIP regression model was fitted using the zeroinfl() function in the pscl package which employs the optim() function to estimate the model parameters in both components simultaneously, by maximum likelihood, using optimization algorithms such as Nelder-Mead or a quasi-Newton method (Zeileis et al., Reference Zeileis, Kleiber and Jackman2008).

Table 3. Summary of the variable selection process of the logistic component of the ZIP regression model

Table 4. List of covariates and the corresponding coefficient parameters in both the components of the ZIP regression model

3.2 Zero-inflated Poisson neural network (ZIPNN) model

A generic feed forward NN comprises an input layer, followed by multiple hidden layers and then the output layer. The structure of a feed forward NN is first discussed using a NN model with an underlying Poisson distribution assumption (model 3 in Table 5). A feature space $\mathcal{X}$ is taken as the input layer with dimension $q_0$ . Assuming a network architecture of $d\in \mathbb{N}$ hidden layers with $q_m\in \mathbb{N}, 1\leq m\leq d$ neurons in each of the layers then a neuron $ z_{j}^{(m)},1 \leq j \leq q_{m},$ in the $m^{th}$ hidden layer $\textbf{z}^{(m)}$ is given by

(8) \begin{equation} z_{j}^{(m)}(\textbf{z})=\psi \left ( \langle \boldsymbol{\beta }_{j}^{(m)},\textbf{z} \rangle,\right ) \end{equation}

with $\textbf{z}^{(m)}$ represented as

(9) \begin{equation} \textbf{z}^{(m)}:\, \mathbb{R}^{q_{m-1}} \rightarrow \mathbb{R}^{q_{m}}, \quad \textbf{z}\mapsto \textbf{z}^{(m)}(\textbf{z})=(1,z_{1}^{(m)}(\textbf{z}),\dots,z_{q_{m}}^{(m)}(\textbf{z}))^\top, \end{equation}

inclusive of the intercept component, where $\boldsymbol{\beta }_{j}^{(m)}= (\beta _{l,j}^{(m)})^\top _{0\leq l \leq q_{m-1}} \in \mathbb{R}^{q_{m-1}+1}$ are the network parameters and $\psi :\, \mathbb{R} \rightarrow \mathbb{R},$ the activation function. The network parameters corresponding to the hidden layer $\textbf{z}^{(m)}$ are the $\left (\boldsymbol{\beta }_{1}^{(m)},\ldots,\boldsymbol{\beta }_{q_{m}}^{(m)}\right ) \in \mathbb{R}^{q_{m}}$ . The comprehensive set of network parameter is given by $\boldsymbol{\beta }=\left (\boldsymbol{\beta }_{1}^{(1)},\ldots,\boldsymbol{\beta }_{q_{d}}^{(d)}, \boldsymbol{\beta }^{(d+1)}\right ) \in \mathbb{R}^{r}$ , where the dimension is $r=\sum _{m=1}^{d} q_{m}\left (q_{m-1}+1\right ) + (q_{d}+1)$ . The predictor of the NN obtained from the output layer is then of the form

(10) \begin{equation} (o_i, \boldsymbol{x}_i) \mapsto \log(\mu _i^{NN})= o_i +\langle \boldsymbol{\beta }^{(d+1)}, (\textbf{z}^{(d)} \circ \dots \circ \textbf{z}^{(1)}) (\boldsymbol{x}_{i})\rangle, \end{equation}

for $i=1,\ldots,n$ , where $\boldsymbol{\beta }^{(d+1)} \in \mathbb{R}^{q_{d}+1}$ are the weights associated with the output layer which connects the last hidden layer $\textbf{z}^{d}$ to the output layer $\mathbb{R}_{+}$ .

Table 5. Predictive performance of Poisson and ZIP regression models and network-based models with and without an additional layer for interpretation based on NLL. The empirical mean of the observed data is 0.0027

For a ZIPNN model, the feature space $\mathcal{X}$ with dimension $q$ is taken as the input layer, i.e., $q_0=q$ . The construction of ZIPNN models involves incorporating zero-inflation considerations through the establishment of a distribution layer within the NN. This distribution layer, referred to as a lambda layerFootnote 2 in machine learning terminology, is a component in a NN that facilitates the integration of custom operations by defining a function applied to input data. In particular, this functionality allows for the inclusion of specific computations not covered by standard layers. In the specific context of zero-inflated NNs, the distribution lambda layer is employed to integrate a zero-inflated distribution assumption into the NN, enabling the model to capture complex relationships in the data. Also, we take $q_d=2$ with $\psi :\, f(x)=x$ . In other words, a dense layer with two neurons without activation is defined before the distribution layer.

The output of these two neurons is then used to calculate the probability $(\pi )$ and rate $(\lambda )$ parameters of the zero-inflated Poisson distribution constructed using the distribution layer that follows. The distribution layer thus created has an underlying zero-inflated Poisson mixture distribution. The model output is then derived as the mean of the distribution given by

(11) \begin{equation} E(Y_i) = (1-\pi _i^{zipnn})\mu ^{zipnn}_i \end{equation}

where $\pi _i^{zipnn}$ and $\mu ^{zipnn}_i$ are the rate and probability parameters of the underlying zero-inflated Poisson distribution. As in the case of a ZIP model, the exposure is taken as an offset term in the count component of the ZIPNN model. An exponential transformation is applied to the sum of the offset term and the output from the neuron associated with $\lambda$ . A sigmoid function is applied to the output from the neuron associated with $\pi$ to restrict the value of $\pi$ to $[0,1]$ interval. The $\mu _i^{zipnn}$ is then given by

(12) \begin{equation} \log (\mu _i^{zipnn}) = z_{1}^{(d)}(\textbf{z})=\psi \left ( \langle \boldsymbol{\beta }_{1}^{(d)},\textbf{z} \rangle \right ) +o_i = \langle \boldsymbol{\beta }_{1}^{(d)},\textbf{z} \rangle + o_i \end{equation}

or equivalently

(13) \begin{equation} \log (\mu _i^{zipnn}) = o_i +\langle \boldsymbol{\beta }_{1}^{(d)}, (\textbf{z}^{(d-1)} \circ \dots \circ \textbf{z}^{(1)}) (\boldsymbol{x}_{i})\rangle. \end{equation}

Similarly, $\pi _i^{zipnn}$ is of the form

(14) \begin{equation} logit(\pi _i^{zipnn}) = z_{2}^{(d)}(\textbf{z})=\psi \left ( \langle \boldsymbol{\beta }_{2}^{(d)},\textbf{z} \rangle \right ) = \langle \boldsymbol{\beta }_{2}^{(d)},\textbf{z} \rangle. \end{equation}

$\boldsymbol{\beta }_{1}^{(d)}$ and $ \boldsymbol{\beta }_{2}^{(d)}$ represent the network weights associated with the two neurons in the last hidden layer corresponding to probability and rate parameters. In essence the Equations (6) and (7) are replaced by Equations (13) and (14). Alternatively, we could also define

(15) \begin{equation} logit(1-\pi _i^{zipnn}) = z_{2}^{(d)}(\textbf{z})=\psi \left ( \langle \boldsymbol{\beta }_{2}^{(d)},\textbf{z} \rangle \right ) = \langle \boldsymbol{\beta }_{2}^{(d)},\textbf{z} \rangle. \end{equation}

Since the zero-inflated distribution is defined as a mixture distribution using the distribution layer, the models could also be constructed using the $p_i= (1-\pi _i)$ parameter with ease. This means that instead of $\pi$ , the probability of zero component, the distribution layer considers p, which is the probability of count component or, in other words, the probability of having a non-zero admission. The only difference is how the mixture distribution is defined. For ease of interpretation, while using the attention layers, we have adopted the latter approach. A schematic representation of a ZIPNN with 20,15,10,2 neurons (model 5 in Table 5) in each dense layer is shown in Fig. 1. Code for implementing the same is given in Listing 3 and additional details regarding the approach can be found in Dürr et al. (Reference Dürr, Sick and Murina2020). The structure of the ZIPNN model in terms of the connection between layers, the shape of input and out of each layer, and the number of parameters is shown in Listing 1. The final layer of the model is the distribution layer, which follows a zero-inflated Poisson distribution. The mean of the distribution is taken as the output or the response, which, in our case, is the admission count.

Figure 1 An illustration of a sample ZIPNN model with three hidden layers and 20,15,10,2 neurons in each layer.

Listing A1. Structure of the ZIPCANN model.

3.3 Zero-inflated Poisson combined actuarial neural network (ZIPCANN) model

The ZIPCANN model contains additional skip connections, as in the case of a CANN model. In a CANN model (model 4 in Table 5), a skip connection connects the input features directly to the output layer. In effect, the predictor of a CANN model contains an additional regression function compared to the predictor of a NN and is of the form

(16) \begin{equation} (o_i, \boldsymbol{x}_i) \mapsto \log(\mu ^{CANN}_i)= o_i + \langle \boldsymbol{\beta }^{reg}, \boldsymbol{x}_i \rangle + \langle \boldsymbol{\beta }^{(d+1)}, (\textbf{z}^{(d)} \circ \dots \circ \textbf{z}^{(1)}) (\boldsymbol{x}_{i})\rangle, \end{equation}

with the parameters from the regression function or skip connection represented by $\boldsymbol{\beta }^{reg}$ (Schelldorfer & Wuthrich, Reference Schelldorfer and Wuthrich2019; Jose et al., Reference Jose, Macdonald, Tzougas and Streftaris2022). As detailed in Schelldorfer & Wuthrich (Reference Schelldorfer and Wuthrich2019), various versions of CANNs exist depending on whether the weights in the regression part are updated or not whilst training the model. More specifically, one can either estimate the weights in the skip connection during training, or alternatively, keep the weights of the regression component fixed as the iterated weighted least squares estimates from the corresponding regression model. In this work, we adopt the former approach for all models containing skip connections.

For the ZIPCANN model, two skip connections are used; one for the rate parameter and one for the probability parameter. Instead of the output layer, the skip connection is made to the two neurons in the last hidden layer. In order to be consistent with the ZIP model, only the age variable was used in the skip connection associated with $\pi$ , whereas all the features were used in the skip connection for $\mu$ . Consequently, the predictor for the ZIPCANN model is of the form

(17) \begin{equation} E(Y_i) = (1-\pi _i^{zipcann})\mu ^{zipcann}_i \end{equation}

where $\mu ^{zipcann}_i$ is given by

(18) \begin{equation} \log (\mu _i^{zipcann}) = \langle \boldsymbol{\beta }^{reg}, \boldsymbol{x}_i \rangle + \langle \boldsymbol{\beta }_{1}^{(d)},\textbf{z} \rangle + o_i \end{equation}

and $\pi _i^{zipcann}$ is given by

(19) \begin{equation} logit(\pi _i^{zipcann}) = \gamma _{age} x^{age}_i+ \langle \boldsymbol{\beta }_{2}^{(d)},\textbf{z} \rangle \end{equation}

or alternatively,

(20) \begin{equation} logit(1-\pi _i^{zipcann}) = \gamma _{age} x^{age}_i+ \langle \boldsymbol{\beta }_{2}^{(d)},\textbf{z} \rangle. \end{equation}

Similarly to the ZIPNN model, the final output layer of the ZIPCANN model is also a distribution layer. The code for fitting the ZIPCANN model is given in Listing 4, and a diagrammatic representation of a sample model with one-hot encoding, skip connections, and 20,15,10,2 neurons in each of the layers (model 6 in Table 5), is shown in Fig. 2. The structure of the ZIPCANN model is given in Listing 5.

Figure 2 An illustration of a sample ZIPCANN model with skip connections, one-hot encoding, and 20,15,10,2 neurons in each of the layers.

4.1 Hyper-parameters

The fitting of NN models requires determining a number of associated hyper-parameters, including the depth of the network, choice of activation function, loss function, and gradient descent method. The hyper-parameter options identified and adopted in Jose et al. (Reference Jose, Macdonald, Tzougas and Streftaris2022) have also been used here, with details given below.

• • Gradient descent method (GDM): the Nesterov-accelerated adaptive moment estimation (Nadam) method was used as the choice of the Gradient descent optimization algorithm for estimating the model weights.

• • Network architecture: for NN and CANN models, a network architecture with three hidden layers with (20,15,10) neurons in each layer was considered. As mentioned earlier, the ZIPNN and ZIPCANN models have an additional layer with two neurons. In addition to the layers mentioned above, an attention layer is used to interpret the results from NN models (see Section 5). Also, for the ZIPNN and ZIPCANN models, alternate architectures with different numbers of neurons in the first three hidden layers were considered. The results, as given in Table A.1, indicate that the choice of (20,15,10) architecture has the best predictive performance compared to other architectures in the case of the ZIPNN model. Moreover, for the ZIPCANN model, the (35,25,20) architecture only minimally outperforms the (20,15,10) architecture. Hence, the combination of (20,15,10) neurons in the initial three hidden layers is adopted for both models.

• • Batch size and epochs: a combination of 30,000 batch size and 500 epochs were used for fitting the NN models. The 30,000 batch size was identified as a reasonable choice in earlier work by Jose et al. (Reference Jose, Macdonald, Tzougas and Streftaris2022), and the large epoch was used since the early stopping approach implemented using the callback method will restrict the model from overfitting.

• • Validation split: a further 80:20 split of the learning set was used as the training set, $\mathcal{D}^{(-)}$ , and a validation data set $\mathcal{V}.$

• • Loss function: the negative log-likelihood (NLL) was used as the objective function, which the GDM algorithm minimizes for estimating the model weights. The log-likelihood of ZIP mixture distribution is given by

  (21) \begin{equation} \begin{aligned} l(\pi,\mu ;y) = \sum _{i=1}^n &u_i \ln \big (\pi _i + (1-\pi _i) \exp (-\mu _i)\big )\\ &+ (1-u_i)\big [ \ln (1-\pi _i) - \mu _i + y_i\ln (\mu _i)- \ln{(y_i!)} \big ], \end{aligned} \end{equation}

where $u_i = I(y_i = 0)$ , $\mu _i$ is the mean of the Poisson component, $\pi _i$ the probability parameter and $y_i$ is the response variable. It is worth noting that in the case of the zero-inflated models, the saturated model reduces to the saturated version of the count component (Martin & Hall, Reference Martin and Hall2016). This follows from the fact that a saturated model has the mean responses equivalent to the data, i.e., $E(y_i)= y_i$ . For the zero-inflated models, this is true when we set $\mu _i = y_i$ and replace $\pi _i$ with $u_i= I(y_i=0)$ since $E(y_i)=(1-\pi _i)\mu _i$ . Hence, the saturated model for ZIP mixture distribution will be a saturated Poisson distribution.

As previously mentioned, for the ZIPNN and ZIPCANN models, the output layer is a distribution layer (see Listing 1 and A5) with underlying Zero-inflated Poisson distribution and the mean of the distribution is taken as the model outcome or response which is the admission number. More precisely the distribution layer is considered as the model output. The NLL is the negative sum of the logarithm of the likelihood assigned to the observed count by the conditional probability distribution ( $-\sum log(P(y|x,w))$ ). As the model’s output is the distribution layer, we can obtain this by defining the loss function as shown in Listing 2.

Listing A2. Negative log-likelihood loss function used for training zero-inflated neural network models.

Additionally, the early stopping and dropout model improvement approaches, as detailed in Jose et al. (Reference Jose, Macdonald, Tzougas and Streftaris2022), were used to avoid overfitting for all the network-based models.

5.1 Interpreting the ZIPNN model

In order to interpret the results obtained from ZIPNN model, we extend the additive decomposition assumption to the $\mu$ and $p=(1-\pi )$ parameters of the underlying ZIP mixture distribution as shown in Equation (4). The regression weights corresponding to each feature are calculated separately for $\mu$ and $p$ . Hence, the dimension of the layer analogous to the LocalGLM layer with the component-wise product of regression weights and features is $2\times q_0$ . Among these $2\times q_0$ weights, $q_0$ weights are associated with $\mu$ and the rest with $p$ .

5.1.1 Regression attention for ZIPNN

The ZIPNN model, as defined in Section 3.2, consists of an input layer, multiple hidden layers, with the last hidden layer containing two neurons, followed by the distribution layer and the output layer. To create the interpretable ZIPNN model, the attention layer is introduced before the distribution layer, effectively replacing the last hidden layer of the ZIPNN model with two neurons (see Fig. 6).

Figure 6 An illustration of an interpretable ZIPNN model with regression attention, one-hot encoding, and 20,15,10 neurons in each of the layers.

Hence, the rate and probability parameters defined in Equations (13) and (14) are replaced by

(24) \begin{equation} \log (\mu _i^{zipnn}) = \beta _{0\mu } + \langle \boldsymbol{\beta _{(q_0+1:2q_0)}(x}_i),\boldsymbol{x_i}\rangle + o_i \end{equation}

and

(25) \begin{equation} logit(p_i^{zipnn}) = \beta _{0\pi } + \langle \boldsymbol{\beta _{(1:q_0)}(x}_i),\boldsymbol{x_i}\rangle \end{equation}

with $\boldsymbol{x_i}\mapsto \boldsymbol{\beta (x_i) = z^{(d:1)}(x_i) = \big (z^{(1)}\circ \ldots .\circ z^{(d)}\big )(x_i)}$ . The attention weights $\boldsymbol{\beta _{(1:q_0)}(x})$ and $\boldsymbol{\beta _{(q_0+1:2q_0)}(x})$ are associated with parameters $p$ and $\mu$ respectively. The interpretations or the impact of the different features on parameters $\mu$ and $p$ are derived by analyzing the covariate contributions: $\beta _j(\boldsymbol{x_i})x_{ij}$ for $p_i$ and $\beta _{q_0+j}(\boldsymbol{x_i})x_{ij}$ for $\mu _i$ , where $j=1,\ldots,q_0$ . Sample code for creating an interpretable ZIPNN model (model 8 in Table 5) is given in Listing 7. The structure of the interpretable ZIPNN model is given in Listing 8.

Although it is possible to create interpretable versions of both CANN and ZIPCANN models similar to those for the NN and ZIPNN models, the introduction of the attention layer nullifies one of the primary purposes of skip connections, which is to distinguish the main effects from the complex effects. This is because the attention layer will integrate the main effect represented by the skip connection and with the complex effect from the network. Hence, in this work, we do not discuss interpretable versions of the CANN and ZIPCANN models.

6.1 Predictive performance

The NLL was used to compare the predictive performance of the implemented models. While other comparison measures are available, such as the mean absolute error and the mean square error, these focus on the deviation of the predicted value from the actual observed value and do not fully account for the considered probability distribution when using a probabilistic model. Also, other alternative measures, such as the predicted probability of a zero count, fall short of providing a comprehensive comparison as they only consider the probability of admission as a dichotomous event, and fail to consider the admission rate. The NLL of a model has the form:

(26) \begin{equation} \mathcal{L}_{\mathcal{A}} (\boldsymbol{\beta }) = -\log \Big (L_{f}\Big ) \end{equation}

where $\mathcal{A}$ represents the data set and $L_{f}$ gives the likelihood of the fitted model (e.g. Equation (21) for the log-likelihood of a ZIP mixture distribution). Performance was compared based on the NLL $\mathcal{L}_{\mathcal{T}}$ of different models for the testing set $\mathcal{T}$ . Values of the NLL for both the testing and training data sets, for all considered models, are shown in Table 5. Table 5 also contains the average fitted mean $\hat{\mu }$ for each model for the full data set. The difference between the average fitted mean of a model and that of the empirical mean of observed data indicates population-level bias.

The results indicate that the ZIPNN model has the lowest testing loss compared to other models. In general, the models with an underlying zero-inflated mixture distribution assumption performed better than those with a basic Poisson distribution assumption. The ZIPNN and ZIPCANN models performed better than the traditional ZIP regression model. Similarly, the NN and CANN models performed better than the conventional Poisson regression model. As the difference between the testing loss of ZIPNN and ZIPCANN models is small, it is not possible to conclusively state the supremacy of one over the other due to the inherent randomness in NN results. The different attributes leading to the randomness in NN model results have been extensively discussed by Richman & Wüthrich (Reference Richman and Wüthrich2020). Approaches such as nagging predictor and k-fold validation could be adopted to address this randomness, as done in previous work by Jose et al. (Reference Jose, Macdonald, Tzougas and Streftaris2022). More than the inherent randomness, the underlying architecture sways the model performance on a larger scale. Hence we analyzed this impact by considering different architectures of varying complexity. As shown in Table A.1, the combination of 20,15,10 neurons in the first three hidden layers is a suitable choice for both ZIPNN and ZIPCANN models.

6.2 Model interpretability

In order to compare the interpretations derived from the model results, we consider the coefficient estimates from the ZIP regression model and the covariate contributions from the ZIPNN model. The coefficient estimates $(\boldsymbol{\beta _{reg},\gamma })$ from the ZIP regression model are given in Table 6. For the ease of interpretation of the results, a sum-to-zero constraint was applied to the coefficient estimates of different levels of all the non-binary categorical variables.

Table 6. Coefficient estimates based on the ZIP regression model with the significance codes ("***,""**,""*," ".,"" ") indicating the level of significance of the estimates at levels (0, 0.001,0.01, 0.05, 0.1, 1 )

The interpretation of coefficient estimates from a ZIP distributional assumption is more complex compared to a conventional regression model, because it refers to a mixture distribution. For the ZIP regression model considered here, the count component considers all features, whereas the zero component only considers the age variable. Additionally, in the case of the ZIPNN model, an attention layer is used for both $\lambda$ and $p$ parameters, which means that the covariate contributions of different variables for $\lambda$ need to be considered in conjunction with those for the $p$ parameter. This makes interpreting the results and comparison with ZIP regression model less straightforward. The comparison was carried out between the most significant coefficient estimates in the ZIP regression model and the corresponding covariate contribution using the ZIPNN model.

The plots showing the covariate contributions were produced using the test data. Figs. 7 and 7 (continued) show the covariate contributions for different variables associated with the rate ( $\lambda$ ) parameter, while Figs. 8 and 8 (continued) show the covariate contributions for the $p$ parameter. The horizontal yellow line in the plots was added as a reference line at levels -0.25 and 0.25. This can help to compare the magnitude of the covariate contributions and how they are associated with the parameters.

Figure 7 Graphical representation of covariate contributions for parameter $\lambda$ in the ZIPNN model: (a) AGE; (b) SEX; (c) UR; (d) REGION; (e) EECLASS; (f) EESTATU.

Figure 7 (continued). Graphical representation of covariate contributions for parameter $\lambda$ in the ZIPNN model: (a) EMPREL; (b) PLANTYP.

The analysis of the covariate contributions using the ZIPNN model for the AGE variable for both the probability and rate parameters (Figs. 7 and 8) indicates an increasing trend over age. This suggests a positive relationship showing that as age increases, both the probability of having an admission and the rate of admission increase. The spline fit (blue line) indicates an almost linear trend, with the greatest variability at younger and older ages. In both cases, this suggests that the impact of the age variable on the parameters varies at different ages. The variation appears to be lowest around the age of 47 before increasing for older ages as well. The highly significant positive coefficient estimate for the AGE variable in the count component of the ZIP regression model (see Table 6) indicates that the rate of admission increases with age. The negative estimate in the zero component associated with the $\pi$ parameter also suggests that the probability of having excess zero decreases with age. In other words, the probability of having non-zero admissions increases with age.

The covariate contributions under the ZIPNN model for the SEX variable indicate a higher probability and rate of admission for females (SEX = 2) compared to males (Figs. 7 and 8). The significant positive coefficient estimate in the count component of the ZIP regression model also implies the same. This is in line with the pattern displayed by the crude rates (see Fig. 5). The low variability in both cases suggests that gender has a consistent and explicit impact on admissions related to respiratory diseases.

Regarding the UR variable, the ZIP regression model suggests a lower admission rate for individuals in the urban area (UR = 2), based on its significant negative coefficient estimate (Table 6). This contradicts the pattern observed under the ZIPNN model. Although not considerably different, the covariate contributions suggest that people in urban areas are more prone to respiratory disease-related admissions (Figs. 7 and 8). This difference could be arising due to the fact that the network architecture captures complex interactions between different features, which is accounted for while estimating the non-linear regression weights, $\boldsymbol{\beta (x)}$ . For instance, there may exist an interaction between the age and UR variables, due to the possibility that retired or older individuals are more likely to reside in rural areas.

Figure 8 Graphical representation of covariate contributions for parameter $p$ in the ZIPNN model: (a) AGE; (b) SEX; (c) UR; (d) REGION; (e) EECLASS; (f) EESTATU.

Figure 8 (continued). Graphical representation of covariate contributions for parameter $p$ in the ZIPNN model: (a) EMPREL; (b) PLANTYP.

The coefficient estimates based on the ZIP regression model for the REGION variable $(\beta _{region_z}, z=1,\ldots,5)$ , given in Table 6, indicate a higher rate of admission when the REGION variable is unknown (REGION = 5). However, the estimate is not significant. The high variability of covariate contribution under the ZIPNN model for level 5 of the region variable is also pointing toward this uncertainty (see Fig. 7). On the other hand, the “west” region (REGION = 4) has a significant negative coefficient estimate under the ZIP regression model, implying a lower admission rate for individuals from the western region. The covariate contribution in the ZIPNN model corresponding to the $\lambda$ parameter for the REGION variable also suggests the same (Fig. 7).

In the case of the EECLASS variable, a significantly higher rate of admission is inferred from the coefficient estimates of the ZIP regression model for “hourly non-union” and “hourly union” (levels 4 and 5) and significantly low rates for individuals with unknown (level 9) employee classification (Table 6). A comparable trend is observed when the covariate contributions for both rate and probability parameters under the ZIPNN model are considered together, which is essential as interpreting the results based on a single parameter could be unclear. For example, the covariate contributions for the “salary union” level (level 2) of the EECLASS variable related to the rate parameter (see Fig. 7) suggest that the admission rate is high for individuals belonging to that employee classification cohort. However, the covariate contribution associated with the probability parameter (Fig. 8) indicates that the probability of having an admission is very low for individuals in that group. In contrast, for individuals with unknown employee classification, the covariate contributions not only indicate a lower probability of admission but also imply a lower admission rate.

Similarly, for the EESTATU variable, the combined analysis of covariate contributions of both parameters (Figs. 7 and 8) points toward a higher propensity for admission related to respiratory diseases for individuals with employment status categorized as “long-term disability” (level 7), whereas it is low for individuals with employment status “Active Part Time or Seasonal” (level 2). The same pattern is inferred from the associated coefficient estimates in the ZIP regression model. It is reasonable to conclude that the higher admission rates for people with long-term disability are as anticipated. Drawing further inferences regarding the relationship between employment details of individuals and hospital admission for respiratory diseases is difficult in the current context of admission data. The data provider has compiled the data set by assigning the employment-related information of the primary beneficiary/employee to their dependents as well. In other words, the employment information of those individuals classified as dependents (EMPREL as “spouse” (level = 2) or “child/other” (level = 3)) in the data set does not accurately represent their actual employment status. For instance, an unemployed dependent of an employee/primary beneficiary will also have the same employment-related information as the employee. Due to these inaccuracies within the data set, it becomes challenging to draw conclusive associations between employment information and the rate of admission due to respiratory diseases.

For the EMPREL variable, the pattern is slightly more evident. The coefficient estimates based on the ZIP regression model indicate a significant difference in the admission rates for different levels of EMPREL covariate (Table 6). The individuals with EMPREL as “child/other” (level = 3) have the highest rate of admission, followed by “spouse”(level = 2) and “employee”(level = 1). Although the same pattern is observed in the covariate contribution related to the $\lambda$ parameter (Fig. 7 (continued)), the covariate contribution associated with $p$ (Fig. 8 (continued)) indicates a lower probability of admission for individuals with EMPREL categorized as "spouse" compared to "employee."

The coefficient estimates for the different levels of the PLANTYP variable in the ZIP regression model indicate that the admission rates are significantly lower for the High-Deductible Health Plan (HDHP) (level = 9) than for other types of plans. This trend is also supported by the ZIPNN model. The covariate contributions not only indicate a lower probability of admission (Fig. 8 (continued)) but also suggest lower admission rates (Fig. 7 (continued)) for individuals under the "HDHP" plan type. This suggests that individuals who are covered under HDHP plans may have a lower likelihood of getting hospitalized for respiratory diseases, and if they do, it occurs at a lower rate. This might be due to the fact that the HDHP plan incentivizes individuals to be more proactive in managing their health and seek care only when essential or have more restrictive coverage.