By representing documents as mixtures of topics, topic modeling has allowed the successful analysis of datasets across a wide spectrum of applications ranging from ecology to genetics. An important body of recent work has demonstrated the computational and statistical efficiency of probabilistic Latent Semantic Indexing (pLSI)-- a type of topic modeling -- in estimating both the topic matrix (corresponding to distributions over word frequencies), and the topic assignment matrix. However, these methods are not easily extendable to the incorporation of additional temporal, spatial, or document-specific information, thereby potentially neglecting useful information in the analysis of spatial or longitudinal datasets that can be represented as tensors. Consequently, in this paper, we propose using a modified higher-order singular value decomposition (HOSVD) to estimate topic models based on a Tucker decomposition, thus accommodating the complexity of tensor data. Our method exploits the strength of tensor decomposition in reducing data to lower-dimensional spaces and successfully recovers lower-rank topic and cluster structures, as well as a core tensor that highlights interactions among latent factors. We further characterize explicitly the convergence rate of our method in entry-wise $\ell_1$ norm. Experiments on synthetic data demonstrate the statistical efficiency of our method and its ability to better capture patterns across multiple dimensions. Additionally, our approach also performs well when applied to large datasets of research abstracts and in the analysis of vaginal microbiome data.