We aim to promote the use of the modified profile likelihood function for estimating the variance parameters of a GLMM in analogy to the REML criterion for linear mixed models. Our approach is based on both quasi-Monte Carlo integration and numerical quadrature, obtaining in either case simulation-free inferential results. We will illustrate our idea by applying it to regression models with binary responses or count data and independent clusters, covering also the case of two-part models. Two real data examples and three simulation studies support the use of the proposed solution as a natural extension of REML for GLMMs. An R package implementing the methodology is available online