multiverse_analysis.Rmd

---
title: "Part 3: Bayesian Meta-analysis of BBS phenotypes"
abstract: 'This is a supplementary file for the paper "Meta-analysis of genotype-phenotype associations in Bardet-Biedl Syndrome uncovers differences among causative genes". Here, we compute how the conclusions of the paper hold under multiple different models. The complete source code for the analysis can be found at https://github.com/martinmodrak/bbs-metaanalysis-bayes or Zenodo, DOI: 10.5281/zenodo.3243264'
output:
  pdf_document: default
---

```{r setup_multiverse, echo=FALSE, message = FALSE, warning=FALSE}
knitr::opts_chunk$set(echo=FALSE, cache = TRUE)
library(rstan)
library(brms)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
library(skimr)
library(readxl)
library(here)
library(mice)
library(tidyverse)
library(tidybayes)
library(bayesplot)
library(cowplot)

source(here("data_processing.R"))
source(here("models.R"))
source(here("models_funcs.R"))
source(here("plots.R"))

include_validation = FALSE
```

In this part we show how conclusions we discuss in the main paper hold under different model choice`r if(include_validation){" and in the validation dataset" }`. A wide variety of models was tested. Those are briefly described in the comparison. Their exact formulations can be found in Part 2. `r if(include_validation){ "For a discussion of the validation dataset and the overall ability of the model to predict validation data, see Part 3." }`


*Note: For historical reasons the feature of "certain loss of function" (cLOF) as discussed in the data is called just "lof" in most analysis code. This part will thus use "lof" and "cLOF" interchangeably.*




```{r}
data <- read_main_data()
genes_to_show <- genes_to_show_from_data(data)
```


```{r load_all_fits}
all_fits <- all_models %>% map(function(def) {
  stored_fit_file <- paste0(here(stored_fits_dir,def$name), ".rds")
  if(!file.exists(stored_fit_file)) {
    stop(paste0("Computed fit for model '", def$name,"' cannot be found at ", stored_fit_file, ".\nYou probably need to run alternative_models.Rmd to compute all fits  (or download the fits from Zenodo)"))
  }
  readRDS(stored_fit_file)
})
```

## Defining Precise Criteria

```{r}
min_interesting_effect <- 2
min_phenotypes_for_diffs <- 3
min_bbsome_pairwise_for_diffs <- 5

evaluators <- list()
```

There is some flexibility in defining the exact model configurations consistent with a given conclusion. In general, when discussing whether there is or isn't a difference, the Bayesian model will always say "yes there is a (possibly small) difference" - the posterior probability of a difference of exactly 0 is 0. Instead we choose a threshold of a "clinically relevant effect", which for us is odds ratio outside of (`r 1/min_interesting_effect`,`r min_interesting_effect`) and any effect outside of these bounds is counted as "different".

Similarly, other plain-English statements need to be transformed to an exact predicate that can be considered true or false for a given assignment of numerical values to model coefficients, to let us evaluate its posterior probability.  The exact definition of the individual tested statements follows.


```{r}
n_genes_lof_different <- 5
```

  *	The severity of BBS is worse in patients with LOF mutations than in patients with other mutations. 
    * Ignored for models that do not include LOF.
    * Measured as posterior probability, that the LOF effect is positive (OR > 1) for at least `r n_genes_lof_different` BBSome genes in `r min_phenotypes_for_diffs` phenotypes.
```{r}

evaluator_lof <- function(def, fit) {
 if("lof" %in% def$additional_components || "lof per gene" %in% def$additional_components)  {
   
   genes_to_check <- bbsome_genes
   data_for_prediction_lof <- data_for_prediction_base_model(def, genes_to_check, phenotypes_to_show = phenotypes_to_show, age_transform = age_transform_from_age(data$age_numbers_groups_guessed))
   samples_tidy <- get_samples_lof_diff(fit, data_for_prediction_lof)
   
   samples_tidy %>% 
     group_by(sample, phenotype) %>% summarise(enough_genes = sum(odds_ratio > 1) >= n_genes_lof_different) %>%
     group_by(sample) %>% summarise(enough_phenotypes = sum(enough_genes) >= min_phenotypes_for_diffs) %>%
     ungroup() %>% summarise(`cLOF mutations more severe` = mean(enough_phenotypes))
 } else {
   tibble(`cLOF mutations more severe` = NA_real_)
 }
}
```
  
  *	The data suggest that mutations in BBS3 have lower severity than mutations in different functional groups of genes.
    * Measured as posterior probability, that odds ratio is < 1 for at least `r min_phenotypes_for_diffs` phenotypes for at least 1/2 of pairwise gene comparisons.
  
