fixed_factors.R (2.2 KB)
Dear friends - I simulated some data to represent actual experimental data on 3 groups of rats on which double exponential decay of a tracer is modelled. With your help I got a fine fit with Rhat values = 1 as long as I only modled one population. Doing all three together as shown below was more problematic. Rhat still being 1 pp_check had an unreproduced large central peak. And for the three groups fitted values didn’t match the original data too well except for the controls. I noticed that it didn’t matter if I had specified as below or as A + B + g1 + g2 ~ TRT. So I wonder how I happened to miss-specify the fixed effects of TRT
fit4 <- brm(bf(YY ~ A * exp( - g1 * TT/1000)+B * exp( -g2 * TT/1000),
A + B + g1 + g2 ~ 1+TRT,nl=TRUE, sigma ~ TRT),data=dats,cores=4,prior = prior1,chains=4,control = list(adapt_delta = 0.8,max_treedepth = 10),iter = 2000)
I uploaded the simulation and the model code for plots
I’m just going back now and picking up unanswered questions. This one might be better tagged a brms question (I assume that’s what the brm is doing). I think they have their own message board, but Paul sometimes answers questions here.
Thanks a lot - it appeared that the model was very sensititive to the prior specification - I still struggle to understand how one specification produces “perfect” reproduction of measured values and posterior checks but the increasing the support in gaussian priors makes the model run wildely astray - so even in spite of bayesianR2 > 0.998 I’m uncertain I have really captured the whole relevant space. Will continue in the suggested forum - even if the issue may be general. Thanks a lot!
I usually try to have all model parameters in a similar scale to ease sampling, such as avoiding one parameter to go in the 100s and another one close to 0. So try to re-scale your equations to that.
I am a bit puzzled as to how the model can deal with having the same variables (TT) appearing twice in the regression formula, I never worked with double decay exponential model before but it looks strange to me, maybe re-parametrizing the model to have TT appearing only once on the right-hand side of the equation would help.
The prior you set on some parameters are super tight, like normal(275,2), this is a very constrained one. Have a look at curve(dnorm(x,275,2),from=250,to=300) to get an idea.
In the model you set the residual variation to depend on TRT, but in your simulation process the residual variation depend on the non-linear predictor, line 19 of your code is:
Thank you very much!!! - yes I can see I have to rescale, and also to put much broader priors to make sure to have the fitting running - but these are intracate individually and in combination even more. The model is over time - I should have made it clear that TT was time - the double exponential decay is old school model in pharmacokinetics to model first order distribution and elimination - in medical school we were taught “curve peeling” in which in a semilogarithmic plot the individual components could be separated if they had sufficiently different rate constants. I just hoped that using brute force in brms/stan I could write up directly the relationship and identify the parameters. Not so easy. As for the residual variation I can see it was not a good idea and have made gaussian errors.
With quite a bit of fiddling around I have managed to fit the actual experimental data very well in a hierarchical mixed model with treatments as fixed effects and random effects corresponding to each of the four parameters corresponding to each rat - but it appears to be quite fragile.