Subject-specific Deep Neural Networks for Count Data with
  High-cardinality Categorical Features
There is a growing interest in subject-specific predictions using deep neural
networks (DNNs) because real-world data often exhibit correlations, which has
been typically overlooked in traditional DNN frameworks. In this paper, we
propose a novel hierarchical likelihood learning framework for introducing
gamma random effects into the Poisson DNN, so as to improve the prediction
performance by capturing both nonlinear effects of input variables and
subject-specific cluster effects. The proposed method simultaneously yields
maximum likelihood estimators for fixed parameters and best unbiased predictors
for random effects by optimizing a single objective function. This approach
enables a fast end-to-end algorithm for handling clustered count data, which
often involve high-cardinality categorical features. Furthermore,
state-of-the-art network architectures can be easily implemented into the
proposed h-likelihood framework. As an example, we introduce multi-head
attention layer and a sparsemax function, which allows feature selection in
high-dimensional settings. To enhance practical performance and learning
efficiency, we present an adjustment procedure for prediction of random
parameters and a method-of-moments estimator for pretraining of variance
component. Various experiential studies and real data analyses confirm the
advantages of our proposed methods.
Hangbin Lee, Il Do Ha, Changha Hwang, Youngjo Lee
# Subject-specific Deep Neural Networks for Count Data with High-cardinality Categorical Features

###### Abstract

There is a growing interest in subject-specific predictions using deep neural networks (DNNs) because real-world data often exhibit correlations, which has been typically overlooked in traditional DNN frameworks. In this paper, we propose a novel hierarchical likelihood learning framework for introducing gamma random effects into the Poisson DNN, so as to improve the prediction performance by capturing both nonlinear effects of input variables and subject-specific cluster effects. The proposed method simultaneously yields maximum likelihood estimators for fixed parameters and best unbiased predictors for random effects by optimizing a single objective function. This approach enables a fast end-to-end algorithm for handling clustered count data, which often involve high-cardinality categorical features. Furthermore, state-of-the-art network architectures can be easily implemented into the proposed h-likelihood framework. As an example, we introduce multi-head attention layer and a sparsemax function, which allows feature selection in high-dimensional settings. To enhance practical performance and learning efficiency, we present an adjustment procedure for prediction of random parameters and a method-of-moments estimator for pretraining of variance component. Various experiential studies and real data analyses confirm the advantages of our proposed methods.

## 1 Introduction

Deep neural networks (DNNs), which have been proposed to capture the nonlinear relationship between input and output variables (LeCun et al., 2015; Goodfellow et al., 2016), provide outstanding marginal predictions for independent outputs. However, in practical applications, it is common to encounter correlated data with high-cardinality categorical features, which can pose challenges for DNNs. While the traditional DNN framework overlooks such correlation, random effect models have emerged in statistics to make subject-specific predictions for correlated data. Lee and Nelder (1996) proposed hierarchical generalized linear models (HGLMs), which allow the incorporation of random effects from an arbitrary conjugate distribution of generalized linear model (GLM) family.

Both DNNs and random effect models have been successful in improving prediction accuracy of linear models but in different ways. Recently, there has been a rising interest in combining these two extensions. Simchoni and Rosset (2021, 2023) proposed the linear mixed model neural network for continuous (Gaussian) outputs with Gaussian random effects, which allow explicit expressions for likelihoods. Lee and Lee (2023) introduced the hierarchical likelihood (h-likelihood) approach, as an extension of classical likelihood for Gaussian outputs, which provides an efficient likelihood-based procedure. For non-Gaussian (discrete) outputs, Tran et al. (2020) proposed a Bayesian approach for DNNs with normal random effects using the variational approximation method (Bishop and Nasrabadi, 2006; Blei
et al., 2017). Mandel et al. (2023) used a quasi-likelihood approach (Breslow and Clayton, 1993) for DNNs, but the quasi-likelihood method has been criticized for its prediction accuracy. Lee and Nelder (2001) proposed the use of Laplace approximation to have approximate maximum likelihood estimators (MLEs). Although Mandel et al. (2023) also applied Laplace approximation for DNNs, their method ignored many terms in the second derivatives due to computational expense, which could lead to inconsistent estimations (Lee et al., 2017). Therefore, a new approach is desired for non-Gaussian DNNs to obtain the exact MLEs for fixed parameters.

Clustered count data are widely encountered in various fields (Roulin and Bersier, 2007; Henderson and Shimakura, 2003; Thall and Vail, 1990; Henry et al., 1998), which often involve high-cardinality categorical features, i.e., categorical variables with a large number of unique levels or categories, such as subject ID or cluster name. However, to the best of our knowledge, there appears to be no available source code for the subject-specific Poisson DNN models. In this paper, we introduce Poisson-gamma DNN for the clustered count data and derive the h-likelihood that simultaneously provides MLEs of fixed parameters and best unbiased predictors (BUPs) of random effects. In contrast to the ordinary HGLM and DNN framework, we found that the local minima can cause poor prediction when the DNN model contains subject-specific random effects. To resolve this issue, we propose an adjustment to the random effect prediction that prevents from violation of the constraints for identifiability. Additionally, we introduce the method-of-moments estimators for pretraining the variance component. It is worth emphasizing that incorporating state-of-the-art network architectures into the proposed h-likelihood framework is straightforward. As an example, we implement a feature selection method with multi-head attention.