```{r}
evaluator_functional_groups_lower <- function(group1, group2) {
  function(samples_diff, samples_diff_bbsome, samples_prediction_bbsome) {
    target_column_name = paste(bbs_labeller(group1),"less severe than",bbs_labeller(group2))
    samples_diff %>% 
      filter(functional_group.x == group1, functional_group.y == group2) %>%
      group_by(sample, phenotype) %>%
      summarise(enough_pairwise_relevant = mean(odds_ratio < 1) >= 1/2) %>%
      group_by(sample) %>%
      summarise(enough_phenotype_relevant = sum(enough_pairwise_relevant) >= min_phenotypes_for_diffs) %>%
      ungroup() %>%
      summarise(!!target_column_name := mean(enough_phenotype_relevant))  
  }
}

evaluators[["bbs3_vs_bbsome"]] <- evaluator_functional_groups_lower("BBS03","BBSome")
evaluators[["bbs3_vs_chaperonins"]] <- evaluator_functional_groups_lower("BBS03","Chaperonins")
```
  *	The data suggest a difference between the severity of BBS in patients with mutations in different BBSome subunits.
    * Measured as posterior probability, that there is clinically relevant effect for at least `r min_phenotypes_for_diffs` phenotypes for at least `r min_bbsome_pairwise_for_diffs` pairwise comparisons. 
```{r}
evaluators[["bbsome_genes_differ"]] <- function(samples_diff, samples_diff_bbsome, samples_prediction_bbsome) {
samples_diff_bbsome %>% 
  filter(as.character(gene.x) > as.character(gene.y)) %>%
  group_by(sample, phenotype) %>%
  summarise(enough_pairwise_relevant = sum(is_relevant) >= min_bbsome_pairwise_for_diffs) %>%
  group_by(sample) %>%
  summarise(enough_phenotype_relevant = sum(enough_pairwise_relevant) >= min_phenotypes_for_diffs)  %>%
  ungroup() %>%
  summarise(`BBSome genes differ` = mean(enough_phenotype_relevant))
}
```

```{r}
    # def <- models_base$gene_source
    # fit <- all_fits[[def$name]]
    # data_for_prediction <- data_for_prediction_base_model(def, genes_to_show, phenotypes_to_show, age_transform = age_transform_from_age(data$age_numbers_groups_guessed))
    # 
    # data_for_prediction_bbsome <- data_for_prediction %>% filter_for_BBSome()
    # 
    # samples <- get_tidy_samples(fit, data_for_prediction)
    # samples_diff <- get_samples_pair♦wise_diff(def, samples) %>% mutate(is_relevant = odds_ratio >= min_interesting_effect | odds_ratio <= 1 /min_interesting_effect)
    # 
    # samples_diff_bbsome <- get_tidy_samples(fit, data_for_prediction_bbsome) %>% get_samples_pairwise_diff(def, .) %>% mutate(is_relevant = odds_ratio >= min_interesting_effect | odds_ratio <= 1 /min_interesting_effect)
    # 
    # samples_prediction_bbsome = get_tidy_samples_prediction(fit, data_for_prediction_bbsome)
    # 
    # samples_linear_bbsome <- get_tidy_samples(fit, data_for_prediction_bbsome, scale = "linear")

```

  *	BBS4 phenotype is the most severe of all BBSome-encoding genes.
    * Probability that odds for BBS4 are among the top 3 odds for at least 5 phenotypes.
    * Probability that a patient with BBS4 has the highest total number of phenotypes present.
