[R] some irritation with heteroskedasticity testing

[John Fox]

Dear matt,

> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On
> Behalf Of Bunny, lautloscrew.com
> Sent: September-18-09 6:18 AM
> To: r-help at r-project.org
> Subject: [R] some irritation with heteroskedasticity testing
> 
> Dear all,
> 
> Trying to test for heteroskedasticity I tried several test from the
> car package respectively lmtest. Now that they produce rather
> different results i am somewhat clueless how to deal with it.
> Here is what I did:
> 
> 1.  I plotted fitted.values vs residuals and somewhat intuitively
> believe, it isn't really increasing...

It's often difficult to judge whether the spread of the residuals changes
with the fitted values; you might alternatively try to visualize the fit
with the spread.level.plot() function in the car package.

> 
> 
> 2. further I ran the following tests
> bptest (studentized and non-studentized), gqtest, ncv.test with the
> following results:
> 
> ncv:
> Non-constant Variance Score Test
> Variance formula: ~ fitted.values
> Chisquare = 13.87429    Df = 1     p = 0.00194580
> 
> Goldfeld-Quandt test
> data:  reg
> GQ = 1.7092, df1 = 327, df2 = 327, p-value = 7.93e-07
> 

You don't show the function call to gqtest(). If you specified
order.by=fitted(reg), then the two tests are similar, except that gqtest()
by default uses a one-sided alternative hypothesis. If you used the default
for order.by, then the test is probably nonsense. Both tests suggest
heteroscedasticity.

> 
> studentized Breusch-Pagan test
> data:  reg
> BP = 15.8291, df = 23, p-value = 0.92
> 
> 
> Breusch-Pagan test
> data:  reg
> BP = 377.5604, df = 23, p-value < 2.8e-18

bp.test() performs the same score test as ncv.test(), except that the
default alternative hypothesis is different -- in bp.test() that the error
variance is a function of a linear combination of the regressors and in
ncv.test() that the error variance is a function of the fitted values (i.e.,
a *particular* linear combination of regressors). Testing against the fitted
values with 1 df will have greater power if this is the real pattern of
heteroscedasticity. That the chisquare is much greater for the general BP
test suggests that the pattern is more complex. The studentized BP test is
robust against non-normal errors. That you would get such a dramatically
different result is surprising and would lead me to guess that the residual
distribution is very heavy-tailed, perhaps with some bad outliers. Without
the data, it's not possible to know.

Regards,
 John

> 
> 
> bptest and gq.test sport pretty straight forward examples saying the
> H0 = homoskedasticity. The ncv.test clarifies the same in its
> description. Thus the studentized bptest appears to be the only one to
> support my first intuition from the graphical solution. I am really
> disturbed what to do know and how to interpret my results... Can
> someone lead the way ?
> 
> thx in advance
> 
> matt
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.