05-sim.Rmd

---
title: Simulation"
author: "Paloma Cárcamo"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output:
  html_document:
    toc: true
    toc_float: true
    df_print: paged
---

```{r setup, echo = FALSE}
if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, rjags, tidybayes)
knitr::opts_chunk$set(cache = TRUE)
```

```{r}
pcv7_st <- c("4","6B","9V","14","18C","19F","23F")
pcv13_st <- c("1","3","4","5","6A","6B","7F","9V","14","18C","19A","19F","23F")
pcv15_st <- c("1","3","4","5","6A","6B","7F","9V","14","18C","19A","19F","22F","23F","33F")
pcv20_st <- c("1","3","4","5","6A","6B","7F","8","9V","10A","11A","12F","14","15B","18C","19A","19F","22F","23F","33F")

data_raw <- read_rds("data/ABCs_st_1998_2021.rds")

data <- data_raw |> 
  rename(agec = "Age.Group..years.",
         year = Year,
         st = IPD.Serotype,
         N_IPD = Frequency.Count)  |> 
  mutate(st = if_else(st == '16', '16F', st))  |> 
  group_by(st, year) |> 
  summarize(N_IPD = sum(N_IPD))  |> 
  ungroup() |> 
  complete(year, st, fill = list(N_IPD = 0)) |> 
  mutate(pcv7 = if_else(st %in% pcv7_st, 1, 0),
         pcv13 = if_else(st %in% pcv13_st, 1, 0),
         pcv15 = if_else(st %in% pcv15_st, 1, 0),
         pcv20 = if_else(st %in% pcv20_st, 1, 0),
         firstvax = case_when(pcv7 == 1 ~ "pcv7",
                              pcv13 == 1 ~ "pcv13",
                              pcv15 == 1 ~ "pcv15",
                              pcv20 == 1 ~ "pcv20",
                              .default = ""))
```

```{r}
mod9_df <- data |> 
  complete(st, year, fill = list(N_IPD = 0))  |> 
  arrange(st, year) |> 
  filter(st %in% pcv13_st)  |> 
  group_by(st) |> 
  arrange(year) |> 
  mutate(year_n = row_number()) |> 
  ungroup()

mod9_mat <- mod9_df |> 
  reshape2::dcast(year ~ st, value.var = 'N_IPD') |> 
  dplyr::select(-year)

st_mapping9 <- setNames(colnames(mod9_mat), 1:length(pcv13_st))
year_mapping9 <- setNames(unique(mod9_df$year), 1:length(unique(mod9_df$year)))
```

```{r}
inits1 = list(".RNG.seed"=c(123), ".RNG.name"='base::Wichmann-Hill')
inits2 = list(".RNG.seed"=c(456), ".RNG.name"='base::Wichmann-Hill')
inits3 = list(".RNG.seed"=c(789), ".RNG.name"='base::Wichmann-Hill')
```

### Create simulated dataset

