  We perform a comprehensive analysis of neutrinoless double beta decay and its
interplay with low-energy flavor observables in a radiative neutrino mass model
with scalar leptoquarks $S_1(\bar{3},1,1/3)$ and $\widetilde{R}_2(3,2,1/6)$. We
carve out the parameter region consistent with constraints from neutrino mass
and mixing, from collider searches, as well as from measurements of several
flavor observables, including muon and electron anomalous magnetic moments,
charged lepton flavor violation and rare (semi)leptonic kaon and $B$-meson
decays. We perform a global analysis to all existing constraints and show the
(anti)correlations between all relevant Yukawa couplings satisfying these
restrictions. We find that the most stringent constraint on the parameter space
comes from $\mu \to e$ conversion in nuclei. We also point out a tension
between the $R_{D^{(*)}}$ and muon $(g-2)$ anomalies in this context. Taking
benchmark points from these combined allowed regions, we study the implications
for neutrinoless double beta decay including the canonical light neutrino and
the leptoquark contributions. We find that for normal ordering of neutrino
masses, the leptoquark contribution removes the cancellation region that occurs
for the canonical case. The effective mass in presence of leptoquarks can lie
in the desert region between the standard normal and inverted ordering cases,
and this can be probed in future ton-scale experiments like LEGEND-1000 and
nEXO. Taking suitable representative values of Yukawa couplings, we also obtain
stringent constraints on leptoquark masses from different observables.