In Section 2 and 3, we present the Poisson-gamma DNN and derive its h-likelihood, respectively. In Section 4, we present the algorithm for online learning, which includes an adjustment of random effect predictors, pretraining variance components, and feature selection method using multi-head attention layer. In Section 5, we provide experimental studies to compare the proposed method with various existing methods. The results of the experimental studies clearly show that the proposed method improves predictions from existing methods. In Section 6, real data analyses demonstrate that the proposed method has the best prediction accuracy in various clustered count data. Proofs for theoretical results are derived in Appendix. Source codes are included in Supplementary Materials.

## 2 Model Descriptions

### Poisson DNN

Let \(y_{ij}\) denote a count output and \(\mathbf{x}_{ij}\) denote a \(p\)-dimensional vector of input features, where the subscript \((i,j)\) indicates the \(j\)th outcome of the \(i\)th subject (or cluster) for \(i=1,...,n\) and \(j=1,...,q_{i}\). For count outputs, Poisson DNN (Rodrigo and Tsokos, 2020) gives the marginal predictor,

\[\eta_{ij}^{m}=\log\mu_{ij}^{m}=\text{NN}(\mathbf{x}_{ij};\mathbf{w},\mathbf{ \beta})=\sum_{k=1}^{p_{L}}g_{k}(\mathbf{x}_{ij};\mathbf{w})\beta_{k}+\beta_{0}, \tag{1}\]

where \(\mu_{ij}^{m}=\text{E}(Y_{ij}|\mathbf{x}_{ij})\) is the marginal mean, \(\text{NN}(\mathbf{x}_{ij};\mathbf{w},\mathbf{\beta})\) is the neural network predictor, \(\mathbf{\beta}=(\beta_{0},\beta_{1},...,\beta_{p_{L}})^{T}\) is the vector of weights and bias between the last hidden layer and the output layer, \(g_{k}(\mathbf{x}_{ij};\mathbf{w})\) denotes the \(k\)-th node of the last hidden layer, and \(\mathbf{w}\) is the vector of all the weights and biases before the last hidden layer. Here the inverse of the log function, \(\exp(\cdot)\), becomes the activation function of the output layer. Poisson DNNs allow highly nonlinear relationship between input and output variables, but only provide the marginal predictions for \(\mu_{ij}^{m}\). Thus, Poisson DNN can be viewed as an extension of Poisson GLM with \(\eta_{ij}^{m}=\mathbf{x}_{ij}^{T}\beta\).

### Poisson-gamma DNN

To allow subject-specific prediction into the model (1), we propose the Poisson-gamma DNN,

\[\eta_{ij}^{c}=\log\mu_{ij}^{c}=\text{NN}(\mathbf{x}_{ij};\mathbf{w},\mathbf{\beta })+\mathbf{z}_{ij}^{T}\mathbf{v}, \tag{2}\]

where \(\mu_{ij}^{c}=\text{E}(Y_{ij}|\mathbf{x}_{ij},v_{i})\) is the conditional mean, \(\text{NN}(\mathbf{x}_{ij};\mathbf{w},\mathbf{\beta})\) is the marginal predictor of the Poisson DNN (1), \(\mathbf{v}=(v_{1},...,v_{n})^{T}\) is the vector of random effects from the log-gamma distribution, and \(\mathbf{z}_{ij}\) is a vector from the model matrix for random effects, representing the high-cardinality categorical features. The conditional mean \(\mu_{ij}^{c}\) can be formulated as

\[\mu_{ij}^{c}=\exp\left\{\text{NN}(\mathbf{x}_{ij};\mathbf{w},\mathbf{\beta}) \right\}\cdot u_{i},\]

where \(u_{i}=\exp(v_{i})\) is the gamma random effect.
Note here that, for any \(\epsilon\in\mathbb{R}\), the model (2) can be expressed as

\[\log\mu^{c}_{ij}=\sum_{k=1}^{p_{L}}g_{k}(\mathbf{x}_{ij};\mathbf{w})\beta_{k}+ \beta_{0}+v_{i}=\sum_{k=1}^{p_{L}}g_{k}(\mathbf{x}_{ij};\mathbf{w})\beta_{k}+( \beta_{0}-\epsilon)+(v_{i}+\epsilon),\]

or equivalently, for any \(\delta=\exp(\epsilon)>0\),

\[\mu^{c}_{ij}=\exp\left\{\text{NN}(\mathbf{x}_{ij};\mathbf{w},\mathbf{\beta}) \right\}\cdot u_{i}=\exp\left\{\text{NN}(\mathbf{x}_{ij};\mathbf{w},\mathbf{\beta} )-\log\delta\right\}\cdot(\delta u_{i}),\]