```{r}
N_years <- max(mod9_df$year_n)
N_sts <- ncol(mod9_mat)

set.seed(100)

# Priors
pi <- 0.5
alpha1 <- 3 # global intercept
tau_global <- 0.7 # st-specific variance

rho_beta1 <- 0.9 # autocorrelation for beta1
rho_beta2 <- 0.9 # autocorrelation for beta2
rho_eps <- 0.9 # autocorrelation for st-specific ar1s

tau_beta1 <- 0.4 # variance for beta1
tau_beta2 <- 0.2 # variance for beta2
tau_eps <- 0.3 # variance for st-specific ar1s

delta1 <- numeric(N_sts)
r <- numeric(N_sts)
grp <- numeric(N_sts)

for(j in 1:N_sts){
  delta1[j] <- 4
  r[j] <- runif(1, 0, 250)  # Serotype dispersion parameter
  grp[j] <- rbinom(1, 1, pi)
}

# Serotype-specific AR(1)
epsilon1 <- matrix(nrow = N_years, ncol = N_sts)

for(j in 1:N_sts){
  epsilon1[1,j] <- rnorm(1, 0, (1 - rho_eps^2) * tau_eps)
  
  for(i in 2:N_years){
    epsilon1[i,j] <- rnorm(1, rho_eps * epsilon1[i-1, j], tau_eps)
  }
}

# Global AR(1) - beta 1

beta1 <- numeric(N_years)

beta1[1] <- 0.2
for(i in 2:N_years){
  beta1[i] <- rnorm(1, rho_beta1 * beta1[i-1], tau_beta1)
}

# Global AR(1) - beta2

beta2 <- numeric(N_years)

beta2[1] <- 2
for(i in 2:N_years){
  beta2[i] <- rnorm(1, rho_beta2 * beta2[i-1], tau_beta2)
}

# Case counts

lambda <- matrix(nrow = N_years, ncol = N_sts)
prob <- matrix(nrow = N_years, ncol = N_sts)
N_IPD <- matrix(nrow = N_years, ncol = N_sts)

for(i in 1:N_years){
  
  for(j in 1:N_sts){
    lambda[i,j] <- exp(alpha1 + delta1[j] + epsilon1[i,j] + grp[j]*beta1[i] + ((1-grp[j])*beta2[i]))
    
    N_IPD[i,j] <- rpois(n = 1, lambda = lambda[i,j])
    
    }
}

pre_simulated_data <- data.frame(N_IPD)
colnames(pre_simulated_data) <- 1:N_sts
rownames(pre_simulated_data) <- 1:N_years

simulated_data <- pre_simulated_data |>
  rownames_to_column("year") |>
  pivot_longer(cols = 2:14, names_to = "st", values_to = "N_IPD") |>
  mutate(year = as.numeric(year),
         N_IPD = as.numeric(N_IPD),
         group = rep(grp, N_years),
         st = as.numeric(st))
```

```{r, echo = FALSE, fig.width = 11}
# simulated_data |> 
#   dplyr::select(group, st) |> 
#   distinct() |> 
#   mutate(group = if_else(group < 0.5, "Group 1 (beta1)", "Group 2 (beta2)")) |> 
#   group_by(group) |> 
#   summarise(serotypes = paste0(st, collapse = " ,"))

simulated_data |> 
  ggplot(aes(x = year, y = N_IPD, color = as.factor(st))) +
  geom_line() +
  labs(color = "serotype") +
  facet_wrap(group~., scales = "free_y") +
  theme_bw()

g0 <- simulated_data |> 
  filter(group == 0) |> 
  ggplot(aes(x = year, y = N_IPD, color = as.factor(st))) +
  geom_line(show.legend = FALSE) +
  labs(color = "serotype", title = "Group = 0") +
  facet_wrap(~st, scales = "free_y") +
  theme_bw()

g1 <- simulated_data |> 
  filter(group == 1) |> 
  ggplot(aes(x = year, y = N_IPD, color = as.factor(st))) +
  geom_line(show.legend = FALSE) +
  labs(color = "serotype", title = "Group = 1") +
  facet_wrap(~st, scales = "free_y") +
  theme_bw()

gridExtra::grid.arrange(g0, g1, ncol = 2)
```

```{r, echo = FALSE}
mod14_df <- simulated_data |> 
  ungroup() |> 
  complete(st, year, fill = list(N_IPD = 0))  |> 
  arrange(st, year) |> 
  group_by(st) |> 
  arrange(year) |> 
  mutate(year_n = row_number()) |> 
  ungroup()

mod14_mat <- mod14_df |> 
  reshape2::dcast(year ~ st, value.var = 'N_IPD') |> 
  dplyr::select(-year)

st_mapping14 <- setNames(colnames(mod14_mat), 1:length(unique(mod14_df$st)))
year_mapping14 <- setNames(unique(mod14_df$year), 1:length(unique(mod14_df$year)))
```

### Define model string

