# Load pre-installed, required packages
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(plotly))
suppressPackageStartupMessages(library(scales))
suppressPackageStartupMessages(library(readr))
suppressPackageStartupMessages(library(DT))

# Install and load the 'ggthemes' package
suppressWarnings(suppressMessages(install.packages("ggthemes", quite=TRUE)))
suppressPackageStartupMessages(library(ggthemes))

# Read the 'Predict Hair Fall' dataset
data <- read_csv("data/Predict Hair Fall.csv", show_col_types=FALSE)

# Encode data types and factors
data <- data %>%
	mutate(Id = as.integer(Id),
		   Age = as.integer(Age),
		   across(where(is.character), as.factor),
		   Stress = factor(Stress, levels=c("Low", "Moderate", "High"), ordered=TRUE),
		   `Hair Loss` = as.factor(`Hair Loss`)
		   )

# Backup original dataset
data1 <- read_csv("data/Predict Hair Fall.csv", show_col_types=FALSE)

# Function to install missing packages
install_if_missing <- function(p) {
  if (!require(p, character.only = TRUE)) {
    suppressWarnings(suppressMessages(install.packages(p, dependencies = TRUE)))
    suppressPackageStartupMessages(library(p, character.only = TRUE))
  }
}

# Function to tabularize summary statistics of a given variable in the dataset
summarize_variable <- function(data, variable=NULL, caption=NULL, dom="t", searching=FALSE, rownames=FALSE, pageLength=10) {
  if (is.null(variable)) {
    datatable(data, 
                  class="cell-border stripe", 
                  options = list(dom=dom, searching=searching, pageLength=pageLength),
                  rownames=rownames,
                  caption=if (!is.null(caption)) htmltools::tags$caption(
                    style = 'caption-side: top; text-align: center;',
                    caption
                  ) else NULL)
  } else {
	summarized_data <- data %>%
      filter(!!sym(variable) != "No Data") %>%
      group_by(!!sym(variable)) %>%
      summarise(Count = n(), .groups="drop") %>%
      {
        if (is.ordered(data[[variable]])) {
          arrange(., desc(!!sym(variable)))
        } else {
          arrange(., desc(Count))
        }
      } %>%
      mutate(Percentage = label_percent(accuracy=0.01)(Count/sum(Count)))
    
    datatable(summarized_data %>% rename(!!variable := !!sym(variable)),
                  class="cell-border stripe", 
                  options = list(dom=dom, searching=searching, pageLength=pageLength),
			      rownames=rownames,
                  caption=if (!is.null(caption)) htmltools::tags$caption(
                    style = 'caption-side: top; text-align: center;',
                    caption
                  ) else NULL)
  }
}
			  
# Function to check if a package is installed, and install it if not
install_if_missing <- function(package) {
  if (!require(package, character.only = TRUE)) {
    install.packages(package, dependencies = TRUE)
    library(package, character.only = TRUE)
  }
}
			  
# Function to create a propotional stack bar plot (used for subplots of multiple variables)
create_plotly_stackbar <- function(data = data, var, showlegend=FALSE, yaxis_text=NULL, yaxis_titleSize=11.5, yaxis_tickSize=10, tickangle=0) {
  yaxis_title <- ifelse(is.null(yaxis_text), var, yaxis_text)
  
  plot_data <- as.data.frame(table(data %>% filter(!!sym(var) != "No Data") %>%
								   		select(`Hair Loss`, all_of(var)))) %>%
    rename_with(~ gsub("\\.", " ", .)) %>%
	group_by(!!sym(var)) %>%
    mutate(Percentage = 100*Freq/sum(Freq),
           `Hair Loss` = ifelse(`Hair Loss` == 1, "With baldness", "Without baldness")
          ) %>%
    rename(Count = Freq)
  
  plot <- plot_ly(
    data = plot_data,
    x = ~Percentage,
    y = as.formula(paste0("~`", var, "`")),
    type = "bar",
    color = ~`Hair Loss`,
    colors = c("#FF6347", "#8c8c8c"),
    showlegend = showlegend,
    text = ~paste0("<b>Percentage:</b> ", label_percent(accuracy=0.01)(Percentage/100), "\n<b>Count:</b> ", Count),
    hoverinfo = "text"
  ) %>%
    layout(
      xaxis = list(title = list(text="Percentage (%)", font = list(family="ITC Officina Sans", size=18))),
      yaxis = list(title = list(text = paste("\n<b>", yaxis_title, "</b>"),
								font = list(family="ITC Officina Sans", size = yaxis_titleSize), 
								standoff=43),
                   tickfont = list(family="ITC Officina Sans", size = yaxis_tickSize), tickangle=tickangle
				  ),
      barmode="stack",
      legend = list(orientation="h", xanchor="center", x=0.438, y=1.11, traceorder="normal",
					font = list(family="ITC Officina Sans", size=15)),
	  font = list(size = 14),
	  hoverlabel = list(font = list(family="ITC Officina Sans", size=15))
    )
  
  return(plot)
}			  
			  
