#Workshop Frequencies
No_Workshops <- data$No_Workshops
## Frequency table including unspecified cases
WS_freq_with_na <- table(No_Workshops, useNA = "ifany")
## Proportions excluding unspecified cases
WS_prop_excl_na <- prop.table(table(No_Workshops, useNA = "no"))
## Combined dataframe
WS_Overview <- data.frame(
  Value = names(WS_freq_with_na),
  Frequency = as.vector(WS_freq_with_na),
  Proportion = c(signif(as.vector(WS_prop_excl_na), 2), rep(NA, length(WS_freq_with_na) - length(WS_prop_excl_na)))
)

datatable(WS_Overview, caption = 'Workshop Frequencies', options = list(pageLength = 10, autoWidth=TRUE))

##Mean number of workshops
No_WS_mean <- mean(data$No_Workshops, na.rm= TRUE)
No_WS_mean

##Workshop frequencies histogram
Workshops_Histogram <- ggplot(data, aes(x = No_Workshops)) +
  geom_histogram(binwidth = 2, color = "black", fill = "grey") +
  labs(
    title = "Histogram of Workshop Frequency",
    x = "Number of Workshops",
    y = "Frequency"
  ) +
  theme_minimal(base_family = "Times New Roman") +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),  # Centered, bold
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )

Workshops_Histogram

##Save as image for importing into thesis chapter (hashtagged to prevent automatic saving)
#ggsave("Workshop frequency histogram.png", width = 6, height = 4, dpi = 300)

#Stakeholder Group Size Frequencies
No_Stakeholders <- data$No_Stakeholders
## Frequency table including unspecified cases
SH_freq_with_na <- table(No_Stakeholders, useNA = "ifany")
## Proportions excluding unspecified cases
SH_prop_excl_na <- prop.table(table(No_Stakeholders, useNA = "no"))
## Combined dataframe
SH_Overview <- data.frame(
  Value = names(SH_freq_with_na),
  Frequency = as.vector(SH_freq_with_na),
  Proportion = c(signif(as.vector(SH_prop_excl_na), 2), rep(NA, length(SH_freq_with_na) - length(SH_prop_excl_na)))
)


datatable(SH_Overview, caption = 'Stakeholder Group Sizes', options = list(pageLength = 10, autoWidth=TRUE))

data$No_Stakeholders <- as.numeric(data$No_Stakeholders)
No_SH_mean <- mean(data$No_Stakeholders, na.rm= TRUE)
No_SH_mean

##Stakeholder group size histogram
Stakeholder_Group_Size_Histogram <- ggplot(data, aes(x = No_Stakeholders)) +
  geom_histogram(binwidth = 2, color = "black", fill = "grey") +
    scale_y_continuous(breaks = seq(0, 10, by = 2)) +
  labs(
    title = "Histogram of Stakeholder Group Size",
    x = "Number of Stakeholders Per Workshop",
    y = "Frequency"
  ) +
  theme_minimal(base_family = "Times New Roman") +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),  # Centered, bold
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )

Stakeholder_Group_Size_Histogram

##Save as image for importing into thesis chapter (hashtagged to prevent automatic saving)
##ggsave("Stakeholder Group Size histogram.png", width = 6, height = 4, dpi = 300)

# Frequency of single vs. multiple facilitation methods
Facilitation_Method_Total <- table(data$No_Method_Categories)
Facilitation_Method_Total 

facilitation_method_counts <- sapply(
  c("Facilitated_PS", "Facilitated_DA", "Facilitated_SD"),
  function(var) sum(data[[var]] == 1, na.rm = TRUE)
)

## Convert to data frame
facilitation_methods_table <- data.frame(
  Variable = names(facilitation_method_counts),
  Count_1 = as.vector(facilitation_method_counts)
)

datatable(facilitation_methods_table, caption = 'Faciltation Methods', options = list(pageLength = 10, autoWidth=TRUE))

## Generating dataframe for plot
fac_plot_data <- data[, c("Facilitated_PS", "Facilitated_DA", "Facilitated_SD")]

## Replace NAs with 0 (method not used)
fac_plot_data[is.na(fac_plot_data)] <- 0
fac_plot_data <- as.data.frame(fac_plot_data)

