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?
 A: 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?
A: The fitted means are not constrained to be positive, which is a potential problem for using the identity link function in that circumstance. 
A: 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. 
