Networks of interacting oscillators abound in nature, and one of the prevailing challenges in science is how to characterize and reconstruct them from measured data. We present a method of reconstruction based on dynamical Bayesian inference that is capable of detecting the effective phase connectivity within networks of time-evolving coupled phase oscillators subject to noise. It not only reconstructs pairwise, but also encompasses couplings of higher degree, including triplets and quadruplets of interacting oscillators. Thus inference of a multivariate network enables one to reconstruct the coupling functions that specify possible causal interactions, together with the functional mechanisms that underlie them. The characteristic features of the method are illustrated by the analysis of a numerically generated example: a network of noisy phase oscillators with time-dependent coupling parameters. To demonstrate its potential, the method is also applied to neuronal coupling functions from single- and multi-channel electroencephalograph recordings. The crossfrequency δ, α to α coupling function, and the θ, α, γ to γ triplet are computed, and their coupling strengths, forms of coupling function, and predominant coupling components, are analysed. The results demonstrate the applicability of the method to multivariate networks of oscillators, quite generally.