  We present a Markov chain Monte Carlo pipeline that can be used for robust
and unbiased constraints of $f(R)$ gravity using galaxy cluster number counts.
This pipeline makes use of a detailed modelling of the halo mass function in
$f(R)$ gravity, which is based on the spherical collapse model and calibrated
by simulations, and fully accounts for the effects of the fifth force on the
dynamical mass, the halo concentration and the observable-mass scaling
relations. Using a set of mock cluster catalogues observed through the thermal
Sunyaev-Zel'dovich effect, we demonstrate that this pipeline, which constrains
the present-day background scalar field $f_{R0}$, performs very well for both
$\Lambda$CDM and $f(R)$ fiducial cosmologies. We find that using an incomplete
treatment of the scaling relation, which could deviate from the usual power-law
behaviour in $f(R)$ gravity, can lead to imprecise and biased constraints. We
also find that various degeneracies between the modified gravity, cosmological
and scaling relation parameters can significantly affect the constraints, and
show how this can be rectified by using tighter priors and better knowledge of
the cosmological and scaling relation parameters. Our pipeline can be easily
extended to other modified gravity models, to test gravity on large scales
using galaxy cluster catalogues from ongoing and upcoming surveys.
