We investigate neutrinoless double beta decay (0νββ) in the presence of sterile neutrinos with Majorana mass terms. These gauge-singlet fields are allowed to interact with Standard-Model (SM) fields via renormalizable Yukawa couplings as well as higher-dimensional gauge-invariant operators up to dimension seven in the Standard Model Effective Field Theory extended with sterile neutrinos. At the GeV scale, we use Chiral effective field theory involving sterile neutrinos to connect the operators at the level of quarks and gluons to hadronic interactions involving pions and nucleons. This allows us to derive an expression for 0νββ rates for various isotopes in terms of phase-space factors, hadronic low-energy constants, nuclear matrix elements, the neutrino masses, and the Wilson coefficients of higher-dimensional operators. The required hadronic low-energy constants and nuclear matrix elements depend on the neutrino masses, for which we obtain interpolation formulae grounded in QCD and chiral perturbation theory that improve existing formulae that are only valid in a small regime of neutrino masses. The resulting framework can be used directly to assess the impact of 0νββ experiments on scenarios with light sterile neutrinos and should prove useful in global analyses of sterile-neutrino searches. We per- form several phenomenological studies of 0νββ in the presence of sterile neutrinos with and without higher-dimensional operators. We find that non-standard interactions involving sterile neutrinos have a dramatic impact on 0νββ phenomenology, and next-generation experiments can probe such interactions up to scales of mathcal{O} (100) TeV.