This paper provides a unified algorithm to explicitly calculate the maximum likelihood estimates of parameters in a general setting of generalised linear mixed models (GLMMs) in terms of Gauss-Hermite quadrature approximation. The score function and observed information matrix are expressed explicitly as analytically closed forms so that Newton-Raphson algorithm can be applied straightforwardly. Compared with some existing methods, this approach can produce more accurate estimates of the fixed effects and variance components in GLMMs, and can serve as a basis of assessing existing approximations in GLMMs. A simulation study and practical example analysis are provided to illustrate this point.