# Function to create a data frame for contingency table of a variable			  
create_contingency_data <- function(data, var) {
	
  # Select relevant columns and filter out "No Data" values
  selected_data <- data %>% 
    select(`Hair Loss`, all_of(var)) %>%
    filter(!!sym(var) != "No Data")
  
  # Create contingency table
  contingency_data <- data.frame(matrix(
    with(selected_data, table(selected_data[[var]], `Hair Loss`)),
    nrow = dim(with(selected_data, table(selected_data[[var]], `Hair Loss`)))[1], 
    dimnames = dimnames(with(selected_data, table(selected_data[[var]], `Hair Loss`)))
  ))
  
  return(contingency_data)
}
			  
# Function to create stack bar plot for a variable
create_stackbar_plot1 <- function(data, var, yaxis_title, yaxis_text=NULL, plot_title=NULL, tickangle=45) {
  stackbar_data <- data %>% filter(!!sym(var) != "No Data") %>%
    mutate(!!sym(var) := factor(!!sym(var),
								levels = rownames(create_contingency_data(data, var) %>% 
                                                    mutate(X01 = X0+X1, X0 = X0/X01, X1 = X1/X01) %>%
													arrange(X1, desc(rownames(.)))
												 )
    ))
  
  plot <- create_plotly_stackbar(data = stackbar_data, 
								 var = var, 
								 showlegend=TRUE,
								 yaxis_text = yaxis_text,
								 yaxis_titleSize=20,
								 yaxis_tickSize=15,
								 tickangle = tickangle)
  
  plot %>%
    layout(
      plot_bgcolor = "#D5E4EB",
      paper_bgcolor = "#D5E4EB",
      title = list(text = plot_title, 
                   font = list(family = "ITC Officina Sans", size = 18), 
                   y = 0.93),
      margin = list(t = 95),
      xaxis = list(tickfont = list(family = "ITC Officina Sans", size = 15))
    )
}			  

# Function to create a contingency table of a variable				  
create_contingency_table <- function(data, var, first_colname=var) {
  # Append the proportion to each value in the dataframe
  data <- as.data.frame(
    apply(data %>% arrange(desc(X1)), c(1, 2), function(x) {
      percentage <- label_percent(accuracy=0.01)(x / sum(data))
      sprintf("%d (%s)", x, percentage)
    })
  ) %>%
    rownames_to_column(var=first_colname) 
  
  # Rename rows and columns
  colnames(data) <- c(first_colname, "Without baldness", "With baldness")
  
  # Reorder columns
  data <- data[,c(first_colname, "With baldness", "Without baldness")]								
  
  return(data)
}
			  
# Function to compute Point-Biserial Correlation
compute_point_biserial <- function(data, x, y, conf=0.95) {
  # Ensure the x and y are in the data
  if (!(x %in% colnames(data)) | !(y %in% colnames(data))) {
    stop("The variable x or y not found in the data.")
  }
  
  # Ensure the y contains only the specified values
  unique_y <- unique(data[[y]])
  if (length(unique_y) != 2) {
    stop("The variable y is not binary.")
  }
  
  # Extract values for individuals with and without the binary outcome
  with_group <- data[[x]][data[[y]] == 1]
  without_group <- data[[x]][data[[y]] == 0]
  
  # Calculate means
  X1 <- mean(with_group)
  X0 <- mean(without_group)
  
  # Calculate standard deviation of the entire x data
  sn <- sd(data[[x]])
  
  # Calculate sample sizes
  n1 <- length(with_group)
  n0 <- length(without_group)
  n <- n1 + n0
  
  # Compute the point biserial correlation coefficient
  rpb <- (X1 - X0) / sn * sqrt((n1 * n0) / (n^2))
  
  # Compute confidence intervals using cor.test
  rpb_cortest <- cor.test(data[[x]], data[[y]], conf.level=conf)
  lower_ci <- rpb_cortest$conf.int[1]
  upper_ci <- rpb_cortest$conf.int[2]
  
  rpb_df <- data.frame(x=x, y=y, X1=X1, X0=X0, rpb=rpb, lower_ci=lower_ci, upper_ci=upper_ci
                      ) %>%
    mutate(across(where(is.numeric) & !where(is.integer), ~ ifelse(abs(.) >= 0.01, 
                                          label_number(accuracy = 0.01)(.), 
                                          label_scientific(digits = 3)(.)))
    ) %>%
    rename(`Correlation coefficient` = rpb) %>%
	mutate(!!paste0(conf * 100, "% CI") := paste0("[", lower_ci, ", ", upper_ci, "]"),
		   Method="Point-biserial") %>%
	select(-c("lower_ci", "upper_ci"))
  
  return(rpb_df)
}			