```{r}
mod14_string <- "
  model {
   
    for(i in 1:N_years){
     
      for(j in 1:N_sts){
        N_IPD[i,j] ~ dnegbin(prob[i,j], r[j])
       
        prob[i,j] <- r[j] / (r[j] + lambda[i,j])  # Likelihood
       
        log(lambda[i,j]) <- alpha1 + # Global intercept
                            delta1[j] + # Serotype intercept
                            epsilon1[i,j] + # Serotype-specific AR(1), 0 centered
                            grp[j]*(beta1[i]) + (1-grp[j])*(beta2[i]) # Global AR(1) version 1 or 2, 0 centered

      }
    }
   
    # Global AR(1) effect for beta1
   
    beta1[1] ~ dnorm(0, (1 - rho_beta1^2) * tau_beta1)
    for(i in 2:N_years){
      beta1[i] ~ dnorm(rho_beta1 * beta1[i-1], tau_beta1)
    }
   
    # Global AR(1) effect for beta2
   
    beta2[1] ~ dnorm(0, (1 - rho_beta2^2) * tau_beta2)
    for(i in 2:N_years){
      beta2[i] ~ dnorm( rho_beta2 * beta2[i-1], tau_beta2)
    }
       
    # Serotype-specific AR(1) effect
    for(j in 1:N_sts){
      epsilon1[1,j] ~ dnorm(0, (1 - rho_eps^2) * tau_eps)
     
      for(i in 2:N_years){
        epsilon1[i,j] ~ dnorm(0, tau_eps)
      }
    }
   
    # Priors for group selection
    for(j in 1:N_sts){
      grp[j] ~ dbern(pi)  # Group assignment - 0 or 1
    }
   
    pi ~ dunif(0,1)
       
    # Priors
    alpha1 ~ dnorm(0, 1e-4)
    tau_global ~ dgamma(0.01, 0.01)

    rho_beta1 ~ dunif(-1, 1)   # Uniform prior for rho_beta1 -- global AR(1) for beta1
    rho_beta2 ~ dunif(-1, 1)   # Uniform prior for rho_beta2 -- global AR(1) for beta2
    rho_eps ~ dunif(-1, 1)     # Prior for rho_eps same for all STs
   
    tau_beta1 ~ dgamma(3, 2)   # Tight prior for tau for beta1
    tau_beta2 ~ dgamma(3, 2)   # Tight prior for tau for beta2
    tau_eps ~ dgamma(3, 2)     # Tight prior for tau, shared for all serotypes
   
    for(j in 1:N_sts){
      delta1[j] ~ dnorm(0, tau_global)  # Serotype means centered around 0
      r[j] ~ dunif(0, 250)  # Serotype dispersion parameter
    }
  }
"
```

### Create JAGS model object

```{r, eval = FALSE}
mod14 <- jags.model(textConnection(mod14_string),
                   data = list("N_IPD" = mod14_mat,
                               "N_years" = max(mod14_df$year_n),
                               "N_sts" = ncol(mod14_mat)),
                   inits = list(inits1, inits2, inits3),
                   n.adapt = 10000,
                   n.chains = 3)
```

### Sample posteriors

```{r, eval = FALSE}
params14 <- c("alpha1", "epsilon1", "delta1", "beta1", "beta2", "tau_beta1", "tau_beta2", "lambda", "rho_beta1", "rho_beta2", "rho_eps", "pi", "grp", "r", "tau_global", "tau_eps", "prob", "N_IPD")

mod14_postsamp <- coda.samples(mod14, params14, n.iter = 50000)

# save(mod14_postsamp, file = "data/interim/mod14_postsamp.rda")
```

```{r, echo = FALSE}
load("data/interim/mod14_postsamp.rda")
```

### Plot traces for beta 1 and beta 2 {.tabset}

```{r, echo = FALSE}
beta_vars14 <- grep("^beta", varnames(mod14_postsamp), value = TRUE)
```

```{r, results = 'asis', echo = FALSE}
for (i in seq_along(beta_vars14)) {
  cat(paste0("\n#### ", "beta", substr(beta_vars14[i], 5, 5), " - ",
             gsub(".*\\[(\\d+)\\].*", "\\1", beta_vars14)[i], " \n"))
  cat("\n")
  traceplot(mod14_postsamp[, beta_vars14[i]], main = beta_vars14[i])
  cat("\n")
}
```

### Plot traces for alpha1 + beta + delta1 {.tabset}

#### Beta1 {.tabset}

```{r, echo = FALSE}
alpha1_samples <- mod14_postsamp[, "alpha1"]
delta1_samples <- mod14_postsamp[, grep("^delta1\\[", varnames(mod14_postsamp), value = TRUE)]
beta1_samples <- mod14_postsamp[, grep("^beta1\\[", varnames(mod14_postsamp), value = TRUE)]

N_years <- ncol(beta1_samples[[1]])
N_sts <- ncol(delta1_samples[[1]])

combined_effect14a <- list()
```

