I also plotted a residual plot and I suspected that there was evidence on heteroskedasticity. The Breush-Pagan test confirmed it. Could I be sure test's result? Should I have use other tests?

First, test whether the data fits to Gaussian (Normal) distribution. 

If YES, then Bartlett test is most powerful to detect heteroskedasticity.

If there is MINOR DEVIATION (see the Q-Q plot from test for normality) from normality, then use Levene test for heteroskedasticity.

If there is MAJOR DEVIATION (see the Q-Q plot from test for normality) from normality, then use either Fligner-Killeen test or Brown–Forsythe test for heteroskedasticity.  

Dear Francesca,

We recommend the Koenker test because it is the simplest and does not require assumptions to be checked first. One drawback is it is not easy to carry out in SPSS. However, I'm not aware of any test for heteroskedasticity in SPSS.

Indeed, there are some tests for heteroskedasticity that are more robust to the violation of the normality assumption of the error term in the regression model. So I also made the choice of the test based on this observation as suggested by Somenath.

You can first use the informal tests such as residual plots, qq-plots. Once these show some deviation from constant variance, I would suggest you use Bruesch Pagan (Lagrange Multiplier) test which gives the same results as Cook–Weisberg in Stata. For Levins, well, I have seen those mostly used in studies that use experimental and agronomic data. For social sciences, Bruesch Pagan, and White's test are likely used. But it's recommended that you know its form before deciding which test to use.

Besides the usual problems with interpreting the impact of sample size and perhaps competing sensitivities of various hypothesis tests, the basic problem here may be in assuming that you either do or do not have heteroscedasticity, when actually you are looking for a kind of effect size:  How much heteroscedasticity do you have?  I think an estimate of the regression weight could be more useful to you, though I don't know the exact nature of your application.   

Here's an idea: If you consider the regression weight to be a function of a size measure, then you can find an estimate of that weight function, perhaps experimentally.  (The size measure can be one of the independent variables, or a function of them.  Perhaps predicted-y may be a good choice.  For the simplest case with one regressor and regression through the origin, one just uses x.  But for most cases, I'd use predicted-y.)  The closer the weights are to a constant, the less the heteroscedasticity.  You mentioned using a residual plot.  If you used the absolute values of the estimated residuals, what weight would you need to divide out to have a regression line through those points look "flat?"  As an example, consider going from Figure 5 to Figure 6 in https://www.researchgate.net/publication/313925533_Comparing_Model_Performances_Graphically

A reference you might consult on regression weights:

Carroll and Ruppert(1988), Transformation and Weighting in Regression, Chapman & Hall, Ltd. London, UK.

Also, there is this:

 https://www.researchgate.net/publication/263032446_Weighting_in_Regression_for_Use_in_Survey_Methodology

remembering that predicted-y can be used for a size measure rather than the simple case of x.  (Note that w on page 2 was a bad choice of notation on my part, as it is used there for a value of the coefficient of heteroscedasticity, not the regression weight.  Sorry if that notation causes any confusion.) 

So the idea then would be to estimate the level/degree/amount of heteroscedasticity. 

Method Comparing Model Performances Graphically

Article Weighting in Regression for Use in Survey Methodology

Note: It is not clear to me if you are looking at heteroscedasticity (regarding variance of estimated prediction errors, where you are looking at conditional distributions of y given x, which is what I assumed), or else you may be looking to see if different populations sampled may have a difference in their variances (i.e., heterogeneity of variance).   I think you need to distinguish between heteroscedasticity and heterogeneity of variance, and I just saw something on Stack Exchange where someone seemed to agree with me.

Anyone else want to comment on that? 

Dear Peter Samuels, Ahamad Daryanto has written syntax for supporting BP and Koenker tests in SPSS. Simply download the sps file by the same and install as an extension in SPSS. Make sure you are using latest versions at least 18. If all goes well you will see the heteroskedasticity icon in the list of regression functions.

Francesca -

As I noted above, a test does not get you anywhere.  However, you can estimate the coefficient of heteroscedasticity, or choose a reasonable default value, as noted in the document and Excel file whose URLs are given below, and use it in the regression weight.  (I noted above that the predicted y can be used as a  "size" measure in this.  That should be the best possibility.)

Estimating the Coefficient of Heteroscedasticity:

https://www.researchgate.net/publication/333642828_Estimating_the_Coefficient_of_Heteroscedasticity

https://www.researchgate.net/publication/333659087_Tool_for_estimating_coefficient_of_heteroscedasticityxlsx

which are part of

 https://www.researchgate.net/project/OLS-Regression-Should-Not-Be-a-Default-for-WLS-Regression

Cheers - Jim

The unfortunate part of most of these tests for heteroskedasticity is that they do not show you the source of the problem, e.g., which variable is at fault.

You can look at heteroscedasticity with respect to one independent variable at a time, but because essential heteroscedasticity is due to population member size, the best measure of size is predicted y, not an individual predictor, unless that is all you have (ratio estimators).

Heteroscedasticity is not a problem to be fixed, nor the sole 'fault' of any one independent variable in multiple regression. It is a natural consequence of the various 'sizes' of members of a population (essential heteroscedasticity). Other factors which may need to be addressed may artificially increase, or decrease the level of heteroscedasticity. Those may possibly be 'fixed.'

Well, it depends on the problem at hand... The independent variable can be at fault for heteroskedasticity. An example is when one is modelling family expenditures on household incomes. You would find that families that have low incomes would spend relatively less on vacations, resulting in low variations in expenditures across such families. However, those families with larger income levels, the amount of discretionary income would be higher. Thus, the expected level of expenditure on vacations would be higher - leading to greater variability among such families and therefore resulting in heteroskedasticity... Alternatively, errors may as well increase as values carried by independent variables are extreme in both directions, e.g. a commonly used example are attitudes that range from extremely negative to extremely positive. As most textbooks say: this would produce an hourglass shape - resulting in heteroskedasticity...

It depends on the complete set of predictors providing the predicted y values.

This explains it:

https://www.researchgate.net/publication/320853387_Essential_Heteroscedasticity

This is modified by this:

https://www.researchgate.net/publication/324706010_Nonessential_Heteroscedasticity

If predictors are included that do not belong, or are left out which are needed; if there is overfitting or underfitting; this can make the predicted values less useful as a size measure, so this can impact heteroscedasticity. But the best size measure should be the predicted y values, if we do not have these other problems.

Looking at individual predictors may not help if they are working together in groups.

It is desirable to have an unbiased model with as little variance as possible, and a coefficient of heteroscedasticity between 0.5 and 1, as defined in the following:

https://www.researchgate.net/publication/333642828_Estimating_the_Coefficient_of_Heteroscedasticity

https://www.researchgate.net/publication/333659087_Tool_for_estimating_coefficient_of_heteroscedasticityxlsx

and Brewer, K.R.W.(2002), Combined survey sampling inference: Weighing Basu's elephants, Arnold: London and Oxford University Press, especially pages 87, 111, 130, 137, 142, and 203.

I'm not certain of the application here. It says "I have two data sample. For each sample I have more variables. I should make some regressions." Well, the tests listed after that make me wonder if the point is to look at some linear regressions and consider heteroscedasticity shown in the estimated residuals in each case, OR as in Levene's test, really be more interested in the two data sets mentioned, and the variances of those two datasets. The goals aren't clear.

Note that the papers and spreadsheet I provided are designed for the first of the two situations I just described.