ClinicalTrialsEstimation.Rproj

@@ -1,13 +0,0 @@
-Version: 1.0
-
-RestoreWorkspace: Default
-SaveWorkspace: Default
-AlwaysSaveHistory: Default
-
-EnableCodeIndexing: Yes
-UseSpacesForTab: Yes
-NumSpacesForTab: 2
-Encoding: UTF-8
-
-RnwWeave: Sweave
-LaTeX: pdfLaTeX

r-analysis/EffectsOfEnrollmentDelay.qmd

@@ -28,14 +28,6 @@ options(mc.cores = parallel::detectCores())
 
 #test installation, shouldn't get any errors
 #example(stan_model, package = "rstan", run.dontrun = TRUE)
-
-
-image_root <- "./output/EffectsOfEnrollmentDelay"
-image_dist_diff_analysis <- paste0(image_root,"/dist_diff_analysis")
-image_trial_details <-paste0(image_root,"/trials_details")
-image_parameters_by_groups <-paste0(image_root,"/betas/parameters_by_group")
-image_parameters_across_groups <-paste0(image_root,"/betas/parameter_across_groups")
-
 ```
 
 ```{r}
@@ -481,7 +473,7 @@ ggplot(data=df3, aes(x=count, fill=final_status)) +
     geom_histogram(binwidth=1) +
     ggtitle("Histogram of snapshots per trial (matched trials)") +
     xlab("Snapshots per trial")
