Regression analysis under the assumption of monotonicity is a well-studied statistical problem and has been used in a wide range of applications. However, there remains a lack of a broadly applicable methodology that permits informa-tion borrowing, for efficiency gains, when jointly estimating multiple monotonic regression functions. We fill this gap in the literature and introduce a methodology which can be applied to both fixed and random designs and any number of explanatory variables (regressors). Our framework penalizes pairwise differences in the values of the monotonic function estimates, with the weight of penalty being determined, for instance, based on a statistical test for equivalence of functions at a point. Function estimates are subsequently derived using an iterative optimization routine which updates the individual function estimates in turn until convergence.Simulation studies for normally and binomially distributed response data illustrate that function estimates are improved when similarities between functions exist, and are not over-smoothed otherwise. We further apply our methodology to analyze two public health data sets: neonatal mortality data for Porto Alegre, Brazil, and stroke patient data for North West England.