To ensure agreement between theoretical calculations and experimental data, parameters to selected nuclear physics models are perturbed and fine-tuned in nuclear data evaluations. This approach assumes that the chosen set of models accurately represents the ‘true’ distribution of considered observables. Furthermore, the models are chosen globally, indicating their applicability across the entire energy range of interest. However, this approach overlooks uncertainties inherent in the models themselves. In this work, we propose that instead of selecting globally a winning model set and proceeding with it as if it was the ‘true’ model set, we, instead, take a weighted average over multiple models within a Bayesian model averaging (BMA) framework, each weighted by its posterior probability. The method involves executing a set of TALYS calculations by randomly varying multiple nuclear physics models and their parameters to yield a vector of calculated observables. Next, computed likelihood function values at each incident energy point were then combined with the prior distributions to obtain updated posterior distributions for selected cross sections and the elastic angular distributions. As the cross sections and elastic angular distributions were updated locally on a per-energy-point basis, the approach typically results in discontinuities or “kinks” in the cross section curves, and these were addressed using spline interpolation. The proposed BMA method was applied to the evaluation of proton-induced reactions on ⁵⁸Ni between 1 and 100 MeV. The results demonstrated a favorable comparison with experimental data as well as with the TENDL-2023 evaluation.