Double exponential using brms

Modeling
1

Dear friends - based on published data I have simulated double exponential decay and tried to recover the known parameters, I get divergensies and bad pairs plots. Did I do something very wrong or is it not possible this way? Below are the simulation and model. All best wishes Troels

library(brms)
Time <- seq(1,60,by=1)
A <- 275; B <- 140; g1 <- 0.1105; g2 <- .0161
Y <- Aexp(-g1Time) + Bexp(-g2Time)
Y <- Y + rnorm(60,0,0.5*sqrt(Y))
plot(Time,Y,pch=19,cex=.5)
plot(Time,log(Y),pch=19,cex=.5)
dat1 <- data.frame(Time,Y)

prior1 <- c(prior(normal(275,2),nlpar=“A”),
prior(normal(140,2),nlpar=“B”),
prior(normal(.1,.002),nlpar=“g1”),
prior(normal(01,.005),nlpar=“g2”)
)
fit1 <- brm(bf(Y~Aexp(-g1Time)+Bexp(-g2Time),A+B+g1+g2 ~1,nl=TRUE),
data=dat1,cores=4,prior = prior1,chains=4,control = list(adapt_delta = 0.99,max_treedepth = 20),iter = 8000)

2

You could try brigning A and B to a similar scale as g1 and g2. So for instance, scale time by 100:

Y~A * exp(-g1 * Time / 100)+B * exp(-g2 * Time / 100)

Of course, you should then adjust the priors for g1 and g2 as well.

1 Like
3

This worked perfectly - thanks a lot Paul!!

All best wishes

Troels

4

Inspired by your help I added random effects like
fit2 <- brm(bf(YY~Aexp(-g1TT/1000)+Bexp(-g2TT/1000),
A~1+(1|RAT), B~1+(1|RAT),g1~1+(1|RAT),g2~(1+1|RAT),nl=TRUE),
data=dats,cores=4,prior = prior1,chains=4,control = list(adapt_delta = 0.99,max_treedepth = 20),iter = 8000)
and it worked very well and I recovered the used parameters and variance components handsomely.

But then I wanted to add treatment with 3 levels which was not allowed as I specified treatment as covariate like below since I want to know how treatment affects the recovery of the parameters in the doouble exponential decay:
fit3 <- brm(bf(YY~TRTAexp(-TRTg1TT/1000)+TRTBexp(-TRTg2TT/1000),
A~1+(1|RAT), B~1+(1|RAT),g1~1+(1|RAT),g2~(1+1|RAT),nl=TRUE),
data=dats,cores=4,prior = prior1,chains=4,control = list(adapt_delta = 0.99,max_treedepth = 20),iter = 8000)

and this doesn’t run - I can fit each group in isolation and compare - but this is against the grain. So how can I add treatment ad fixed effect to inform on each parameter? Is that possible?

5

Sorry for the confusion - I think I found the proper formulation -
fit3 <- brm(bf(YY~Aexp(-g1TT/1000)+Bexp(-g2TT/1000),
A~ 1+TRT+(1|RAT), B~ 1 + TRT+(1|RAT),g1~ 1+ TRT+(1|RAT),g2~ 1 + TRT +(1|RAT),nl=TRUE),
data=dats,cores=4,prior = prior1,chains=4,control = list(adapt_delta = 0.99,max_treedepth = 20),iter = 8000)

looking forward to see it end :-)

bw Troels