which leads to an identifiability problem. Thus, it is necessary to place constraints on either the fixed parts \(\text{NN}(\mathbf{x}_{ij};\mathbf{w},\mathbf{\beta})\) or the random parts \(u_{i}\). Lee and Lee (2023) developed subject-specific DNN models with Gaussian outputs, imposing the constraint \(\text{E}(v_{i})=0\), which is common for normal random effects. For Poisson-gamma DNNs, we use the constraints \(\text{E}(u_{i})=\text{E}(\exp(v_{i}))=1\) for subject-specific prediction of count outputs. The use of constraints \(\text{E}(u_{i})=1\) has an advantage that the marginal predictions for multiplicative model can be directly obtained, because

\[\mu^{m}_{ij}=\text{E}[\text{E}(Y_{ij}|\mathbf{x}_{ij},u_{i})]=\text{E}\left[ \exp\left\{\text{NN}(\mathbf{x}_{ij};\mathbf{w},\mathbf{\beta})\right\}\cdot u_{ i}\right]=\exp\left\{\text{NN}(\mathbf{x}_{ij};\mathbf{w},\mathbf{\beta})\right\}.\]

Thus, we employ \(v_{i}=\log u_{i}\) to the model (2), where \(u_{i}\sim\text{Gamma}(\lambda^{-1},\lambda^{-1})\) with \(\text{E}(u_{i})=1\) and \(\text{var}(u_{i})=\lambda\). By allowing two separate output nodes, the Poisson-gamma DNN provides both marginal and subject-specific predictions,

\[\widehat{\mu}^{m}_{ij}=\exp(\widehat{\eta}^{m}_{ij})=\exp\left\{\text{NN}( \mathbf{x}_{ij};\widehat{\mathbf{w}},\widehat{\mathbf{\beta}})\right\}\quad\text {and}\quad\widehat{\mu}^{c}_{ij}=\exp\left\{\text{NN}(\mathbf{x}_{ij};\widehat {\mathbf{w}},\widehat{\mathbf{\beta}})+\mathbf{z}^{T}_{ij}\widehat{\mathbf{v}} \right\},\]

where the hats denote the predicted values. Subject-specific prediction can be achieved by multiplying the marginal mean predictor \(\widehat{\mu}^{m}_{ij}\) and the subject-specific predictor of random effect \(\widehat{u}_{i}=\exp(\widehat{v}_{i})\). Note here that

\[\text{var}(Y|\mathbf{x})=\text{E}(\text{var}(Y|\mathbf{x},\mathbf{v}))+\text {var}(\text{E}(Y|\mathbf{x},\mathbf{v}))\geq\text{E}(\text{var}(Y|\mathbf{x},\mathbf{v})),\]

where \(\text{var}(\text{E}(Y|\mathbf{x},\mathbf{v}))\) represents the between-subject variance and \(\text{E}(\text{var}(Y|\mathbf{x},\mathbf{v}))\) represents the within-subject variance. To enhance the predictions, Poisson DNN improves the marginal predictor \(\mathbf{\mu}^{m}=\text{E}(Y|\mathbf{x})=\)E\(\{\text{E}(Y|\mathbf{x},\mathbf{v})\}\) by allowing highly nonlinear function of \(\mathbf{x}\), whereas Poisson-gamma DNN further uses the conditional predictor \(\mathbf{\mu}^{c}=\text{E}(Y|\mathbf{x},\mathbf{v})\), eliminating between-subject variance. Figure 1 illustrates an example of the proposed model architecture including feature selection in Section 4.

## 3 Construction of h-likelihood

For subject-specific prediction via random effects, it is important to define the objective function for obtaining exact MLEs of fixed parameters \(\mathbf{\theta}=(\mathbf{w},\mathbf{\beta},\lambda)\). In the context of linear mixed models, Henderson et al. (1959) proposed to maximize the joint density with respect to fixed and random parameters. However, it cannot yield MLEs of variance components. There have been various attempts to extend joint maximization schemes with different justifications

Figure 1: An example of the proposed model architecture. The input features are denoted by \(\mathbf{x}_{ij}\) and the high-cardinality categorical features are denoted by \(\mathbf{z}_{ij}\) in the proposed model (2).
(Gilmour et al., 1985; Harville and Mee, 1984; Schall, 1991; Breslow and Clayton, 1993; Wolfinger, 1993), but failed to obtain simultaneously the exact MLEs of all fixed parameters and BUPs of random parameters by optimizing a single objective function. It is worth emphasizing that defining joint density requires careful consideration because of the Jacobian term associated with the random parameters. For \(\mathbf{\theta}\) and \(\mathbf{u}\), an extended likelihood (Lee et al., 2017) can be defined as

\[\ell_{e}(\mathbf{\theta},\mathbf{u})=\sum_{i,j}\log f_{\mathbf{\theta}}(y_{ij}|u_{i})+ \sum_{i}\log f_{\mathbf{\theta}}(u_{i}). \tag{3}\]

