The aim of this paper is to present a regression model for multivariate claim frequency data with dependence structures across the claim count responses, which may be of different sign and range, and overdispersion from the unobserved heterogeneity due to systematic effects in the data. For illustrative purposes, we consider the bivariate Poisson-lognormal regression model with varying dispersion. Maximum likelihood estimation of the model parameters is achieved through a novel Monte Carlo expectation–maximization algorithm, which is shown to have a satisfactory performance when we exemplify our approach to Local Government Property Insurance Fund data from the state of Wisconsin.