-make range of values more constrained in shiny app

-make app color coding fixed and same as in the paper
-fix issue with plots in the demo

Author: ArCourtin

apps & wrappers/BRMS demo.Rmd

@@ -196,22 +196,22 @@ fit_eip_HRDT <-
     prior = empirically_informed_priors_HRDT
     )
 
-fit_wip_HRDT <-
-  brm(
-    bff_HRDT,
-    data=sim_data_HRDT,
-    chains = 4,
-    cores = 4,
-    seed = seed,
-    backend = "cmdstan",
-    warmup=2000,
-    iter=4000,           
-    control=list(
-     adapt_delta=0.9,
-     max_treedepth = 10
-     ),
-    prior = weakly_informative_priors_HRDT
-    )
+# fit_wip_HRDT <-
+#   brm(
+#     bff_HRDT,
+#     data=sim_data_HRDT,
+#     chains = 4,
+#     cores = 4,
+#     seed = seed,
+#     backend = "cmdstan",
+#     warmup=2000,
+#     iter=4000,           
+#     control=list(
+#      adapt_delta=0.9,
+#      max_treedepth = 10
+#      ),
+#     prior = weakly_informative_priors_HRDT
+#     )
 ```
 
 ## Model diagnosing
@@ -296,15 +296,24 @@ sprintf(
 We can then construct plots that show the posterior distribution for the group-level differences in threshold and slope. We will also add a line range materializing the 95% CI.
 
 ```{r}
+summary_stats_HRDT_plot<-
+  summary_stats_HRDT %>% 
+  filter(grepl('diff',summary_stats_HRDT$name)) %>% 
+  mutate(name=factor(name,c('mean_alpha_diff','mean_beta_diff'),c('Threshold','Log(Slope)'))) 
+
+  
 as_draws_df(fit_eip_HRDT) %>% 
   select(b_alpha_conditiontreatment,b_beta_conditiontreatment) %>% 
+  rename(
+    mean_alpha_diff=b_alpha_conditiontreatment,
+    mean_beta_diff=b_beta_conditiontreatment
+  ) %>%
   pivot_longer(everything()) %>% 
-  full_join(summary_stats_HRDT) %>% 
-  mutate(name=factor(name,c('b_alpha_conditiontreatment','b_beta_conditiontreatment'),c('Threshold','Log(Slope)'))) %>% 
+  mutate(name=factor(name,c('mean_alpha_diff','mean_beta_diff'),c('Threshold','Log(Slope)'))) %>% 
   ggplot()+
   geom_histogram(aes(x=value),color='black')+
   geom_vline(aes(xintercept = 0),linetype='dotted',color='grey',linewidth=1)+
-  geom_pointrange(aes(x=m,xmin=LB,xmax=UB,y=-20,group = name),linewidth=1)+
+  geom_pointrange(data=summary_stats_HRDT_plot,aes(x=m,xmin=LB,xmax=UB,y=-20,group = name),linewidth=1)+
   labs(
     title='Group mean differences between conditions',
     subtitle='HRDT',
@@ -645,22 +654,22 @@ fit_eip_RRST <-
     prior = empirically_informed_priors_RRST
     )
 
-fit_wip_RRST <-
-  brm(
-    bff_RRST,
-    data=sim_data_RRST,
-    chains = 4,
-    cores = 4,
-    seed = seed,
-    backend = "cmdstan",
-    warmup=1000,
-    iter=2000,   
-    control=list(
-     adapt_delta=0.99,
-     max_treedepth = 10
-     ),
-    prior = weakly_informative_priors_RRST
-  )
+# fit_wip_RRST <-
+#   brm(
+#     bff_RRST,
+#     data=sim_data_RRST,
+#     chains = 4,
+#     cores = 4,
+#     seed = seed,
+#     backend = "cmdstan",
+#     warmup=1000,
+#     iter=2000,   
+#     control=list(
+#      adapt_delta=0.99,
+#      max_treedepth = 10
+#      ),
+#     prior = weakly_informative_priors_RRST
+#   )
 ```
 
 ## Model diagnosing
@@ -780,6 +789,7 @@ spread_diffs<-
 spread_diffs_stats<-
   spread_diffs %>% 
   summarise(
+    name='spread_difference',
     m=mean(difference),
     LB=quantile(difference,0.025),
     UB=quantile(difference,0.975),
@@ -833,28 +843,25 @@ sprintf(
 We can then construct plots that show the posterior distribution for the group-level differences in threshold and spread. We will also add a linerange materializing the 95% CI.
 
 ```{r}
-spread_differences_long<-
+spread_differences<-
   spread_diffs %>% 
-  select(difference) %>% 
-  rename(value=difference) %>% 
-  mutate(
-    name='spread_difference',
-    m=spread_diffs_stats$m,
-    LB=spread_diffs_stats$LB,
-    UB=spread_diffs_stats$UB,
-    p_larger_than_0=spread_diffs_stats$p_larger_than_0,
-    pseudo_p_value=spread_diffs_stats$pseudo_p_value
-    ) 
+  select(difference)
+summary_RRST_plot<-
+  summary_stats_RRST %>% 
+  filter(name=='mean_alpha_diff') %>% 
+  full_join(spread_diffs_stats) %>% 
+  mutate(name=factor(name,c('mean_alpha_diff','spread_difference'),c('Log(Threshold)','Spread')))
+
 as_draws_df(fit_eip_RRST) %>% 
