Thanks and nice recommendation. Here is what I tried, reduced draw time from 50s (model only took 23s to fit with 4 chians,1000 warm up and 1000 samples) to 3s, if I just use fit$draws
cmdstanfiles <- list()
for (f in fit$output_files()) {
cmdstanfiles[[f]] <- fread(cmd= paste0("grep -v '^#' ", f))
}
var_of_interest <- cmdstanfiles %>% bind_rows() %>% dplyr::select(starts_with("var_of_interest")) %>% as.matrix()