# Function to compute Rank-Biserial Correlation
compute_rank_biserial <- function(data, x, y, o=0, l=1, conf=0.95) {
  # Install and load the 'rcompanion' package if not already installed
  if (!requireNamespace("rcompanion", quietly = TRUE)) {
    suppressWarnings(suppressMessages(install.packages("rcompanion", quiet=TRUE)))
  }
  suppressPackageStartupMessages(library(rcompanion))
  
  # Ensure the x and y are in the data
  if (!(x %in% colnames(data)) | !(y %in% colnames(data))) {
    stop("The variable x or y not found in the data.")
  }
  
  # Ensure the y contains only the specified values
  unique_y <- unique(data[[y]])
  if (length(unique_y) != 2) {
    stop("The variable y is not binary.")
  }
    
  # Rank the predictor variable
  data$ranked_factor <- rank(data[[x]], ties.method = "average") 
  with_group <- data$ranked_factor[data[[y]] == 1]
  without_group <- data$ranked_factor[data[[y]] == 0]	
	
  # Calculate mean ranks for individuals with and without the binary outcome
  R1 <- mean(with_group)
  R0 <- mean(without_group)
 
  # Compute outputs using the wilcoxonRG function
  group <- factor(c(rep("With group", length(with_group)), rep("Without group", length(without_group))))	
  rrb_wilcoxonRG <- wilcoxonRG(x = c(with_group, without_group), g = group, ci=TRUE, conf=conf)
	
  # Extract coefficient and confidence intervals
  rrb <- rrb_wilcoxonRG[[1]]
  lower_ci <- rrb_wilcoxonRG[[2]]
  upper_ci <- rrb_wilcoxonRG[[3]]
  
  rrb_df <- data.frame(x=x, y=y, R1=R1, R0=R0, rrb=rrb, lower_ci=lower_ci, upper_ci=upper_ci
					  ) %>%
    mutate(across(where(is.numeric) & !where(is.integer), ~ ifelse(abs(.) >= 0.01, 
                                          label_number(accuracy = 0.01)(.), 
                                          label_scientific(digits = 3)(.)))
    ) %>%
    rename(`Correlation coefficient` = rrb) %>%
	mutate(!!paste0(conf * 100, "% CI") := paste0("[", lower_ci, ", ", upper_ci, "]"),
		   Method="Rank-biserial") %>%
	select(-c("lower_ci", "upper_ci")) 
  
  return(rrb_df)
}

# Function to compute Phi Coefficient
compute_phi <- function(data, x, y, o=0, l=1, u, t, conf=0.95) {
  # Install and load the 'statpsych' package if not already installed
  if (!requireNamespace("statpsych", quietly = TRUE)) {
    suppressWarnings(suppressMessages(install.packages("statpsych", quiet=TRUE)))
  }
  suppressPackageStartupMessages(library(statpsych))
  
  # Ensure the x and y are in the data
  if (!(x %in% colnames(data)) | !(y %in% colnames(data))) {
    stop("The variable x or y not found in the data.")
  }
  
  # Ensure the x and y contain only binary values
  unique_x <- unique(data[[x]])
  unique_y <- unique(data[[y]])
  if (length(unique_x) != 2 | length(unique_y) != 2) {
    stop("The variables x and y must be binary.")
  }
  
  # Create contingency table
  contingency_table <- table(data[[x]], data[[y]])
  
  # Convert o, l, u, t to character
  o <- as.character(o)
  l <- as.character(l)
  u <- as.character(u)
  t <- as.character(t)
  
  # Extract counts from the contingency table using provided levels
  a <- contingency_table[t, l]
  b <- contingency_table[t, o]
  c <- contingency_table[u, l]
  d <- contingency_table[u, o]
  
  # Compute outputs using the function ci.phi
  ci_phi <- ci.phi(alpha = 1-conf, d, c, b, a)
  
  # Compute confidence interval limits
  phi <- ci_phi[[1]]
  ci_lower <- ci_phi[[3]]
  ci_upper <- ci_phi[[4]]
  
  # Create a dataframe to return the results
  phi_df <- data.frame(x = x, y = y, counts=paste(a,b,c,d, sep=", "), phi = phi, 
                       ci = paste0("[", round(ci_lower, 2), ", ", round(ci_upper, 2), "]")) %>%
    mutate(across(where(is.numeric) & !where(is.integer), ~ ifelse(abs(.) > 0.01, 
                                                                  label_number(accuracy = 0.01)(.), 
                                                                  label_scientific(digits = 3)(.)))
    ) %>%
    rename(`a, b, c, d` = counts,
           `Phi coefficient` = phi,
		   !!paste0(conf * 100, "% CI") := ci)
  
  return(phi_df)
}		

