Bunny, lautloscrew.com
2009-Sep-18 10:17 UTC
[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... 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 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 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
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-18bp.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 guidehttp://www.R-project.org/posting-guide.html> and provide commented, minimal, self-contained, reproducible code.