In Generalized Linear Models, likelihood equations are intractable and do not have explicit solutions; thus, they must be solved by using Newton-type algorithms. Solving these equations by iterations, however, can be problematic: the iterations might converge to wrong values or the iterations might not converge at all. In this study, we derive the modified maximum likelihood estimators for Poisson regression model and study their properties. We also search the robustness of these estimators when there are outliers in the covariates.

Keywords: count data, poisson regression, modified maximum likelihood, newton-type algorithms, dixon’s outlier model

Poisson regression is widely used for modeling count data, especially when there is no over- or under- dispersion.¹ Since the likelihood equations from this model are intractable, solving these equations requires using iterative methods, such as Newton Raphson or Fisher scoring. However, using iterative methods to find maximum likelihood estimators (MLEs) can generally be problematic and time should be spent to investigate the stability of such solutions.^2-4 Specifically the following difficulties can arise: the iterations might converge to wrong values if the likelihood equations have multiple roots, or the iterations might not converge at all. See⁵ and⁶ for a discussion about situations where one encounters these difficulties in solving MLEs. The common software, Stata for example, is known to be very sensitive to numerical iterations. Researchers have reported problems in getting Poisson regression estimates with the “poisson” command, which encounters problems in locating the maximum and does not converge.⁷ Note that in Poisson regression modeling, additional problems might occur. The most common problem is the over- or under- dispersion in data, in which case using a more flexible model such as negative binomial regression is more appropriate. The second problem is analogous to the complete separation or quasi-complete separation problem in binary regression: the MLEs may not exist for certain data configurations, see.^7,8 In this study we do not consider either of these problems and focus only on the case where it is appropriate to model the data with Poisson regression.

Unlike the maximum likelihood (ML) technique, modified maximum likelihood (MML) methodology produces explicit estimates. MML achieves explicit estimates by linearizing intractable functions within the likelihood equations using the ordered statistics.⁹ Asymptotically, MML estimators (MMLEs) are known to be unbiased and have minimum variances, i.e. they are fully efficient. For small sample sizes MMLEs have negligible bias and their variances are only marginally bigger than the minimum variance bounds, i.e. they are highly efficient.^10-15

In the GLM setting, Tiku and Vaughan³ used the MML methodology to extend the techniques of traditional logistic regression to non-logistic density functions. Oral and Gunay¹⁶ and Oral¹⁷ later extended the work in³ to the binary regression model with one stochastic covariate. Oral¹⁸ derived the MMLEs in general GLMs which use canonical link when there is only one risk factor. In this study, we derive the explicit MMLEs of Poisson regression model, generalize the derivations to more than one covariate, and study their robustness properties via simulations.

The Univariate Poisson regression model is given by

$E(Y_i|X_i=x_i)={\mu }_i$

$=\exp (z_i)$  (1)

where $z_i=\alpha +\beta x_i$ , for $1\le i\le n$ , and the outcome $Y$ has the probability distribution

$f_y(y_i)=\frac{\exp (-{\mu }_i){\mu }^{y_i}}{y_i!},y_i=0,1,2,\mathrm{...}$ (2)

Note that in equation (2), $Y$ is presumed to increase with X so that b is a priori greater than zero. For the model given in (1)-(2), the log-likelihood function of the random sample $(y_i,x_i)$  can be written as

$\ln L\propto {\displaystyle \sum _{i=1}^n{y}_i{z}_i}-{\displaystyle \sum _{i=1}^{n}g(z_i)}$ , (3)

where $g(z_i)=\exp (z_i)$ . The likelihood equations for estimating a and b do not have explicit solutions because of the nonlinear function $g(z_i)$ , $1\le i\le n$ . To obtain the MMLEs, we first express the likelihood equations in terms of the ordered variates $z_{(1)}\le z_{(2)}\le \mathrm{...}\le z_{(n)}$ . The likelihood equations can be re-written as

$\frac{\partial \ln L}{\partial \alpha }={\displaystyle \sum _{i=1}^n\left\{y_{[i]}-g(z_{(i)})\right\}}=0$  (4)

and

