library(cmdstanr) library(rstan) library(bayesplot) #set.seed(54) N <- 3000 N_groups <- 7 groups <- rep(1:N_groups, length.out = N) intercept <-0 group_sd <- 1 group_mean <- 1.5 group_r <- rnorm(N_groups, mean = group_mean, sd = group_sd) mu <- intercept + group_r[groups] y <- rnorm(N, mu, 0.25) data_stan <- list(N = N, N_groups = N_groups, groups = groups, y = y) n_chains <- 3 mod <- cmdstan_model("sumtozero_gaussian_QR_grouplevel_mean.stan") fit <- mod$sample( data = data_stan , seed = 12 , chains = n_chains , parallel_chains = n_chains # , init = init_ll , refresh = 10 , iter_warmup = 1000 , iter_sampling = 1000 , max_treedepth = 14 , adapt_delta = 0.8 ) lm_waic <- rstan::read_stan_csv(fit$output_files()) extr<-extract(lm_waic) apply(extr$group_r_hat,2,mean) #mean(extr$intercept) mean(extr$group_sd) mean(extr$group_mean) mean(extr$sigma) group_r #pairs(lm_waic,pars=c("intercept","group_r_hat","group_mean")) pairs(lm_waic,pars=c("group_r_hat","group_mean"))