The analysis of overdispersed counts has been the focus of a large amount of literature, with the general objective of providing reliable parameter estimates in the presence of heterogeneity or dependence among subjects. In this paper we extend the standard variance component models to the analysis of multivariate counts, defining the dependence among counts through a set of correlated random coefficients. Estimation is carried out by numerical integration through an EM algorithm without parametric assumptions upon the random coefficients distribution. The proposed model is computationally parsimonious and, when applied to a real dataset, seems to produce better results than parametric models. A simulation study has been carried out to investigate the behavior of the proposed models in a series of empirical situations.