-ggsave(paste0(image_trial_details,"/HistSnapshots.png"))
+ggsave("./Images/HistSnapshots.png")
 
 #Plot duration for terminated vs completed
 df4 <- dbGetQuery(
@@ -504,7 +496,7 @@ ggplot(data=df4, aes(x=duration,fill=overall_status)) +
     ggtitle("Histogram of trial durations") +
     xlab("duration")+
     facet_wrap(~overall_status)
-ggsave(paste0(image_trial_details,"/HistTrialDurations_Faceted.png"))
+ggsave("./Images/HistTrialDurations_Faceted.png")
 
 df5 <- dbGetQuery(
     con,
@@ -535,7 +527,7 @@ ggplot(data=df5, aes(x=duration,y=snapshot_count,color=overall_status)) +
     ggtitle("Comparison of duration, status, and snapshot_count") +
     xlab("duration") +
     ylab("snapshot count") 
-ggsave(paste0(image_trial_details,"/SnapshotsVsDurationVsTermination.png"))
+ggsave("./Images/SnapshotsVsDurationVsTermination.png")
 
 dbDisconnect(con)
 
@@ -544,14 +536,14 @@ group_trials_by_category <- as.data.frame(aggregate(category_id ~ nct_id, df, ma
 group_trials_by_category <- as.data.frame(group_trials_by_category)
 
 ggplot(data = group_trials_by_category, aes(x=category_id)) +
-    geom_bar(binwidth=1,color="black",fill="lightblue") +
+    geom_bar(binwidth=1,color="black",fill="seagreen") +
     scale_x_continuous(breaks=scales::pretty_breaks(n=22)) + 
     labs(
         title="bar chart of trial categories"
         ,x="Category ID"
         ,y="Count"
     )
-ggsave(paste0(image_trial_details,"/CategoryCounts.png"))
+ggsave("./Images/CategoryCounts.png")
     
 
 
@@ -582,25 +574,25 @@ print(fit)
 gx <- c()
 
 #grab parameters for every category with more than 8 observations
-for (i in category_count$category_id[category_count$n >= 0]) {
+for (i in category_count$category_id[category_count$n >= 8]) {
     print(i)
     
     #Print parameter distributions
     gi <- group_mcmc_areas("beta",beta_list,fit,i) #add way to filter groups
     ggsave(
-        paste0(image_parameters_by_groups,"/group_",i,"_",gi$name,".png")
+        paste0("./Images/DirectEffects/Parameters/group_",i,"_",gi$name,".png")
         ,plot=gi$plot
         )
     gx <- c(gx,gi)
 
     #Get Quantiles and means for parameters
-#    table <- xtable(gi$quantiles,
+    table <- xtable(gi$quantiles,
-#      floating=FALSE
+      floating=FALSE
-#      ,latex.environments = NULL
+      ,latex.environments = NULL
-#      ,booktabs = TRUE
+      ,booktabs = TRUE
-#      ,zap=getOption("digits")
+      ,zap=getOption("digits")
-#      )
+      )
-#    write_lines(table,paste0("./latex_output/DirectEffects/group_",gi$name,".tex"))
+    write_lines(table,paste0("./latex_output/DirectEffects/group_",gi$name,".tex"))
 }
 ```
 
@@ -615,19 +607,19 @@ for (i in c(1,2,3,9,10,11,12)) {
     #Print parameter distributions
     pi <- parameter_mcmc_areas("beta",beta_list,fit,i) #add way to filter groups
     ggsave(
-        paste0(image_parameters_across_groups,"/parameters_",i,"_",pi$name,".png")
+        paste0("./Images/DirectEffects/Parameters/parameters_",i,"_",pi$name,".png")
         ,plot=pi$plot
         )
     px <- c(px,pi)
 
     #Get Quantiles and means for parameters
-#    table <- xtable(pi$quantiles,
+    table <- xtable(pi$quantiles,
-#      floating=FALSE
+      floating=FALSE
-#      ,latex.environments = NULL
+      ,latex.environments = NULL
-#      ,booktabs = TRUE
+      ,booktabs = TRUE
-#      ,zap=getOption("digits")
+      ,zap=getOption("digits")
-#      )
+      )
-#    write_lines(table,paste0("./latex_output/DirectEffects/parameters_",i,"_",pi$name,".tex"))
+    write_lines(table,paste0("./latex_output/DirectEffects/parameters_",i,"_",pi$name,".tex"))
     
 }
 ```
@@ -639,13 +631,8 @@ Note these have 95% outer CI and 80% inner (shaded)
 
 
 ```{r}
-#get the generic and uspdc parameters
 print(px[4]$plot + px[7]$plot)
-ggsave(paste0(image_parameters_across_groups,"2+3_generic_and_uspdc.png"))
+ggsave("./Images/DirectEffects/Parameters/2+3_generic_and_uspdc.png")
-
-#get the parameters associated with duration
-px[16]$plot + px[19]$plot
-ggsave(paste0(image_parameters_across_groups,"11+12_statusREC_and_statusANR.png"))
 ```
 
 
@@ -680,7 +667,7 @@ ggplot(df_ib_p, aes(x=p_prior)) +
         ,x="Probability Domain 'p'"
         ,y="Probability Density"
     )
-ggsave(paste0(image_dist_diff_analysis,"/prior_p.png"))
+ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/prior_p.png")
 
 #p_posterior
 ggplot(df_ib_p, aes(x=p_predicted)) +
@@ -691,7 +678,7 @@ ggplot(df_ib_p, aes(x=p_predicted)) +
         ,x="Probability Domain 'p'"
         ,y="Probability Density"
     )
-ggsave(paste0(image_dist_diff_analysis,"/posterior_p.png"))
+ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/posterior_p.png")
 
 #mu_prior
 ggplot(df_ib_prior) +
@@ -702,7 +689,7 @@ ggplot(df_ib_prior) +
         ,x="Mu"
         ,y="Probability"
     )
-ggsave(paste0(image_dist_diff_analysis,"/prior_mu.png"))
+ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/prior_mu.png")
 
 #sigma_posterior
 ggplot(df_ib_prior) +
@@ -713,7 +700,7 @@ ggplot(df_ib_prior) +
         ,x="Sigma"
         ,y="Probability"
     )
-ggsave(paste0(image_dist_diff_analysis,"/prior_sigma.png"))
+ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/prior_sigma.png")
 ```
 
 
@@ -744,7 +731,7 @@ ggplot(counterfact_predicted_ib, aes(x=p_predicted_default)) +
         ,x="Probability Domain 'p'"
         ,y="Probability Density"
     )
-ggsave(paste0(image_dist_diff_analysis,"p_no_intervention.png"))
+ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/default_p_generic_intervention_base.png")
 
 ggplot(counterfact_predicted_ib, aes(x=p_predicted_intervention)) +
     geom_density() + 
@@ -754,7 +741,7 @@ ggplot(counterfact_predicted_ib, aes(x=p_predicted_intervention)) +
         ,x="Probability Domain 'p'"
         ,y="Probability Density"
     )
-ggsave(paste0(image_dist_diff_analysis,"/p_delay_intervention.png"))
+ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/default_p_generic_intervention_interv.png")
 
 ggplot(counterfact_predicted_ib, aes(x=predicted_difference)) +
     geom_density() + 
@@ -764,7 +751,7 @@ ggplot(counterfact_predicted_ib, aes(x=predicted_difference)) +
         ,x="Difference in 'p' under treatment"
         ,y="Probability Density"
     )
-ggsave(paste0(image_dist_diff_analysis,"/p_delay_intervention_distdiff_1.png"))
+ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/default_p_generic_intervention_distdiff.png")
 ```
 
 
@@ -786,6 +773,7 @@ category_names <- sapply(1:length(beta_list$groups),
 pddf_ib <- data.frame(extract(generated_ib, pars="predicted_difference")$predicted_difference) |>
     pivot_longer(X1:X168) #CHANGE_NOTE: moved from X169 to X168
 
+#TODO: Fix Category names
 pddf_ib["entry_idx"] <- as.numeric(gsub("\\D","",pddf_ib$name))
 pddf_ib["category"] <-  sapply(pddf_ib$entry_idx, function(i) df$category_id[i])
 pddf_ib["category_name"] <- sapply(
@@ -804,7 +792,8 @@ ggplot(pddf_ib, aes(x=value,)) +
         ,y = "Probability Density"
     ) + 
     geom_vline(aes(xintercept = 0), color = "skyblue", linetype="dashed") 
-ggsave(paste0(image_dist_diff_analysis,"/p_delay_intervention_distdiff_styled.png"))
+    #todo: add median, mean, 40/60 quantiles as well as 
+ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/p_delay_intervention_distdiff_styled.png")
 
 ggplot(pddf_ib, aes(x=value,)) +
     geom_density(adjust=1/5) +
@@ -823,7 +812,7 @@ ggplot(pddf_ib, aes(x=value,)) +
     ) + 
     geom_vline(aes(xintercept = 0), color = "skyblue", linetype="dashed") +
     theme(strip.text.x = element_text(size = 8))
-ggsave(paste0(image_dist_diff_analysis,"/p_delay_intervention_distdiff_by_group.png"))
+ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/p_delay_intervention_distdiff_by_group.png")
 
 ggplot(pddf_ib, aes(x=value,)) +
     geom_histogram(bins=300) +
@@ -842,7 +831,7 @@ ggplot(pddf_ib, aes(x=value,)) +
     ) + 
     geom_vline(aes(xintercept = 0), color = "skyblue", linetype="dashed") +
     theme(strip.text.x = element_text(size = 8))