# Function to compute Cramer's V
compute_cramersV <- function(data, x, y, conf=0.95) {
  # Ensure the x and y are in the data
  if (!(x %in% colnames(data)) | !(y %in% colnames(data))) {
    stop("The variable x or y not found in the data.")
  }
  
  # Create contingency table
  contingency_table <- table(data[[x]], data[[y]])
  
  # Compute chi-square statistic
  chi2 <- chisq.test(contingency_table, simulate.p.value = TRUE)$statistic
  
  # Get the number of rows and columns in the contingency table
  n <- sum(contingency_table)
  r <- nrow(contingency_table)
  c <- ncol(contingency_table)
  rc <- paste(r, c, sep=", ")	
  
  # Compute outputs using the cramerV function
  cramersV <- cramerV(contingency_table, ci=TRUE, conf=conf)
  
  # Extract coefficient and confidence intervals
  cramers_v <- cramersV[[1]]
  lower_ci <- cramersV[[2]]
  upper_ci <- cramersV[[3]]
  
  # Create a dataframe to return the results
  cramers_v_df <- data.frame(
    n = n, 
    rc = rc, 
    x = x, 
    y = y, 
    chi2 = chi2, 
    cramers_v = cramers_v, 
    lower_ci = lower_ci, 
    upper_ci = upper_ci
  ) %>%
    mutate(across(
      where(is.numeric) & !where(is.integer), 
      ~ ifelse(abs(.) > 0.01, label_number(accuracy = 0.01)(.), label_scientific(digits = 3)(.))
    )) %>%
    rename(
      `Chi-square statistic` = chi2,
      `r, c` = rc,
      `Cramer's V` = cramers_v
    ) %>%
    mutate(!!paste0(conf * 100, "% CI") := paste0("[", lower_ci, ", ", upper_ci, "]")) %>%
    select(-c("lower_ci", "upper_ci"))
  
  rownames(cramers_v_df) <- c(1:nrow(cramers_v_df))
  
  return(cramers_v_df)
}
			
# Function for evaluating model performance
evaluate_model <- function(model, test, y, model_name) {
  
  # Check if pROC package is installed, if not, install it
  if (!requireNamespace("pROC", quietly = TRUE)) {
    install.packages("pROC")
  }
  suppressPackageStartupMessages(library(pROC))
  
  # Perform predictions using the model and test data
  prob <- predict(model, test, type = "prob")[, 2]
  predicted <- predict(model, test)
  
  # Create a confusion matrix to evaluate the model's performance
  conf_matrix <- confusionMatrix(predicted, test[[y]])
  
  # Compute model metrics
  roc <- roc(test[[y]], prob)
  auc <- auc(roc)
  accuracy <- conf_matrix$overall['Accuracy']
  precision <- conf_matrix$byClass['Pos Pred Value']
  recall <- conf_matrix$byClass['Sensitivity']
  f1 <- 2 * (precision * recall) / (precision + recall)
  
  # Return performance metrics as a dataframe
  results <- data.frame(
    AUC = round(auc, 4),
    Accuracy = round(accuracy, 4),
    Precision = round(precision, 4),
    Recall = round(recall, 4),
    f1 = round(f1, 4)
  ) %>%
  rename(`F1 Score` = f1)
  
  rownames(results) <- model_name
  
  return(results)
}			  
			  		  
# Function for cluster statistics
cstats.table <- function(dist, tree, k) {
clust.assess <- c("cluster.number","n","within.cluster.ss","average.within","average.between",
                  "wb.ratio","dunn2","avg.silwidth")
clust.size <- c("cluster.size")
stats.names <- c()
row.clust <- c()
output.stats <- matrix(ncol = k, nrow = length(clust.assess))
cluster.sizes <- matrix(ncol = k, nrow = k)
for(i in c(1:k)){
  row.clust[i] <- paste("Cluster-", i, " size")
}
for(i in c(2:k)){
  stats.names[i] <- paste("Test", i-1)
  
  for(j in seq_along(clust.assess)){
    output.stats[j, i] <- unlist(cluster.stats(d = dist, clustering = cutree(tree, k = i))[clust.assess])[j]
    
  }
  
  for(d in 1:k) {
    cluster.sizes[d, i] <- unlist(cluster.stats(d = dist, clustering = cutree(tree, k = i))[clust.size])[d]
    dim(cluster.sizes[d, i]) <- c(length(cluster.sizes[i]), 1)
    cluster.sizes[d, i]
    
  }
}
output.stats.df <- data.frame(output.stats)
cluster.sizes <- data.frame(cluster.sizes)
cluster.sizes[is.na(cluster.sizes)] <- 0
rows.all <- c(clust.assess, row.clust)
output <- rbind(output.stats.df, cluster.sizes)[ ,-1]
colnames(output) <- stats.names[2:k]
rownames(output) <- rows.all
is.num <- sapply(output, is.numeric)
output[is.num] <- lapply(output[is.num], round, 2)
output
}			  		  
			  		  