$\frac{\partial \ln L}{\partial \beta }={\displaystyle \sum _{i=1}^n{x}_{(i)}\left\{y_{[i]}-g(z_{(i)})\right\}}=0$  (5)

Where $y_{[i]}$  is the concomitant of $x_{(i)}$   $(1\le i\le n)$  and $g(z_{(i)})=\exp (z_{(i)})$ . Linearizing the intractable function $g(z_i)$  by using the first two terms of its Taylor series expansion around the population quantiles $t_{(i)}=E(z_{(i)})$  ( $1\le i\le n$ ), we find

$g(z_{(i)})\cong a_i+b_i{z}_{(i)}$ , (6)

where $a_i=\exp (t_{(i)})\left\{1-t_{(i)}\right\}$  and $b_i=\exp (t_{(i)})$  for $1\le i\le n$ . In order to calculate $t_{(i)}=E(z_{(i)})$  values, we define a dummy random variable U with the probability density function

$f(u)=\exp (u),u<0$  (7)

and re-write the model (1) as

$E(Y_i|X_i=x_i)={\mu }_i=F\left(z_i\right)$ , (8)

where $F(u)=\exp (u)$ , $u<0$ . Thus, the $t_{(i)}$  values can be obtained from the equation

$t_{(i)}=\ln \left(i/\left(n+1\right)\right)$ , (9)

For $1\le i\le n$ ; asymptotically $z_{(i)}-t_{(i)}\cong 0.$ Alternatively, when n is large one can utilize the standard normal distribution

$t_{(i)}={\Phi }^{-1}\left(i/\left(n+1\right)\right)$ . (10)

Incorporating (6) into (4)-(5) and solving the resulting modified likelihood equations yield the explicit MMLEs below:

$\widehat{\alpha }=\frac{\delta }{m}-\widehat{\beta }{\overline{x}}_a$  and $\widehat{\beta }=\frac{{\displaystyle \sum _{i=1}^n{\delta }_i(x_{(i)}-{\overline{x}}_a})}{{\displaystyle \sum _{i=1}^n{b}_i{(x_{(i)}-{\overline{x}}_a)}^2}}$ , (10)

where

$\delta ={\displaystyle \sum _{i=1}^n{\delta }_i}$ , ${\delta }_i=y_{[i]}-a_i$ , $m={\displaystyle \sum _{i=1}^n{b}_i}$ , (11)

and

${\overline{x}}_a=\left({\displaystyle \sum _{i=1}^{n}b{}_{i}x{}_{(i)}}\right)/m$ . (12)

The MMLEs derived above are asymptotically equivalent to their corresponding MLEs, giving them the same attractive asymptotic properties; however, one can refine the estimates by re-calculating $a_i$  and $b_i$  values by replacing the theoretical population quantiles with their estimated values $t_{(i)}=\widehat{\alpha }+\widehat{\beta }{x}_{(i)}$ $(1\le i\le n)$ . This process might be repeated until a desired convergence is met. The stabilization generally is reached within a few iterations.

Asymptotic variances and co-variances

Vaughan and Tiku¹⁵ proved rigorously that the MMLEs are asymptotically unbiased and their variances and co-variances are exactly the same as those of the MLEs. In the present situation, therefore, the asymptotic variances and the covariance of $\widehat{\alpha }$  and $\widehat{\beta }$  are given by $I^{-1}\left(\alpha ,\beta \right)$ , where I is the Fisher Information matrix consisting of the elements $-E\left({\partial }^2\ln L^{*}/\partial {\alpha }^2\right),$ $-E\left({\partial }^2\ln L^{*}/\partial \alpha \partial \beta \right)$ , and $-E\left({\partial }^2\ln L^{*}/\partial {\beta }^2\right)$ . From the modified likelihood equations

$\frac{\partial \ln L}{\partial \alpha }\propto \frac{\partial \ln L^{*}}{\partial \alpha }={\displaystyle \sum _{i=1}^n\left\{y_i-(a_i+b_i{z}_i)\right\}}=0$ ,