```{r}
evaluators[["bbs04_top"]] <- function(samples_diff, samples_diff_bbsome, samples_prediction_bbsome) {
samples_diff_bbsome %>% 
  filter(gene.x == "BBS04", gene.y != "BBS04") %>%
  group_by(sample, phenotype) %>%
  summarise(bbs04_top = mean(odds_ratio < 1) <= 2/length(bbsome_genes)) %>% #Switched to mean to work well with validation
  group_by(sample) %>%
  summarise(bbs04_top_all = sum(bbs04_top) >= 5 ) %>%
  summarise(`BBS4 top odds` = mean(bbs04_top_all))
}
```
```{r}
evaluators[["bbs04_most_pheno"]] <- function(samples_diff, samples_diff_bbsome, samples_prediction_bbsome) {
  samples_prediction_bbsome %>% 
    group_by(sample, gene) %>%
    summarise(mean_phenotypes = mean(value)) %>% #Using mean instead of sum to play nice with uneven number of reported (for validation data)
    group_by(sample) %>%
    mutate(has_most_pheno = mean_phenotypes >= max(mean_phenotypes) ) %>%
    ungroup() %>%
    filter(gene == "BBS04") %>%
    #group_by(gene) %>%
    summarise(`BBS4 most #phenotypes` = mean(has_most_pheno))
}
```
  *	Mutations in different BBSome subunits predispose to different renal phenotype.
    *  Measured as posterior probability, that there is clinically relevant effect for at least `r min_bbsome_pairwise_for_diffs` pairwise comparisons. 
```{r}
evaluators[["REN_bbsome_genes_differ"]] <- function(samples_diff, samples_diff_bbsome, samples_prediction_bbsome) {
  samples_diff_bbsome %>% 
    filter(phenotype == "REN", as.character(gene.x) > as.character(gene.y)) %>%
    group_by(sample) %>%
    summarise(enough_pairwise_relevant = sum(is_relevant) >= min_bbsome_pairwise_for_diffs) %>%
    ungroup() %>%
    summarise(`REN: BBSome genes differ` = mean(enough_pairwise_relevant))
}
```

  *	Differences between small groups of genes: probability that all pairwise odds ratios are greater/less than 1 (depending on the direction of the comparison). 
    * Cognitive impairment is less frequent in BBS3 patients compared to other patients (all canonical BBS genes).
```{r}
evaluator_prob_group_diff_all <- function(phenotype_name, group1, group2, group1_name = paste0(group1, collapse = "_"),group2_name = paste0(group2, collapse = "_")) {
    function (samples_diff, samples_diff_bbsome, samples_prediction_bbsome) {

      target_column_name = paste0(phenotype_name,": ",bbs_labeller(group1_name)," more frequent than ", bbs_labeller(group2_name))
      samples_diff %>% 
        filter(phenotype == phenotype_name, gene.x %in% group1, gene.y %in% group2) %>%
        group_by(sample) %>%
        summarise(all_higher = all(odds_ratio >= 1)) %>%
        ungroup() %>%
        summarise(!!target_column_name := mean(all_higher))   
    }
}

evaluators[["CI_BBSO3_all"]] <- evaluator_prob_group_diff_all("CI", paste0("BBS0", c(1,2,4,5,7,8,9)), "BBS03", group1_name = "BBSome")
```

    * Cognitive impairment is more frequent in BBS7 patients compared to other patients with mutations in BBSome-encoding genes.
```{r}
evaluators[["CI_BBS07_all"]] <-  evaluator_prob_group_diff_all("CI", "BBS07", paste0("BBS0", c(1,2,4,5,8,9)), group2_name = "others")

```
  
    *	Renal involvement is less frequent in BBS1, BBS4 and BBS8 patients compared to BBS2, BBS7 and BBS9 patients.
```{r}
evaluators[["REN_groups_all"]] <- evaluator_prob_group_diff_all("REN", paste0("BBS0", c(2,7,9)),  paste0("BBS0", c(1,4,8)), group1_name = "BBS2,7,9", group2_name = "BBS1,4,8")
```
    *	Patients with BBS2 mutations are more likely to have heart anomalies compared to patients with other mutations in BBSome-encoding genes.
```{r}
evaluators[["HEART_BBS2_all"]] <-  evaluator_prob_group_diff_all("HEART", "BBS02", paste0("BBS0", c(1,4,5,7,8,9)), group2_name = "others")
```
    *	Patients with BBS5 mutations are more likely to have liver anomalies compared to patients with other mutations in BBSome-encoding genes.
```{r}
evaluators[["LIVER_BBS5_all"]] <-  evaluator_prob_group_diff_all("LIV", "BBS05", paste0("BBS0", c(1,2,4,7,8,9)), group2_name = "others")
```
    *	Patients with BBS4 mutations are less likely to have liver anomalies compared to patients with other BBSome mutations.

```{r}
evaluators[["LIVER_BBS4_all"]] <-  evaluator_prob_group_diff_all("LIV", paste0("BBS0", c(1,2,5,7,8,9)), "BBS04", group1_name = "others")
```  
  
  *	Individual differences in phenotype between two genes: directly the probability that the OR > 1 for patients with cLOF mutation.
    * Polydactyly is more frequent in BBS2 patients compared to BBS1 patients.
```{r}
evaluator_prob_pairwise_diff <- function(phenotype_name, gene_more_frequent, gene_less_frequent) {
  function (samples_diff, samples_diff_bbsome, samples_prediction_bbsome) {
     target_column_name = paste0(phenotype_name,": ", bbs_labeller(gene_more_frequent), " more frequent than ", bbs_labeller(gene_less_frequent))
      samples_diff %>% 
        filter(phenotype == phenotype_name, gene.x == gene_more_frequent, gene.y == gene_less_frequent) %>%
        group_by(sample) %>%
        summarise(all_higher = all(odds_ratio >= 1)) %>%
        ungroup() %>%
        summarise(!!target_column_name := mean(all_higher))   
  }
}