# Tabularize first five rows of the dataset			  
summarize_variable(slice(data, 1:5), dom="t")		  

# Install and load the 'moments' package for calculating skewness
suppressWarnings(suppressMessages(install.packages("moments", quiet=TRUE)))
suppressPackageStartupMessages(library(moments))

# Summary statistics of the age
age_stats <- data.frame(Value = unclass(summary(data$Age))) %>%
  bind_rows(data.frame(Value = c(sd(data$Age), skewness(data$Age)), row.names = c("Std. Dev.", "Sk\u209a"))) %>%
  mutate_if(is.numeric, round, 2)

# Age distribution plot
age_dist_plot <- ggplot(data, aes(x = Age, text = paste0("<b>Count:</b> ", after_stat(count)))) +
  geom_histogram(breaks = hist(data$Age, breaks="FD", plot = FALSE)$breaks) +
  geom_vline(aes(xintercept = mean(Age), text = paste0("<b>Mean: </b>", round(mean(Age), 2))), col="#FFA53E", linewidth=0.5) +
  geom_vline(aes(xintercept = median(Age), text = paste0("<b>Median: </b>", round(median(Age), 2))), col="#4C6DB4", linewidth=0.5) +
  theme_economist() + 
  scale_color_economist() +
  scale_x_continuous(breaks = seq(15, 50, by=5)) +
  ggtitle("\nFIG. 1﹕HISTOGRAM OF THE AGE DISTRIBUTION OF SURVEYED INDIVIDUALS") +
  theme(
    plot.title = element_text(size=13, hjust=0.5),
    axis.title.x = element_text(family="ITC Officina Sans", size=13),
    axis.title.y = element_text(family="ITC Officina Sans", size=13)
  ) +
  labs(y="Count")

# Convert ggplot to plotly
age_dist_plotly <- ggplotly(age_dist_plot, tooltip=c("text")) %>%
  layout(
    hoverlabel = list(font = list(family="ITC Officina Sans", size=15)),
    xaxis = list(tickfont = list(family="ITC Officina Sans", size=13)),
    yaxis = list(tickfont = list(family="ITC Officina Sans", size=13))
  ) %>%
  add_annotations(
    text = paste0(
      "<b>Std. Deviation (<i>\u03c3</i>): </b>", round(sd(data$Age), 2), 
      ";  <b>Skewness (<i>Sk\u209a</i>): </b>", round(skewness(data$Age), 2)
    ),
    x = 0.037, y = 0.97, xref = "paper", yref = "paper",
    showarrow = FALSE, font = list(family="ITC Officina Sans", size=13.5), align="left"
  )

age_dist_plotly

# Tabularize summary statistics of Genetics
summarize_variable(
  data = data %>% rename(`Has a family history of baldness` = Genetics),
  variable = "Has a family history of baldness",
  rownames = FALSE
)

# Tabularize summary statistics of Hormonal Changes
summarize_variable(
  data %>% rename(`Has experienced hormonal changes` = `Hormonal Changes`),
  variable = "Has experienced hormonal changes",
  rownames = FALSE
)

# Tabularize summary statistics of Medical Conditions
summarize_variable(
  data = data %>% rename(`Medical condition that may lead to baldness` = `Medical Conditions`),
  variable = "Medical condition that may lead to baldness",
  dom = "tip",
  pageLength = 5
)

# Tabularize summary statistics of Medications & Treatments
summarize_variable(
  data = data %>% rename(`Medication or treatment that may cause hair loss` = `Medications & Treatments`),
  variable = "Medication or treatment that may cause hair loss",
  dom = "tip",
  pageLength = 5
)

# Tabularize summary statistics of Nutritional Deficiencies
summarize_variable(
  data %>% rename(`Nutritional deficiency that may contribute to hair loss` = `Nutritional Deficiencies`),
  variable = "Nutritional deficiency that may contribute to hair loss",
  dom = "tip",
  pageLength = 5
)

# Tabularize summary statistics of Stress
summarize_variable(
  data %>% 
    mutate(
      Stress = factor(Stress, levels = c("High", "Moderate", "Low"), ordered = TRUE)
    ) %>%
    rename(`Stress level` = Stress),
  variable = "Stress level",
  rownames = FALSE
)

# Tabularize summary statistics of Poor Hair Care Habits
summarize_variable(
  data %>%
    rename(
      `Practices poor hair care habits` = `Poor Hair Care Habits`
    ),
  variable = "Practices poor hair care habits",
  rownames = FALSE
)
			  
# Tabularize summary statistics of Environmental Factors
summarize_variable(
  data %>%
    rename(
      `Exposed to environmental factors that may contribute to hair loss` = `Environmental Factors`
    ),
  variable = "Exposed to environmental factors that may contribute to hair loss",
  rownames = FALSE
)
			  
# Tabularize summary statistics of Smoking
summarize_variable(
  data,
  variable = "Smoking",
  rownames = FALSE
)
			  
