Why is a Breusch-Pagan test returning significant heteroskedasticity when the fitted value chart indicates homoskedasticity?

Question

I was working the results of a regression equation and I wanted to test to see if there was significant heteroskedasticity in the residuals. Checking the results of the graph of fitted values versus residuals using plot() produced a straight line that is almost exactly zero for most of the resulting graph and the data fits rather closely.

However, a studentized Breusch-Pagan test of the residuals using the bptest() function indicated significant heteroskedasticity (BP = 8.4085, df = 1, p = 0.0037). I can't figure out why the Breusch-Pagan test is returning significant heteroskedasticity. The regression doesn't look heteroskedastic, the regression line is pretty consistent overall but there are a few outliers in the right half of the graph. There isn't a clearly non-linear pattern like increasing variance in residuals with increasing fitted values like one would expect from heteroskedastic results. So then why is the Breush-Pagan test telling me that I have significant heteroskedasticity in my regression?

Answer 1

I see some evidence of heteroscedasticity, with seemingly larger variability to the right. More importantly, there is a suggestion of nonlinearity (not related to heteroscedasticity) that can be better visualized by overlaying the loess smooth.

But tests for model specification are mostly irrelevant. If the BP test "passed", it is simply a Type II error, because heteroscedasticity is true in reality when $Y$ and $X$ are related. Rather than test for heteroscedasticity, why not just model it? Then you will have a better idea as to the size of the effect. Use that information, along with subject matter considerations, supplemented with simulation study as needed, to decide whether to ignore it.