$\frac{\partial \ln L}{\partial \beta }\propto \frac{\partial \ln L^{*}}{\partial \beta }={\displaystyle \sum _{i=1}^n{x}_i\left\{y_i-(a_i+b_i{z}_i)\right\}}=0$ ,

the Fisher Information matrix can be easily obtained as

$V=I^{-1}({\gamma }_0,{\gamma }_1)={\left[\begin{array}{cc}{\displaystyle \sum _{i=1}^n{Q}_i}& {\displaystyle \sum _{i=1}^n{Q}_i{x}_i}\\ {\displaystyle \sum _{i=1}^n{Q}_i{x}_i}& {\displaystyle \sum _{i=1}^n{Q}_i{x}_i^2}\end{array}\right]}^{-1}$ , (13)

where $Q_i=\exp (z_i)$ . V is estimated by replacing $Q_i$  with its estimate ${\widehat{Q}}_i=\exp ({\widehat{z}}_i)$ , ${\widehat{z}}_i=\widehat{\alpha }+\widehat{\beta }{x}_i$  ( $1\le i\le n$ ). Since $z_i$  values converge to $t_i$  values as n tends to infinity, $a_i$  and $b_i$  values are treated as constant coefficients for large n, see also.³ Hence, the asymptotic variances can be estimated by

$\text{Var}(\widehat{\alpha })={\displaystyle \sum _{i=1}^n{\widehat{Q}}_i{x}_i^2}/\left\{{\displaystyle \sum _{i=1}^n{\widehat{Q}}_i}{\displaystyle \sum _{i=1}^n{\widehat{Q}}_i{x}_i^2-{\left({\displaystyle \sum _{i=1}^n{\widehat{Q}}_i{x}_i}\right)}^2}\right\}$ , (14)

$\text{Var}(\widehat{\beta })={\displaystyle \sum _{i=1}^n{\widehat{Q}}_i}/\left\{{\displaystyle \sum _{i=1}^n{\widehat{Q}}_i}{\displaystyle \sum _{i=1}^n{\widehat{Q}}_i{x}_i^2-{\left({\displaystyle \sum _{i=1}^n{\widehat{Q}}_i{x}_i}\right)}^2}\right\}$ . (15)

Hypothesis testing

Testing the null hypothesis $H_0:\beta =0$  is of great practical importance in Poisson regression modelling. The likelihood ratio statistic for testing $H_0$  is $LR=-2(L_0-L_1)$ , where $L_0$  and $L_1$  denote the maximized log-modified likelihood functions under the null and alternative hypotheses, respectively. The null distribution of LR is asymptotically a chi-square with 1 degree of freedom. Large values of LR lead the rejection of $H_0$ . Alternatively, the Wald statistic W (the ratio of $\widehat{\beta }$  to its standard error) might be used. Since $\widehat{\beta }$  is asymptotically equivalent to the MLE, the null distribution of W is asymptotically normal N(0,1). Large values of W lead to the rejection of $H_0$ .

Numerical example

To compare ML and MML estimates numerically, we analyzed the data given on page 82 of Agresti.¹⁹ The data is from a study of nesting horseshoe crabs where the response Y is the number of satellites that each female crab has, and the corresponding values of the covariate X is the carapace width of 173 crabs. The study investigates the relationship between Y and X. We calculated the MMLEs from equations (10)-(12) and their approximate standard errors from (14)-(15). The FORTRAN code written to carry out the calculations can be obtained from the author. Our results are completely consistent with those given in,¹⁹ which is expected; see Table 1.

Remark: In solving (10)-(12), Oral¹⁸ proposed to calculate the initial $t_{(i)}$  values from the least squares estimators (LSEs), which is a different approach than using equation (9) (Approach 1) or equation (10) (Approach 2). Since $t_{(i)}=E(\alpha +\beta x_{(i)})$ , $t_{(i)}$  values can be approximated by ${\tilde{t}}_{(i)}=\tilde{\alpha }+\tilde{\beta }{x}_{(i)}$ , where