```{r, results = 'asis', echo = FALSE}
for (i in 1:N_years) {
  
  cat(paste0("\n##### ", year_mapping14[i], " {.tabset}", "\n"))
  cat("\n")
  
  for (j in 1:N_sts) {
    for (chain in 1:3) { 
      combined_effect14a[[chain]] <- alpha1_samples[[chain]] + delta1_samples[[chain]][, j] + beta1_samples[[chain]][, i]
      }
    
    combined_mcmc14a <- mcmc.list(combined_effect14a)
    
    cat(paste0("\n###### ", st_mapping14[j]," \n"))
    cat("\n")
    traceplot(combined_mcmc14a, main = paste("Year", year_mapping14[i], "ST", st_mapping14[j]), 
              xlab = "Iteration", ylab = "Combined Effect")
    cat("\n")
  }
}
```

#### Beta2 {.tabset}

```{r, echo = FALSE}
beta2_samples <- mod14_postsamp[, grep("^beta2\\[", varnames(mod14_postsamp), value = TRUE)]

combined_effect14b <- list()
```

```{r, results = 'asis', echo = FALSE}
for (i in 1:N_years) {
  
  cat(paste0("\n##### ", year_mapping14[i], " {.tabset}", "\n"))
  cat("\n")
  
  for (j in 1:N_sts) {
    for (chain in 1:3) { 
      combined_effect14b[[chain]] <- alpha1_samples[[chain]] + delta1_samples[[chain]][, j] + beta2_samples[[chain]][, i]
      }
    
    combined_mcmc14b <- mcmc.list(combined_effect14b)
    
    cat(paste0("\n###### ", st_mapping14[j]," \n"))
    cat("\n")
    traceplot(combined_mcmc14b, main = paste("Year", year_mapping14[i], "ST", st_mapping14[j]), 
              xlab = "Iteration", ylab = "Combined Effect")
    cat("\n")
  }
}
```

### Plot traces for alpha1 + delta[j] + grp[j]\*(beta1[i] ) + (1-grp[j])\*(beta2[i] ) {.tabset}

```{r, echo = FALSE}
grp_samples <- mod14_postsamp[, grep("^grp\\[", varnames(mod14_postsamp), value = TRUE)]

combined_effect14c <- list()
```

```{r, results = 'asis', echo = FALSE}
for (i in 1:N_years) {
  
  cat(paste0("\n#### ", year_mapping14[i], " {.tabset}", "\n"))
  cat("\n")
  
  for (j in 1:N_sts) {
    for (chain in 1:3) { 
      combined_effect14c[[chain]] <- alpha1_samples[[chain]] + 
        delta1_samples[[chain]][, j] + 
        beta1_samples[[chain]][, i]*grp_samples[[chain]][, j] +
        beta2_samples[[chain]][, j]*(1-grp_samples[[chain]][, j])
      }
    
    combined_mcmc14c <- mcmc.list(combined_effect14c)
    
    cat(paste0("\n##### ", st_mapping14[j]," \n"))
    cat("\n")
    traceplot(combined_mcmc14c, main = paste("Year", year_mapping14[i], "ST", st_mapping14[j]), 
              xlab = "Iteration", ylab = "Combined Effect")
    cat("\n")
  }
}
```

### AR1 groups

```{r, warning = FALSE, echo = FALSE}
groups14 <- gather_draws(mod14_postsamp, grp[j]) |> 
  median_hdci() |> 
  rename(st = j) |> 
  mutate(st = as.numeric(st_mapping14[as.character(st)]),
         group = if_else(.value > 0.5, "Group 1 (beta1)", "Group 2 (beta2)"))

simulated_data |> 
  dplyr::select(group, st) |> 
  distinct() |> 
  mutate(group = if_else(group < 0.5, "Group 1 (beta1)", "Group 2 (beta2)")) |> 
  group_by(group) |> 
  summarise(original_assignment = paste0(st, collapse = ", ")) |> 
  left_join(groups14 |> 
  group_by(group) |> 
  summarise(model_assignment = paste0(st, collapse = ", ")), by = "group")

groups14 |> 
  ggplot(aes(x = st)) +
  geom_pointrange(aes(y = .value, ymin = .lower, ymax = .upper)) +
  scale_x_continuous(breaks = c(1:13)) +
  theme_bw()
```

