Relative risks (RRs) are generally considered preferable to odds ratios in prospective studies. However, unlike logistic regression for odds ratios, the standard log-binomial model for RR regression does not respect the natural parameter constraints and is therefore often subject to numerical instability. In this paper, we develop a reliable and flexible method for fitting log-binomial models. We use an Expectation-Maximization (EM) algorithm where the multiplicative event probability is viewed as the joint probability for a collection of latent binary outcomes. This gives a simple iterative scheme that provides stable convergence to the maximum likelihood estimate. In addition to reliability, the method offers some flexible generalizations, including models with unspecified isotonic regression functions. We examine the method's performance using simulations and data analyses of the age-specific RR of mortality following heart attack. These analyses demonstrate the potential for numerical instability in RR regression and show how this can be overcome using the proposed approach. Source code to implement the method in R is provided as supplementary material available at Biostatistics online.