$\tilde{\alpha }=\overline{y}-\tilde{\beta }\overline{x}$  and $\tilde{\beta }={\displaystyle \sum _{i=1}^n\left(x_i-\overline{x}\right)y_i}/{\displaystyle \sum _{i=1}^n{\left(x_i-\overline{x}\right)}^2}$

are the LSEs; see also³ and.²⁰ When using this approach (say, Approach 3), the $t_{(i)}$  values need to be revised after the first iteration with their estimated values $t_{(i)}=\widehat{\alpha }+\widehat{\beta }{x}_{(i)}$ $(1\le i\le n)$  as described above. Estimating population quantiles from the LSEs changes neither the derivations/solutions (10)-(12) nor the results. However, the total revision number needed for stabilization under different approaches is not the same, see also³ page 889. To compare the performance of these three approaches, we conducted a simulation study where we calculated the bias values and variances of the resulting MLEs as well as the coverage probabilities. We also provided average revision numbers needed for stabilization from each approach. We set $\alpha$  to zero and considered various values for $\beta$  and sample size n. Our results from 10,000 Monte Carlo runs are given in Table 2. 

As can be seen from Table 2, all approaches provide same biases and variances after stabilization, which is expected. The resulting coverage probabilities from all approaches are close to 0.95, and as sample size increases, both bias values and variances decrease, also as expected. For a given $(\alpha ,\beta )$  value, it can be seen that the fastest stabilization is achieved by Approach 2 (i.e. equation (10)). Thus, although all three approaches yield the same results, we suggest to calculate initial $t_{(i)}$  values $(1\le i\le n)$  from equation (10) because the stabilization from this approach is the fastest one.

Generalization to multivariable case

Now consider k $(k\ge 2)$  covariates and assume all of them take positive values without loss of generality. The Poisson regression model with k covariates can be written as

$E(Y_i|x_{i1},x_{i2},...,x_{ik})=F\left(z_i\right)$  (16)

where $F\left(z_i\right)=\exp (z_i)$  and

$z_i={\beta }_0+{\displaystyle \sum _{j=1}^k{\beta }_j{x}_{ij}}$ , (17)

for $1\le i\le n.$ , $1\le j\le k$ . Following the same lines of,³ in order to rank the z-values we can assume that all covariates are equally effective in increasing the response Y, i.e. we initially take ${\beta }_j$ ’s all equal, and order the z-values that would correspond to the ordered x-values, where $x_i=x_{i1}+x_{i2}+\mathrm{...}+x_{ik}$ $(1\le i\le n)$ . In other words, the ordered z-values become

$z_{(i)}={\beta }_0+{\beta }_1{x}_{i1}^{*}+...+{\beta }_k{x}_{ik}^{*}$ , (18)

where the vector ${\left[\begin{array}{cccccc}1& x_{i1}^{*}& x_{i2}^{*}& .& .& x_{ik}^{*}\end{array}\right]}^{}$  is the ith row of the matrix

$X^{*}=\left[\begin{array}{ccccc}1& x_{11}^{*}& x_{12}^{*}& \mathrm{...}& x_{1k}^{*}\\ 1& x_{21}^{*}& x_{22}^{*}& \mathrm{...}& x_{2k}^{*}\\ .& .& .& \mathrm{...}& .\\ .& .& .& \mathrm{...}& .\\ 1& x_{n1}^{*}& x_{n2}^{*}& \mathrm{...}& x_{nk}^{*}\end{array}\right]$ , (19)

which is constructed by arranging the rows of the X matrix

$X=\left[\begin{array}{ccccc}1& x_{11}& x_{12}& \mathrm{...}& x_{1k}\\ 1& x_{21}& x_{22}& \mathrm{...}& x_{2k}\\ .& .& .& \mathrm{...}& .\\ .& .& .& \mathrm{...}& .\\ 1& x_{n1}& x_{n2}& \mathrm{...}& x_{nk}\end{array}\right]$ ,

so as to correspond to the ordered $x_{(i)}$  value $(1\le i\le n)$ . The MMLEs can be obtained along the same lines as in the Univariate case:

$\widehat{\Gamma }={\left(X^{*}{}^{T}M{X}^{*}\right)}^{-1}{X}^{*}{}^T\Delta$  (20)

