A general phase reduction method for a network of coupled dynamical elements exhibiting collective oscillations, which is applicable to arbitrary networks of heterogeneous dynamical elements, is developed. A set of coupled adjoint equations for phase sensitivity functions, which characterize phase response of the collective oscillation to small perturbations applied to individual elements, is derived. Using the phase sensitivity functions, collective oscillation of the network under weak perturbation can be described approximately by a one-dimensional phase equation. As an example, mutual synchronization between a pair of collectively oscillating networks of excitable and oscillatory FitzHugh-Nagumo elements with random coupling is studied.Networks of coupled dynamical elements exhibiting collective oscillations often play important functional roles in real-world systems. Here, a method for dimensionality reduction of such networks is proposed by extending the classical phase reduction method for nonlinear oscillators. By projecting the network state to a single phase variable, a simple one-dimensional phase equation describing the collective oscillation is derived. As an example, synchronization between collectively oscillating random networks of neural oscillators is studied. The derived phase equation is general and will have wide applicability in control and optimization of collectively oscillating networks.