In this report, we develop a procedure to analyze the relationship between the observed multi-dimensional counts and a set of explanatory variables. The counts follow a multivariate Poisson distribution or a multivariate zero-inflated Poisson distribution. Maximum likelihood estimates (MLE) for the model parameters are obtained by the Newton-Raphson (NR) iteration and the expectation-maximization (EM) algorithm, respectively. In Newton-Raphson method, the first and second derivatives of the log-likelihood function are derived to carry out the numerical evaluation. Formulas using EM algorithm are also introduced. A comparison of the estimation performance is made from simulation studies.