where $\Delta ={\left[\begin{array}{ccccc}{\delta }_1& {\delta }_2& .& .& {\delta }_n\end{array}\right]}^T$  , ${\delta }_i$  is given by (11) and M is the nxn diagonal matrix

$M=\left[\begin{array}{cccc}{b}_1& 0& \mathrm{...}& 0\\ 0& b_2& \mathrm{...}& 0\\ \mathrm{...}& \mathrm{..}& \mathrm{...}& \mathrm{..}\\ 0& 0& \mathrm{...}& b_n\end{array}\right]$ .

The asymptotic variance-covariance matrix V of the estimators can be derived from the Fisher information matrix $V=I^{-1}\left({\beta }_0,{\beta }_1,\mathrm{...},{\beta }_k\right)$  as given below

$V={\left[\begin{array}{llll}{\displaystyle \sum Q_i}\hfill & {\displaystyle \sum Q_i}{x}_{1i}\hfill & \mathrm{...}\hfill & {\displaystyle \sum Q_i}{x}_{ki}\hfill \\ {\displaystyle \sum Q_i}{x}_{1i}\hfill & {\displaystyle \sum Q_i}{x}^2{}_{1i}\hfill & \mathrm{...}\hfill & {\displaystyle \sum Q_i}{x}_{ki}{x}_{1i}\hfill \\ \mathrm{...}\hfill & \mathrm{...}\hfill & \mathrm{...}\hfill & \mathrm{...}\hfill \\ {\displaystyle \sum Q_i}{x}_{ki}\hfill & {\displaystyle \sum Q_i}{x}_{ki}{x}_{1i}\hfill & \mathrm{...}\hfill & {\displaystyle \sum Q_i}{x}^2{}_{ki}\hfill \end{array}\right]}^{-1}$  (21)

where $Q_i=\exp (z_i)$ . V is estimated by replacing $Q_i$  by

${\widehat{Q}}_i=\exp ({\widehat{z}}_i)$ , ${\widehat{z}}_{(i)}={\widehat{\beta }}_0+{\widehat{\beta }}_1{x}_{i1}^{*}+...+{\widehat{\beta }}_k{x}_{ik}^{*}$   $(1\le i\le n)$ .

Robustness

Measures of influence considered in linear regression models, such as high leverage values, are analogous in the GLM framework. Large leverage values typically mean that there are outliers in covariates. When outliers present in the data, inferences based on MLEs becomes unreliable. In fact, it has been showed that MLEs are not robust in GLMs.²¹ In Poisson regression setting, if there are outliers in the continuous covariates, the estimates can be influenced. Thus, we also studied the robustness properties of the derived MMLEs under several outlier models. We considered the Univariate model given by (1) for simplicity. We assumed that $\alpha =0$  and performed a Monte-Carlo study for different values of $\beta$ , where (n-r) of the observations $X_1,X_2,...,X_n$  (we don’t know which) come from the Standard Normal Distribution with $\sigma =1$  and the remaining r observations come from the Normal distribution with a scale $c\sigma$  where c is a positive constant. We calculated the value of r from the equation $r=\left[0.1n+0.5\right]$  (Dixon’s outlier model). The outlier models considered are:

(a) (n-r) come from N(0,1) and r come from N(0,1) (No outlier situation),

(b) (n-r) come from N(0,1) and r come from N(0,1.5),

(c) (n-r) come from N(0,1) and r come from N(0,2),

(d) (n-r) come from N(0,1) and r come from N(0,4).

Note that the model (a) above does not involve outliers and is given for the sake of comparisons. In order to be able to make direct comparisons, after generating the X values we divided them by the standard deviation of the distribution for each model. After generating the X values, we calculated $z_i=\alpha +\beta x_i$ and ${\mu }_i=\exp (z_i)$  for $1\le i\le n$  to generate $Y_i$  values from Poisson distribution with mean ${\mu }_i$  ( $1\le i\le n$ ). The values obtained from 5000 runs are given in (Table 3).

