  This study presents a comparative analysis of three predictive models with an
increasing degree of flexibility: hidden dynamic geostatistical models (HDGM),
generalised additive mixed models (GAMM), and the random forest spatiotemporal
kriging models (RFSTK). These models are evaluated for their effectiveness in
predicting PM$_{2.5}$ concentrations in Lombardy (North Italy) from 2016 to
2020. Despite differing methodologies, all models demonstrate proficient
capture of spatiotemporal patterns within air pollution data with similar
out-of-sample performance. Furthermore, the study delves into station-specific
analyses, revealing variable model performance contingent on localised
conditions. Model interpretation, facilitated by parametric coefficient
analysis and partial dependence plots, unveils consistent associations between
predictor variables and PM$_{2.5}$ concentrations. Despite nuanced variations
in modelling spatiotemporal correlations, all models effectively accounted for
the underlying dependence. In summary, this study underscores the efficacy of
conventional techniques in modelling correlated spatiotemporal data,
concurrently highlighting the complementary potential of Machine Learning and
classical statistical approaches.