evaluators[["PD_BBS02_BBS01"]] <- evaluator_prob_pairwise_diff("PD","BBS02","BBS01")
```

    * Polydactyly is more frequent in BBS10 patients compared to BBS1 patients.
```{r}
evaluators[["PD_BBS10_BBS01"]] <- evaluator_prob_pairwise_diff("PD","BBS10","BBS01")

```

    *	Renal involvement is less frequent in BBS1 patients compared to BBS2 patients.
```{r}
evaluators[["REN_BBS1_BBS2"]] <- evaluator_prob_pairwise_diff("REN","BBS02","BBS01")
```
    *	Renal involvement is less frequent in BBS1 patients compared to BBS10 patients.
```{r}
evaluators[["REN_BBS1_BBS10"]] <- evaluator_prob_pairwise_diff("REN","BBS10","BBS01")
```
    *	Liver involvement is less frequent in BBS1 patients compared to BBS2 patients.
```{r}
evaluators[["LIV_BBS1_BBS2"]] <- evaluator_prob_pairwise_diff("LIV","BBS02","BBS01")
```
    *	Liver involvement is less frequent in BBS1 patients compared to BBS10 patients.
```{r}
evaluators[["LIV_BBS1_BBS10"]] <- evaluator_prob_pairwise_diff("LIV","BBS10","BBS01")
```
    


## Bayesian Comparison

```{r compute_evaluation}
#models_for_test <- all_models[c("gene_source_lof","gene_source_filtered_lof","gene_imputed_age_sex")]
#evaluation <- models_for_test %>%
evaluation <- all_models %>% 
 map_df(function(def) {
    fit <- all_fits[[def$name]]
    data_for_prediction <- data_for_prediction_base_model(def, genes_to_show, phenotypes_to_show, age_transform = age_transform_from_age(data$age_numbers_groups_guessed))
    
    data_for_prediction_bbsome <- data_for_prediction %>% filter_for_BBSome()
    
    samples <- get_tidy_samples(fit, data_for_prediction)
    samples_diff <- get_samples_pairwise_diff(def, samples) %>% mutate(is_relevant = odds_ratio >= min_interesting_effect | odds_ratio <= 1 /min_interesting_effect)
    
    samples_diff_bbsome <- get_tidy_samples(fit, data_for_prediction_bbsome) %>% get_samples_pairwise_diff(def, .) %>% mutate(is_relevant = odds_ratio >= min_interesting_effect | odds_ratio <= 1 /min_interesting_effect)
    
    samples_prediction_bbsome = get_tidy_samples_prediction(fit, data_for_prediction_bbsome)
    
    #samples_linear_bbsome <- get_tidy_samples(fit, data_for_prediction_bbsome, scale = "linear")
    
    lof_evaluation <- evaluator_lof(def, fit)
    
    evaluators %>% map(function(evaluator) {
      evaluator(samples_diff, samples_diff_bbsome, samples_prediction_bbsome)
    }) %>% 
      do.call(cbind, .) %>% cbind(lof_evaluation, .) %>%
      mutate(model_name = def$name)
})

```



```{r load_freqentist}
#Freq. conclusions contain order 
freq_conclusions <- readr::read_csv(here("data","BBS - Frequentist conclusions - p.csv"), col_types = cols(
  Conclusion = col_character(),
  `p-value-all` = col_double(),
  `p-value-cLOF` = col_double(),
  order = col_integer(),
  include_in_main = col_logical()
))


freq_conclusions_order <- order(freq_conclusions$order)
freq_conclusions <- freq_conclusions %>% 
  mutate(Conclusion = factor(Conclusion, levels = rev(Conclusion[freq_conclusions_order]))) %>%
  gather("model_name","p", `p-value-all`, `p-value-cLOF`) %>%
  mutate(model_nice = paste0("frequentist ", gsub("p-value-","", model_name)), Category = "Freq.")
```


```{r, fig.width = 11.5, fig.height = 6.5}
nice_model_name <- function(x)  {
  x %>% gsub("only","none", .) %>% 
    gsub("lof_per_gene","cLOF (by gene)", .) %>% 
    gsub("ethnic_group","ethnic group", .) %>% 
    gsub("lof","cLOF", .) %>% 
    gsub("gene_","",.) %>% gsub("filtered_","filtered ", .) %>%
    gsub("very_narrow","very narrow", .) %>% 
    gsub("cor"," corr.", .) %>% 
    gsub("imputed_","imputed ", .) %>%  gsub("_"," + ", .)
}

