We construct higher order spectral shift functions, extending the perturbation theory results of M.G. Krein [M.G. Krein, On a trace formula in perturbation theory, Mat. Sb. 33 (1953) 597-626 (in Russian)] and L.S. Koplienko [L.S. Koplienko, Trace formula for perturbations of nonnuclear type, Sibirsk. Mat. Zh. 25 (1984) 62-71 (in Russian); translation in: Trace formula for nontrace-class perturbations, Siberian Math. J. 25 (1984) 735-743] on representations for the remainders of the first and second order Taylor-type approximations of operator functions. The higher order spectral shift functions represent the remainders of higher order Taylor-type approximations; they can be expressed recursively via the lower order (in particular, Krein's and Koplienko's) ones. We also obtain higher order spectral averaging formulas generalizing the Birman-Solomyak spectral averaging formula. The results are obtained in the semi-finite von Neumann algebra setting, with the perturbation taken in the Hilbert-Schmidt class of the algebra.