-ggsave(paste0(image_dist_diff_analysis,"p_delay_intervention_histdiff_by_group.png"))
+ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/p_delay_intervention_histdiff_by_group.png")
 ```
 
 
@@ -858,8 +847,6 @@ p3 <- ggplot(pddf_ib, aes(x=value,)) +
     ) + 
     geom_vline(aes(xintercept = 0), color = "skyblue", linetype="dashed") 
 
-density_hight <- max(density(pddf_ib$value)$y)
-
 stats <- list(
  p5 = quantile(pddf_ib$value, 0.05),
  p10 = quantile(pddf_ib$value, 0.10),
@@ -870,8 +857,7 @@ stats <- list(
  p90 = quantile(pddf_ib$value, 0.90),
  p95 = quantile(pddf_ib$value, 0.95),
  max_height = max(ggplot_build(p3)$data[[1]]$count),
- y_offset = -max(ggplot_build(p3)$data[[1]]$count) * 0.05,
+ y_offset = -max(ggplot_build(p3)$data[[1]]$count) * 0.05
- y_offset_density = -density_hight * 0.15
 )
 
 p3 + 
@@ -921,65 +907,7 @@ p3 +
    x = stats$mean,
    y = stats$y_offset * 1.5
  ), aes(x = x, y = y))
-ggsave(paste0(image_dist_diff_analysis,"/p_delay_intervention_histdiff_boxplot.png"))
+ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/p_delay_intervention_histdiff_boxplot.png")
-
-
-p4 <- ggplot(pddf_ib, aes(x = value)) +
-  geom_density() +
-    labs(
-        title = "Distribution of predicted differences"
-        ,subtitle = "Intervention: Delay close of enrollment"
-        ,x = "Difference in probability due to intervention"
-        ,y = "Probability Density"
-        ,caption = "Vertical marks: 5/10/25/50/75/90/95th percentiles. Dot shows mean."
-    ) + 
-    geom_vline(aes(xintercept = 0), color = "skyblue", linetype="dashed") +
- geom_segment(data = data.frame(
-   x = c(stats$q1, stats$q3, stats$med),
-   xend = c(stats$q1, stats$q3, stats$med),
-   y = rep(stats$y_offset_density, 3), 
-   yend = rep(stats$y_offset_density * 2, 3)
- ), aes(x = x, xend = xend, y = y, yend = yend)) +
- geom_segment(data = data.frame(
-   x = rep(stats$q1, 2),
-   xend = rep(stats$q3, 2),
-   y = c(stats$y_offset_density, stats$y_offset_density * 2),
-   yend = c(stats$y_offset_density, stats$y_offset_density * 2)
- ), aes(x = x, xend = xend, y = y, yend = yend)) +
- # Inner whiskers (Q1->P10, Q3->P90)
- geom_segment(data = data.frame(
-   x = c(stats$q1, stats$q3),
-   xend = c(stats$p10, stats$p90),
-   y = rep(stats$y_offset_density * 1.5, 2),
-   yend = rep(stats$y_offset_density * 1.5, 2)
- ), aes(x = x, xend = xend, y = y, yend = yend)) +
- # Crossbars at P10/P90
- geom_segment(data = data.frame(
-   x = c(stats$p10, stats$p90),
-   xend = c(stats$p10, stats$p90),
-   y = stats$y_offset_density * 1.3,
-   yend = stats$y_offset_density * 1.7
- ), aes(x = x, xend = xend, y = y, yend = yend)) +
- # Outer whiskers (P10->P5, P90->P95)
- geom_segment(data = data.frame(
-   x = c(stats$p10, stats$p90),
-   xend = c(stats$p5, stats$p95),
-   y = rep(stats$y_offset_density * 1.5, 2),
-   yend = rep(stats$y_offset_density * 1.5, 2)
- ), aes(x = x, xend = xend, y = y, yend = yend)) +
- # Crossbars at P5/P95
- geom_segment(data = data.frame(
-   x = c(stats$p5, stats$p95),
-   xend = c(stats$p5, stats$p95),
-   y = stats$y_offset_density * 1.3,
-   yend = stats$y_offset_density * 1.7
- ), aes(x = x, xend = xend, y = y, yend = yend)) +
- # Mean dot
- geom_point(data = data.frame(
-   x = stats$mean,
-   y = stats$y_offset_density * 1.5
- ), aes(x = x, y = y))
-ggsave(paste0(image_dist_diff_analysis,"/p_delay_intervention_distdiff_boxplot.png"))
 ```
 
 ```{r}
@@ -992,7 +920,7 @@ ggsave(paste0(image_dist_diff_analysis,"/p_delay_intervention_distdiff_boxplot.p
         y = "Cumulative Proportion"
     ) 
 
-ggsave(paste0(image_dist_diff_analysis,"/p_delay_intervention_cumulative_distdiff.png"))
+ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/p_delay_intervention_cumulative_distdiff.png")
 ```
 
 
