Layered queueing networks (LQNs) are a class of performance models for software systems in which multiple distributed resources may be possessed simultaneously by a job. Estimating response times in a layered system is an essential but challenging analysis dimension in Quality of Service (QoS) assessment. Current analytic methods are capable of providing accurate estimates of mean response times. However, accurately approximating response time distributions used in service-level objective analysis is a demanding task. This paper proposes a novel hybrid framework that leverages phase-type (PH) distributions and neural networks to provide accurate density estimates of response times in layered queueing networks. The core step of this framework is to recursively obtain response time distributions in the submodels that are used to analyse the network by means of decomposition. We describe these response time distributions as a mixture of density functions for which we learn the parameters through a Mixture Density Network (MDN). The approach recursively propagates MDN predictions across software layers using PH distributions and performs repeated moment-matching based refitting to efficiently estimate end-to-end response time densities. Extensive numerical experiment results show that our scheme significantly improves density estimations compared to the state-of-the-art.