library(brms)
options(mc.cores = parallel::detectCores())
library(ggplot2)

set.seed(1234)

id <- rep(1:30, each=4)
n <- length(id)
time <- rep(1:4, times=30)
a <- 10
b <- 5
z <- rnorm(30,0,1)
y <- (a + z[id]) + b*time
sd <- 1
outcome <- rnorm(n, y, sd)

id2 <- rep(31:90, each=4)
n2 <- length(id2)
time2 <- rep(1:4, times=60)
a2 <- 10
b2 <- 10
z2 <- rnorm(90,0,1)   
y2 <- (a2 + z2[id2]) + b2*time2
sd2 <- 1
outcome2 <- rnorm(n2, y2, sd2)

ID <- c(id, id2)
TIME <- c(time, time2)
OUTCOME <- c(outcome, outcome2)

data1 <- cbind(ID,TIME,OUTCOME)
data1 <- data.frame(data1)
data1$ID <- factor(data1$ID)

str(data1)

p <- ggplot(data1, aes(x=TIME, y=OUTCOME, colour=ID)) +
    geom_point() +
    geom_line()
p


#fit mixture models
mix <- mixture(gaussian, gaussian)
prior1 <- c(
  prior(normal(8, 2), Intercept, dpar = mu1),
  prior(normal(8, 2), Intercept, dpar = mu2),
  prior(normal(8, 2.5), b, dpar = mu1),
  prior(normal(8, 2.5), b, dpar = mu2),
  prior(normal(0, 2), sigma1),
  prior(normal(0, 2), sigma2),
  prior(normal(0, 1), sd, dpar = mu1),
  prior(normal(0, 1), sd, dpar = mu2)
)

m1 <- brm(bf(OUTCOME ~ 1 + TIME + (1|ID)), data = data1, family = mix,
            prior = prior1, iter = 2000, chains = 4, control = list(adapt_delta = 0.9999, stepsize=0.5, max_treedepth=13)) 
m1

prior2 <- c(
  prior(normal(8, 2), Intercept, dpar = mu1),
  prior(normal(12, 2), Intercept, dpar = mu2),
  prior(normal(8, 2.5), b, dpar = mu1),
  prior(normal(8, 2.5), b, dpar = mu2),
  prior(normal(0, 2), sigma1),
  prior(normal(0, 2), sigma2)
)

m2 <- brm(bf(OUTCOME ~ 1 + TIME), data = data1, family = mix,
            prior = prior2, iter = 2000, chains = 4, control = list(adapt_delta = 0.9999, stepsize=0.5, max_treedepth=13)) 
m2