@@ -1058,10 +986,8 @@ The simulation above shows that this results in a percentage-point increase of a
 ```{r}
 #| eval: true
 #trace plots
-image_diagnostics <- paste0(image_root,"/diagnostics")
-
 plot(fit, pars=c("mu"), plotfun="trace")
-ggsave(paste0(image_diagnostics,"/trace_plot_mu.png"))
+ggsave("./EffectsOfEnrollmentDelay/Images/diagnostics/trace_plot_mu.png")
 
 
 for (i in 1:3) {
@@ -1079,7 +1005,7 @@ for (i in 1:3) {
              scale_colour_ghibli_d("KikiMedium")
     )
     mu_range <- paste0(4*i-3,"-",4*i)
-    filename <- paste0(image_diagnostics,"/trace_rank_plot_mu_",mu_range,".png")
+    filename <- paste0("./EffectsOfEnrollmentDelay/Images/diagnostics/trace_rank_plot_mu_",mu_range,".png")
     ggsave(filename)
 }
 ```
@@ -1087,7 +1013,7 @@ for (i in 1:3) {
 ```{r}
 #| eval: true
 plot(fit, pars=c("sigma"), plotfun="trace")
-ggsave(paste0(image_diagnostics,"/traceplot_sigma.png"))
+ggsave("./EffectsOfEnrollmentDelay/Images/diagnostics/traceplot_sigma.png")
 
 for (i in 1:3) {
     print(
@@ -1104,7 +1030,7 @@ for (i in 1:3) {
              scale_colour_ghibli_d("KikiMedium")
     )
     sigma_range <- paste0(4*i-3,"-",4*i)
-    filename <- paste0(image_diagnostics,"/trace_rank_plot_sigma_",sigma_range,".png")
+    filename <- paste0("./EffectsOfEnrollmentDelay/Images/diagnostics/trace_rank_plot_sigma_",sigma_range,".png")
     ggsave(filename)
 }
 ```
@@ -1123,7 +1049,7 @@ posterior <- as.array(fit)
 color_scheme_set("darkgray")
 div_style <- parcoord_style_np(div_color = "green", div_size = 0.05, div_alpha = 0.4)
 mcmc_parcoord(posterior, regex_pars = "mu", np=nuts_prmts, np_style = div_style, alpha = 0.05)
-ggsave(paste0(image_diagnostics,"/parcoord_mu.png"))
+ggsave("./EffectsOfEnrollmentDelay/Images/diagnostics/parcoord_mu.png")
 ```
 
 ```{r}
@@ -1142,7 +1068,7 @@ for (i in 1:3) {
         )
     )
     mu_range <- paste0(4*i-3,"-",4*i)