However, a nonlinear transformation \(v_{i}=v(u_{i})\) of random effects \(u_{i}\) leads to different extended likelihood due to the Jacobian terms:

\[\ell_{e}(\mathbf{\theta},\mathbf{v}) =\sum_{i,j}\log f_{\mathbf{\theta}}(y_{ij}|v_{i})+\sum_{i}\log f_{\bm {\theta}}(v_{i})\] \[=\sum_{i,j}\log f_{\mathbf{\theta}}(y_{ij}|u_{i})+\sum_{i}\log f_{ \mathbf{\theta}}(u_{i})+\sum_{i}\log\left|\frac{du_{i}}{dv_{i}}\right|\neq\ell_{e} (\mathbf{\theta},\mathbf{u}).\]

The two extended likelihoods \(\ell_{e}(\mathbf{\theta},\mathbf{u})\) and \(\ell_{e}(\mathbf{\theta},\mathbf{v})\) lead to different estimates, raising the question on how to obtain the true MLEs. In Poisson-gamma HGLMs, Lee and Nelder (1996) proposed the use of \(\ell_{e}(\mathbf{\theta},\mathbf{v})\) that can give MLEs for \(\mathbf{\beta}\) and BUPs for \(\mathbf{u}\) by the joint maximization. However, it could not yield MLE for the variance component \(\lambda\). In this paper, we derive the new h-likelihood whose joint maximization simultaneously yields MLEs of the whole fixed parameters including the variance component \(\lambda\), BUPs of the random effects \(\mathbf{u}\), and conditional expectations \(\mathbf{\mu}^{c}\).

Suppose that \(\mathbf{v}^{*}=(v_{1}^{*},...,v_{n}^{*})^{T}\) is a transformation of \(\mathbf{v}\) such that

\[v_{i}^{*}=v_{i}\cdot\exp\{-c_{i}(\mathbf{\theta};\mathbf{y}_{i})\},\]

where \(c_{i}(\mathbf{\theta};\mathbf{y}_{i})\) is a function of \(\mathbf{\theta}\) and \(\mathbf{y}_{i}=(y_{i1},...,y_{iq_{i}})^{T}\) for \(i=1,2,...,n\). Then we define the h-likelihood as

\[h(\mathbf{\theta},\mathbf{v})\equiv\log f_{\mathbf{\theta}}(\mathbf{v}^{*}|\mathbf{y} )+\log f_{\mathbf{\theta}}(\mathbf{y})=\ell_{e}(\mathbf{\theta},\mathbf{v})+\sum_{i=1 }^{n}c_{i}(\mathbf{\theta},\mathbf{y}_{i}), \tag{4}\]

if the joint maximization of \(h(\mathbf{\theta},\mathbf{v})\) leads to MLEs of all the fixed parameters and BUPs of the random parameters. A sufficient condition for \(h(\mathbf{\theta},\mathbf{v})\) to yield exact MLEs of all the fixed parameters in \(\mathbf{\theta}\) is that \(f_{\mathbf{\theta}}(\widetilde{\mathbf{v}}^{*}|\mathbf{y})\) is independent of \(\mathbf{\theta}\), where \(\widetilde{\mathbf{v}}^{*}\) is the mode,

\[\widetilde{\mathbf{v}}^{*}=\operatorname*{arg\,max}_{\mathbf{v}^{*}}h(\mathbf{ \theta},\mathbf{v}^{*})=\operatorname*{arg\,max}_{\mathbf{v}^{*}}\log f_{\mathbf{ \theta}}(\mathbf{v}^{*}|\mathbf{y}).\]

For the proposed model, Poisson-gamma DNN, we found that the following function

\[c_{i}(\mathbf{\theta};\mathbf{y}_{i})=c_{i}(\lambda;y_{i+})=(y_{i+}+\lambda^{-1}) +\log\Gamma(y_{i+}+\lambda^{-1})-(y_{i+}+\lambda^{-1})\log(y_{i+}+\lambda^{-1})\]

satisfies the sufficient condition,

\[\log f(\widetilde{\mathbf{v}}^{*}|\mathbf{y})=\sum_{i=1}^{n}\log f_{\mathbf{ \theta}}(\widetilde{v}_{i}^{*}|\mathbf{y})=\sum_{i=1}^{n}\left\{\log f_{\mathbf{ \theta}}(\widetilde{v}_{i}|\mathbf{y})+c_{i}(\mathbf{\theta};\mathbf{y}_{i}) \right\}=0,\]

where \(y_{i+}=\sum_{j=1}^{q_{i}}y_{ij}\) is the sum of outputs in \(\mathbf{y}_{i}\) and \(\widetilde{v}_{i}=\widetilde{v}_{i}^{*}\cdot\exp\{c_{i}(\mathbf{\theta};\mathbf{y}_ {i})\}\). Then, the h-likelihood at mode \(h(\mathbf{\theta},\widetilde{\mathbf{v}})\) becomes the classical (marginal) log-likelihood,