conclusion_names <- names(evaluation)[names(evaluation) != "model_name"]

if(!identical(sort(levels(freq_conclusions$Conclusion)), sort(conclusion_names))) {
  print(setdiff(levels(freq_conclusions$Conclusion), conclusion_names))
  print(setdiff(conclusion_names, levels(freq_conclusions$Conclusion)))
  stop("Conclusions mismatch")
}

min_prob = 1e-3
max_prob = 1-1e-3

evaluation_processed <- evaluation %>% 
  arrange(grepl("imputed", model_name), grepl("filtered", model_name), 
          desc(grepl("source", model_name)), grepl("none", model_name), grepl("family", model_name), 
          grepl("lof", model_name)) %>%
  mutate(model_nice = factor(nice_model_name(model_name), levels = c(
   nice_model_name(model_name)
  ))) %>%
  gather("Conclusion","Probability", -model_name, -model_nice) %>%
  mutate(Probability = pmin(pmax(Probability, min_prob),max_prob), #Avoid probabilities of 0 and 1
         Conclusion = factor(Conclusion, levels = levels(freq_conclusions$Conclusion)),
         Category = case_when(model_name %in% c("gene_source", "gene_source_lof", "gene_source_filtered_lof") ~ "Main",
                              !grepl("source", model_name) ~ "Problematic fits",
                              model_name %in% c("gene_source_very_narrow") ~ "Problematic fits", 
                              TRUE ~ "Secondary") %>% factor(levels = c("Main","Secondary","Problematic fits"))) 

plots_na_color <- "#606060"

build_conclusion_plot <- function(evaluation_processed, categories = TRUE) {
  if(categories) {
    facet <- facet_grid(.~Category, scales = "free_x", space = "free_x") 
  } else {
    facet <- NULL
  }
  breaks = c(0.01,0.1,0.5,0.9,0.99)
  logit <- scales::logit_trans()$transform
  transformed_values <- c(0, (logit(c(0.01,0.25,0.75,0.99)) - logit(min_prob)) / (logit(max_prob) - log(min_prob)), 1)

  evaluation_processed %>%
    ggplot(aes(x = model_nice, y = Conclusion, fill = Probability, colour = "")) + geom_tile() + 
      theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust =0.5)) + 
      scale_fill_gradientn("Posterior prob.", 
                           #colours = c("#1b7837","#f0f0f0", "#f0f0f0", "#762a83"),
                           # values = c(0, 0.25,0.75,1),  
                           # limits = c(0,1), na.value = plots_na_color) +
                           colours = c("#1b7837","#1b7837","#f0f0f0", "#f0f0f0", "#762a83", "#762a83"),
                           values = transformed_values,
                           limits = c(min_prob,max_prob),
                           breaks = breaks,
                           na.value = plots_na_color, trans = "logit") +
      scale_x_discrete("Model components/modifications besides gene") +
      scale_colour_manual(values=NA) +              
      #Using color scale to create legend for NA
      (if(any(is.na(evaluation_processed$Probability))) {
        guides(fill = guide_colorbar(order = 10), colour=guide_legend("ND", override.aes=list(colour=plots_na_color, fill = plots_na_color))) 
      } else { 
        guides(colour = FALSE)
      }
      ) +
      facet + base_theme
}

conclusions_plot <- build_conclusion_plot(evaluation_processed)
conclusions_plot