#Producing plot
upset(
  fac_plot_data,
  sets = c("Facilitated_PS", "Facilitated_DA", "Facilitated_SD"),
  order.by = "freq",
  main.bar.color = "grey30",
  sets.bar.color = "grey60",
  text.scale = c(1.6, 1.2, 1.2, 1.2, 1.4, 1.2),
  point.size = 4,
  line.size = 0.6,
  keep.order = TRUE,
  number.angles = 0,
  sets.x.label = "Number of Cases per Method",
  mainbar.y.label = "Cases per Combination"
)

##Upset has known issues with exporting graphs which prevents font changes. Manual export required. 

#Create list of all variables starting with "Evaluation_" and set NAs to 0s (Method not used)
eval_vars <- grep("^Evaluation_", names(data), value = TRUE)
eval_matrix <- data[eval_vars]
eval_matrix[is.na(eval_matrix)] <- 0 

#Count number of facilitation method categories used for each case. Then create variables for whether methods were used alone (method_counts= 1) or in combination (method_counts= 2 or 3)
method_counts <- rowSums(eval_matrix)
evaluation_method_summary_table <- sapply(eval_vars, function(var) {
  used <- eval_matrix[[var]] == 1
  only_used <- used & (method_counts == 1)
  with_others <- used & (method_counts > 1)
  
  #Create dataframe showing ns and %s of cases using each method alone and in combination with others 
  total_used <- sum(used)
  only_count <- sum(only_used)
  with_others_count <- sum(with_others)
  proportion_used <- total_used / nrow(data)
  
  c(
    Only_Method = only_count,
    With_Others = with_others_count,
    Total_Used = total_used,
    Proportion_Used = round(proportion_used, 3)
  )
})
evaluation_method_summary_table <- as.data.frame(t(evaluation_method_summary_table))
evaluation_method_summary_table$Variable <- rownames(evaluation_method_summary_table)
rownames(evaluation_method_summary_table) <- NULL
evaluation_method_summary_table <- evaluation_method_summary_table[, c("Variable", "Only_Method", "With_Others", "Total_Used", "Proportion_Used")]


datatable(evaluation_method_summary_table, caption = 'Evaluation Methods', options = list(pageLength = 10, autoWidth=TRUE))

# Creating graph
#Regenerating plot data (to ensure compatible structure) that only retains configurations with more than 3 cases
eval_plot_data <- data %>% select(starts_with("Evaluation_"))
eval_plot_data[is.na(eval_plot_data)] <- 0
eval_plot_data <- as.data.frame(eval_plot_data)
colnames(eval_plot_data)

# Create a new column with a combination key (e.g., "10110")
eval_plot_data$combo <- apply(eval_plot_data, 1, paste0, collapse = "")

# Count combinations
combo_freq <- eval_plot_data %>%
  count(combo) %>%
  filter(n >= 3)

# Keep only rows with a combo that appears at least 3 times
eval_plot_filtered <- eval_plot_data %>%
  filter(combo %in% combo_freq$combo) %>%
  select(-combo)  # remove combo column if you want


upset(
  eval_plot_filtered,
  sets = colnames(eval_plot_filtered),
  order.by = "freq",
  keep.order = TRUE,
  main.bar.color = "grey30",
  sets.bar.color = "grey60",
  text.scale = c(1.6, 1.2, 1.2, 1.2, 1.4, 1.2),
  point.size = 3,
  line.size = 1,
  sets.x.label = "Cases using each method",
  mainbar.y.label = "Method combinations",
)

# Create dataframe for qca
qca_IN <- data.frame(
  PJ = data$CH_PJ_rescale,
  AP = data$CH_SH_Engagement,
  DIV= data$CH_Diverse_SH_Recruitment,
  IN= data$Outcome_Info_Share
)

# Recalibrate PJ and DIV into binary variables due to low ns. 
  ## Stricter recalibration of DIV retains only cases of high diversity as successful outcome. 
  ## More lenient callibration of PJ is necessary here due to low frequency of cases with PJ > 1. 

