This paper presents a two-step pseudo likelihood estimation technique for the generalized linear mixed models (GLMM) with random effects being correlated (possibly between subjects). Due to the use of the two-step estimation technique the proposed algorithm outperforms the conventional pseudo likelihood algorithms, e.g. Wolfinger and O’Connell (1993), in terms of computational time. Moreover, it does not require any reparametarisation of the model such as Lindstrom and Bates (1989). Multivariate Taylor’s approximation has been used to approximate the intractable integrals in the likelihood function of the GLMM. Based on the analytical expression for the estimator of the covariance matrix of the random effects, a condition has been presented as to when such a covariance matrix can be estimated through the estimates of the random effects. An application of the estimation technique with a binary response variable is presented using a real data set on credit defaults.