# How are fitted Poisson means constrained to be positive when the identity link is used in Poisson regression?

I apologise for the imprecise notation that follows, but hopefully I have conveyed the idea sufficiently. In Poisson regression of $Y \sim \mathbf{x}$ , the canonical link function $\ln$ constrains the fitted mean parameter $\mu_i$ of the $i$-th observation. to be positive, as:

$\mu_i = \exp(\sum_j \beta_j x_{ji}$)

However, when the identity link is chosen:

$\mu_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \dots$

Some combinations of $\beta$'s obviously make the RHS non-positive, (e.g. $\beta_0 = -1, \beta_i = 0, \forall i > 0$), however this produces a nonsensical Poisson mean. My question is: when fitting Poisson models with the identity link, how in practice are the fitted means constrained to be positive?

• I am not sure whether I understand this question: isn't a model with the identity link (instead of the canonical $ln$ link) just a regular linear regression model (and not a Poisson model)?
– IWS
Feb 28 '17 at 10:52
• No it is not: it still has a poisson likelihood function Feb 28 '17 at 11:17
• If it's properly implemented, fitted means (the mean for the fitted model at an observation) cannot go negative (recall we're maximizing likelihood, so think about at the likelihood function; likelihood is 0 there, log-likelihood is $-\infty$). Indeed even when the likelihood is maximized by a fitted value at 0, a well implemented fitting algorithm should (and most programs generally will) warn you if $\mu_i=0$ at any data point. For example in R you can get warning messages like: "algorithm stopped at boundary value" and "fitted rates numerically 0 occurred" -- ...ctd Feb 28 '17 at 12:14
• ctd... if you actually ignore those warnings and act as if all is okay, your results may mean little. Careful assessment would be needed in that situation. However, R"s algorithm for GLMs is less than ideal (see the package glm2 and its associated paper, for example, which discusses some of the problems. That package can achieve convergence in some cases where glm fails). But even with a perfect algorithm, predicted values not at observations can easily be negative. Feb 28 '17 at 12:34
• I guess that if all the beta's could be constrained to be positive, then the predictions would also remain positive also with the identity link. But the standard glm function doesn't support bounds/nonnegativity. glmnet does, but then that one doesn't support poisson with an identity link. ordinis does though (github.com/jaredhuling/ordinis) - you could have a try with that? Aug 15 '18 at 10:29

The standard glm() implementation in R with family=poisson(link=identity)) tries to fit this model using iteratively reweighted least squares and that is known to be unstable and often fail, because it does not constrain coefficients to be positive, which is the only way to ensure that the overall predicted values would remain positive and be sensible for ANY covariate value.

There is the nnpois function in the addreg package, however, which does fit nonnegative identity-link Poisson models using an EM-algorithm. The computational details can be found in Marschner, I. C. (2010). Stable computation of maximum likelihood estimates in identity link Poisson regression. Journal of Computational and Graphical Statistics 19(3): 666--683. Only downside of that algorithm is that it constrains ALL coefficients to be positive, and does not allow some to be zero. Hence, the solution that it returns is not sparse.

As an alternative, there is the nnlm function in the NNLM package, where you can specify loss="mkl" corresponding to Kullback-Leibler divergence, which then amounts to fitting a nonnegative identity-link Poisson model. In that case, the model is fit using a coordinate descent algorithm, and starting values are required (e.g. nnls estimates).

There is also the glm.cons function in the zetadiv package, which allows you to constrain some or all coefficients of your GLM fit to be positive, see https://www.rdocumentation.org/packages/zetadiv/versions/1.1.1/topics/glm.cons. This algorithm uses the IRLS algorithm but replaces the standard weighted least squares fit step by a weighted nonnegative least squares fit one (nnpls function) to produce only positive coefficients. In contrast to addreg this does return a sparse solution.

Finally, there is also the restriktor package, which can fit GLMs with particular constraints on individual coefficients. Instead of nnpls it uses quadratic programming in the IRLS algorithm to get the right constraints and like glm.cons does produce a sparse solution.

Out of these options only the restriktor package provides correct inference out of the box (using Monte Carlo approaches). The glm.cons function is a bit deceptive in that the summary of your fit provides standard errors and p values, but these are not correct, as they are calculated from the standard Fisher information matrix, but that's not a correct thing to do when many of the coefficients are on the constraint boundary. For all methods above you could use nonparametric bootstrapping though to get correct confidence intervals and p values on your coefficients, e.g. via the boot package.

Hope this helps?

The fitted means are not constrained to be positive, which is a potential problem for using the identity link function in that circumstance.

• I understand that depending on particular combinations of covariates, particularly when interpolating and extrapolating the model, one my obtain means that are not positive. What I did not understand was how the fitting process itself constrained means to be non-negative.
– Alex
Feb 28 '17 at 22:52
• It does not do so. When you get negative predictions, then the model typically will have trouble converging. Mar 1 '17 at 8:38
• The npois function from the addreg package and nnlm from the NNLM package can stably fit nonnegative identity-link Poisson regression models (the standard GLM algo in R indeed often fails for identity link Poisson models because it cannot constrain coefficients to be positive). I've added this as an answer below! Apr 25 '19 at 15:28

One way you could (I won't say should) constrain all the observed sets of predictions to be non-negative is to use Sequential Quadratic Programming (SQP).

To explain how SQP works, we can backstep to thinking about how Newton's method works. One way to view it is we make a quadratic approximation of the function given the current derivatives. Quadratic functions have a closed form solution for maximization, so we use that solution as our next step. We then retake the derivatives and repeat.

Now, quadratic programming (QP) is a method for finding a global max/min of a quadratic function with linear equality or inequality constraints. We can make the constraints be all the different combinations of the covariates in our dataset and constrain them to be non-negative. Then, we can repeatedly make a quadratic approximation of our function given the current derviatives, maximize this quadratic function conditional on the constraint that all predictors are non-negative via quadratic programming, update the derivatives and repeat. This is SQP and, under concavity of the log-likelihood function, is provable to find the global max under the linear constraints.

Note that this just guarantees that none of the observed set of predictor result in a negative estimated mean, but does not have such guarantees about a new set of predictors.

Or we can use the log-link.