An infectious disease spreads not only over a single population or community but also across multiple and heterogeneous communities. Moreover, its transmissibility varies over time because of various factors such as seasonality and epidemic control, which results in strongly nonstationary behavior. In conventional methods for assessing transmissibility trends or changes, univariate time-varying reproduction numbers are calculated without taking into account transmission across multiple communities. In this paper, we propose a multivariate-count time series model for epidemics. We also propose a statistical method for estimating the transmission of infections across multiple communities and the time-varying reproduction numbers of each community simultaneously from a multivariate time series of case counts. We apply our method to incidence data for the novel coronavirus disease 2019 (COVID-19) pandemic to reveal the spatiotemporal heterogeneity of the epidemic process.