I am currently trying to fit a Bayesian model on repeated measures data, but seem to be facing an obscure problem.

For a start, I did prior predictive check and my model seemed to be behaving well (Rhats = 1 and 0 divergent transitions), hence I decided to proceed with model fitting. I run 4 chains of 6000 iterations (with warmup 1000), which seemed to be enough - my full model, again, converged with no indications of problems.

I then proceeded to testing (on my training data) whether the model would predict it well, however, the results I got were way off (with order of magnitude around 15 for most predictions, compared to the real data). I then tried to fit the same model with lme4 and the output I got was reasonable.

**So I am very curious where this discrepancy might be coming from? And if, for example, the problem is with my priors, why was the predictive prior check OK and didn’t throw any problems?**

My code is given below:

```
> library(brms)
> library(bayesplot)
> library(rstanarm)
> # Read source data in .csvformat
> ldlData <- read.csv("sorted_with_category.csv", header = TRUE)
>
> # Rename columns to make faster / more intuitive use
> colnames (ldlData) <- c ("Sample", "Time", "LDL","Cat" )
>
> #Turning string column into a column of numbers (if necessary)
> ldlData$LDL <-as.numeric(as.character(ldlData$LDL))
> ldlData$Time <-as.numeric(as.character(ldlData$Time))
> View(ldlData)
> quantile(ldlData[,"LDL"], na.rm=T, probs=seq(0, 1, by=.05))
> plot (LDL ~ Time, data = ldlData, pch = 16)
>
> # defining the prior
> "prior <- brms::get_prior(
> brms::bf(LDL ~ 15 + (alpha-15) / (1 + exp ((gamma-Time) / delta)),
> alpha ~ 1 + (1|Cat)+(1|Sample),
> gamma ~ 1 + (1|Cat)+(1|Sample),
> delta ~ 1,
> nl = TRUE),
> data = ldlData,
> family = lognormal())
> prior"
> fit<- brm(
> bf(LDL ~ 15 + (alpha-15) / (1 + exp ((gamma-Time) / delta)),
> alpha ~ 1 + (1|Cat)+(1|Sample),
> gamma ~ 1 + (1|Cat)+(1|Sample),
> delta ~ 1,
> nl = TRUE),
> prior = c(
> prior(lognormal(8.33, 0.15), class="b", lb=3800, nlpar = "alpha"),
> prior(gamma(42,15), class="b", lb=0,ub=5, nlpar = "gamma"),
> prior(normal(0.4,0.02), class="b", lb=0, ub=1, nlpar = "delta"),
> prior(student_t(3, 0, 10), class="sigma"),
> prior(student_t(3, 0, 10), class="sd", nlpar="alpha"),
> prior(student_t(3, 0, 10), class="sd", nlpar="gamma"),
> prior(normal(-0.002420144,0.482032688 ), class="sd", group="Sample", coef="Intercept", nlpar="gamma"),
> prior(normal(0.01815935,0.52960660 ), class="sd", group="Sample", coef="Intercept", nlpar="alpha"),
> prior(normal(-0.002420144,0.482032688 ), class="sd", group="Cat", coef="Intercept", nlpar="gamma"),
> prior(normal(0.01815935,0.52960660 ), class="sd", group="Cat", coef="Intercept", nlpar="alpha")),
> data = ldlData, family = lognormal,
> chains = 1,
> iter=3000,
> warmup = 1000,
> cores = getOption("mc.cores",1L),
> thin = 5,
> control = list(adapt_delta = 0.99999, stepsize=0.000001, max_treedepth=10)
> )
>
> predict(fit)
```