# Tabularize summary statistics of Weight Loss
summarize_variable(
  data %>%
    rename(
      `Has experienced significant weight loss` = `Weight Loss`
    ),
  variable = "Has experienced significant weight loss",
  rownames = FALSE
)	
					  
# Tabularize summary statistics of Hair Loss
summarize_variable(
  data %>%
    mutate(
      `Hair Loss` = ifelse(`Hair Loss` == 1, "Yes", "No"),
      `Hair Loss` = factor(`Hair Loss`, levels = c("No", "Yes"), ordered = TRUE)
    ) %>%
    rename(`With baldness` = `Hair Loss`),
  variable = "With baldness",
  rownames = FALSE
)

# Data frame for the pie graph, by Hair Loss
hair_loss <- data %>%
  select(`Hair Loss`) %>%
  group_by(`Hair Loss`) %>%
  summarise(Count = n(), .groups="drop") %>% 
  mutate(Proportion = Count/sum(Count)) %>%
  arrange(desc(`Hair Loss`), desc(Count))

# Factorize the Hair Loss variable
hair_loss$`Hair Loss` <- factor(hair_loss$`Hair Loss`, levels=c(0, 1))

# Create a pie graph for hair loss experience
hair_loss_pie <- plot_ly(
  type="pie",
  labels=paste0("<b>",ifelse(hair_loss$`Hair Loss` == 1, "With baldness", "Without baldness"),"</b>"),
  values = hair_loss$Count,
  textinfo="label",
  textposition="inside",
  hoverinfo="text",
  hovertext = paste0(
    "<b>Percentage:</b> ", label_percent(accuracy=0.01)(hair_loss$Proportion),
    "<b>\nCount:</b> ", hair_loss$Count
  ),
  marker = list(
    colors = case_when(
      hair_loss$`Hair Loss`==1 ~ "#FF6347",
      hair_loss$`Hair Loss`==0 ~ "#8c8c8c"
    )
  ),
  textfont = list(
    family="ITC Officina Sans",
    size=18,
    color="white"
  )
) %>%
  layout(
    title = list(
      text="𝗙𝗜𝗚. 𝟮﹕𝗣𝗜𝗘 𝗚𝗥𝗔𝗣𝗛 𝗢𝗙 𝗧𝗛𝗘 𝗦𝗨𝗥𝗩𝗘𝗬𝗘𝗗 𝗜𝗡𝗗𝗜𝗩𝗜𝗗𝗨𝗔𝗟𝗦\n𝗕𝗬 𝗛𝗔𝗜𝗥 𝗟𝗢𝗦𝗦 𝗢𝗨𝗧𝗖𝗢𝗠𝗘",
      font = list(size=18)
    ),
    paper_bgcolor="#D5E4EB",
    margin = list(l=0, r=0, b=20, t=80),
    showlegend = FALSE  # Hide the legend
  ) %>%
  style(
    hoverlabel = list(
      font = list(
        family="ITC Officina Sans",
        size=15
      )
    )
  )

hair_loss_pie

# Create dataframe for Age Group
age_data <- data %>%
  select(Age, `Hair Loss`) %>%
  mutate(
    `Age Group` = case_when(
      Age >= 17 & Age <= 30 ~ "Young (17-30)",
      Age >= 31 & Age <= 45 ~ "Middle-aged (31-45)",
      Age > 45 ~ "Old (>45)",
      TRUE ~ NA_character_
    ),
    `Age Group` = factor(`Age Group`, levels = c("Young (17-30)", "Middle-aged (31-45)", "Old (>45)"), ordered = TRUE)
  )

# Create contingency table by Age Group and Hair Loss outcome
age_table <- create_contingency_table(create_contingency_data(age_data, "Age Group"), "Age Group")

# Create a boxplot by Age Group and Hair Loss outcome
age_hair_loss_boxplot <- plot_ly(
  data = age_data,
  x = ~factor(`Hair Loss`, labels = c("Without Baldness", "With Baldness")),
  y = ~Age,
  type = "box",
  boxpoints = "outliers", # Show only outliers
  jitter = 0.3, # Add some jitter to the points
  color = ~factor(`Hair Loss`, labels = c("Without Baldness", "With Baldness")), # Add color by Hair Loss
  colors = c('#8c8c8c', '#FF6347') # Define colors
) %>%
  layout(
    title = list(
      text = "\n𝗙𝗜𝗚. 𝟯﹕𝗕𝗢𝗫𝗣𝗟𝗢𝗧 𝗢𝗙 𝗧𝗛𝗘 𝗔𝗚𝗘 𝗗𝗜𝗦𝗧𝗥𝗜𝗕𝗨𝗧𝗜𝗢𝗡 𝗕𝗬 𝗛𝗔𝗜𝗥 𝗟𝗢𝗦𝗦 𝗢𝗨𝗧𝗖𝗢𝗠𝗘",
      font = list(size = 18),
      x = 0.5,
      xanchor = "center",
      y = 0.95,
      yanchor = "top"
    ),
    margin = list(t = 100), # Add margin between title and plot
    xaxis = list(
      title = "Hair Loss Outcome",
      tickfont = list(family = "ITC Officina Sans"),
      titlefont = list(family = "ITC Officina Sans", size = 18),
      showticklabels = FALSE # Do not show x ticks texts
    ),
    yaxis = list(
      title = "Age",
      tickfont = list(family = "ITC Officina Sans"),
      titlefont = list(family = "ITC Officina Sans", size = 18)
    ),
    paper_bgcolor = "#D5E4EB",
    plot_bgcolor = "#D5E4EB",
    legend = list(
      orientation = "h",
      x = 0.5,
      xanchor = "center",
      y = 1.1,
      yanchor = "top"
    )
  ) %>%
  style(
    hoverlabel = list(
      font = list(
        family = "ITC Officina Sans",
        size = 15
      )
    )
  )