As can be seen from the table, the biases in the estimates are negligible for all models. The variances $\text{Var}(\widehat{\beta })$  (hence the Wald statistics W) are almost the same for a given n for the models (a), (b), (c) and (d), which means that the MMLEs are robust to outliers in the covariate. Note that the MML methodology achieves robustness through the $t_{(i)}$  ( $1\le i\le n$ ) values.

Poisson regression serves as a useful technique to model count data. The MLEs in Poisson regression are obtained via Newton-type algorithms; however these algorithms might not converge or converge to inaccurate values. In this study we derived the explicit MMLEs for Poisson regression. We also considered the case where there are outliers in the (continuous) covariate, which generally is the case in real life applications, and searched the properties of the derived MMLEs under several data violations. Although the scope of the simulations reported here is limited, we can conclude that MML methodology provides robust estimation in Poisson regression.

None.

1. Cameron AC, Trivedi PK. Regression Analysis of Count Data. Cambridge Univ. Press, NY, USA. 1998.
2. Casella G, Berger R. Statistical Inference. (2nd edn), Thomson Learning, Pacific Grove, CA, USA. 2002.
3. Tiku ML, Vaughan DC. Logistic and nonlogistic density functions in binary regression with nonstochastic covariates. Biometrical Journal. 1997;39(8):883‒898.
4.  Vaughan DC. The generalized secant hyperbolic distribution and its properties. Communications in Statistics ‒ Theory and Methods. 2002;31(2):219‒238.
5. Barnett VD. Evaluation of the maximum likelihood estimator when the likelihood equation has multiple roots. Biometrika. 1966;53(1):151‒165.
6. Lee KR, Kapadia CH, Dwight BB. On estimating the scale parameter of Rayleigh distribution from censored samples. Statistische Hefte. 1980;21(1):14‒20.
7. Silva JMC, Tenreyro S. Poisson: Some convergence Issues. The Stata Journal. 2011;11(2):207‒212.
8. Silva JMC, Tenreyro S. On the Existence of the Maximum Likelihood Estimates in Poisson Regression. Economics Letters. 2010;107:310‒312.
9. Tiku ML. Estimating the mean and standard deviation from a censored normal sample. Biometrika. 1967;54(1):155‒165.
10. Bhattacharyya GK. The asymptotics of maximum likelihood and related estimators based on type II censored data. Journal of the American Statistical Association. 1985;80:398‒404.
11. Tan WY. On Tiku’s robust procedure‒a Bayesian insight. Journal of Statistical Planning and Inference. 1985;11(3):329‒340.
12. Tiku ML, Tan WY, Balakrishnan N. Robust Inference. Marcel Dekker, New York, USA. 1986.
13. Tiku ML, Suresh RP. A new method of estimation for location and scale parameters. Journal of Statitical Planning and Inference. 1992;30:281‒292.
14. Vaughan DC. On the Tiku‒Suresh method of estimation. Communications in Statistics‒Theory and Methods. 1992;21(2):451‒469.
15. Vaughan DC, Tiku ML. Estimation and hypothesis testing for nonnormal bivariate distribution with applications. Mathematical and Computer Modelling. 2000;32:53‒67.
16. Oral E, Gunay S. Stochastic Covariates in Binary Regression. Hacettepe Journal of Mathematics and Statistics. 2004;33:97‒109.
17. Oral E. Binary Regression with Stochastic Covariates. Communications in Statistics Theory and Methods. 2006;35:1429‒1447.
18. Oral E. Parameter Estimation in Generalized Linear Models through Modified Maximum Likelihood. International Statistical Institute. Proceedings of the 58th World Statistical Congress. 2011.
19. Agresti A. Categorical Data Analysis. John Wiley and Sons, New York, USA. 1996.
20. Kwan R. Lee, Kapadia C H , Dwight B. et al. On Estimating the Scale Parameter of the Rayleigh distribution from Doubly Censored Samples. Statistische Hefte. 1980;21(1):14‒29.
21. Cantoni E, Ronchetti E. Robust Inference for Generalized Linear Models. Journal of the American Statistical Association. 2001;96(455):1022‒1030.