library(mrbrt, lib.loc = "/ihme/code/mscm/R/packages/")
## Loading required package: reticulate
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(RColorBrewer)



set.seed(1)
k_studies <- 20
n_per_study <- 5
tau_1 <- 8
sigma_1 <- 2
tau_2 <- 0.6
sigma_2 <- 0.2

df_sim_study <- data.frame(study_id = as.factor(1:k_studies)) %>%
  mutate(
    study_effect1 = rnorm(n = k_studies, mean = 0, sd = tau_1),
    study_effect2 = rnorm(n = k_studies, mean = 0, sd = tau_2)
    # study_colors = brewer.pal(n = k_studies, "Spectral")
  )

df_sim <- do.call("rbind", lapply(1:nrow(df_sim_study), function(i) {
  df_sim_study[rep(i, n_per_study), ] })) %>%
  mutate(
    x1 = runif(n = nrow(.), min = 0, max = 10),
    y1 = 1 + 2*x1 + study_effect1 + rnorm(nrow(.), mean = 0, sd = sigma_1),
    y2_true = 2 * sin(0.43*x1-2.9),
    y2 = y2_true - min(y2_true) + study_effect2 + rnorm(nrow(.), mean = 0, sd = sigma_2),
    is_outlier = FALSE) %>%
  arrange(x1)

df_sim2 <- df_sim %>%
  rbind(., .[(nrow(.)-7):(nrow(.)-4), ] %>% mutate(y1=y1-18, y2=y2-4, is_outlier = TRUE))
  
# with(df_sim, plot(x1, y, xlim = c(0, 10), ylim = c(0, max(df_sim$y))))
with(df_sim2, plot(x1, y1, xlim = c(0, 10)))

with(df_sim, plot(x1, y1, xlim = c(0,10), type = "n"))

with(df_sim, plot(x1, y1, xlim = c(0,10)))
for (x in unique(df_sim$study_id)) with(filter(df_sim, study_id == x), lines(x1, y1))

with(df_sim2, plot(x1, y1, xlim = c(0,10)))
# with(df_sim2, plot(x1, y1, xlim = c(0,10), type = "n"))
for (x in unique(df_sim2$study_id)) with(filter(df_sim2, study_id == x), lines(x1, y1))

## Loading required package: Matrix

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.