Merge branch 'exploring_fixed_effects'

ClinicalTrialsEstimation.Rproj

@@ -0,0 +1,13 @@
+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,6 +28,14 @@ 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}
@@ -473,7 +481,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("./Images/HistSnapshots.png")
+ggsave(paste0(image_trial_details,"/HistSnapshots.png"))
 
 #Plot duration for terminated vs completed
 df4 <- dbGetQuery(
@@ -496,7 +504,7 @@ ggplot(data=df4, aes(x=duration,fill=overall_status)) +
     ggtitle("Histogram of trial durations") +
     xlab("duration")+
     facet_wrap(~overall_status)
-ggsave("./Images/HistTrialDurations_Faceted.png")
+ggsave(paste0(image_trial_details,"/HistTrialDurations_Faceted.png"))
 
 df5 <- dbGetQuery(
     con,
@@ -516,7 +524,7 @@ df5 <- dbGetQuery(
     )
     select a.nct_id, a.overall_status, a.duration,b.snapshot_count
     from cte1 as a
-        join cte2 as b
+        join cte2 as b  
             on a.nct_id=b.nct_id
     ;"
     )
@@ -527,7 +535,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("./Images/SnapshotsVsDurationVsTermination.png")
+ggsave(paste0(image_trial_details,"/SnapshotsVsDurationVsTermination.png"))
 
 dbDisconnect(con)
 
@@ -536,14 +544,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="seagreen") +
+    geom_bar(binwidth=1,color="black",fill="lightblue") +
     scale_x_continuous(breaks=scales::pretty_breaks(n=22)) + 
     labs(
         title="bar chart of trial categories"
         ,x="Category ID"
         ,y="Count"
     )
-ggsave("./Images/CategoryCounts.png")
+ggsave(paste0(image_trial_details,"/CategoryCounts.png"))
     
 
 
@@ -574,25 +582,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 >= 8]) {
+for (i in category_count$category_id[category_count$n >= 0]) {
     print(i)
     
     #Print parameter distributions
     gi <- group_mcmc_areas("beta",beta_list,fit,i) #add way to filter groups
     ggsave(
-        paste0("./Images/DirectEffects/Parameters/group_",i,"_",gi$name,".png")
+        paste0(image_parameters_by_groups,"/group_",i,"_",gi$name,".png")
         ,plot=gi$plot
         )
     gx <- c(gx,gi)
 
     #Get Quantiles and means for parameters
-    table <- xtable(gi$quantiles,
-      floating=FALSE
-      ,latex.environments = NULL
-      ,booktabs = TRUE
-      ,zap=getOption("digits")
-      )
-    write_lines(table,paste0("./latex_output/DirectEffects/group_",gi$name,".tex"))
+#    table <- xtable(gi$quantiles,
+#      floating=FALSE
+#      ,latex.environments = NULL
+#      ,booktabs = TRUE
+#      ,zap=getOption("digits")
+#      )
+#    write_lines(table,paste0("./latex_output/DirectEffects/group_",gi$name,".tex"))
 }
 ```
 
@@ -607,19 +615,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("./Images/DirectEffects/Parameters/parameters_",i,"_",pi$name,".png")
+        paste0(image_parameters_across_groups,"/parameters_",i,"_",pi$name,".png")
         ,plot=pi$plot
         )
     px <- c(px,pi)
 
     #Get Quantiles and means for parameters
-    table <- xtable(pi$quantiles,
-      floating=FALSE
-      ,latex.environments = NULL
-      ,booktabs = TRUE
-      ,zap=getOption("digits")
-      )
-    write_lines(table,paste0("./latex_output/DirectEffects/parameters_",i,"_",pi$name,".tex"))
+#    table <- xtable(pi$quantiles,
+#      floating=FALSE
+#      ,latex.environments = NULL
+#      ,booktabs = TRUE
+#      ,zap=getOption("digits")
+#      )
+#    write_lines(table,paste0("./latex_output/DirectEffects/parameters_",i,"_",pi$name,".tex"))
     
 }
 ```
@@ -631,8 +639,13 @@ 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("./Images/DirectEffects/Parameters/2+3_generic_and_uspdc.png")
+ggsave(paste0(image_parameters_across_groups,"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"))
 ```
 
 
@@ -667,7 +680,7 @@ ggplot(df_ib_p, aes(x=p_prior)) +
         ,x="Probability Domain 'p'"
         ,y="Probability Density"
     )
-ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/prior_p.png")
+ggsave(paste0(image_dist_diff_analysis,"/prior_p.png"))
 
 #p_posterior
 ggplot(df_ib_p, aes(x=p_predicted)) +
@@ -678,7 +691,7 @@ ggplot(df_ib_p, aes(x=p_predicted)) +
         ,x="Probability Domain 'p'"
         ,y="Probability Density"
     )
-ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/posterior_p.png")
+ggsave(paste0(image_dist_diff_analysis,"/posterior_p.png"))
 
 #mu_prior
 ggplot(df_ib_prior) +
@@ -689,7 +702,7 @@ ggplot(df_ib_prior) +
         ,x="Mu"
         ,y="Probability"
     )
-ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/prior_mu.png")
+ggsave(paste0(image_dist_diff_analysis,"/prior_mu.png"))
 
 #sigma_posterior
 ggplot(df_ib_prior) +
@@ -700,7 +713,7 @@ ggplot(df_ib_prior) +
         ,x="Sigma"
         ,y="Probability"
     )
-ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/prior_sigma.png")
+ggsave(paste0(image_dist_diff_analysis,"/prior_sigma.png"))
 ```
 
 
@@ -731,7 +744,7 @@ ggplot(counterfact_predicted_ib, aes(x=p_predicted_default)) +
         ,x="Probability Domain 'p'"
         ,y="Probability Density"
     )
-ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/default_p_generic_intervention_base.png")
+ggsave(paste0(image_dist_diff_analysis,"p_no_intervention.png"))
 
 ggplot(counterfact_predicted_ib, aes(x=p_predicted_intervention)) +
     geom_density() + 
@@ -741,7 +754,7 @@ ggplot(counterfact_predicted_ib, aes(x=p_predicted_intervention)) +
         ,x="Probability Domain 'p'"
         ,y="Probability Density"
     )
-ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/default_p_generic_intervention_interv.png")
+ggsave(paste0(image_dist_diff_analysis,"/p_delay_intervention.png"))
 
 ggplot(counterfact_predicted_ib, aes(x=predicted_difference)) +
     geom_density() + 
@@ -751,7 +764,7 @@ ggplot(counterfact_predicted_ib, aes(x=predicted_difference)) +
         ,x="Difference in 'p' under treatment"
         ,y="Probability Density"
     )
-ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/default_p_generic_intervention_distdiff.png")
+ggsave(paste0(image_dist_diff_analysis,"/p_delay_intervention_distdiff_1.png"))
 ```
 
 
@@ -773,7 +786,6 @@ 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(
@@ -792,8 +804,7 @@ ggplot(pddf_ib, aes(x=value,)) +
         ,y = "Probability Density"
     ) + 
     geom_vline(aes(xintercept = 0), color = "skyblue", linetype="dashed") 
-    #todo: add median, mean, 40/60 quantiles as well as 
-ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/p_delay_intervention_distdiff_styled.png")
+ggsave(paste0(image_dist_diff_analysis,"/p_delay_intervention_distdiff_styled.png"))
 
 ggplot(pddf_ib, aes(x=value,)) +
     geom_density(adjust=1/5) +
@@ -812,7 +823,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("./EffectsOfEnrollmentDelay/Images/DirectEffects/p_delay_intervention_distdiff_by_group.png")
+ggsave(paste0(image_dist_diff_analysis,"/p_delay_intervention_distdiff_by_group.png"))
 
 ggplot(pddf_ib, aes(x=value,)) +
     geom_histogram(bins=300) +
@@ -831,7 +842,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("./EffectsOfEnrollmentDelay/Images/DirectEffects/p_delay_intervention_histdiff_by_group.png")
+ggsave(paste0(image_dist_diff_analysis,"p_delay_intervention_histdiff_by_group.png"))
 ```
 
 
@@ -847,6 +858,8 @@ 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),
@@ -857,7 +870,8 @@ 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 + 
@@ -907,7 +921,65 @@ p3 +
    x = stats$mean,
    y = stats$y_offset * 1.5
  ), aes(x = x, y = y))
-ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/p_delay_intervention_histdiff_boxplot.png")
+ggsave(paste0(image_dist_diff_analysis,"/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}
@@ -920,7 +992,7 @@ ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/p_delay_intervention_his
         y = "Cumulative Proportion"
     ) 
 
-ggsave("./EffectsOfEnrollmentDelay/Images/DirectEffects/p_delay_intervention_cumulative_distdiff.png")
+ggsave(paste0(image_dist_diff_analysis,"/p_delay_intervention_cumulative_distdiff.png"))
 ```
 
 
@@ -986,8 +1058,10 @@ 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("./EffectsOfEnrollmentDelay/Images/diagnostics/trace_plot_mu.png")
+ggsave(paste0(image_diagnostics,"/trace_plot_mu.png"))
 
 
 for (i in 1:3) {
@@ -1005,7 +1079,7 @@ for (i in 1:3) {
              scale_colour_ghibli_d("KikiMedium")
     )
     mu_range <- paste0(4*i-3,"-",4*i)
-    filename <- paste0("./EffectsOfEnrollmentDelay/Images/diagnostics/trace_rank_plot_mu_",mu_range,".png")
+    filename <- paste0(image_diagnostics,"/trace_rank_plot_mu_",mu_range,".png")
     ggsave(filename)
 }
 ```
@@ -1013,7 +1087,7 @@ for (i in 1:3) {
 ```{r}
 #| eval: true
 plot(fit, pars=c("sigma"), plotfun="trace")
-ggsave("./EffectsOfEnrollmentDelay/Images/diagnostics/traceplot_sigma.png")
+ggsave(paste0(image_diagnostics,"/traceplot_sigma.png"))
 
 for (i in 1:3) {
     print(
@@ -1030,7 +1104,7 @@ for (i in 1:3) {
              scale_colour_ghibli_d("KikiMedium")
     )
     sigma_range <- paste0(4*i-3,"-",4*i)
-    filename <- paste0("./EffectsOfEnrollmentDelay/Images/diagnostics/trace_rank_plot_sigma_",sigma_range,".png")
+    filename <- paste0(image_diagnostics,"/trace_rank_plot_sigma_",sigma_range,".png")
     ggsave(filename)
 }
 ```
@@ -1049,7 +1123,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("./EffectsOfEnrollmentDelay/Images/diagnostics/parcoord_mu.png")
+ggsave(paste0(image_diagnostics,"/parcoord_mu.png"))
 ```
 
 ```{r}
@@ -1068,7 +1142,7 @@ for (i in 1:3) {
         )
     )
     mu_range <- paste0(4*i-3,"-",4*i)
-    filename <- paste0("./EffectsOfEnrollmentDelay/Images/diagnostics/correlation_plot_mu_",mu_range,".png")
+    filename <- paste0(image_diagnostics,"/correlation_plot_mu_",mu_range,".png")
     ggsave(filename)
 }
 
@@ -1079,7 +1153,7 @@ for (i in 1:3) {
 ```{r}
 #| eval: true
 mcmc_parcoord(posterior,regex_pars = "sigma", np=nuts_prmts, alpha=0.05)
-ggsave("./EffectsOfEnrollmentDelay/Images/diagnostics/parcoord_sigma.png")
+ggsave(paste0(image_diagnostics,"/parcoord_sigma.png"))
 ```
 
 ```{r}
@@ -1099,7 +1173,7 @@ for (i in 1:3) {
         )
     )
     sigma_range <- paste0(4*i-3,"-",4*i)
-    filename <- paste0("./EffectsOfEnrollmentDelay/Images/diagnostics/correlation_plot_sigma_",sigma_range,".png")
+    filename <- paste0(image_diagnostics,"/correlation_plot_sigma_",sigma_range,".png")
     ggsave(filename)
 }
 ```
@@ -1123,7 +1197,7 @@ for (i in 1:3) {
     )
     
     beta_range <- paste0("k_",k,"_i_",4*i-3,"-",4*i)
-    filename <- paste0("./EffectsOfEnrollmentDelay/Images/diagnostics/correlation_plot_beta_",beta_range,".png")
+    filename <- paste0(image_diagnostics,"/correlation_plot_beta_",beta_range,".png")
     ggsave(filename)
 }}
 ```
@@ -1133,11 +1207,4 @@ 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,6 +26,8 @@ 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;
@@ -57,6 +59,13 @@ generated quantities {
     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
     {

r-analysis/Hierarchal_Logistic_with_status_diff.stan

@@ -0,0 +1,102 @@
+//
+// 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];
+    }
+}