In recent years, the Poisson lognormal mixed model has been frequently used in modeling count data because it can accommodate both the over-dispersion of the data and the existence of within-subject correlation. Since the likelihood function of this model is expressed in terms of an intractable integral, estimating the parameters and obtaining inference for the parameters are challenging problems. Some approximation procedures have been proposed in the literature; however, they are computationally intensive. Moreover, the existing studies of approximate parameter inference using the Gaussian variational approximation method are usually restricted to models with only one predictor. In this paper, we consider the Poisson lognormal mixed model with more than one predictor. By extending the Gaussian variational approximation method, we derive explicit forms for the estimators of the parameters and examine their properties, including the asymptotic distributions of the estimators of the parameters. Accurate inference for the parameters is also obtained. A real-life example demonstrates the applicability of the proposed method, and simulation studies illustrate the accuracy of the proposed method.