[R-sig-ME] Gamma transformation

[Ben Bolker]

Iker Vaquero Alba <karraspito at ...> writes:

> 
> 
>    Dear list: 
> 
>    After implementing a model, I tested it for homocedasticity, 
> and I was not actually very happy, as the
> variance seems to increase for higher values, according to 
> the fitted values vs. residuals plot
> (attached - "Unmodified"). I tried to log-transform the
>  response variable, but the result was similar.
> After looking for possible transformations of the response variable, 
> I found that "gamma(x)" worked
> really well, and stabilized the variance of the model for all 
> the range when plotting the fitted values vs
> residuals (attached - "Gamma"). 

  I won't say it's impossible, and maybe someone will chime in,
but transforming by the gamma distribution seems pretty weird to me.
Power transformations are much more common.  If gamma-transforming
does work, though, it is consistent with log transformation not
being satisfactory, because the gamma is *accelerating* for x>1,
while the log is monotonically *decelerating* -- it suggests that
if you are going to use a power transformation you should use
x^b with b > 1 ...

> 
>    I am happy with that, but it's just that I would like to make
>  sure wether this is a correct step, and if it's
> possible, know what is the reason why this happens when 
> applying that transformation. I am not an expert
> but I would really like to try to learn something more about 
> this. I have tried to find some information in
> the web, but I'm not very happy with it. Obviously, only if 
> you have the time and the willingness to explain
> it, or help me find some information about it.
>  I appreciate it very much in any case.

  A few things to think about:

1. boxcox() (in MASS) and powerTransform (in car) don't work with
lme-type objects, but they would work on a simpler lm-approximation
of your model (either leaving out the random effects or treating
them as fixed), and should give you at least a rough idea of the
appropriate transformation.

2. you can also try using a varStruct (e.g.
weights=varPower() ; by default this will fit
a model with the residuals proportional to 
a power of the fitted value)

  Ben Bolker

>