```

A heatmap of posterior probability of statements characterizing individual conclusions. Note that the probability is on logit scale. All Bayesian models include gene as covariate, but may also include additional covariates: source, age, sex, ethnicity or ethnic group, family and certain loss of function (cLOF) - either as a global covariate or by gene. Since age and sex are not available for all data, we can either fit the model only to patients where those are reported (filtered) or impute missing data (imputed). Instead of using cLOF as a covariate, we can fit the model using only patients with cLOF mutations (filtered cLOF). For most models we include a correlation structure across phenotypes (e.g., that two phenotypes occur frequently together across all genes), but this structure may be absent (no corr.) or replaced with a correlation structure across genes (gene corr. - e.g., that two genes have similar pattern of effects across all phenotypes). We also tried modifying the width of prior distributions (wide, narrow, very narrow). See Part 2 of this supplement for a detailed description of all models and the imputation procedure. Dark grey indicates that the question could not be evaluated for the given model (currently only asking for cLOF differences in models that exclude cLOF). The "Problematic fits" category is reserved for models we know do not capture some important variability in the dataset, as discussed in Part 2.

Note: posterior probabilities are clamped to be at least $0.001$ and at most $0.999$ as the sampling scheme used does not let us be very confident in the tails of the distribution. It is possible that with more computational resources some of the posterior probabilities will be more extreme.

Some patterns to notice:

* When fitting filtered datasets, there is less certainty implying less strong evidence in both directions.
* Using very narrow priors on gene coefficients ($N(0,0.1)$, i.e., that almost all odds ratios should be less than $\sim 1.2$) unsurprisingly results in little evidence for directional differences between genes.
* Other than noted above, the conclusions are not sensitive to model choice.

## Frequentist results 

```{r, fig.width=5.7, fig.height=5.1}
build_freq_conclusions_plot <- function(freq_conclusions, categories = TRUE) {

  significance_pos <- 1 - (log10(0.05) / log10(min(freq_conclusions$p, na.rm = TRUE)))
  
  if(categories) {
    facet <- facet_grid(.~Category, scales = "free_x", space = "free_x") 
  } else {
    facet <- NULL
  }
  
  freq_conclusions %>% ggplot(aes(x = model_nice, y = Conclusion, fill = p)) + geom_tile() + 
      theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust =0.5)) + 
      scale_fill_gradientn("p-value", 
                           colours = c("#7f0000", "#fdd49e","#f0f0f0", "#f0f0f0"),
                           values = c(0, significance_pos - 0.01, significance_pos, 1), 
                           trans = "log10",
                           breaks = c(1e-8,1e-4,0.05), labels = c("1e-8","1e-4","0.05"),
                           na.value = plots_na_color) +
      scale_x_discrete("Freq. type") +
      guides(fill = guide_colorbar(reverse = TRUE)) +
      expand_limits(fill = 1) +
      facet + base_theme
}

freq_conclusions_plot <-  build_freq_conclusions_plot(freq_conclusions)

freq_conclusions_plot
```

Heatmap of p-values from frequentist analysis for the same conclusions (exact computations described in the main body of the paper and in the analysis code and report on Zenodo, DOI: 10.5281/zenodo.3243400).

`r if(include_validation) {"## Validation dataset"}`


```{r evaluation_validation}
if(include_validation) {
  validation_data <- read_validation_data()
  validation_data_long <- data_long_from_data(validation_data)
  
  #Turn validation data into "samples" that the evaluators can manage
  samples_validation <- validation_data_long %>%
    group_by(source, phenotype, loss_of_function_certain, gene, functional_group) %>%
    summarise(value = mean(phenotype_value)) %>%
    ungroup() %>%
    mutate(sample = source, odds = value / (1 - value)) 
  
  
  samples_diff_validation <- get_samples_pairwise_diff(main_model_def, samples_validation) %>%
    mutate(odds_ratio = if_else(odds.x == odds.y, 1, odds_ratio)) %>% 
    mutate(is_relevant = odds_ratio >= min_interesting_effect | odds_ratio <= 1 /min_interesting_effect)
  
  if(any(is.na(samples_diff_validation$odds_ratio))) {
    stop("Problem")
  }
  
  
  samples_diff_bbsome_validation <- samples_diff_validation %>% 
    filter(gene.x %in% bbsome_genes, gene.y %in% bbsome_genes)
  
  samples_prediction_bbsome_validation <- validation_data_long %>%
    filter(gene %in% bbsome_genes) %>%
    mutate(value = phenotype_value, sample = source)
  
  
  samples_lof_diff_validation <- samples_validation %>%
    filter(loss_of_function_certain == 1) %>%
      inner_join(samples_validation %>% filter(loss_of_function_certain == 0), 
                 by = c("phenotype" = "phenotype", "gene" = "gene", "source",
                        "functional_group" = "functional_group")) %>%
      mutate(odds_ratio = odds.x / odds.y) %>%
      mutate(odds_ratio = if_else(odds.x == odds.y, 1, odds_ratio)) 
  
  if(any(is.na(samples_lof_diff_validation$odds_ratio))) {
    stop("Problem")
  }
  
  
  
  #Note: this slightly differs from lof evaluation for model as I need to renormalize to have the same retio of genes/phenotypes required to be different. I also take equal odds of 0/Inf as odds_ratio > 1, i.e. equal odds are counted as negative evidence only when
  lof_evaluation_validation <-  samples_lof_diff_validation %>% 
       group_by(source, phenotype) %>% summarise(enough_genes = mean(odds_ratio > 1 | (odds_ratio == 1 & (odds.x == 0 | odds.x == Inf))) >= n_genes_lof_different / length(bbsome_genes)) %>%
       group_by(source) %>% summarise(enough_phenotypes = mean(enough_genes) >= min_phenotypes_for_diffs / 9) %>%
       ungroup() %>% summarise(`cLOF mutations more severe` = mean(enough_phenotypes))
   
  #Apply the evaluators
  evaluation_validation <- evaluators %>% map( ~ .x(samples_diff_validation, samples_diff_bbsome_validation, samples_prediction_bbsome_validation)) %>% 
        do.call(cbind, .) %>% cbind(lof_evaluation_validation, .) %>%
    mutate(model_name = "Validation dataset", model_nice = model_name, Category = "V")
  
  evaluation_validation_processed <- evaluation_validation %>% 
    gather("Conclusion","Result", -model_name, -model_nice, -Category) %>%
    mutate(Conclusion = factor(Conclusion, levels = levels(freq_conclusions$Conclusion)))

}         

