[R] using hglm to fit a gamma GLMM with nested random effects?

[Benjamin Caldwell]

Apologies for continuing to ask about this but . .  in my quest to fit a
gamma GLMM model to my data (see partial copy of thread below), I'm
exploring using hglm today. The question of the day has to do with the
errors I'm currently getting from the hglm package. Can hglm handle a model
with nested random effects? I don't see an example of one of those in the
package documentation. If it can, can anyone tell me what these errors are
trying to tell me?

If no, I promise I'll let this rest and just take B. Boker's advice to go
with a nice safe modified log transform of the data.

Best


test.gamma<-hglm(fixed=post.f.crwn.length~lg.shigo.av+dbh+leaf.area+
bark.thick.bh+ht.any, random=~1|site/transect/plot, family=Gamma(link=log),
data=rws30.BL)
Error in `contrasts<-`(`*tmp*`, value = "contr.treatment") :
  contrasts can be applied only to factors with 2 or more levels
In addition: Warning messages:
1: In Ops.factor(site, transect) : / not meaningful for factors
2: In Ops.factor(site/transect, plot) : / not meaningful for factors
> test.gamma<-hglm(fixed=post.f.crwn.length~lg.shigo.av+dbh+leaf.area+
bark.thick.bh+ht.any, random=~1|site, family=Gamma(link=log), data=rws30.BL)
Error in hglm.default(X = X, y = Y, Z = z, family = family, rand.family =
rand.family,  :
  Length of X and Z differ.
>




* Dennis Murphy <djmuser at gmail.com> Tue, May 17, 2011 at 6:18 PM
To: Benjamin Caldwell <btcaldwell at berkeley.edu>
 Hi:

Someone else (Wayne Zhang, CNA) asked a similar question re
hierarchical Gamma models on R-help today and responded to suggestions
as follows:

"Hglm does the work! Thanks!

Also, I find that the developing version of lme4, called lme4a, has
the capability to fit Gamma models. And both lme4a and hglm produce
results consistent with the published ones.

Problems solved!"

Perhaps you might find success following his lead?

Dennis
Ben Bolker <bbolker at gmail.com>
* *Tue, May 17, 2011 at 4:50

PM*
To: Benjamin Caldwell <btcaldwell at berkeley.edu>
Cc: "r-sig-mixed-models at r-project.org" <r-sig-mixed-models at r-project.org>
  [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<https://mail.google.com/mail/?ui=2&ik=938097cb0f&view=pt&search=inbox&msg=130005e6516c0ad7&dsqt=1>
 <mailto:btcaldwell at berkeley.edu<https://mail.google.com/mail/?ui=2&ik=938097cb0f&view=pt&search=inbox&msg=130005e6516c0ad7&dsqt=1>>>
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
>
>
>


*Ben Caldwell*

PhD Candidate
University of California, Berkeley
-------------- next part --------------
A non-text attachment was scrubbed...
Name: residuals.png
Type: image/png
Size: 9189 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20110518/110df56d/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: transformed response var, and subsets.png
Type: image/png
Size: 10643 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20110518/110df56d/attachment-0001.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Distribution response var.png
Type: image/png
Size: 8376 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20110518/110df56d/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: qqnorm plot of gaussian hierarchical fit.png
Type: image/png
Size: 7515 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20110518/110df56d/attachment-0003.png>