library(brms) library(tidyr) library(dplyr) library(DirichletReg) library(ggplot2) options(mc.cores = parallel::detectCores()) set.seed(63234) make_AR = function(df) { # df$y1[1] = 1 df$y_d2[1] = 1 df$y_d3[1] = 1 for (t in 2:N) { df$mu2[t] = df$r_i_G2[t] + df$r_i_T[t] + 0 * df$x[t] + 0 * df$cat[t] + 0.8 * df$y_d2[t-1] df$mu3[t] = df$r_i_G3[t] + df$r_i_T[t] + 0 * df$x[t] + 0 * df$cat[t] + 0.8 * df$y_d3[t-1] df$y1[t] = exp(df$mu1[t]) / (exp(df$mu1[t]) + exp(df$mu2[t]) + exp(df$mu3[t])) df$y2[t] = exp(df$mu2[t]) / (exp(df$mu1[t]) + exp(df$mu2[t]) + exp(df$mu3[t])) df$y3[t] = exp(df$mu3[t]) / (exp(df$mu1[t]) + exp(df$mu2[t]) + exp(df$mu3[t])) df[t,c("y_d1", "y_d2", "y_d3")] = rdirichlet(1, cbind(df$y1[t], df$y2[t], df$y3[t])*exp(phi)) } return(df) } bind <- function(...) cbind(...) phi = 1 N = 20 G = 10 df <- data.frame( id = rep(1:G, each=N), year = rep(1:N, G), r_i_G2 = rep(rnorm(G,0, 1), each=N), r_i_G3 = rep(rnorm(G,0, 1), each=N), r_i_T = rep(rnorm(N,0, 0.5), G), y1 = 0, y2 = NA, y3 = NA, x = rnorm(N*G, 0,2), cat = rep(sample(0:1,G, replace=T), each=N), y_d1 = NA, y_d2 = NA, y_d3 = NA, mu1 = 0, mu2 = NA, mu3 = NA ) AR_df = df %>% group_by(id) %>% do(make_AR(.)) # Serial Correlation is present: AR_df %>% select(id, year, x, y_d1, y_d2, y_d3) %>% mutate_all(funs(round (., 3))) %>% filter(id == 9) %>% pivot_longer(c(y_d1,y_d2,y_d3) ) %>% ggplot(aes(x=year, y=value, col=name)) + geom_line(size=1) my_formula = bf( bind(y_d1, y_d2, y_d3) ~ 1, muyd2 ~ 1 + (1|id) + (1|year), muyd3 ~ 1 + (1|id) + (1|year), autocor = cor_ar(form = ~ year | id, 1, cov=T) ) get_prior(my_formula, AR_df, family=dirichlet()) priors <- c(set_prior("normal(0,100)", class = "Intercept"), set_prior("cauchy(0,5)", class = "sd", dpar = "muyd2"), set_prior("cauchy(0,5)", class = "sd", dpar = "muyd3") ) fit <- brm(my_formula, AR_df, family=dirichlet(), prior = priors, chains=4) summary(fit)