Traditional multivariate methods often lead to overfitting when the structure of the data is not considered. In this paper, we propose a family of regression models that takes tensor responses and covariates of any number of dimensions. We deal with the curse of dimensionality in two ways: first, by imposing various low-rank tensor structures on one of the regression parameters (the systematic component), and second, by modeling the errors (the random component) using a family of elliptically contoured distributions with Kronecker separable covariance. We provide iterative algorithms for obtaining the maximum likelihood estimators (MLEs) of our parameters under both tensor normal and tensor t distributed errors, as well for obtaining identifiable MLEs of the parameters involved in the tensor normal and tensor t distributions. We study the exact and asymptotic distributions of our parameters. The methodology is applied to a two factor tensor ANOVA of functional magnetic resonance imaging, where the first factor identifies the control vs suicidal patients and the second factor corresponds to stimulus concepts that are either positive, negative or related to suicide.