Why is there heteroskedasticity even though no relationship seems present in the residual plot?

Question

I estimated a random effects panel model and performed the Breusch-Pagan (BP) test for heteroskedasticity. The test is significant, meaning that there is heteroskedasticity. However, the residual plot seems not to have a clear fanning-in our fanning-out pattern.

Also, even after removing outliers with very low or high residuals, the BP test remains significant. How should I interpret the significance of the test that is in contrast with the residual plot?

Residual plot before removing outliers (sig BP test):

Residual plot after removing outliers (also sig BP test):

The model I use is the following:

plm(ln(gini_eurostat) ~ 
    ln(intraEU_trade_bymemberstate_pct) + 
    FDI_in_pctGDP + FDI_out_pctGDP +
    unempl_pct + 
    lnpop + 
    marketcap_pctGDP + 
    naturalresources_pctGDP + 
    ln(socialbenefits) + 
    ln(tech_exports_pctexports) +
    ln(GDP_percap_constantLCU) + 
    ln(intraEU_trade_bymemberstate_pct) * 
    lowGDP_percap_currentUSD,
    data=pdf, model="random")

Answer 1

First, focus onthe first plot, without outlier removal: To me it is clear that the spread of the points around the horizontal zero line is increasing with increasing fitted values! If you calculate variance of residuals by intervals --- or better, plot a windowed smooth of the variance above your plot --- that would be clear.

That you do not see this so clearly might be because the spread of points horizontally is somewhat clumped, with interacts with the view of the vertical spread. You could investigate this with further plots, to see if the variance heterogeneity is associated with some specific predictor.

Edit I digitised the values from the first plot in the post (without outlier removal), using the R (CRAN) package metaDigitise. Below is the reconstructed plot, wit added smooths (using R function lowess(locally weighted scatterplot smoother), of the residuals, and of the residuals squared, the last giving a local smooth estimator of variance. In the plot its square root:

The non-constancy is rather clear. Below the used R code:

library(metaDigitise)
rawdata <- metaDigitise(dir="./Digitise/",
                        summary=FALSE) 
    # and a lot of clicking ...
plotdata <- with(rawdata[[1]][[1]], 
    data.frame(Residuals=y, Fitted_values=x))
smoothvar <- with(plotdata, lowess(Fitted_values, 
     Residuals^2, iter=0))

smoothvar.1 <- within(smoothvar,  y <-sqrt(y))
smoothvar.2 <- within(smoothvar,  y <- -sqrt(y))

with(plotdata, plot(Fitted_values, Residuals,
    main=paste("Residuals with smooth,  locally", 
    "smoothed standard dev (blue)", sep="\n")))
with(plotdata, lines(lowess(Fitted_values, 
     Residuals), col="red"))
lines(smoothvar.1, col="blue")
lines(smoothvar.2, col="blue")