  Modified theories of gravity encompass a class of $f(R)$-models that seek to
elucidate the observed late time accelerated expansion of the universe. In this
study, we examine a set of viable $f(R)$ models (Hu-Sawicki: two cases,
Satrobinsky, Tsujikawa, exponential and arcTanh models) in metric formalism,
using recent cosmological data sets: type Ia supernovae data, cosmic
chronometer observations, baryonic acoustic oscillations data, data from
H\textsc{ii} starburst galaxies, and local measurements of the Hubble parameter
$H_0$. The model parameters are constrained using a Bayesian analysis with the
Monte Carlo Markov Chain method. We employ statistical tools such as the Akaike
Information Criterion, Bayesian Information Criterion, and reduced chi-square
statistics to conduct a comparative investigation of these models. We determine
the transition redshift, the evolution of total equation-of-state (EoS)
parameter, and the EoS for the component responsible for current accelerated
expansion to characterize the expansion's evolution. Taking into account the
``Hubble tension," we perform the study with and without a Gaussian prior for
$H_0$ from local measurements. Our findings are as follows: (i) in many cases
the $f(R)$ models are strongly favored over the standard $\Lambda$CDM model,
(ii) the deviation parameter ($b$) significantly deviates from zero in several
cases, (iii) the inclusion of local $H_0$ not only increases the fitted value
of $H_0$ (as expected) but also affects the gap between predictions of $f(R)$
models and the $\Lambda$CDM model, and (iv) the relevant quantities
characterizing the (accelerated) expansion of the universe obtained in our
models are consistent with those obtained in a model-independent way by others.
Our investigation and results present a compelling case for pursuing further
research on $f(R)$ models with future observations to come.
