# Regression modelling with unequal variance

I would like to fit a linear model (lm) where the residuals variance is clearly dependent on the explanatory variable.

The way I know to do this is by using glm with the Gamma family to model the variance, and then put its inverse into the weights in the lm function (example: http://nitro.biosci.arizona.edu/r/chapter31.pdf)

I was wondering:

• Is this the only technique?
• What other approaches are relevant?
• What R packages/functions relevant to this type of modelling? (other then glm, lm)
• Where do they use glm() then lm() in the chapter you link to. Seems to me the glm() is all that is required and used there, but I may have missed something. You can try generalised least squares (gls() in nlme) which allows weights to be estimated to control for the type of heteroscedasticity you mention; see ?varFunc and follow the links from there. IIRC varFixed() will do what you want. Aug 14, 2012 at 21:35
• In 'proc mixed', 'subject=option' produces a block-diagonal structure in the variance-covariance matrix of the residuals. Have you thus considered a general linear mixed model to alter the homoscedasticity hypothesis? Aug 15, 2012 at 12:16
• Thanks Gavin, I've looked a bit at these functions. Two questions: 1) Do you recommend any tutorials? (I suspect that the MASS book my be a good start, but was wondering if you had a thought on it). 2) Since the model that I am fitting is a simple OLS, how different will the estimation be for when using the gls function? (If I rememeber correctly - not much, since it should be working on some iterative first degree approximation, but I am not at all sure about this). Ocram - thanks, but I do not use SAS. Aug 15, 2012 at 15:47
• Here in Section 2 it is explained how to do this in STATA for quasipoisson regression: stata.com/meeting/fnasug08/gutierrez.pdf. If somebody could suggest a way to recode this in R, I would be very grateful. Jul 25, 2013 at 12:10

Pills against the "megaphone effect" include (among others):

1. Use log or square root transform $Y$. This is not exact but sometimes it tames the widening.
2. Use weighted least square regression. In this approach, each observation is given its own variance factor. This answer shows how to use WLSR in R (for instance if the variance of the residuals is proportional to the means, you can provide as weights the inverse of the fitted value in the unweighted model).
3. Use robust regression. The funciton rlm() in the MASS package of R does M-estimation, which is supposed to be robust to inequality of variances.

July 2017 edit: It seems that generalized least squares, as suggested in the answer of Greg Snow, is one of the best options.

• I built on this answer for a Stack Overflow question. Jan 14, 2017 at 6:23
• Could be worth pointing out the generalized least squares option as well, using gls with the weights option set to varFixed() - to me this would seem like one of the more elegant options... Jul 13, 2017 at 9:01
• @TomWenseleers I agree. Notice that this is the answer of Greg Snow. Jul 13, 2017 at 16:28

The gls function in the nlme package for R can estimate the regression and the relationship with the variance at the same time. See the weights argument and the 2nd example on the help page.

With the gamlss package you can model the error distribution of the response as a linear, a non-linear, or a smooth function of the explanatory variables. This seems to be a quite powerful approach (I learned a lot about all the possibilities that might arise during the model selection process) and everything is nicely explained in several publications (including books) that is referenced at the link above.