## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----duration-simulation------------------------------------------------------
library(gp3bayes)

duration_simulation <- simulate_hierarchical_duration_data(
  n_participants = 16,
  trials_per_participant = 10,
  n_items = 8,
  baseline_median = 500,
  random_slope_sd = 0.15,
  residual_sd = 0.40,
  outcome_unit = "milliseconds",
  seed = 2030
)

duration_simulation
head(duration_simulation$data)

## ----duration-contract--------------------------------------------------------
duration_contract <- create_model_contract(
  family = "duration",
  outcome_col = "duration",
  participant_col = "participant_id",
  item_col = "item_id",
  trial_col = "trial_id",
  condition_col = "condition",
  predictors = c(
    "participant_covariate",
    "trial_covariate"
  ),
  interaction = c(
    "condition",
    "participant_covariate"
  ),
  random_slope = TRUE,
  outcome_unit = "milliseconds"
)

duration_contract

## ----duration-preparation-----------------------------------------------------
duration_prepared <- prepare_hierarchical_duration_data(
  duration_simulation$data,
  duration_contract,
  condition_levels = c(
    "control",
    "treatment"
  )
)

duration_prepared
duration_prepared$decision_log

## ----duration-unit-conversion-------------------------------------------------
duration_seconds <- prepare_hierarchical_duration_data(
  duration_simulation$data,
  duration_contract,
  condition_levels = c(
    "control",
    "treatment"
  ),
  outcome_multiplier = 0.001,
  converted_unit = "seconds"
)

duration_seconds$transformations$outcome

## ----duration-specification---------------------------------------------------
duration_specification <- specify_duration_model(
  duration_prepared,
  baseline = 500,
  intercept_scale = 1,
  coefficient_scale = 0.5,
  group_sd_scale = 1,
  residual_scale = 1,
  correlation_eta = 2,
  student_df = 3
)

duration_specification
duration_specification$priors$table

## ----duration-prior-predictive------------------------------------------------
duration_prior_predictive <- check_duration_prior_predictive(
  duration_specification,
  draws = 100,
  seed = 2031
)

duration_prior_predictive
duration_prior_predictive$checks

## ----duration-fit, eval=FALSE-------------------------------------------------
# duration_translation <- translate_duration_model_to_brms(
#   duration_specification
# )
# 
# duration_fit <- fit_duration_model(
#   duration_specification,
#   chains = 2,
#   iter = 2000,
#   warmup = 1000,
#   cores = 2,
#   seed = 2032,
#   adapt_delta = 0.95,
#   max_treedepth = 12,
#   refresh = 100
# )

## ----duration-validation, eval=FALSE------------------------------------------
# duration_diagnostics <- diagnose_duration_fit(
#   duration_fit
# )
# 
# duration_posterior <- summarise_duration_posterior(
#   duration_fit
# )
# 
# duration_predictive <- check_duration_posterior_predictive(
#   duration_fit,
#   draws = 500,
#   seed = 2033
# )

## ----duration-sensitivity-recovery, eval=FALSE--------------------------------
# duration_sensitivity <- assess_duration_prior_sensitivity(
#   duration_fit,
#   scale_multipliers = c(
#     tighter = 0.5,
#     wider = 2
#   )
# )
# 
# duration_recovery <- run_duration_recovery(
#   repetitions = 20,
#   baseline_median = 500,
#   outcome_unit = "milliseconds",
#   seed = 4001
# )
# 
# duration_report_file <- tempfile(fileext = ".md")
# 
# duration_report <- create_duration_model_report(
#   duration_fit,
#   diagnostics = duration_diagnostics,
#   posterior_summary = duration_posterior,
#   posterior_predictive = duration_predictive,
#   prior_sensitivity = duration_sensitivity,
#   recovery = duration_recovery,
#   file = duration_report_file
# )
# 
# duration_report
# unlink(duration_report_file)