age_hair_loss_boxplot

# Data frame for the treemap, by Age Group and by Hair Loss outcome
age_grp_and_hair_loss <- data %>%
  select(Age, `Hair Loss`) %>%
  mutate(
    `Age Group` = case_when(
      Age >= 17 & Age <= 30 ~ "Young (17-30)",
      Age >= 31 & Age <= 45 ~ "Middle-aged (31-45)",
      Age > 45 ~ "Old (>45)",
      TRUE ~ NA_character_
    )
  ) %>%
  group_by(`Hair Loss`, `Age Group`) %>%
  summarise(Count = n(), .groups = "drop") %>%
  mutate(Proportion = Count / sum(Count)) %>%
  arrange(desc(`Hair Loss`), desc(Count))

# Factorize Age Group
age_grp_and_hair_loss$`Age Group` <- factor(
  age_grp_and_hair_loss$`Age Group`,
  levels = c("Young (17-30)", "Middle-aged (31-45)", "Old (>45)"),
  ordered = TRUE
)

# Create a treemap for by Age Group and Hair Loss outcome
colors <- c("#175952", "#762127", "#0b385b", "#01877c", "#ae1726", "#005289")
age_hair_loss_treemap <- plot_ly(
  type = "treemap",
  labels = paste0(
    "<b>", stringr::str_to_sentence(age_grp_and_hair_loss$`Age Group`),
    "\n", ifelse(age_grp_and_hair_loss$`Hair Loss` == 1, "with baldness", "without baldness"), "</b>"
  ),
  parents = rep("", nrow(age_grp_and_hair_loss)),
  values = age_grp_and_hair_loss$Count,
  domain = list(column = 0),
  name = " ",
  hoverinfo = "text",
  hovertext = paste0(
    "<b>Percentage:</b> ", label_percent(accuracy = 0.01)(age_grp_and_hair_loss$Proportion),
    "<b>\nCount:</b> ", age_grp_and_hair_loss$Count
  ),
  marker = list(
    colors = case_when(
      age_grp_and_hair_loss$`Age Group` == "Young (17-30)" & age_grp_and_hair_loss$`Hair Loss` == 1 ~ colors[1],
      age_grp_and_hair_loss$`Age Group` == "Middle-aged (31-45)" & age_grp_and_hair_loss$`Hair Loss` == 1 ~ colors[2],
      age_grp_and_hair_loss$`Age Group` == "Old (>45)" & age_grp_and_hair_loss$`Hair Loss` == 1 ~ colors[3],
      age_grp_and_hair_loss$`Age Group` == "Young (17-30)" & age_grp_and_hair_loss$`Hair Loss` == 0 ~ colors[4],
      age_grp_and_hair_loss$`Age Group` == "Middle-aged (31-45)" & age_grp_and_hair_loss$`Hair Loss` == 0 ~ colors[5],
      age_grp_and_hair_loss$`Age Group` == "Old (>45)" & age_grp_and_hair_loss$`Hair Loss` == 0 ~ colors[6]
    )
  ),
  textfont = list(
    family = "ITC Officina Sans",
    size = 15,
    color = "white"
  ),
  textinfo = "label+text",
  textposition = "middle center"
) %>%
  layout(
    title = list(
      text = "\n𝗙𝗜𝗚. 𝟰﹕𝗧𝗥𝗘𝗘 𝗠𝗔𝗣 𝗢𝗙 𝗧𝗛𝗘 𝗦𝗨𝗥𝗩𝗘𝗬𝗘𝗗 𝗜𝗡𝗗𝗜𝗩𝗜𝗗𝗨𝗔𝗟𝗦 \n𝗕𝗬 𝗔𝗚𝗘 𝗚𝗥𝗢𝗨𝗣 𝗔𝗡𝗗 𝗛𝗔𝗜𝗥 𝗟𝗢𝗦𝗦 𝗢𝗨𝗧𝗖𝗢𝗠𝗘",
      font = list(size = 18)
    ),
    paper_bgcolor = "#D5E4EB",
    margin = list(l = 30, r = 30, b = 25, t = 45)
  ) %>%
  style(
    hoverlabel = list(
      font = list(
        family = "ITC Officina Sans",
        size = 15
      )
    )
  )

