## 3 Statistical Inference: The EM-Type Algorithm
Firstly, taking the necessary derivatives of the Q-function with respect to \(\mathbf{\beta}_{1}\), we obtain the following results \[h_{1}(\mathbf{\beta}_{1})=\frac{\partial Q\Big{(}\mathbf{\theta};\mathbf{\theta}^{(r)}\Big{)} }{\partial\beta_{1,j}}=\sum_{i=1}^{n}\Bigl{(}-\mu_{i}^{(r)}w_{1_{i}}+k_{i} \Bigr{)}x_{1,ij},\] (16) and \[H_{1}(\mathbf{\beta}_{1})=\frac{\partial^{2}Q\Big{(}\mathbf{\theta};\mathbf{\theta}^{(r)} \Big{)}}{\partial\beta_{1,j}\partial\beta_{1,j}^{\top}}=\sum_{i=1}^{n}\Bigl{(} -\mu_{i}^{(r)}w_{1_{i}}\Bigr{)}x_{1,ij}x_{1,ij}^{\top}=\mathbf{X}_{1}^{\top}\mathbf{W} _{1}\mathbf{X}_{1},\] (17) for \(i=1,\ldots,n\) and \(j=1,\ldots,p_{1}\), and where \(\mathbf{W}_{1}=\text{diag}\biggl{\{}-\frac{k_{i}}{\mu_{i}^{(r)}}w_{1_{i}}\biggr{\}}\). Then, the iterative procedure for the Newton-Raphson algorithm for \(\mathbf{\beta}_{1}\) goes as follows \[\mathbf{\beta}_{1}^{(r+1)}\equiv\mathbf{\beta}_{1}^{(r)}-\Bigl{[}H_{1}\Big{(}\mathbf{ \beta}_{1}^{(r)}\Big{)}\Bigr{]}^{-1}h_{1}\Bigl{(}\mathbf{\beta}_{1}^{(r)}\Bigr{)}.\] (18)
* Secondly, differentiating the Q-function with respect to \(\mathbf{\beta}_{2}\) gives \[h_{2}(\mathbf{\beta}_{2})=\frac{\partial Q(\mathbf{\theta};\mathbf{\theta}^{(r)})}{ \partial\beta_{2,l}}=\frac{\partial\mathbb{E}_{z_{i}}[\log g(z_{i};\mathbf{\phi}_{ i}^{(r)})]}{\partial\beta_{2,l}}\] (19) and \[H_{2}(\mathbf{\beta}_{2})=\frac{\partial^{2}Q(\mathbf{\theta};\mathbf{\theta}^{(r)})}{ \partial\beta_{2,l}\partial\beta_{2,l}^{\top}}=\frac{\partial\mathbb{E}_{z_{i }}[\log g(z_{i};\mathbf{\phi}_{i}^{(r)}]}{\partial\beta_{2,l}\partial\beta_{2,l}^ {\top}},\] (20) where for computing \(h_{1}(\mathbf{\beta}_{2})\) and \(H_{2}(\mathbf{\beta}_{2})\), we need to use the pseudo-values \(w_{k_{i}}\) for \(i=1,\ldots,n\) and \(k=1,\ldots,\
u\), because in this case, the maximization of the Q-function reduces to the maximization of the conditional expectation of the log-likelihood of the mixing distribution \(g(z_{i};\phi_{i})\). Then, the Newton-Raphson iterative algorithm for \(\mathbf{\beta}_{2}\) is as follows \[l_{c}(\boldsymbol{\theta})=\sum_{i=1}^{n}\Big{[}-\mu_{i}z_{i}+k_{i}\log(\mu_{i })+k_{i}\log(z_{i})-\log(k_{i}!)\Big{]}+\sum_{i=1}^{n}\text{log}(g(z_{i}; \phi_{i})) \tag{12}\]0 (21) for \(i=1,\ldots,n\) and \(l=1,\ldots,p_{2}\). * Finally, iterate between the E- and the M-Steps until some convergence criterion is satisfied, for instance \[l_{c}(\boldsymbol{\theta})=\sum_{i=1}^{n}\Big{[}-\mu_{i}z_{i}+k_{i}\log(\mu_{i })+k_{i}\log(z_{i})-\log(k_{i}!)\Big{]}+\sum_{i=1}^{n}\text{log}(g(z_{i}; \phi_{i})) \tag{12}\]1 (22) where \(l^{(r)}\) is the value of the log-likelihood after the \(r\)-th iteration, and where \(tol\) is a small number usually of the form \(10^{-m}\), where \(m\in\mathbb{Z}^{+}\). The stopping criterion refers to the progress of the likelihood function (i.e., its convergence). If the stopping criterion is satisfied, the EM algorithm stops iterating, and the estimate of \(\mathbf{\theta}\) is \(\mathbf{\theta}^{(r+1)}\). Otherwise, \(\mathbf{\theta}\) is updated by \(\mathbf{\theta}^{(r+1)}\), and the algorithm returns to the E-step. EM Estimation for the PLN Regression Model with Varying Dispersion_  
In this section, we implement the EM algorithm for finding the ML estimates of the parameters of the PLN regression model with varying dispersion (Algorithm 1). The complete data log-likelihood of the model is given by  
\[l_{c}(\boldsymbol{\theta})=\sum_{i=1}^{n}\Big{[}-\mu_{i}z_{i}+k_{i}\log(\mu_{i })+k_{i}\log(z_{i})-\log(k_{i}!)\Big{]}+\sum_{i=1}^{n}\text{log}(g(z_{i}; \phi_{i})) \tag{12}\]2 \[l_{c}(\boldsymbol{\theta})=\sum_{i=1}^{n}\Big{[}-\mu_{i}z_{i}+k_{i}\log(\mu_{i })+k_{i}\log(z_{i})-\log(k_{i}!)\Big{]}+\sum_{i=1}^{n}\text{log}(g(z_{i}; \phi_{i})) \tag{12}\]3  
for \(i=1,\ldots,n\).