Difficulty with specification of random (linear) effects in nonlinear models.

[EmmanuelCharpentier]

I have trouble unterstanding how to add a nonlinear effect to an (already working) mixed linear model.

Working example : this compiles and gives sensible results :

brm(log(Duree)~Classe+Arcade+ReTt+NoTf+NoWO+(1|Etudiant), data=brmsData)

First attempt :

brm(log(Duree)~Classe+Arcade+ReTt+NoTf+NoWO+(1|Etudiant)+(-log(km+(1-km)*exp(-betaExp*Exp))),
data=brmsData, nonlinear=list(km~1,betaExp~1),prior=c(set_prior("cauchy(0,3)",nlpar="km"),
set_prior("cauchy(0,3)",nlpar="betaExp")))
Erreur : Random effects in non-linear models should be specified in the 'nonlinear' argument.

Well... second attempt :

brm(log(Duree)~Classe+Arcade+ReTt+NoTf+NoWO+(1|Etudiant)+(-log(km+(1-km)*exp(-betaExp*Exp))),
data=brmsData, nonlinear=list(km~1,betaExp~1,1~Etudiant),
prior=c(set_prior("cauchy(0,3)",nlpar="km"),set_prior("cauchy(0,3)",nlpar="betaExp"),
set_prior("cauchy(0,10)",nlpar="1")))
Erreur : Random effects in non-linear models should be specified in the 'nonlinear' argument.

Ah.

Added an "Intercept" column of ones to the dataset. Then :

brm(log(Duree)~0+Intercept+Classe+Arcade+ReTt+NoTf+NoWO+log(km+(1-km)*exp(-betaExp*Exp)),
data=brmsData, nonlinear=list(km~1,betaExp~1,Intercept~Etudiant),
prior=c(set_prior("cauchy(0,3)",nlpar="km"),set_prior("cauchy(0,3)",nlpar="betaExp"),
set_prior("cauchy(5,3)", nlpar="Intercept")))
SYNTAX ERROR, MESSAGE(S) FROM PARSER:


