Adam -

Do you mean that "bptest" assumes heteroscedasticity is due to a single independent variable? Well, you can see from https://www.researchgate.net/publication/320853387_Essential_Heteroscedasticity that (essential) heteroscedasticity is caused by the change in prediction values, which are considered size measurements in that paper, which discusses a concept mid-page 111 in Brewer, K.R.W.(2002), Combined Survey Sampling Inference: Weighing Basu's Elephants, Arnold: London and Oxford University Press. Heteroscedasticity can be effectively reduced or enlarged in impact by model issues, such as omitted variables, or data/population groupings, or data quality issues, but basic (essential) heteroscedasticity is a function of prediction size and the coefficient of heteroscedasticity.

Using z (sub_i) for a size variable, where the best one is predicted y (sub_i), and gamma is the coefficient of heteroscedasticity here, we have for the estimated residuals, e (sub_i),

e = (e0)(z^gamma),

where e0 is the random factor of a given estimated residual. The regression weight is then

w = z^(-2gamma).

One of the updates to https://www.researchgate.net/project/OLS-Regression-Should-Not-Be-a-Default-for-WLS-Regression is with regard to there being no reason to do an hypothesis test for heteroscedasticity in regression. The question is not is there or is there not heteroscedasticity, but how much is there, which could include none, although that could indicate you have some other problem. If you do a hypothesis test and find "none," does that really mean there is none? If you find some, how much is there? What you really want to know is what is a reasonable amount of heteroscedasticity for your data, and then see how that impacts your results. The possible impact on your results is what matters here. (Impact is mostly with regard to variance, but in some cases predictions might change substantially.) Using the coefficient of heteroscedasticity smooths 'rough spots' in your calculations caused by natural fluctuations, and might help more with undetected or unproven outliers. This is covered, along with examples, in https://www.researchgate.net/publication/333642828_Estimating_the_Coefficient_of_Heteroscedasticity. Please particularly see the summary there. You can apply this to your own data using https://www.researchgate.net/publication/333659087_Tool_for_estimating_coefficient_of_heteroscedasticityxlsx.

So if I understand your question correctly Adam, "fitted" is what you want. Looking at just one regressor variable may not give you a good picture. But what effect size do you then have? What next? How much difference does it make? In which cases? - I suggest that you use your data to see what coefficient of heteroscedasticity is reasonable. Also see comments in the "tool" spreadsheet regarding default values. The theory is discussed in the papers above, and examples are given. After deciding on a coefficient of heteroscedasticity, gamma, use that in the regression weight. Hopefully the software you use will allow you to enter the regression weight. Then coefficients which minimize the sum of the squared random factors of the estimated residuals will be used in your regression of form y = y* + e.

Cheers - Jim

There is a third option for the Breusch-Pagan test in R language: function ols_test_breusch_pagan from the library olsrr.

Regarding whether "model terms" or "fitted values" are the basis for heteroscedasticity in regression:

I'll try to explain this a bit more explicitly, having thought about this question and being more certain of what is being asked. As can be seen in https://www.researchgate.net/publication/320853387_Essential_Heteroscedasticity, which discusses a concept from mid-page 111 in Brewer, K.R.W.(2002), Combined Survey Sampling Inference: Weighing Basu's Elephants, Arnold: London and Oxford University Press, heteroscedasticity in regression occurs naturally due to differences in size of members of a population, with regard to variable y. The size associated with each y is best measured by the best predicted y. There can be problems with obtaining the best predicted y (i.e., 'fitted') values, omitted variables for example, and there are other influences due to data quality and stratification, but basic ("essential") naturally occurring heteroscedasticity is due to the differing values of predicted y. Model and data problems may either mask or inflate this, but homoscedasticity for estimated residuals is not the ideal situation. Many examples of OLS regression people see when learning about regression use artificial data contrived to be homoscedastic, or are actually heteroscedastic, and this was not noticed. (See real data examples in https://www.researchgate.net/publication/333642828_Estimating_the_Coefficient_of_Heteroscedasticity, for regressions of form y = y* + e, for predictions using finite populations.) But under ideal circumstances, we can expect heteroscedasticity as long as members of a population vary in size. (This applies separately for every dependent variable, so that if your population is farms in a given region, and say you have data for crop yields for several different crops, and you are estimating total crop yield for the region for each of those crops, then some farms are large for some yields and small or even zero for others. Each one is looked at separately. As a further aside, that means that unequal probability sampling that works well for one crop would be terrible for another. However, using prediction, the best weighted least squares regression would be available for each crop separately, using a sample that is a compromise including at least the largest members for each crop.)

But back to "fitted values," one can expect that a larger such value will have a larger sigma for associated estimated residuals, intuitively, in addition to the mathematical/statistical justification in the first references above.

Anyway, the answer to the question at hand is "fitted values." One independent variable might often seem 'responsible' for most of this, but the best answer is "fitted values."