Something to be aware of is that Abaqus uses a C0 formulation (at least for most cases). This means that while the solution variable (displacements in your case) are continuous at the element boundaries (i.e. nodes) -- meaning that there is one and only one solution value at each node -- that stress being a derived value (i.e. proportional to derivatives of the solution) is discontinuous (C-1) at the element boundaries -- meaning that each element that shares a node provides its own unique stress value at the node. This is why Abaqus uses a volume-weighted average in its contour plots to make stress appear continuous, when in fact it isn't. This is one of the reasons why it is typically recommended to query stresses from quadrature points. As a side note, from theory it's clear to see that derived values will converge at slower rates than solution values - if you would expect displacement to converge at a rate of hp+1 then stress would converge at a rate of hp -- so you need a lot more elements to converge to stress than are needed to converge to displacement.
Anyways, I've never heard of the advice to "average the centroidal stresses of the neighboring elements" because "you can't trust results from a single element". The truth is, for a well-defined FEM-problem you absolutely can trust the results from a single element as you can show that the best-approximation property holds for derived values. The caveat to this is contained within well-defined as stress-singularities, elements at a boundary-condition, locking, etc. can have issues to be sure. What I typically advise is to look at a quilt plot of stresses and if you see a region that has either oscillatory stresses or a stress that qualitatively looks out-of-place compared to its neighbors, to then explore further to see whether you have a singularity, locking, or haven't yet converged to accurate stress results.