################################################################################
### setup
################################################################################

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

# packages
# https://github.com/stan-dev/cmdstanr
# install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
require(pacman)
pacman::p_load(tidyverse, tidybayes, tidytext, brms, cmdstanr, rstan, coda,  
               bayesplot, posterior, performance, r2mlm, purrr, scales,    
               emmeans, broom, broom.mixed, patchwork, ggpomological, magick,
               marginaleffects, ggdist, ggokabeito, gghalves, ggbeeswarm,
               GGally, bayestestR, modelr, ggridges, ggdist, distributional,
               grid, gridExtra, ggh4x, flextable, janitor)

# global options
options(scipen=20)

# global options for marginal effects
options("marginaleffects_posterior_interval" = "hdi")
options("marginaleffects_posterior_center" = "median")


################################################################################
### utility functions
################################################################################

# source functions file
source("0a_functions.R")

  
################################################################################
### data
################################################################################

# load and clean data
read_csv("../data/foot_trait_and_physical_activity_data.csv", na = "") |> 
  janitor::clean_names(case="snake") |>
  filter(!if_all(everything(), is.na)) |> 
  mutate(sex = factor(sex, levels=c("female", "male")),
         daily_step_count_1000 = daily_step_count / 1000,
         mvpa_time_per_day_15 = mvpa_time_per_day / 15,
         static_arch_stiffness_index_1000 = static_arch_stiffness_index / 1000,
         dynamic_arch_stiffness_index_1000 = dynamic_arch_stiffness_index / 1000,
         abductor_hallucis_csa_10 = abductor_hallucis_csa / 10,
         flexor_digitorum_brevis_csa_10 = flexor_digitorum_brevis_csa / 10,
         abductor_digiti_minimi_csa_10 = abductor_digiti_minimi_csa / 10) ->             
dat


################################################################################
### data inspection
################################################################################

glimpse(dat)

# outcomes
# arch_height_index                         
# static_arch_stiffness_index         
# dynamic_arch_stiffnes_index               
# abductor_hallucis_csa       
# flexor_digitorum_brevis_csa 
# abductor_digiti_minimi_csa

# covariates
# walking_speed_m_sec
# truncated_foot_length
# age
# sex

# predictors
# mvpa_time_per_day
# daily_step_count

# SUMMARIES
summary(dat$arch_height_index)
summary(dat$static_arch_stiffness_index)
summary(dat$dynamic_arch_stiffness_index)
summary(dat$abductor_hallucis_csa)
summary(dat$flexor_digitorum_brevis_csa)
summary(dat$abductor_digiti_minimi_csa)
summary(dat$daily_step_count)
summary(dat$mvpa_time_per_day)
summary(dat$age)
summary(dat$sex)
summary(dat$walking_speed)
summary(dat$truncated_foot_length)

# DENSITIES
plot(density(na.omit(dat$arch_height_index))) # positive below 1
plot(density(na.omit(dat$static_arch_stiffness_index))) # positive in thousands
plot(density(na.omit(dat$dynamic_arch_stiffness_index))) # positive in thousands
plot(density(na.omit(dat$abductor_hallucis_csa))) # positive in single digits
plot(density(na.omit(dat$flexor_digitorum_brevis_csa))) # positive in single digits
plot(density(na.omit(dat$abductor_digiti_minimi_csa))) # positive in single digits
plot(density(na.omit(dat$daily_step_count))) # positive in thousands
plot(density(na.omit(dat$mvpa_time_per_day))) # positive in hundreds



################################################################################
### priors
################################################################################

# student_t prior visualization
expand.grid(df = c(3, 5, 10, 30),
            scale = c(0.5, 1, 2.5, 5)) |>
  ggplot(aes(y = 0, dist="student_t", arg1=df, arg2=0, 
             arg3=scale, color=ordered(df))) +
  stat_slab(p_limits = c(.01, .99), fill = NA, linewidth=0.5) +
  scale_y_continuous(breaks = NULL) +
  facet_grid( ~ scale) +
  labs(title = "dstudent_t(x, df, 0, sigma)",
       subtitle = "Scale (sigma)", x=NULL, y=NULL) +
  theme_ggdist() ->
priors_student_t


# weak priors
tibble(sigma = c(dist_normal(0, 3),
                 dist_normal(0, 1), 
                 dist_student_t(5, 0, 1),
                 dist_student_t(5, 0, 2.5),
                 dist_student_t(3, 0, 2.5))) |>
  ggplot(aes(xdist=sigma, y=as.character(sigma))) +
  stat_slab() +
  xlim(-10, 10) +
  theme_pomol() ->
priors_weak

# more informative priors
tibble(sigma = c(dist_normal(0, 0.1),
                 dist_normal(0, 0.3), 
                 dist_student_t(3, 0, 0.1),
                 dist_student_t(5, 0, 0.1),
                 dist_student_t(5, 0, 0.3))) |>
  ggplot(aes(xdist=sigma, y=as.character(sigma))) +
  stat_slab() +
  xlim(-1, 1) +
  theme_pomol() ->
priors_inform

# variance priors
tibble(sigma = c(exp(dist_normal(0, 1)), 
                 exp(dist_student_t(5, 0, 2.5)),
                 exp(dist_student_t(3, 0, 2.5)))) |>
  ggplot(aes(xdist=sigma, y=as.character(sigma))) +
  stat_slab(p_limits = c(0, 0.8)) +
  xlim(0, 5) +
  theme_pomol() ->
priors_variance