\[\ell(\mathbf{\theta};\mathbf{y})=\log f_{\mathbf{\theta}}(\mathbf{y})=\log\int\exp\left\{ \ell_{e}(\mathbf{\theta},\mathbf{v})\right\}d\mathbf{v}. \tag{5}\]

Thus, joint maximization of the h-likelihood (4) provides exact MLEs for the fixed parameters \(\mathbf{\theta}\), including the variance component \(\lambda\). BUPs of \(\mathbf{u}\) and \(\mathbf{\mu}^{c}\) can be also obtained from our h-likelihood,

\[\widetilde{\mathbf{u}}=\exp(\widetilde{\mathbf{v}})=\mathrm{E}(\mathbf{u}| \mathbf{y})\quad\text{and}\quad\widetilde{\mathbf{\mu}}^{c}=\exp(\widetilde{ \mathbf{v}})\cdot\exp\left\{\text{NN}(\mathbf{X};\mathbf{w},\mathbf{\beta})\right\}= \mathrm{E}(\mathbf{\mu}^{c}|\mathbf{y}).\]

The proof and technical details for the theoretical results are derived in Appendix A.1.

## 4 Learning algorithm for Poisson-gamma DNN models

In this section, we introduce the h-likelihood learning framework for handling the count data with high-cardinality categorical features. We decompose the negative h-likelihood loss for online learning The entire learning algorithm of the proposed method is briefly described in Algorithm 1.
### Loss function for online learning

The proposed Poisson-gamma DNN can be trained by optimizing the negative h-likelihood loss,

\[\text{Loss}=-h(\mathbf{\theta},\mathbf{v})=-\log f_{\mathbf{\theta}}(\mathbf{y}|\mathbf{v })-\log f_{\mathbf{\theta}}(\mathbf{v})-c(\mathbf{\theta};\mathbf{y}),\]

which is a function of the two separate output nodes \(\mu_{ij}^{m}=\text{NN}(\mathbf{x}_{ij};\mathbf{w},\mathbf{\beta})\) and \(v_{i}=\mathbf{z}_{ij}^{T}\mathbf{v}\). To apply online stochastic optimization methods, the proposed loss function is expressed as

\[\text{Loss}=\sum_{i,j}\left[-y_{ij}\left(\log\mu_{ij}^{m}+v_{i}\right)+e^{v_{i} }\mu_{ij}^{m}-\frac{v_{i}-e^{v_{i}}}{q_{i}\lambda}+a_{i}(\lambda;\mathbf{y}_{i })\right], \tag{6}\]

where \(a_{i}(\lambda;\mathbf{y}_{i})=q_{i}^{-1}\left\{\lambda^{-1}\log\lambda+\log \Gamma(\lambda^{-1})-c_{i}(\lambda,y_{i+})\right\}.\)

### Random Effect Adjustment

While DNNs often encounter local minima, Dauphin et al. (2014) claimed that in ordinary DNNs, local minima may not necessarily result in poor predictions. In contrast to HGLM and DNN, we observed that the local minima can lead to poor prediction when the network reflects subject-specific random effects. In Poisson-gamma DNNs, we impose the constraint \(\mathds{E}(u_{i})=1\) for identifiability, because for any \(\delta>0\),

\[\mu_{ij}^{c}=\exp\left\{\text{NN}(\mathbf{x}_{ij};\mathbf{w},\mathbf{\beta})\right\} \cdot u_{i}=\left[\exp\left\{\text{NN}(\mathbf{x}_{ij};\mathbf{w},\mathbf{\beta}) -\log\delta\right\}\right]\cdot\left(\delta u_{i}\right).\]

However, in practice, Poisson-gamma DNNs often end with local minima that violate the constraint. To prevent poor prediction due to local minima, we introduce an adjustment to the predictors of \(u_{i}\),

\[\widehat{u}_{i}\leftarrow\frac{\widehat{u}_{i}}{\frac{1}{n}\sum_{i=1}^{n} \widehat{u}_{i}}\quad\text{and}\quad\widehat{\beta}_{0}\leftarrow\widehat{ \beta}_{0}+\log\left(\frac{1}{n}\sum_{i=1}^{n}\widehat{u}_{i}\right) \tag{7}\]

to satisfy \(\sum_{i=1}^{n}\widehat{u}_{i}/n=1\). The following theorem shows that the proposed adjustment improves the local h-likelihood prediction. The proof is given in Appendix A.2.

**Theorem 1**.: _In Poisson-gamma DNNs, suppose that \(\widehat{\beta}_{0}\) and \(\widehat{u}_{i}\) are estimates of \(\beta_{0}\) and \(u_{i}\) such that \(\sum_{i=1}^{n}\widehat{u}_{i}/n=1+\epsilon\) for some \(\epsilon\in\mathbb{R}\). Let \(\widehat{u}_{i}^{*}\) and \(\widehat{\beta}_{0}^{*}\) be the adjusted estimators in (7). Then,_