qca_IN <- qca_IN %>%
  mutate(
    PJ = case_when(
      PJ <= 0.67 ~ 0,
      PJ == 1 ~ 1,
      PJ >= 1.33 ~ 1
    ),
    DIV = case_when(
      DIV == -1 ~ 0,
      DIV == 0 ~ 0,
      DIV == 1 ~ 1
    )
  )

#Remove NAs
qca_IN_exclusions <- qca_IN%>% 
  drop_na(AP, PJ, IN, DIV)

#Generate truth table
tt_IN <- truthTable(qca_IN_exclusions, outcome = "IN", conditions = c("AP", "PJ", "DIV"), type = "mv", incl.cut = 0.8, n.cut= 3, show.cases = TRUE)
tt_IN

solution_IN <- minimize(tt_IN, details= TRUE)
solution_IN

# Generate dataframe for qca
qca_AP_PJ_IN <- data.frame(
  PJ = data$CH_PJ_rescale,
  IN= data$Outcome_Info_Share,
  AP = data$CH_SH_Engagement
)

#Recalibrate PJ. 
  ## More lenient, binary calibration required due to low frequency of cases where PJ > 1.
qca_AP_PJ_IN  <- qca_AP_PJ_IN  %>%
  mutate(
    PJ = case_when(
      PJ <= 0.67 ~ 0,
      PJ == 1 ~ 1,
      PJ >= 1.33 ~ 1
    )
  )

# Drop nas
qca_AP_PJ_IN_exclusions <- qca_AP_PJ_IN%>% 
  drop_na(PJ, IN, AP)

#Generate truth table
  ##Note: n.cut not used here due to low variance of IN
tt_AP_PJ_IN <- truthTable(qca_AP_PJ_IN_exclusions, outcome = "IN", conditions = c("AP", "PJ"), type = "mv", incl.cut = 0.8, show.cases = TRUE)
tt_AP_PJ_IN

solution_AP_PJ_IN <- minimize(tt_AP_PJ_IN, details= TRUE)
solution_AP_PJ_IN 

# Create dataframe for qca
qca_PE_DIV_IN <- data.frame(
  PE = data$CH_PE_rescale,
  DIV = data$CH_Diverse_SH_Recruitment,
  IN= data$Outcome_Info_Share
)

#Recalibrate PE and DIV.
  ## Again, strict recalibration of DIV means only cases of high diversity are coded as 1.
qca_PE_DIV_IN <- qca_PE_DIV_IN %>%
    mutate(
    PE = case_when(
      PE <= 0.5 ~ 0,
      PE == 1 ~ 1,
      PE >= 1.5 ~ 2
    ),
    DIV = case_when(
      DIV == -1 ~ 0,
      DIV == 0 ~ 0,
      DIV == 1 ~ 1
    )
  )

#Remove NAs
qca_PE_DIV_IN_exclusions <- qca_PE_DIV_IN%>% 
  drop_na(PE, DIV, IN)

#Generate truth table
  ##Note: n.cut not used here due to low frequency of IN=0. 
tt_PE_DIV_IN <- truthTable(qca_PE_DIV_IN_exclusions, outcome = "IN", conditions = c("PE", "DIV"), incl.cut = 0.8, show.cases = TRUE)
tt_PE_DIV_IN

solution_PE_DIV_IN <- minimize(tt_PE_DIV_IN, details= TRUE)
solution_PE_DIV_IN 

# Create dataframe for qca
qca_CMB_CLEARN_UND <- data.frame(
  CMB = data$CH_Int_Model_Dev,
  UND = data$Outcome_Shared_Understanding,
  CLEARN = data$Outcome_Learning_Collective
)

#Remove NAs
qca_CMB_CLEARN_UND_exclusions <- qca_CMB_CLEARN_UND%>% 
  drop_na(CMB, UND, CLEARN)

#Recalibrate MB
qca_CMB_CLEARN_UND_exclusions <- qca_CMB_CLEARN_UND_exclusions %>%
    mutate(
    CMB = case_when(
      CMB == 1 ~ 2,
      CMB == 0 ~ 1,
      CMB == -1 ~ 0
    )
 )

