Identity link Poisson regression is useful when the mean of a count variable depends additively on a collection of predictor variables. It is particularly important in epidemiology, for modeling absolute differences in disease incidence rates as a function of covariates. A complication of such models is that standard computational methods for maximum likelihood estimation can be numerically unstable due to the nonnegativity constraints on the Poisson means. Here we present a straightforward and flexible method that provides stable maximization of the likelihood function over the constrained parameter space. This is achieved by conducting a sequence of maximizations within subsets of the parameter space, after which the global maximum is identified from among the subset maxima. The method adapts and extends EM algorithms that are useful in specialized applications involving Poisson deconvolution, but which do not apply in more general regression contexts. As well as allowing categorical and continuous covariates, the method has the flexibility to accommodate covariates with an unspecified isotonic form. Its computational reliability makes it particularly useful in bootstrap analyses, which may require stable convergence for thousands of implementations. Computations are illustrated using epidemiological data on occupational mortality, and biological data on crab population counts. This article has supplementary material online.