+  select(b_alpha_conditiontreatment) %>% 
+  rename(mean_alpha_diff=b_alpha_conditiontreatment) %>% 
+  mutate(spread_difference=spread_differences$difference) %>% 
   pivot_longer(everything()) %>% 
-  full_join(summary_stats_RRST) %>% 
-  filter(name=='b_alpha_conditiontreatment') %>% 
-  full_join(spread_differences_long) %>% 
-  mutate(name=factor(name,c('b_alpha_conditiontreatment','spread_difference'),c('Log(Threshold)','Spread'))) %>% 
+  mutate(name=factor(name,c('mean_alpha_diff','spread_difference'),c('Log(Threshold)','Spread'))) %>% 
   ggplot()+
   geom_histogram(aes(x=value),color='black')+
   geom_vline(aes(xintercept = 0),linetype='dotted',color='grey',linewidth=1)+
-  geom_pointrange(aes(x=m,xmin=LB,xmax=UB,y=-10,group = name))+
+  geom_pointrange(data=summary_RRST_plot,aes(x=m,xmin=LB,xmax=UB,y=-10,group = name))+
   labs(
     title='Group mean differences between conditions',
     subtitle='RRST',

apps & wrappers/shiny app.R

@@ -22,8 +22,8 @@ ui <- navbarPage("Bayesian Power Simulation Explorer",
                             sidebarPanel(
                               selectInput("parameter", "Parameter:", choices = c("Threshold", "Slope")),
                               numericInput("es", "Effect size (d):", value = 0.5, min = 0, max = 1, step = 0.01),
-                              sliderInput("subjects", "N subjects:", min = 5, max = 150, value = c(10, 60), step = 1),
-                              sliderInput("trials", "N trials:", min = 10, max = 150, value = c(10, 60), step = 1),
+                              sliderInput("subjects", "N subjects:", min = 15, max = 120, value = c(10, 60), step = 1),
+                              sliderInput("trials", "N trials:", min = 30, max = 90, value = c(10, 60), step = 1),
                               radioButtons(
                                 inputId = "powerlevel",
                                 label = "Desired power level:",
@@ -33,7 +33,8 @@ ui <- navbarPage("Bayesian Power Simulation Explorer",
                               )
                             ),
                             mainPanel(
-                              plotOutput("gridPlot")
+                              plotOutput("gridPlot"),
+                              p(em('Disclaimer: power lines may fall out of the displayed range. If none is visible, try changing the values'))
                             )
                           )
                  ),
@@ -43,8 +44,8 @@ ui <- navbarPage("Bayesian Power Simulation Explorer",
                           sidebarLayout(
                             sidebarPanel(
                               selectInput("parameter2", "Parameter:", choices = c("Threshold", "Slope")),
-                              sliderInput("subjects2", "N subjects:", min = 5, max = 150, value = 30, step = 1),
-                              sliderInput("trials2", "N trials:", min = 10, max = 150, value = 30, step = 1),
+                              sliderInput("subjects2", "N subjects:", min = 15, max = 120, value = 30, step = 1),
+                              sliderInput("trials2", "N trials:", min = 30, max = 90, value = 30, step = 1),
                               radioButtons(
                                 inputId = "powerlevel2",
                                 label = "Desired power level:",
@@ -63,8 +64,8 @@ ui <- navbarPage("Bayesian Power Simulation Explorer",
                           sidebarLayout(
                             sidebarPanel(
                               selectInput("parameter3", "Parameter:", choices = c("Threshold", "Slope")),
-                              numericInput("manual_subjects", "Number of subjects:", value = 30, min = 1, max = 150, step = 1),
-                              numericInput("manual_trials", "Number of trials:", value = 40, min = 1, max = 150, step = 1),
+                              numericInput("manual_subjects", "Number of subjects:", value = 30, min = 15, max = 120, step = 1),
+                              numericInput("manual_trials", "Number of trials:", value = 40, min = 30, max = 90, step = 1),
                               numericInput("manual_es", "Effect size (d):", value = 0.5, min = 0, max = 1, step = 0.01),
                               actionButton("manual_submit", "Submit")
                             ),
@@ -97,6 +98,7 @@ server <- function(input, output, session) {
                       ylim = c(input$subjects[1], input$subjects[2])) +
       scale_x_continuous(breaks = scales::pretty_breaks(n = 3)) +
       scale_y_continuous(breaks = scales::pretty_breaks(n = 3)) +
+      scale_color_manual(labels=c('Hierarchical','Simple t-test','Uncertainty prop. t-test'),values =c('#D55E00','#CC79A7','#0072B2'))+
       theme(text = element_text(size = 16))
   })
 
@@ -274,6 +276,8 @@ server <- function(input, output, session) {
       labs(x = "Effect size (d)", y = "Power",colour='Test type',fill='Test type') +
       scale_y_continuous(limits = c(0, 1)) +
       scale_x_continuous(limits = c(0, 1)) +
+      scale_color_manual(labels=c('Hierarchical','Simple t-test','Uncertainty prop. t-test'),values =c('#D55E00','#CC79A7','#0072B2'))+
+      scale_fill_manual(labels=c('Hierarchical','Simple t-test','Uncertainty prop. t-test'),values =c('#D55E00','#CC79A7','#0072B2'))+
       theme(text = element_text(size = 16))+
       geom_hline(yintercept = 0.05, linetype = 2)+
       geom_hline(yintercept = as.numeric(desired_power), linetype = 2)