30
$\begingroup$

I have a plot of residual values of a linear model in function of the fitted values where the heteroscedasticity is very clear. However I'm not sure how I should proceed now because as far as I understand this heteroscedasticity makes my linear model invalid. (Is that right?)

  1. Use robust linear fitting using the rlm() function of the MASS package because it's apparently robust to heteroscedasticity.

  2. As the standard errors of my coefficients are wrong because of the heteroscedasticity, I can just adjust the standard errors to be robust to the heteroscedasticity? Using the method posted on Stack Overflow here: Regression with Heteroskedasticity Corrected Standard Errors

Which would be the best method to use to deal with my problem? If I use solution 2 is my predicting capability of my model completely useless?

The Breusch-Pagan test confirmed that the variance is not constant.

My residuals in function of the fitted values looks like this:

https://i.gyazo.com/9407a829a168492b31dfa3d1dd33a21d.png

(larger version)

Improve this question
$\endgroup$
3
  • $\begingroup$ Do you mean 'stackoverflow' rather than 'stackexchange'? (you're still on stackexchange here.) If it was SO, it's generally better to migrate the question rather than posting a second copy (the help asks not to post the same Q multiple times but choose one best place). $\endgroup$
    – Glen_b
    Commented Apr 19, 2015 at 2:04
  • $\begingroup$ The variation in spread isn't so much that the impact will be severe (that is, while it will bias your standard errors and so impact inference, it's probably not going to make a huge difference). I'd be inclined to consider whether spread was related to mean, and perhaps look at a GLM or possibly transformation (it sure looks related to fitted). What's the y-variable? $\endgroup$
    – Glen_b
    Commented Apr 19, 2015 at 2:07
  • 2
    $\begingroup$ Another possibility is to model the heteroscedasticity, e.g., using gls and one of the variance structures from package nlme. $\endgroup$
    – Roland
    Commented Apr 20, 2015 at 7:03

4 Answers 4

Reset to default
35
$\begingroup$

It's a good question, but I think it's the wrong question. Your figure makes it clear that you have a more fundamental problem than heteroscedasticity, i.e. your model has a nonlinearity that you haven't accounted for. Many of the potential problems that a model can have (nonlinearity, interactions, outliers, heteroscedasticity, non-Normality) can masquerade as each other. I don't think there's a hard and fast rule, but in general I would suggest dealing with problems in the order

outliers > nonlinearity > heteroscedasticity > non-normality

(e.g., don't worry about nonlinearity before checking whether there are weird observations that are skewing the fit; don't worry about normality before you worry about heteroscedasticity).

In this particular case, I would fit a quadratic model y ~ poly(x,2) (or poly(x,2,raw=TRUE) or y ~ x + I(x^2) and see if it makes the problem go away.

Improve this answer
$\endgroup$
4
  • $\begingroup$ The plot is small & the axes aren't labeled. I don't know if it's a residuals vs fitted plot. I assumed the OP included a squared term, eg. If not, you're clearly right. $\endgroup$
    – gung - Reinstate Monica
    Commented Jan 3, 2016 at 23:15
  • 1
    $\begingroup$ in my browser I can see that the y-axis range goes from -4 to 3, which seems to suggest a residuals vs. fitted plot / rule out a scale-location plot ... $\endgroup$
    – Ben Bolker
    Commented Jan 3, 2016 at 23:18
  • 2
    $\begingroup$ Hi Ben, love what you do. Can you expand on the idea that "outliers" are the biggest issue? Do you include single high-leverage points as "outliers" even if they have a small residual? I deal with extreme value observations all the time in my line of work (environmental statistics), and I find that some people (the EPA in particular) tend to blow outliers way out of proportion (pardon any unintentional pun) and are way to eager to exclude them. I tend to adopt a tolerant attitude to outliers if I can't find good evidence that they are clearly the result of data (collection, entry) error. $\endgroup$
    – Dalton Hance
    Commented Jan 4, 2016 at 20:48
  • 1
    $\begingroup$ @DaltonHance: we're probably pretty much on the same page. My point is just that if you have outliers (by whatever definition) and they're not taken into account by whatever statistical model/approach you're using (mixture models, robust statistics, fat-tailed distributions, etc.), then it will tend to screw up all the rest of your diagnostics -- it will make residuals look nonlinear/heteroscedastic/non-Normal. I certainly agree that you shouldn't just thoughtlessly/reflexively throw them out. $\endgroup$
    – Ben Bolker
    Commented Jan 4, 2016 at 21:00
13
$\begingroup$

I list a number of methods of dealing with heteroscedasticity (with R examples) here: Alternatives to one-way ANOVA for heteroskedastic data. Many of those recommendations would be less ideal because you have a single continuous variable, rather than a multi-level categorical variable, but it might be nice to read through as an overview anyway.

For your situation, weighted least squares (perhaps combined with robust regression if you suspect there may be some outliers) would be a reasonable choice. Using the Huber-White sandwich errors would also be good.

Here are some answers to your specific questions:

  1. Robust regression is a viable option, but would be better if paired with weights in my opinion. If you aren't worried that the heteroscedasticity is due to outliers, you could just use regular linear regression with weights. Be aware that the variance can be very sensitive to outliers, and your results can be sensitive to inappropriate weights, so what might be more important than using robust regression for the final model would be using a robust measure of dispersion to estimate the weights. In the linked thread, I use 1/IQR, for example.
  2. The standard errors are wrong because of the heteroscedasticity. You can adjust the standard errors with the Huber-White sandwich estimator. That is what @GavinSimpson is doing in the linked SO thread.

The heteroscedasticity does not make your linear model totally invalid. It primarily affects the standard errors. If you don't have outliers, least squares methods should remain unbiased. Therefore the predictive accuracy of point predictions should be unaffected. The coverage of interval predictions would be affected if you didn't model the variance as a function of $X$ and use that to adjust the width of your prediction intervals conditional on $X$.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ using robust regression from the lmrob package would automatically infer some weights, why not use those instead in #1? $\endgroup$
    – tool.ish
    Commented Oct 17, 2016 at 12:55
1
$\begingroup$

Load the sandwich package and compute the var-cov matrix of your regression with var_cov<-vcovHC(regression_result, type = "HC4") (read the manual of sandwich). Now with the lmtest package use the coeftest function:

coeftest(regression_result, df = Inf, var_cov)
Improve this answer
$\endgroup$
-1
$\begingroup$

How does the distribution of your data looks like? Does it look like a bell curve at all? From the subject matter, can it be normally distributed at all? Duration of a phone call may not be negative, for example. So in that specific case of calls a gamma distribution describes it well. And with gamma you can use generalized linear model (glm in R)

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.