### Plot trajectory of global AR(1)s {.tabset}

#### AR1 - beta1

```{r, echo = FALSE}
mod14_beta1 <- gather_draws(mod14_postsamp, beta1[i]) |> 
  median_hdci() |> 
  mutate(year = year_mapping14[as.character(i)])

mod14_beta1 |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_line() +
  labs(x = "", y = "Global AR(1) - Beta 1") +
  theme_minimal()
```

#### AR1 - beta2

```{r, echo = FALSE}
mod14_beta2 <- gather_draws(mod14_postsamp, beta2[i]) |> 
  median_hdci() |> 
  mutate(year = year_mapping14[as.character(i)])

mod14_beta2 |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_line() +
  labs(x = "", y = "Global AR(1) - Beta 2") +
  theme_minimal()
```

### Plot obs vs pred {.tabset}

```{r, echo = FALSE}
mod14_summary <- gather_draws(mod14_postsamp, lambda[i,j]) |> 
  median_hdci() |> 
  rename(year = i,
         st = j) |> 
  mutate(st = as.numeric(st_mapping14[as.character(st)]),
         year = year_mapping14[as.character(year)]) |> 
  left_join(mod14_df, by = c("year", "st"))
```

#### Standard

```{r, fig.width = 11, echo = FALSE}
mod14_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()
```

#### Log-transformed y-axis

```{r, fig.width = 11, echo = FALSE}
mod14_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal() +
  coord_trans(y = "log1p")
```

#### Free y-axis

```{r, fig.width = 11, echo = FALSE}
mod14_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st, scales = "free_y") +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()
```

### Compare real vs model-estimated parameters {.tabset}

#### Single params

```{r, echo = FALSE, fig.width = 11, fig.height = 8}
single_pars <- gather_draws(mod14_postsamp, c(pi, alpha1, tau_global, rho_beta1, rho_beta2, rho_eps, tau_beta1, tau_beta2, tau_eps)) |> 
  median_hdci() |>
  mutate(obs = case_when(.variable == "pi" ~ pi,
                         .variable == "alpha1" ~ alpha1,
                         .variable == "tau_global" ~ tau_global,
                         .variable == "rho_beta1" ~ rho_beta1,
                         .variable == "rho_beta2" ~ rho_beta2,
                         .variable == "rho_eps" ~ rho_eps,
                         .variable == "tau_beta1" ~ tau_beta1,
                         .variable == "tau_beta2" ~ tau_beta2,
                         .variable == "tau_eps" ~ tau_eps))

single_pars |> 
  ggplot(aes(x = 1)) +
  geom_pointrange(aes(ymin = .lower, ymax = .upper, y = .value, color = "Estimated")) +
  geom_point(aes(y = obs, color = "Real")) +
  labs(x = "", y = "", color = "") +
  facet_wrap(~.variable, scales = "free_y") +
  theme_bw() +
  theme(axis.text.x = element_blank())
```

#### delta1

```{r, echo = FALSE}
delta_est <- gather_draws(mod14_postsamp, delta1[j]) |> 
  median_hdci() |> 
  rename(st = j) |> 
  mutate(obs = delta1)

delta_est |> 
  ggplot(aes(x = st)) +
  geom_pointrange(aes(ymin = .lower, ymax = .upper, y = .value, color = "Estimated")) +
  geom_point(aes(y = obs, color = "Real")) +
  labs(x = "Serotype", y = "", color = "", title = "delta1") +
  theme_bw() +
  scale_x_continuous(breaks = c(1:13))
```

#### r

```{r, echo = FALSE}
r_est <- gather_draws(mod14_postsamp, r[j]) |> 
  median_hdci() |> 
  rename(st = j) |> 
  mutate(obs = r)

r_est |> 
  ggplot(aes(x = st)) +
  geom_pointrange(aes(ymin = .lower, ymax = .upper, y = .value, color = "Estimated")) +
  geom_point(aes(y = obs, color = "Real")) +
  labs(x = "Serotype", y = "", color = "", title = "r") +
  theme_bw() +
  scale_x_continuous(breaks = c(1:13))
```