age_hair_loss_treemap

# Compute point-biserial correlation coefficient between Age and Hair Loss outcome
age_rpb <- compute_point_biserial(
  data1, 
  x = "Age", 
  y = "Hair Loss"
)

# Tabularize results 
summarize_variable(
  age_rpb[, c(1, 3:4, 7, 5, 6)] %>% column_to_rownames("x"), 
  rownames = TRUE
)

# Create contingency table for Stress by Hair Loss outcome
stress_table <- create_contingency_table(
  create_contingency_data(data, "Stress"), 
  "Stress"
)

# Data frame for the treemap, by Stress and Hair Loss
stress_lvl_and_hair_loss <- data %>%
	select(Stress, `Hair Loss`) %>%
	group_by(`Hair Loss`, `Stress`) %>%
	summarise(Count = n(), .groups="drop") %>%
	mutate(Proportion = Count/sum(Count)) %>%
	arrange(desc(`Hair Loss`), desc(Count))

# Factorize the Stress variable
stress_lvl_and_hair_loss$Stress <- factor(stress_lvl_and_hair_loss$Stress, levels=c("Low", "Moderate", "High"),
										  ordered=TRUE)

# Create a treemap for by stress level and hair loss experience
colors1 <- c("#567c07", "#996f00", "#a73735", "#77b200", "#e19c00", "#ec3330") 
stress_hair_loss_treemap <- plot_ly(
  type="treemap",
  labels=paste0("<b>",stringr::str_to_sentence(stress_lvl_and_hair_loss$`Stress`)," stress level",
				 "\n",ifelse(stress_lvl_and_hair_loss$`Hair Loss`==1, "with baldness", "without baldness"),"</b>"),
  parents=rep("", nrow(stress_lvl_and_hair_loss)),
  values=stress_lvl_and_hair_loss$Count,
  domain = list(column=0),
  name=" ",
  hoverinfo="text",
  hovertext=paste0(
    "<b>Percentage:</b> ", label_percent(accuracy=0.01)(stress_lvl_and_hair_loss$Count/sum(stress_lvl_and_hair_loss$Count)),
    "<b>\nCount: </b>", stress_lvl_and_hair_loss$Count
  ),
  marker = list(
    colors=case_when(
      stress_lvl_and_hair_loss$`Stress`=="Low" & stress_lvl_and_hair_loss$`Hair Loss`==1 ~ colors1[1],
      stress_lvl_and_hair_loss$`Stress`=="Moderate" & stress_lvl_and_hair_loss$`Hair Loss`==1 ~ colors1[2],
      stress_lvl_and_hair_loss$`Stress`=="High" & stress_lvl_and_hair_loss$`Hair Loss`==1 ~ colors1[3],
      stress_lvl_and_hair_loss$`Stress`=="Low" & stress_lvl_and_hair_loss$`Hair Loss`==0 ~ colors1[4],
      stress_lvl_and_hair_loss$`Stress`=="Moderate" & stress_lvl_and_hair_loss$`Hair Loss`==0 ~ colors1[5],
      stress_lvl_and_hair_loss$`Stress`=="High" & stress_lvl_and_hair_loss$`Hair Loss`==0 ~ colors1[6]
    )
  ),
  textfont = list(
    family="ITC Officina Sans",
    size=15,
    color="white"
  ),
  textinfo="label+text",
  textposition="middle center"
) %>%
  layout(
    title = list(
      text="\n𝗙𝗜𝗚. 𝟱﹕𝗧𝗥𝗘𝗘 𝗠𝗔𝗣 𝗢𝗙 𝗧𝗛𝗘 𝗦𝗨𝗥𝗩𝗘𝗬𝗘𝗗 𝗜𝗡𝗗𝗜𝗩𝗜𝗗𝗨𝗔𝗟𝗦\n𝗕𝗬 𝗦𝗧𝗥𝗘𝗦𝗦 𝗟𝗘𝗩𝗘𝗟 𝗔𝗡𝗗 𝗛𝗔𝗜𝗥 𝗟𝗢𝗦𝗦 𝗢𝗨𝗧𝗖𝗢𝗠𝗘",
      font = list(size=18)
    ),
    paper_bgcolor="#D5E4EB",
    margin = list(l=30, r=30, b=25, t=45)
  ) %>%
  style(
    hoverlabel = list(
      font = list(
        family="ITC Officina Sans",
		size=15
      )
    )
  )

stress_hair_loss_treemap