## GGPLOT -> BOXPLOTS (Figures 2-4)

library(readODS)
measurements_cor <- read_ods("raw_data.ods", sheet = "Sheet1")

library(dplyr)
library(tidyr)
library(ggplot2)

measurements_cor %>%
  gather(key = "Measur.", value = "Distance", -c(1:6)) %>%
  mutate(across(WUSTL, ~factor(., levels = c("W023", "W031", "W096", "W188", "W264", "W281", "W297", "W343", "W451", "W481", "W489", "W535", "W614", "W620", "W650", "W728", "W781", "W805", "W914", "W955")))) %>%
  ggplot(aes(x = Skull, y = Distance, fill = Skull)) + 
  geom_boxplot(aes(color =  Skull), alpha = 0.5, outliers = F) +
  geom_point(aes(color = Skull, alpha = 0.75), size = 3) +
  facet_wrap(~WUSTL, scales = "free") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
        legend.position="none") +
  ylab("Basion-Prosthion distance (mm)") + xlab("Skull model")


# PAIRED T-TESTS (Table 1)

library(ggpubr)
library(rstatix)
stat.test_cor <- measurements_cor %>%
  gather(key = "Measur.", value = "Distance", -c(1:6)) %>%
  #filter(!Skull == "CT_Joganic") %>%
  filter(!Measurer == "Joganic") %>%
  group_by(WUSTL) %>%
  t_test(Distance ~ Skull,
         paired = TRUE) %>%
  adjust_pvalue(method = "bonferroni") %>%
  add_significance("p.adj")

library(kableExtra)
library(knitr)
kable(stat.test_cor) %>%
  kable_classic(full_width = T, 
                html_font = "Cambria", 
                font_size = 10)

writexl::write_xlsx(stat.test_cor, "stat_test_cor.xlsx")

# PREDICTION INTERVALS (Tables 2-3)

joganic_function_0.75 <- function(skull, level) {
  test <- measurements_cor %>%
    filter(WUSTL == skull) %>%
    select(3, 7:16) %>%
    t() %>% as.data.frame() 
  
  colnames(test) <- test[1,]
  test <- test %>%
    slice(2:nrow(test)) %>%
    mutate_all(as.numeric)
  
  joganic <- test[1, ncol(test)]
  
  p1 <- as.data.frame(predict(lm(test$CT_0.75  ~ 1), interval = "predict", level = level))
  p1$Sample <- "CT_0.75"
  p1$Min <- min(test$CT_0.75)
  p1$Max <- max(test$CT_0.75)
  p2 <- as.data.frame(predict(lm(test$Original ~ 1), interval = "predict", level = level))
  p2$Sample <- "Original"
  p2$Min <- min(test$Original)
  p2$Max <- max(test$Original)
  
  df_pr <- rbind(p1[1,], p2[1,])
  df_pr$Joganic <- joganic
  
  df_pr_final <- df_pr %>%
    mutate(Inside_Prediction = case_when(
      joganic >= lwr  & joganic < upr ~ "TRUE" ,
      joganic < lwr ~ "FALSE",
      joganic > upr ~ "FALSE")) %>%
    mutate(Skull = skull) %>%
    select(Skull, Sample, Min, Max, Joganic, lwr, upr, Inside_Prediction)
  
  df_pr_final
}

final_table_0.75 <- rbind(
  joganic_function_0.75("W023", level = 0.99),
  joganic_function_0.75("W031", level = 0.99),
  joganic_function_0.75("W096", level = 0.99),
  joganic_function_0.75("W188", level = 0.99),
  joganic_function_0.75("W264", level = 0.99),
  joganic_function_0.75("W281", level = 0.99),
  joganic_function_0.75("W297", level = 0.99),
  joganic_function_0.75("W343", level = 0.99),
  joganic_function_0.75("W451", level = 0.99),
  joganic_function_0.75("W481", level = 0.99)
)

joganic_function_0.60 <- function(skull, level) {
  test <- measurements_cor %>%
    filter(WUSTL == skull) %>%
    select(3, 7:16) %>%
    t() %>% as.data.frame() 
  
  colnames(test) <- test[1,]
  test <- test %>%
    slice(2:nrow(test)) %>%
    mutate_all(as.numeric)
  
  joganic <- test[1, ncol(test)]
  
  p1 <- as.data.frame(predict(lm(test$CT_0.60  ~ 1), interval = "predict", level = level))
  p1$Sample <- "CT_0.60"
  p1$Min <- min(test$CT_0.60)
  p1$Max <- max(test$CT_0.60)
  p2 <- as.data.frame(predict(lm(test$Original ~ 1), interval = "predict", level = level))
  p2$Sample <- "Original"
  p2$Min <- min(test$Original)
  p2$Max <- max(test$Original)
  p3 <- as.data.frame(predict(lm(test$CT_0.58  ~ 1), interval = "predict", level = level))
  p3$Sample <- "CT_0.58"
  p3$Min <- min(test$CT_0.58)
  p3$Max <- max(test$CT_0.58)
  
  df_pr <- rbind(p1[1,], p2[1,], p3[1,])
  df_pr$Joganic <- joganic
  
  df_pr_final <- df_pr %>%
    mutate(Inside_Prediction = case_when(
      joganic >= lwr  & joganic < upr ~ "TRUE" ,
      joganic < lwr ~ "FALSE",
      joganic > upr ~ "FALSE")) %>%
    mutate(Skull = skull) %>%
    select(Skull, Sample, Min, Max, Joganic, lwr, upr, Inside_Prediction)
  
  df_pr_final
}

final_table_0.60 <- rbind(
  joganic_function_0.60("W489", level = 0.99),
  joganic_function_0.60("W535", level = 0.99),
  joganic_function_0.60("W614", level = 0.99),
  joganic_function_0.60("W620", level = 0.99),
  joganic_function_0.60("W650", level = 0.99),
  joganic_function_0.60("W728", level = 0.99),
  joganic_function_0.60("W781", level = 0.99),
  joganic_function_0.60("W805", level = 0.99),
  joganic_function_0.60("W914", level = 0.99),
  joganic_function_0.60("W955", level = 0.99)
)

final_table <- rbind(final_table_0.75, final_table_0.60)