In this paper we propose a method for estimating the covariance matrix of random effects in mixed effect models. We adopt a stochastic approach to approximating the gradient of the log marginal likelihood function. This approach uses a Metropolis-Hastings algorithm to sample random effects from distributions conditional on observed information. In particular, to make the algorithm efficient, this approach uses a data-driven proposal for sampling random effects from target distributions. As a result of that, this approach can achieve high accuracy in approximating the log marginal likelihood. We apply the approach to estimate the covariance matrix of the random effects in various models, including generalized linear mixed models and the Gaussian graphical model.