We provide a method for fitting monotone polynomials to data with both fixed and random effects. In pursuit of such a method, a novel approach to least squares regression is proposed for models with functional constraints. The new method is able to fit models with constrained parameter spaces that are closed and convex, and is used in conjunction with an expectation-maximisation algorithm to fit monotone polynomials with mixed effects. The resulting mixed effects models have constrained mean curves and have the flexibility to include either unconstrained or constrained subject-specific curves. This new methodology is demonstrated on real-world repeated measures data with an application from sleep science. Code to fit the methods described in this paper is available online.