This paper presents a multivariate functional data statistical approach, for spatiotemporal prediction of COVID-19 mortality counts. Specifically, spatial heterogeneous nonlinear parametric functional regression trend model fitting is first implemented. Classical and Bayesian infinite-dimensional log-Gaussian linear residual correlation analysis is then applied. The nonlinear regression predictor of the mortality risk is combined with the plug-in predictor of the multiplicative error term. An empirical model ranking, based on random K-fold validation, is established for COVID-19 mortality risk forecasting and assessment, involving Machine Learning (ML) models, and the adopted Classical and Bayesian semilinear estimation approach. This empirical analysis also determines the ML models favored by the spatial multivariate Functional Data Analysis (FDA) framework. The results could be extrapolated to other countries.