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.
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?