-    filename <- paste0(image_diagnostics,"/correlation_plot_mu_",mu_range,".png")
+    filename <- paste0("./EffectsOfEnrollmentDelay/Images/diagnostics/correlation_plot_mu_",mu_range,".png")
     ggsave(filename)
 }
 
@@ -1153,7 +1079,7 @@ for (i in 1:3) {
 ```{r}
 #| eval: true
 mcmc_parcoord(posterior,regex_pars = "sigma", np=nuts_prmts, alpha=0.05)
-ggsave(paste0(image_diagnostics,"/parcoord_sigma.png"))
+ggsave("./EffectsOfEnrollmentDelay/Images/diagnostics/parcoord_sigma.png")
 ```
 
 ```{r}
@@ -1173,7 +1099,7 @@ for (i in 1:3) {
         )
     )
     sigma_range <- paste0(4*i-3,"-",4*i)
-    filename <- paste0(image_diagnostics,"/correlation_plot_sigma_",sigma_range,".png")
+    filename <- paste0("./EffectsOfEnrollmentDelay/Images/diagnostics/correlation_plot_sigma_",sigma_range,".png")
     ggsave(filename)
 }
 ```
@@ -1197,7 +1123,7 @@ for (i in 1:3) {
     )
     
     beta_range <- paste0("k_",k,"_i_",4*i-3,"-",4*i)
-    filename <- paste0(image_diagnostics,"/correlation_plot_beta_",beta_range,".png")
+    filename <- paste0("./EffectsOfEnrollmentDelay/Images/diagnostics/correlation_plot_beta_",beta_range,".png")
     ggsave(filename)
 }}
 ```
@@ -1207,4 +1133,11 @@ for (i in 1:3) {
 
 
 
-  
+# TODO
+- [ ] Double check data flow. (Write summary of this in human readable form)
+    - Is it the data we want from the database 
+        - Training
+        - Counterfactual Evaluation
+            - choose a single snapshot per trial.
+    - Is the model in STAN well specified.
+- [ ] work on LOO validation of model

r-analysis/Hierarchal_Logistic.stan

@@ -26,8 +26,6 @@ data {
   array[Nx] int<lower=1, upper=L> llx;//vec of categories
   array[Nx] row_vector[D] counterfact_x_tilde; // Posterior Prediction intervention
   array[Nx] row_vector[D] counterfact_x; // Posterior Prediction intervention
-  //the two statuses to catch relative difference
-  array[2] int status_indexes; //the two status indexes to compare (order is x[1] - x[2])
 }
 parameters {
   array[D] real mu;
@@ -60,13 +58,6 @@ generated quantities {
     real p_predicted_intervention[Nx]; //predicted p_values
     real predicted_difference[Nx]; //difference in predicted values
 
-    //collect array of relative status differences between 
-    array[L] real status_diff;
-    //GENERATE RELATIVE DIFFERENCES BETWEEN STATUSES
-    for (l in 1:L) {
-      status_diff[l] = beta[l,status_indexes[1]] - beta[l,status_indexes[2]];
-    }
-
     //GENERATE PRIOR PREDICTIONS
     {
         vector[D] beta_prior[L];//local var

r-analysis/Hierarchal_Logistic_with_status_diff.stan

@@ -1,102 +0,0 @@
-//
-// This Stan program defines a simple model, with a
-// vector of values 'y' modeled as normally distributed
-// with mean 'mu' and standard deviation 'sigma'.
-//
-// Learn more about model development with Stan at:
-//
-//    http://mc-stan.org/users/interfaces/rstan.html
-//    https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started
-//
-
-// The input data is a vector 'y' of length 'N'.
-data {
-  int<lower=1> D; // number of features
-  int<lower=1> N; // number of observations
-  int<lower=1> L; // number of categories
-  array[N] int<lower=0, upper=1> y; // vector of dependent variables
-  array[N] int<lower=1, upper=L> ll; // vector of categories
-  array[N] row_vector[D] x; // independent variables
-  real mu_mean; //hyperprior
-  real<lower=0> mu_stdev; //hyperprior
-  real sigma_shape; //hyperprior
-  real sigma_rate; //hyperprior
-  //counterfactuals
-  int<lower=0> Nx;
-  array[Nx] int<lower=1, upper=L> llx;//vec of categories
-  array[Nx] row_vector[D] counterfact_x_tilde; // Posterior Prediction intervention
-  array[Nx] row_vector[D] counterfact_x; // Posterior Prediction intervention
-  //the two statuses to catch relative difference
-  array[2] int status_indexes; //the two status indexes to compare (order is x[1] - x[2])
-}
-parameters {
-  array[D] real mu;
-  array[D] real<lower=0> sigma;
-  array[L] vector[D] beta;
-}
-model {
-  sigma ~ gamma(sigma_shape,sigma_rate); //hyperprior for stdev: shape and inverse scale
-  mu ~ normal(mu_mean, mu_stdev); //hyperprior for mean //TODO: convert to mvnormal
-  for (l in 1:L) {
-    beta[l] ~ normal(mu, sigma);
-  }
-  {
-    vector[N] x_beta_ll;
-    for (n in 1:N) {
-      x_beta_ll[n] = x[n] * beta[ll[n]];
-    }
-    y ~ bernoulli_logit(x_beta_ll);
-  }
-}
-generated quantities {
-    //SETUP PRIOR PREDICTION
-    //preallocate
-    real mu_prior[D];
-    real sigma_prior[D];
-    real p_prior[N]; // what I have priors about
-    real p_predicted[N]; //predicted p_values
-    //intervention
-    real p_predicted_default[Nx]; //predicted p_values
-    real p_predicted_intervention[Nx]; //predicted p_values
-    real predicted_difference[Nx]; //difference in predicted values
-    
-    //collect array of relative status differences between 
-    array[L] real status_diff;
-    //GENERATE RELATIVE DIFFERENCES BETWEEN STATUSES
-    for (l in 1:L) {
-      status_diff[l] = beta[l,status_indexes[1]] - beta[l,status_indexes[2]];
-    }
-
-    //GENERATE PRIOR PREDICTIONS
-    {
-        vector[D] beta_prior[L];//local var
-        //sample parameters
-        for (d in 1:D) {
-            mu_prior[d] = normal_rng(mu_mean,mu_stdev);
-            sigma_prior[d] = gamma_rng(sigma_shape,sigma_rate);
-        }
-        for (l in 1:L) {
-            for (d in 1:D) {
-                beta_prior[l,d] = normal_rng(mu_prior[d],sigma_prior[d]);
-            }
-        }
-        //generate probabilities
-        vector[D] b_prior[N];//local var
-        for (n in 1:N){
-            b_prior[n] = beta_prior[ll[n]];
-            p_prior[n] = inv_logit(  x[n] * b_prior[n] );
-        }
-    }
-    //GENERATE POSTERIOR PREDICTIONS
-    for (n in 1:N) {
-        p_predicted[n] = inv_logit( x[n] * beta[ll[n]] );
-    }
-    //GENERATE POSTERIOR  DISTRIBUTION OF DIFFERENCES
-    for (n in 1:Nx) {
-        p_predicted_default[n] = inv_logit(  counterfact_x[n] * beta[llx[n]] );
-        p_predicted_intervention[n] = inv_logit(  counterfact_x_tilde[n] * beta[llx[n]] ); //intervention
-
-        //intervention - base case
-        predicted_difference[n] = p_predicted_intervention[n] - p_predicted_default[n];
-    }
-}