I look at which point is furthest from the generated fit and remove that.
Then I recalculate the fit and repeat until a decent regression coefficient is achieved and/or errors are below a threshold and sufficient stars remain.
Its not an ideal approach and has trouble when there are few stars in the scene, but the results I got suggested it was doing quite a good job – certainly much better than not culling the outliers. It may not make a huge difference generally, but its easy to code up and test so possibly worth a look as a potential refinement.
Had expected stars with the most extreme colour indices would be causing the data points far from the curve, but while they were not great they were not always the worst outliers.
Could probably dig out the code if you wanted it. I was using Gaia DR2 and cooled sensor data captured using a normal Silicon CCD.