Skip to content

Instantly share code, notes, and snippets.

@dlebauer
Created September 15, 2021 21:44
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save dlebauer/11796f7e55ac51fa45b0091513d14a1c to your computer and use it in GitHub Desktop.
Save dlebauer/11796f7e55ac51fa45b0091513d14a1c to your computer and use it in GitHub Desktop.
sum_out <- broom.mixed::tidyMCMC(coda_out, conf.int = TRUE, conf.level = 0.95) %>%
rename(mean = estimate, sd = std.error, pc2.5 = conf.low, pc97.5 = conf.high) %>%
mutate(sig = if_else(pc2.5 * pc97.5 > 0, TRUE, FALSE)) # test if UCL and LCL both have the same sign
years <- data.frame(year_index = 1:length(unique(x$year)), year = unique(x$year))
states <- data.frame(state_index = 1:length(unique(x$state)), state = unique(x$state))
sum_rep <- broom.mixed::tidyMCMC(coda_rep, conf.int = TRUE, conf.level = 0.95) %>%
rename(mean = estimate, sd = std.error, pc2.5 = conf.low, pc97.5 = conf.high) %>%
separate(term, sep = '\\[|\\]|,', into = c('term', 'i', 'j', 'other'), convert = TRUE) %>%
mutate(var = case_when(j == 1 ~ "area~(kha)",
j == 2 ~ "number~(k)",
j == 3 ~ "sales~(MM)"),
year_index = i %% 6 + 1,
state_index = i %% 50 + 1) %>%
left_join(years) %>%
left_join(states)
sum_rep2 <- postjags::coda.fast(chains = 3, burn.in = 0, thin = 1, coda = coda_rep)
sum_rep$year <- rep(dat$year + 2014, each = 3)
sum_rep$Var <- rep(c("area~(kha)", "number~(k)", "sales~(MM)"), 300)
sum_rep$state <- rep(dat$state, each = 3)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment