[R-sig-ME] hierarchical gamma model in lme4

[Ben Bolker]

  [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
> 
> 
>