Maximum likelihood estimation in random effects models for non-Gaussian data is a computationally challenging task that currently receives much attention. This article shows that the estimation process can be facilitated by the use of automatic differentiation, which is a technique for exact numerical differentiation of functions represented as computer programs. Automatic differentiation is applied to an approximation of the likelihood function, obtained by using either Laplace's method of integration or importance sampling. The approach is applied to generalized linear mixed models. The computational speed is high compared to the Monte Carlo EM algorithm and the Monte Carlo Newton-Raphson method.