##Generate truth table
##Note: n.cut not used here due to low variance of shared understanding
tt_CMB_CLEARN_UND <- truthTable(qca_CMB_CLEARN_UND_exclusions,
  outcome = "UND",
  conditions = c("CMB", "CLEARN"),
  incl.cut = 0.95,
  show.cases = TRUE
)

tt_CMB_CLEARN_UND

solution_CMB_CLEARN_UND <- minimize (tt_CMB_CLEARN_UND, details= TRUE)
solution_CMB_CLEARN_UND

# generate dataframe for qca
qca_CMB_CLEARN <- data.frame(
  CLEARN = data$Outcome_Learning_Collective,
  CMB = data$CH_Int_Model_Dev
)

#Remove NAs
qca_CMB_CLEARN_exclusions <- qca_CMB_CLEARN%>% 
  drop_na(CLEARN, CMB)

#Recalibrate model building variable
qca_CMB_CLEARN_exclusions <- qca_CMB_CLEARN_exclusions %>%
    mutate(
    CMB = case_when(
      CMB == 1 ~ 2,
      CMB == 0 ~ 1,
      CMB == -1 ~ 0
    )
 )

#Generate truth table. Low variance of collective learning means n.cut not used and analysis limited to descriptive only.
tt_CMB_CLEARN <- truthTable(
  qca_CMB_CLEARN_exclusions,
  outcome = "CLEARN",
  conditions = c("CMB"),
  incl.cut = 0.95,
  show.cases = TRUE
)
tt_CMB_CLEARN

# Create qca df. 
qca_COM_IMP <- data.frame(
  IMP = data$Composite_AP_Implementation,
  COM = data$Composite_AP_Commitment
)

#Remove NAs
qca_COM_IMP_exclusions <- qca_COM_IMP %>% 
  drop_na(COM, IMP)

#Recalibrate COM AND IMP. 
  ## Low ns necessitate binary calibration. Strict calibration means partial values = 0
qca_COM_IMP_exclusions <- qca_COM_IMP_exclusions %>%
    mutate(
    IMP = case_when(
      IMP == -1 ~ 0,
      IMP == 0 ~ 0,
      IMP == 1 ~ 1
    ),
    COM = case_when(
      COM == -1 ~ 0,
      COM == 0 ~ 0,
      COM == 1 ~ 1
    )
  )

#Generate truth table
  ## Note: n.cut not used here due to low frequencies of 0s. Configurational analysis not conducted.
tt_COM_IMP <- truthTable(
  qca_COM_IMP_exclusions,
  outcome = "IMP",
  conditions = c("COM"),
  incl.cut = 0.8,
  type= "cs",
  show.cases = TRUE
)
tt_COM_IMP

# Create qca df. 
qca_DEC_IMP <- data.frame(
  IMP = data$Composite_AP_Implementation,
  DEC = data$Outcome_Dec_Quality,
  POS = data$Outcome_Positive_Org_Change
)

# Remove NAs
qca_DEC_IMP_exclusions <- qca_DEC_IMP %>% 
  drop_na(IMP, POS)

# Recalibrate DEC, IMP and POS
## EXPLORATORY: Retain cases not mentioning decision quality due to low n and variance
qca_DEC_IMP_exclusions <- qca_DEC_IMP_exclusions %>%
    mutate(
    DEC = case_when(
      is.na(qca_DEC_IMP_exclusions$DEC) ~ 0,
      DEC == 1 ~ 1
    ),
    IMP = case_when(
      IMP == -1 ~ 0,
      IMP == 0 ~ 0,
      IMP == 1 ~ 1
    ),
    POS = case_when(
      POS == -1 ~ 0,
      POS == 0 ~ 0,
      POS == 1 ~ 1
    )
  )

#Generate truth table
  ## Note: n.cut not used here due to low frequencies of 0s
tt_DEC_IMP <- truthTable(
  qca_DEC_IMP_exclusions,
  outcome = "POS",
  conditions = c("DEC", "IMP"),
  incl.cut = 0.8,
  type= "cs",
  show.cases = TRUE
)
tt_DEC_IMP

solution_DEC_IMP <- minimize(tt_DEC_IMP, details= TRUE, method= "QMC", negate= TRUE)
solution_DEC_IMP