\[h(\widehat{\mathbf{\theta}}^{*},\widehat{\mathbf{v}}^{*})\geq h(\widehat{\mathbf{ \theta}},\widehat{\mathbf{v}}),\]

_and the equality holds if and only if \(\epsilon=0\), where \(\widehat{\mathbf{\theta}}\) and \(\widehat{\mathbf{\theta}}^{*}\) are vectors of the same fixed parameter estimates but with different \(\widehat{\beta}_{0}\) and \(\widehat{\beta}_{0}^{*}\) for \(\beta_{0}\), respectively._

Theorem 1 shows that the adjustment (7) improves the random effect prediction. According to our experience, even though limited, this adjustment becomes important, especially when the cluster size is large. Figure 2 is the plot of \(\widehat{u}_{i}\) against the true \(u_{i}\) under \((n,q)=(100,100)\) and \(\lambda=1\). Figure 2 (a) shows that the use of fixed effects for subject-specific effects (PF-NN) produces poor prediction of \(u_{i}\). Figure 2 (b) and (c) show that the use of random effects for subject-specific effects (PG-NN) improves the subject-specific prediction, and the proposed adjustment improves it further.
### Pretraining variance components

We found that the MLE for variance component \(\lambda=\text{var}(u_{i})\) could be sensitive to the choice of initial value, giving a slow convergence. We propose the use of method-of-moments estimator (MME) for pretraining \(\lambda\),

\[\widehat{\lambda}=\left[\frac{1}{n}\sum_{i=1}^{n}(\widehat{u}_{i}-1)^{2}\right] \left[\frac{1}{2}+\sqrt{\frac{1}{4}+\frac{n\sum_{i}^{n}\widehat{\mu}_{i+}^{-1} (\widehat{u}_{i}-1)^{2}}{\left\{\sum_{i}^{n}(\widehat{u}_{i}-1)^{2}\right\}^{2 }}}\right], \tag{8}\]

where \(\widehat{\mu}_{i+}=\sum_{j=1}^{q_{i}}\widehat{\mu}_{ij}^{m}\). Convergence of the MME (8) is shown in Appendix A.3. Figure 3 shows that the proposed pretraining accelerates the convergence in various settings. In Appendix A.4, we demonstrate an additional experiments for verifying the consistency of \(\widehat{\lambda}\) of the proposed method.

### Feature Selection in High Dimensional Settings

Feature selection methods can be easily implemented to the proposed PG-NN. As an example, we implemented feature selection using the multi-head attention layer with sparsemax function (Martins and Astudillo, 2016; Skrlj et al., 2020; Arik and Pfister, 2021),

\[\text{sparsemax}(\mathbf{z})=\operatorname*{arg\,min}_{p\in\Delta^{K-1}}|| \mathbf{p}-\mathbf{z}||^{2},\]

where \(\Delta^{K-1}=\{\mathbf{p}\in\mathbb{R}^{K}:\mathbf{1}^{T}\mathbf{p}=1,\mathbf{ p}\geq 0\}\). As a high dimensional setting, we generate input features \(x_{kij}\) from \(N(0,1)\) for \(k=1,...,100\), including 10 genuine features \((k\leq 10)\) and 90 irrelevant features \((k>10)\). The output

Figure 3: Learning curve for the variance component \(\lambda\) when (a) \(\lambda=0\), (b) \(\lambda=0.5\), and (c) \(\lambda=1\).

Figure 2: Predicted values of \(u_{i}\) from two replications (marked as o and x for each) when \(u_{i}\) is generated from the Gamma distribution with \(\lambda=1\), \(n=100\), \(q=100\).
is generated from \(\text{Poi}(\mu_{ij}^{c})\) with the mean model,

\[\mu_{ij}^{c}=u_{i}\cdot\exp\left[0.2\left\{1+\cos x_{1ij}+\cdots+\cos x_{6ij}+(x _{7ij}^{2}+1)^{-1}+\cdots+(x_{10ij}^{2}+1)^{-1}\right\}\right],\]

where \(u_{i}=e^{v_{i}}\) is generated from \(\text{Gamma}(2,2)\). The number of subjects is \(n=10^{4}\), i.e. the cardinality of the categorical feature is \(10^{4}\). The number of repeated measures (cluster size) is set to be \(q=20\), which is smaller than number of features \(p=100\). We employed a multi-layer perceptron with 20-10-10 number of nodes and three-head attention layer for feature selection. Other details are derived in Section 5. Figure 4 shows the average attention scores of 50 repetitions. It is evident that all the genuine features are ranked at the top 10.

## 5 Experimental Studies

To investigate the performance of the Poisson-gamma DNN, we conducted experimental studies. The five input variables \(\mathbf{x}_{ij}=(x_{1ij},...,x_{5ij})^{T}\) are generated from the AR(1) process with autocorrelation \(\rho=0.5\) for each \(i=1,...,n\) and \(j=1,...,q\). The random effects are generated from either \(u_{i}\sim\text{Gamma}(\lambda^{-1},\lambda^{-1})\) or \(v_{i}\sim\text{N}(0,\lambda)\) where \(v_{i}=\log u_{i}\). When \(\lambda=0\), the conditional mean \(\mu_{ij}^{c}\) is identical to the marginal mean \(\mu_{ij}^{m}\). The output variable \(y_{ij}\) is generated from \(\text{Poisson}(\mu_{ij}^{c})\) with

