hierarchical gamma model in lme4 - R-SIG-mixed-models

Tue, May 17, 2011 4:50 PM #

[forwarding to r-sig-mixed-models list ...]

  As of today, Gamma models are (still) not feasible in lme4 -- they are
somewhat more numerically challenging than the other families, so Doug
Bates is having to do some re-engineering.
  There is a *possibility* that I can get Gamma fitting to work in the
'alpha'/bleeding-edge development version of glmmADMB, but it will
definitely be bleeding-edge ... if you are interested in trying that,
please contact me off-list.

   In the meantime, my standing advice is to try a LMM on the
log-transformed data (zero values in the response are problematic, but
they would be problematic in a Gamma GLMM in any case if the shape
parameter is ever < 1 ...)

  Ben Bolker

On 11-05-17 07:32 PM, Benjamin Caldwell wrote:

Addendum: I tried a gamma fit in glmmPQL and got the same errors.
*Ben Caldwell*

PhD Candidate
University of California, Berkeley




On Tue, May 17, 2011 at 3:51 PM, Benjamin Caldwell
<btcaldwell at berkeley.edu <mailto:btcaldwell at berkeley.edu>> wrote:

    Hello
    After seeing this
    (https://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q1/005213.html) email
    I thought I would check the issue with a gamma family in lme4 hadn't
    been fixed; can I fit a hierarchical gamma model in lme4 at this
    time? There doesn't seem to be another package capable of it at this
    time.

    My thought process: 
    1. took a look at the response variable and some subsets to see what
    it looked like, ("bppfcl" and "transformed response var"), attached
    2. took a look at a gamma and gaussian fit to the response variable.
    3. ran hierarchical gaussian model in nlme to look at residuals
    (more familiar with graphs from that package) ("qqnorm" and "residuals")

    Given the residual output for the gaussian model it looks like I
    could remove the values at the end of the distribution and get a
    decent fit. I'd still like to try a gamma model though, if that's
    possible. Is it possible in lme4 or another package I don't know about?

    ---This is the code I'm running---

    rws30.BL$site <- factor(rws30.BL$site)
    rws30.BL$transect <- interaction(rws30.BL$site, rws30.BL$transect,
    drop = TRUE)
    rws30.BL$plot <- interaction(rws30.BL$site, rws30.BL$transect,
    rws30.BL$plot, drop = TRUE)
    hist(rws30.BL$post.f.crwn.length)
    rws30.BL$gpost.f.crwn.length

    library("nlme")
    burnedmodel1.3<-lme(post.f.crwn.length~lg.shigo.av+dbh+leaf.area+bark.thick.bh
    <http://bark.thick.bh>+ht.any+ht.alive,
    random=(~1|site/transect/plot),na.action=na.omit, data=rws30.BL)
    Error: no valid set of coefficients has been found: please supply
    starting values
    In addition: Warning message:
    In log(ifelse(y == 0, 1, y/mu)) : NaNs produced

    --- I thought the problem might be a convergence error, and so tried
    a reduced model ----
     glmer(gpost.f.crwn.length~dbh+leaf.area+(1|site/transect/plot),
    family=Gamma, na.action=na.omit, data=rws30.BL)
    Error in mer_finalize(ans) : 
      mu[i] must be positive: mu = -0.00780625, i = 3

    Any clarity I could get would be much appreciated.

    Best

    *Ben Caldwell*

    PhD Candidate
    University of California, Berkeley