[R-sig-eco] Fixing heteroscedasticity in mixed-effects model?

[Malin Pinsky]

Hi all,

I'm having problems fitting a mixed-effects model for an ecological
meta-analysis, and I'm curious if anyone has advice. In particular,
it's pretty clear that the variance in the residuals increases with
the predicted mean, but my normal fixes don't seem to be working. The
model is:

mod1 <- lmer(logCd ~ logRe + Hab + logRe:Hab + (logRe|Study), data=temp)

where Cd is a drag coefficient (>0 before log-transformation), Re is a
physical quantity called a Reynolds number (also >0 before
transformation), Hab is a categorical variable for habitat, and Study
is a categorical variable for the study the data came from. I know
from fluid dynamics theory that logCd and logRe can be linearly
related, but I expect that the slope and intercept vary between
habitat types and between studies.

A plot of

plot(fitted(mod1), resid(mod1))

looks like a nice cone with the highest spread on the right-hand side
of the graph. My initial thought was to use a different transformation
for Cd, but I couldn't find a power that made any difference for the
problem of increasing variance. I then tried a Gamma error, thinking
that this corrects for problems of increasing variance. I had to
change the transformation of Cd slightly, now using log(Cd+1) so that
there were no negative values:

mod1 <- lmer(logCd ~ logRe + Hab + logRe:Hab + (logRe|Study),
data=temp, family=Gamma)

but the model won't converge. In particular, I get "Warning in
mer_finalize(ans) : singular convergence (7)." My dataset has 686
observations across 4 habitat levels and 24 studies.

Is there something obvious I'm missing, or are the other avenues to
try? Any advice would be welcome. I'm using lmer() from the lme4
package in R 2.14.1 on Mac OS X 10.6.8.

And, if this belongs on the R-sig-ME list, let me know.

Thanks,
Malin

--
Malin Pinsky
David H. Smith Conservation Research Fellow
Department of Ecology and Evolutionary Biology
Princeton University
malin.pinsky at gmail.com