\[\mu_{ij}^{c}=u_{i}\cdot\exp\left[0.2\left\{1+\cos x_{1ij}+\cos x_{2ij}+\cos x_ {3ij}+(x_{4ij}^{2}+1)^{-1}+(x_{5ij}^{2}+1)^{-1}\right\}\right],\]

Results are based on the 100 sets of simulated data. The data consist of \(q=10\) observations for \(n=1000\) subjects. For each subject, 6 observations are assigned to the training set, 2 are assigned to the validation set, and the remaining 2 are assigned to the test set.

For comparison, we consider the following models.

* **P-GLM** Classic Poisson GLM for count outputs using R.
* **N-NN** Conventional DNN for continuous outputs.
* **P-NN** Poisson DNN for count outputs.
* **PN-GLM** Poisson-normal HGLM using lme4 (Bates et al., 2015) package in R.
* **PG-GLM** Poisson-gamma HGLM using the proposed method.
* **NF-NN** Conventional DNN with fixed subject-specific effects for continuous outputs.
* **NN-NN** DNN with normal random effects for continuous outputs (Lee and Lee, 2023).
* **PF-NN** Conventional Poisson DNN with fixed subject-specific effects for count outputs.
* **PG-NN** The proposed Poisson-gamma DNN for count outputs.

To evaluate the prediction performances, we consider the root mean squared Pearson error (RMSPE)

\[\text{RMSPE}=\sqrt{\frac{1}{N}\sum_{i,j}\frac{(y_{ij}-\widehat{\mu}_{ij})^{2}} {V(\widehat{\mu}_{ij})}},\]

Figure 4: Average attention scores. First 10 features are genuine and the others are irrelevant.
where \(\text{Var}(y_{ij}|u_{i})=\phi V(\mu_{ij})\) and \(\phi\) is a dispersion parameter of GLM family. For Gaussian outputs, the RMSPE is identical to the ordinary root mean squared error, since \(\phi=\sigma^{2}\) and \(V(\mu_{ij})=1\). For Poisson outputs, since \(\phi=1\) and \(V(\mu_{ij})=\mu_{ij}\), the RMSPE for test set is given by

\[\text{RMSPE}=\sqrt{\frac{1}{N_{\text{test }}}\sum_{(i,j)\in\text{test}}\frac{(y_{ij}- \widehat{\mu}_{ij})^{2}}{\widehat{\mu}_{ij}}}.\]

P-GLM, N-NN, and P-NN give marginal predictions \(\widehat{\mu}_{ij}=\widehat{\mu}_{ij}^{m}\), while the others give subject-specific predictions \(\widehat{\mu}_{ij}=\widehat{\mu}_{ij}^{c}\). N-NN, NF-NN, and NN-NN are models for continuous outputs, while the others are models for count outputs. For NF-NN and PF-NN, predictions are made by maximizing the conditional likelihood \(\sum_{i,j}\log f_{\mathbf{\theta}}(y_{ij}|v_{i})\). On the other hand, for PN-GLM, PG-GLM, NN-NN, and PG-NN, subject-specific predictions are made by maximizing the h-likelihood.

PN-GLM is the generalized linear mixed model with random effects \(v_{i}\sim N(0,\lambda)\). Current statistical software for PN-GLM and PG-GLM (1me4 and dhglm) provide approximate MLEs using Laplace approximation. The proposed learning algorithm can yield exact MLEs for PG-GLM, using solely the input and output layers while excluding the hidden layers. Among various methods for NN-NN (Tran et al., 2020; Mandel et al., 2023; Simchoni and Rosset, 2021, 2023; Lee and Lee, 2023), we applied the state-of-the-art method proposed by Lee and Lee (2023). All the DNNs and PG-GLMs were implemented in Python using Keras (Chollet et al., 2015) and TensorFlow (Abadi et al., 2015). For all DNNs, we employed a standard multi-layer perceptron (MLP) consisting of 3 hidden layers with 10 neurons and leaky ReLU activation function. We applied the Adam optimizer with a learning rate of 0.001 and an early stopping process based on the validation loss while training the DNNs. NVIDIA Quadro RTX 6000 were used for computations.

Table 1 shows the mean and standard deviation of test RMSPEs from the experimental studies. When the true model does not have random effects (G(0) and N(0)), the PG-NN is comparable to the P-NN without random effects, which should perform the best (marked by the bold face) in terms of RMSPE. N-NN (P-NN) without random effects is also better than NF-NN and NN-NN (PF-NN and PG-NN) with random effects. When the distribution of random effects is correctly specified (G(0.5) and G(1)), the PG-NN performs the best in terms of RMSPE. Even when the distribution of random effects is misspecified (N(0.5), N(1)), the PG-NN still performs the best. This result is in accordance with the simulation results of McCulloch and Neuhaus (2011), namely, in GLMMs, the prediction accuracy is little affected for violations of the distributional assumption for random effects: see similar performances of PN-GLM and PG-GLM.