```


```{r}

if(include_validation) {
  build_validation_conclusions_plot <- function(eval_conclusions, categories = TRUE) {
  
    if(categories) {
      facet <- facet_grid(.~Category, scales = "free_x", space = "free_x") 
    } else {
      facet <- NULL
    }
    
    eval_conclusions %>% ggplot(aes(x = model_nice, y = Conclusion, fill = factor(Result))) + geom_tile() + 
        theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust =0.5)) + 
        scale_fill_manual("Validated", 
                             values = c("0" = "#b2182b","0.5" = "#f7f7f7", "1" = "#2166ac", "NaN" = plots_na_color),
                             #trans = "log10",
                             breaks = c(0,0.5,1), labels = c("No","Mixed","Yes"),
                             na.value = plots_na_color) +
        scale_x_discrete("") +
        guides(fill = guide_legend(reverse = TRUE)) +
        facet + base_theme
  }
  
  validation_conclusions_plot <-  build_validation_conclusions_plot(evaluation_validation_processed)
  
  print(validation_conclusions_plot)
} else {
  validation_conclusions_plot <- NULL
}
```

`r if(include_validation) { '
Heatmap of evaluation on validation dataset, using exactly the same criteria as for the evaluation in Bayesian model, taking each source in the validation data as a single "sample". "Mixed" means that the relationship holds in one source in the validation data but not the other. Claims that could be tested only in one source because the other had no data on the mutations/phenotypes involved do not count as "Mixed". The claim on BBS5 could not be validated as there are no patients with BBS5 mutation in the validation dataset. '} `

## Combining all results

```{r, fig.width= 12, fig.height=7}


build_combined_conclusions_plot <- function(freq_conclusions_plot, conclusions_plot, validation_plot, rel_widths_plots, rel_widths_legend, rel_heights_legend) {
  hide_guide_theme <- theme(legend.position = 'none')
  hide_y_axis_theme <- theme(axis.title.y = element_blank(), 
                             axis.text.y = element_blank(), 
                             axis.line.y = element_blank(),
                             axis.ticks.y = element_blank())
  
  if(!is.null(validation_plot)) {
    validation_plot_processed <- validation_plot + hide_y_axis_theme + hide_guide_theme + theme(plot.margin = margin(7,7,7, 0, "pt"))
    validation_legend <- cowplot::get_legend(validation_plot)
  } else {
    validation_plot_processed <- NULL
    validation_legend <- NULL
  }
  
  plot_grid(
    plot_grid(freq_conclusions_plot + hide_guide_theme + theme(plot.margin = margin(7,0,7,7, "pt")), 
              conclusions_plot + hide_y_axis_theme + hide_guide_theme + theme(plot.margin = margin(7,0,7, 0, "pt")), 
              validation_plot_processed, 
              align = "h", nrow = 1, rel_widths = rel_widths_plots),
    plot_grid(
      cowplot::get_legend(freq_conclusions_plot), 
      validation_legend,
      cowplot::get_legend(conclusions_plot),
      ncol = 1, nrow = 3, rel_heights = rel_heights_legend),
    rel_widths = rel_widths_legend
  )
}

evaluation_for_main <- evaluation_processed %>% 
  inner_join(freq_conclusions %>% select(Conclusion, include_in_main), by = c("Conclusion" = "Conclusion")) %>% filter(include_in_main)

freq_conclusions_for_main <- freq_conclusions %>% filter(include_in_main)

if(include_validation) {
  rel_widths_plots <- c(1.0, 1.1,0.1)
  rel_widths_legend <-  c(4,1)
  rel_heights_legend <-  c(1,0.9,1.5)
  
  evaluation_validation_for_main  <- evaluation_validation_processed %>% 
    inner_join(freq_conclusions %>% select(Conclusion, include_in_main), by = c("Conclusion" = "Conclusion")) %>% filter(include_in_main)
  
} else {
  rel_widths_plots <- c(1, 1.28, 0.01)
  rel_widths_legend <-  c(7,1)
  rel_heights_legend <-  c(1.0,0.01,1.1)
}

build_combined_conclusions_plot(freq_conclusions_plot, conclusions_plot, validation_conclusions_plot, rel_widths_plots = rel_widths_plots, rel_widths_legend = rel_widths_legend, rel_heights_legend = rel_heights_legend) 

#TODO: check accesibility for those with impaired color perception.
```

```{r, fig.width= 12, fig.height=6}
#Version with filtered conclusions for main paper
if(include_validation) {
  validation_plot <- build_validation_conclusions_plot(evaluation_validation_for_main)
} else {
  validation_plot <- NULL
}

combined_conclusions_plot <- build_combined_conclusions_plot(
  build_freq_conclusions_plot(freq_conclusions_for_main), 
  build_conclusion_plot(evaluation_for_main),
  validation_plot, 
  rel_widths_plots = rel_widths_plots, 
  rel_widths_legend = rel_widths_legend, 
  rel_heights_legend = rel_heights_legend) 

combined_conclusions_plot %>% ggsave(here("tmp_pics","conclusions.pdf"), ., width = 12, height = 6)
```


This is a combined plot for both frequentist and Bayesian analysis, once again showing mostly consistent results even when frequentist analyses are taken into account. The only big difference is that some of the HEART and LIV conclusions are not supported by most Bayesian analyses, or only those ignoring between-study variability. `r if(include_validation) { "The same conclusions also fail to hold in the validation dataset, indicating that between-study variability is indeed important."} ` Also, Part 1 discusses some of those conclusions and shows how between-study variability is likely important for them.

`r if(include_validation) { '
Note that the conclusions included here and their exact criteria were set before we had access to the validation dataset.

Let us investigate the inconsistencies between validation data and the main model. First, the claim that BBS4 is the most severe, which we would actually expect NOT to validate, as it is likely false according to our model, but the "BBS4 top odds" conclusion was validated. However, the BBS4 is only present in 5 patients "Elizabeth" source and although it has one of highest odds (taking all phenotypes together for simplicity), it is very similar to the odds for others, so this is comparable to noise level. Also the "among top 3 odds" clause we used is easier to satisfy when not all BBSome genes are present in the validation data: ' }`

```{r}
if(include_validation) {
validation_data_long %>%
  filter(gene %in% bbsome_genes) %>%
  group_by(source) %>%
  filter("BBS04" %in% gene) %>%
  group_by(source,gene, ID) %>%
  summarise(mean_pheno_1 = mean(phenotype_value), n_pheno = n()) %>%
  group_by(source,gene) %>%
  summarise(mean_pheno = mean(mean_pheno_1), n_patients = n()) %>%
  arrange(source, desc(mean_pheno)) %>% 
  kable()
}
```

Finally, a brief summary of the main analyses as shown in the main paper.

```{r fig.width=6.25, fig.height=3.5}
validation_in_main = include_validation

hide_titles_theme <- theme(axis.title.x = element_blank())


if(validation_in_main) {

  validation_plot_main <- build_validation_conclusions_plot(evaluation_validation_for_main, categories = FALSE) + hide_titles_theme
  
  rel_widths_plots <-  c(14.1,1.1,1.5)
  rel_heights_legend <-  c(1,0.6,1.1)
  
} else {
  validation_plot_main <- NULL
  rel_widths_plots <-  c(13.1,1.0,0.2)
  rel_heights_legend <-  c(1,0.01,1.1)
}

combined_conclusions_main_plot <- build_combined_conclusions_plot(
  build_freq_conclusions_plot(freq_conclusions_for_main %>% filter(model_nice == "frequentist all") %>% mutate(model_nice = "frequentist")
                              , categories = FALSE) + hide_titles_theme,
  build_conclusion_plot(evaluation_for_main %>% filter(model_name == "gene_source_lof") %>% mutate(model_nice = "Bayesian"),
                        categories = FALSE) + hide_titles_theme,
  validation_plot_main,
  rel_widths_plots = rel_widths_plots,
  rel_widths_legend = c(3.2,1),
  rel_heights_legend = rel_heights_legend
)

combined_conclusions_main_plot

combined_conclusions_main_plot %>% ggsave(here("tmp_pics","conclusions_main.pdf"), ., width = 6.25, height = 3.5)
```