We consider a general monotone regression estimation where we allow for independent and dependent regressors. We propose a modification of the classical isotonic least squares estimator and establish its rate of convergence for the integrated $L_1$-loss function. The methodology captures the shape of the data without assuming additivity or a parametric form for the regression function. Furthermore, the degree of smoothing is chosen automatically and no auxiliary tuning is required for the theoretical analysis. Some simulations and two real data illustrations complement the study of the proposed estimator.