#### grp

```{r, echo = FALSE}
grp_est <- gather_draws(mod14_postsamp, grp[j]) |> 
  median_hdci() |> 
  rename(st = j) |> 
  mutate(obs = grp,
         obs = if_else(grp == 0, 1, 0))

grp_est |> 
  ggplot(aes(x = st)) +
  geom_pointrange(aes(ymin = .lower, ymax = .upper, y = .value, color = "Estimated"), size = 1.1) +
  geom_point(aes(y = obs, color = "Real")) +
  labs(x = "Serotype", y = "", color = "", title = "Groups") +
  theme_bw() +
  scale_x_continuous(breaks = c(1:13))
```

#### beta1

```{r, echo = FALSE}
beta1_est <- gather_draws(mod14_postsamp, beta1[i]) |> 
  median_hdci() |> 
  rename(year = i) |> 
  mutate(obs = beta1)

beta1_est |> 
  ggplot(aes(x = year)) +
  geom_ribbon(aes(ymin = .lower, ymax = .upper, fill = "Estimated"), alpha = 0.3, show.legend = FALSE) +
  geom_line(aes(y = .value, color = "Estimated")) +
  geom_line(aes(y = obs, color = "Real")) +
  labs(x = "Year", y = "", color = "", title = "beta1") +
  theme_bw()
```

#### beta2

```{r, echo = FALSE}
beta2_est <- gather_draws(mod14_postsamp, beta2[i]) |> 
  median_hdci() |> 
  rename(year = i) |> 
  mutate(obs = beta2)

beta2_est |> 
  ggplot(aes(x = year)) +
  geom_ribbon(aes(ymin = .lower, ymax = .upper, fill = "Estimated"), alpha = 0.3, show.legend = FALSE) +
  geom_line(aes(y = .value, color = "Estimated")) +
  geom_line(aes(y = obs, color = "Real")) +
  labs(x = "Year", y = "", color = "", title = "beta2") +
  theme_bw()
```

#### epsilon1

```{r, echo = FALSE, fig.width = 11, fig.height = 8}
epsilon_est <- gather_draws(mod14_postsamp, epsilon1[i,j]) |> 
  median_hdci() |> 
  rename(st = j,
         year = i) |> 
  left_join(data.frame(epsilon1) |> 
              rownames_to_column("year") |> 
              pivot_longer(cols = c(2:14), names_to = "st", values_to = "obs", names_prefix = "X") |> 
              mutate(year = as.numeric(year),
                     st = as.numeric(st)), 
            by = c("year", "st"))

epsilon_est |> 
  ggplot(aes(x = year)) +
  geom_ribbon(aes(ymin = .lower, ymax = .upper, fill = "Estimated"), alpha = 0.3, show.legend = FALSE) +
  geom_line(aes(y = .value, color = "Estimated")) +
  geom_line(aes(y = obs, color = "Real")) +
  labs(x = "Year", y = "", color = "", title = "epsilon1") +
  facet_wrap(~st) +
  theme_bw()
```

#### lambda

```{r, echo = FALSE, fig.width = 11, fig.height = 8}
lambda_est <- gather_draws(mod14_postsamp, lambda[i,j]) |> 
  median_hdci() |> 
  rename(st = j,
         year = i) |> 
  left_join(data.frame(lambda) |> 
              rownames_to_column("year") |> 
              pivot_longer(cols = c(2:14), names_to = "st", values_to = "obs", names_prefix = "X") |> 
              mutate(year = as.numeric(year),
                     st = as.numeric(st)), 
            by = c("year", "st"))

lambda_est |> 
  ggplot(aes(x = year)) +
  geom_ribbon(aes(ymin = .lower, ymax = .upper, fill = "Estimated"), alpha = 0.3, show.legend = FALSE) +
  geom_line(aes(y = .value, color = "Estimated")) +
  geom_line(aes(y = obs, color = "Real"), lty = 2) +
  labs(x = "Year", y = "", color = "", title = "lambda") +
  facet_wrap(~st, scales = "free_y") +
  theme_bw()
```