It has been known that handling the high-cardinality categorical features as random effects has advantages over handling them as fixed effects (Lee et al., 2017), especially when the cardinality of categorical feature is close to the sample size, i.e., the number of observations in each category (cluster size \(q\)) is relatively small. Thus, to emphasize the advantages of PG-NN over PF-NN in high-cardinality categorical features, we consider two additional scenarios for experimental study with cluster size \(q_{\text{train}}=3\) and \(q_{\text{train}}=1\) where \(\lambda=0.2\) and \(n=1000\). Mean and standard deviation of RMSPE of PF-NN are 1.269 (0.038) and 1.629 (0.086) for \(q_{\text{train}}=3\) and \(q_{\text{train}}=1\), respectively. Those of PG-NN are 1.124 (0.028) and 1.284 (0.049) for each scenario. Therefore, the proposed method enhances subject-specific predictions as the cardinality of categorical features becomes high.

\begin{table}
\begin{tabular}{c c c c c c} \hline \hline  & \multicolumn{4}{c}{Distribution of random effects (\(\lambda\))} \\ \cline{2-5} Model & G(0) \& N(0) & G(0.5) & G(1) & N(0.5) & N(1) \\ \hline P-GLM & 1.046 (0.029) & 1.501 (0.055) & 1.845 (0.085) & 1.745 (0.113) & 2.818 (0.467) \\ N-NN & 1.013 (0.018) & 1.473 (0.042) & 1.816 (0.074) & 1.713 (0.097) & 1.143 (0.432) \\ P-NN & **1.011 (0.018)** & 1.470 (0.042) & 1.812 (0.066) & 1.711 (0.099) & 1.161 (0.440) \\ PN-GLM & 1.048 (0.029) & 1.112 (0.033) & 1.115 (0.035) & 1.124 (0.030) & 1.152 (0.034) \\ PG-GLM & 1.048 (0.020) & 1.123 (0.027) & 1.106 (0.023) & 1.139 (0.026) & 1.161 (0.028) \\ NF-NN & 1.152 (0.029) & 1.301 (0.584) & 1.136 (0.311) & 1.241 (1.272) & 1.402 (0.298) \\ NN-NN & 1.020 (0.020) & 1.121 (0.026) & 1.209 (0.067) & 1.256 (0.097) & 2.773 (0.384) \\ PF-NN & 1.147 (0.025) & 1.135 (0.029) & 1.128 (0.027) & 1.129 (0.024) & 1.128 (0.027) \\ PG-NN & 1.016 (0.019) & **1.079 (0.024)** & **1.084 (0.023)** & **1.061 (0.022)** & **1.085 (0.026)** \\ \hline \hline \end{tabular}
\end{table}
Table 1: Mean and standard deviation of test RMSPEs of simulation studies over 100 replications. G(0) implies the absence of random effects, i.e., \(v_{i}=0\) for all \(i\). Bold numbers indicate the minimum.
## 6 Real Data Analysis

To investigate the prediction performance of clustered count outputs in practice, we examined the following five real datasets:

* **Epilepsy data**(Thall and Vail, 1990)
* **CD4 data**(Henry et al., 1998)
* **Bolus data**(Henderson and Shimakura, 2003)
* **Owls data**(Roulin and Bersier, 2007)
* **Fruits data**(Santa et al., 2010)

For all the DNNs, a standard MLP with one hidden layer of 10 neurons and a sigmoid activation function were employed. For longitudinal data (Epilepsy, CD4, Bolus), the last observation for each patient was used as the test set. For clustered data (Owls, Fruits), an observation was randomly selected as the test set from each cluster. RMSPEs are reported in Table 2, which shows that the use of subject-specific models (PG-NN and PG-GLM) for count data are the best. Throughout the datasets, P-GLM performs better than P-NN, implying that non-linear model does not improve the linear model in the absence of subject-specific random effects. Meanwhile, in the presence of subject-specific random effects, PG-NN is always preferred to PG-GLM except for Fruits data. The results imply that introducing subject-specific random effects in DNNs can help to identify the nonlinear effects of the input variables. Therefore, while DNNs are widely recognized for improving predictions in independent datasets, introducing subject-specific random effects could be necessary for DNNs to improve their predictions in correlated datasets with high-cardinality categorical features.

## 7 Concluding Remarks

When the data contains high-cardinality categorical features, introducing random effects into DNNs is advantageous. We develop subject-specific Poisson-gamma DNN for clustered count data. The h-likelihood enables a fast end-to-end learning algorithm using the single objective function. By introducing subject-specific random effects, DNNs can effectively identify the nonlinear effects of the input variables. Various state-of-the-art network architectures can be easily implemented into the h-likelihood framework, as we demonstrate with the feature selection based on multi-head attention.
