Initialization problems using brms with Gamma distribution and identity link

Interfaces brms
1

Hi

Short summary of the problem :

I try to fit tree growth height as a function of tree growth circumference with a non-linear hierarchical model in brms.

acc_h ~ exp(a)(1-exp(-exp(b)*acc_c) with :
acc_h = tree growth height
acc_c = tree growth circumference
a and b the parameter to b estimate

I would like to use a Gamma distribution (growth can’t be negative) and an identity link in order to not change the form of the relation between height growth and circumference growth.

The expectation of mu should always be positive because acc_c>0.

The sampler has trouble finding initial values for the scale parameter and throws the following error multiple times:

Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[39] is inf, but must be finite! (in ‘model1eff7e31495d_93ef33a66d298bf1ece66700bd599138’ at line 85)

nlform_C <- bf(acc_h ~ exp(alpha)*(1-exp(-exp(beta)*acc_c)),
               alpha ~ 1+ (1|flib_placette/fperiode)   , beta ~  1+(1|flib_placette/fperiode)  , nl = TRUE)

nlprior_C <- c(prior(normal(4, 2), nlpar = "alpha"),
               prior(normal(0, 2), nlpar = "beta"))

fit_C <- brm(formula = nlform_C, data = arbre_periode3_brassy, prior = nlprior_C, family = Gamma(link="identity"),
             warmup=1000,iter=1200,chains = 1)

Is there a way to fix this problem ?

thanx

Please also provide the following information in addition to your question:

  • Operating System: Ubuntu 20
  • brms Version: 2.16
2

Hi @matthias56
You can set init values using the inits option as described here, maybe that helps.

3

Hi @scholz

Thanks. I had already try with fixed init (ex inits = 0.1 or init_r =0.1) but it doesn’t work. Maybe I should set init only for some parameter.

Du to my formula acc_h ~ exp(a)(1-exp(-exp(b)*acc_c) :

  • acc_h = 0 for acc_c=0
  • acc_h <0 for acc_c <0

So i modified my formula by multiply it by step(acc_c) and add a small constant :

nlform_C <- bf(acc_h ~ 0.0001+step(acc_c)*exp(alpha)*(1-exp(-exp(beta)*acc_c)),
               alpha ~ 1+ (1|flib_placette/fperiode)   , beta ~  1+(1|flib_placette/fperiode)  , nl = TRUE)

With this new form acc_h is always strictely positive and there is no more initialization issue. It work’s even without the step(acc_c). I don’t understand why because without it acc_h is <0 for acc_c < 0.

Is this solution (add 0.0001+step(acc_c)* to my formula) correct ?
Is there a more elegant solution ? like specify i want to model acc_h only for acc_c > 0 ?

4

I don’t think I can really help you with that as I have no experience with that kind of advanced formulas.