Best way to deal with heteroscedasticity? 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?)


*

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

*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:  
 
(larger version)

 A: 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.
A: 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:  


*

*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.  

*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$.  
A: 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)

A: 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)
