In this talk I am going to present a recent proposed regression modelling framework to deal with multivariate count data. In this framework the covariance structure for each response variable is defined in terms of a covariance link function combined with a matrix linear predictor involving known matrices. To specify the joint covariance matrix for the multivariate count response vector, we employ the generalized Kronecker product. The count nature of the data is taken into account by means of the power dispersion function associated with the Poisson-Tweedie distribution. This specification provides a flexible and efficient multivariate regression methodology for a comprehensive family of count models including multivariate analagous to the Hermite, Neyman Type A, Pólya–Aeppli, negative binomial and Poisson-inverse Gaussian. We discuss extensions of the orthodox multivariate analysis of variance (MANOVA) for count data. Furthermore, we present the computational implementation in R through the package mcglm. Illustrations include a five response variables regression model in the context of health services care and a bivariate longitudinal data in the context of bushmeat in Pico Basilé, Bioko Island, Equatorial Guinea.