Robust and fast in vivo treatment verification is expected to increase the clinical efficacy of proton therapy. The combined detection of prompt gamma rays and neutrons has recently been proposed for this purpose and shown to increase the monitoring accuracy. However, the potential of this technique is not fully exploited yet since the proton range reconstruction relies only on a simple landmark of the particle production distributions. Here, we apply machine learning based feature selection and multivariate modelling to improve the range reconstruction accuracy of the system in an exemplary lung cancer case in silico. We show that the mean reconstruction error of this technique is reduced by 30%–50% to a root mean squared error per spot of 0.4, 1.0, and 1.9 mm for pencil beam scanning spot intensities of 10⁸, 10⁷, and 10⁶ initial protons, respectively. The best model performance is reached when combining distribution features of both gamma rays and neutrons. This confirms the advantage of hybrid gamma/neutron imaging over a single-particle approach in the presented setup and increases the potential of this system to be applied clinically for proton therapy treatment verification.