Refinements on the exact method to solve the numerical difficulties in fitting the log binomial regression model for estimating relative risk

Fitting a log binomial regression model using standard software can result in numerical difficulties because its functional form does not restrict the estimated probabilities to values not exceeding unity. The common approaches to resolve the issue are to introduce a constraint to limit the results of each iteration. However, if the ML solution lies on the boundary of the allowable parameter space, some fitted probabilities are equal to unity (named as boundary vector). The common approaches have trouble dealing with those boundary vectors and can only reach somewhere before the ML solution. Previously a remedy has been proposed, but without the details necessary to implement it. Here we provide these details, including formulae for estimating the covariances essential to implement the method, an explanation of inter-dependency between coefficient estimates, and a proof that the method can be implemented in general. Code written for R implements choice of fitting algorithms, finding appropriate starting values, identifying the covariate vector(s) with a fitted probability of unity, and strategies for covariate ordering to handle the issues caused by the common values in two distinct boundary vectors. The model-fitting results are compared with two alternative methods by simulation and example data.