ERROR at line 35

 33:      for (n in 1:N) { 
 34:        // compute non-linear predictor 
 35:        eta[n] <- (0 + eta_Intercept[n] + C[n, 1] + C[n, 2] + C[n, 3] + C[n, 4] + C[n, 5] + log(eta_km[n] +); 
                                                                                                               ^
 36:      } 

PARSER EXPECTED: <expression>
Error in stanc(model_code = paste(program, collapse = "\n"), model_name = model_cppname,  : 
  failed to parse Stan model 'file1f1447d26961' due to the above error.

And I'm stuck... The error seems unrelated to the previous errors, but I'm not sure.

Any ideas ?

[paul-buerkner]

The last error you encountered appears to be a bug in the formula parsing. I will try to fix this later on today.

The moment you use the nonlinear argument, the formula argument will be interpreted completely different than it usually is in the sense that it now describes a non-linear model. That is you have to explicitely write down every parameter in the predictor part of the model. All variables on the right hand side of ~ that are not labled as parameters in the nonlinear argument are considered as data. There is no need to put an intercept column into your data, because you defined the Intercept as a parameter (not as data).

My suggestion would be to write the model as follows:

formula = log(Duree) ~ b0 + b1*Classe + b2*Arcade + b3*ReTt + b4*NoTf + b5*NoWO + log(km+(1-km)*exp(-betaExp*Exp))
nonlinear = list(b0 ~ (1|Etudiant), b1 ~ 1, b2 ~ 1, b3 ~ 1, b4 ~ 1,  b5 ~ 1, km ~ 1, betaExp ~ 1)

Don't forget to set priors on all non-linear parameters. However, this won't work until I fixed the formula parsing. I will write again as soon as I have fixed it

[EmmanuelCharpentier]

Le dimanche 21 février 2016 à 06:34 -0800, Paul-Christian Bürkner a
écrit :

The last error you encountered appears to be a bug in the formula
parsing. I will try to fix this later on today.
The moment you use the nonlinear argument, the formula argument will be
interpreted completely different than it usually is in the sense that
it now describes a non-linear model. That is you have to explicitely
write down every parameter in the predictor part of the model.
This is not readily apparent at the reading of the vignette or the
documentation of the "brm" function...
 All variables on the right hand side of ~ that are not labled as
parameters in the nonlinear argument are considered as data. There is
no need to put an intercept column into your data, because you defined
the Intercept as a parameter (not as data).

My suggestion would be to write the model as follows:

formula> > => log(> Duree> ) > ~> > b0> > +> > b1> *> Classe> > +> > b2> *> Arcade> > +> > b3> *> ReTt> > +> > b4> *> NoTf> > +> > b5> *> NoWO> > +> log(> km> +> (> 1> -> km> )> *<
/span>exp(> -> betaExp> *> Exp> ))

nonlinear> > => > list> (> b0> > ~> (> 1> |> Etudiant> ), > b1> > ~> > 1> , > b2> > ~> > 1> , > b3> > ~> > 1> , > b4> > ~> > 1> , > b5> > ~> > 1> , > km> > ~> > 1> , > betaExp> > ~> > 1> )
Real close to writing the code in Stan, IIUYC. Almost another way to
specify a model (a. k. a. a new modeling (mini-)language)...
Will this code manage to recognize that Classe is a 3-level factor ?
And recode it as a base level confounded with the intercept and two
boolean variables (hence a two-dimensional b1) ?
Should I declare the standard deviations ("Etudiant" random effect and
residual) ? If so, how ?

Don't forget to set priors on all non-linear parameters. However, this won't work until I fixed the formula parsing. I will write again as soon as I have fixed it
Thanks a lot ! This level of responsiveness is unusual (to say the
least...).

Sincerely,

Emmanuel Charpentier

—
Reply to this email directly or view it on GitHub.

[paul-buerkner]

You are right, the documentation needs some improvements here.

All formulas in nonlinear are treated as ordinary formulas so no need to declare the standard deviation of Etudiant. Also the residual SD does not need to be defined as it is not part of any formula anyway but a "property" of the gaussian family.

The way brms handles non-linear formulas is very similar to what the nlme package does. A few minutes ago, Alex Forrence posted some examples in #37, maybe this can help you getting more familiar with the syntax.

I wouldn't call it a mini-language, but rather an explicit way of writing a model. You cannot treat parts of the formula as usual and the other as "non-linear". If you look at the Stan code via stancode(model) after fitting your model (or via directly calling make_stancode), you will see that there is more going on other than writing the predictor part of the model.

Currently, a 3-level factor will not be allowed in a non-linear formula (you will recieve an error message before the model is compiled), but it's allowed to put it in nonlinear.

In fact, as I think about it, you can greatly simplify your model by writing:

formula = log(Duree) ~ b + log(km+(1-km)*exp(-betaExp*Exp))
nonlinear = list(b ~ Classe+Arcade+ReTt+NoTf+NoWO+(1|Etudiant), km ~ 1, betaExp ~ 1)

Your linear predictor will be plugged in for b in the non-linear formula and the actual non-linear part log(km+(1-km)*exp(-betaExp*Exp)) is added.

[paul-buerkner]

The fix for the formula parsing issue is now on github. You can install the latest dev version of brms via

library(devtools)
install_github("paul-buerkner/brms")

Hopefully, this will allow your model to be fitted.

[EmmanuelCharpentier]

[ Sorry for tjhe (non-maskable) interruption : RealLife(TM) happens...
]
Yes, that works : I have been able to recompile my test case and to get
sensible results (close, but not identical, to a hand-made rstan fit).
I'm still a bit confused about the expression of the models and the
priors. But that's a documentation problem...

Again, thank you very much !

Emmanuel Charpentier
Le dimanche 21 février 2016 à 07:28 -0800, Paul-Christian Bürkner a
écrit :

The fix for the formula parsing issue is now on github. You can
install the latest dev version of brms via
library(devtools)
install_github("paul-buerkner/brms")
Hopefully, this will allow your model to be fitted.
—
Reply to this email directly or view it on GitHub.

[paul-buerkner]

Ok great! The documentation will be improved with the next version. Also, I am currently writing a vignette about non-linear models (see #37).

Thanks for opening this issue and helping me to improve brms!