In this article we present and analyse new multilevel adaptations of stochastic approximation algorithms for the computation of a zero of a function $f\colon D \to \mathbb R^d$ defined on a convex domain $D\subset \mathbb R^d$, which is given as a parameterised family of expectations. Our approach is universal in the sense that having multilevel implementations for a particular application at hand it is straightforward to implement the corresponding stochastic approximation algorithm. Moreover, previous research on multilevel Monte Carlo can be incorporated in a natural way. This is due to the fact that the analysis of the error and the computational cost of our method is based on similar assumptions as used in Giles (2008) for the computation of a single expectation. Additionally, we essentially only require that $f$ satisfies a classical contraction property from stochastic approximation theory. Under these assumptions we establish error bounds in $p$-th mean for our multilevel Robbins-Monro and Polyak-Ruppert schemes that decay in the computational time as fast as the classical error bounds for multilevel Monte Carlo approximations of single expectations known from Giles (2008).