27

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)
kjetil b halvorsen
  • 77,844
Tal Galili
  • 21,541
  • 5
    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. – Gavin Simpson Aug 14 '12 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? – ocram Aug 15 '12 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.

     – Tal Galili Aug 15 '12 at 15:47
  • Here in Section 2 it is explained how to do this in STATA for quasipoisson regression: http://www.stata.com/meeting/fnasug08/gutierrez.pdf. If somebody could suggest a way to recode this in R, I would be very grateful. – a11msp Jul 25 '13 at 12:10

3 Answers3

20

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.

gui11aume
  • 14,703
14

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.

Greg Snow
  • 51,722
11

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.

kjetil b halvorsen
  • 77,844
ils
  • 561