We analyze neutrinoless double beta decay (0νββ) within the framework of the Standard Model Effective Field Theory. Apart from the dimension-five Weinberg operator, the first contributions appear at dimension seven. We classify the operators and evolve them to the electroweak scale, where we match them to effective dimension-six, -seven, and -nine operators. In the next step, after renormalization group evolution to the QCD scale, we construct the chiral Lagrangian arising from these operators. We develop a power-counting scheme and derive the two-nucleon 0νββ currents up to leading order in the power counting for each lepton-number-violating operator. We argue that the leading-order contribution to the decay rate depends on a relatively small number of nuclear matrix elements. We test our power counting by comparing nuclear matrix elements obtained by various methods and by different groups. We find that the power counting works well for nuclear matrix elements calculated from a specific method, while, as in the case of light Majorana neutrino exchange, the overall magnitude of the matrix elements can differ by factors of two to three between methods. We calculate the constraints that can be set on dimension-seven lepton-number-violating operators from 0νββ experiments and study the interplay between dimension-five and -seven operators, discussing how dimension-seven contributions affect the interpretation of 0νββ in terms of the effective Majorana mass mββ.