We present an extension of the Kuramoto-Sakaguchi model for networks, deriving the second-order phase approximation for a paradigmatic model of oscillatory networks-an ensemble of nonidentical Stuart-Landau oscillators coupled pairwisely via an arbitrary coupling matrix. We explicitly demonstrate how this matrix translates into the coupling structure in the phase equations. To illustrate the power of our approach and the crucial importance of high-order phase reduction, we tackle a trendy setup of nonlocally coupled oscillators exhibiting a chimera state. We reveal that our second-order phase model reproduces the dependence of the chimera shape on the coupling strength that is not captured by the typically used first-order Kuramoto-like model. Our derivation contributes to a better understanding of complex networks' dynamics, establishing a relation between the coupling matrix and multibody interaction terms in the high-order phase model.