# Load necessary libraries
library(readxl)
## Warning: package 'readxl' was built under R version 4.2.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.2.3
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(plotly)
## Warning: package 'plotly' was built under R version 4.2.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.2.3
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
# Define file path
file_path <- "C:/Users/szhuo/Downloads/PhD/final_CP/final_data/HR_CP_b_u_t50_type1_error_10000.xlsx"

# Read in the data
combined_data <- read_excel(file_path, sheet = "Overall")

# Define the function to create, display, and print surface plots
create_surface_plot <- function(x, y, z_matrix, plot_name, colorscale, z_title, treatment_rate) {
  # Calculate the z-axis range and extend it by 10% on both sides
  z_range <- range(z_matrix, na.rm = TRUE)
  z_margin <- 0.1 * (z_range[2] - z_range[1])
  z_axis_range <- c(z_range[1] - z_margin, z_range[2] + z_margin)
  
  # Calculate dynamic dtick for the z-axis based on the data range (e.g., 10% of the range)
  z_dtick <- (z_range[2] - z_range[1]) / 10
  
  # Find the maximum value and its coordinates
  max_value <- max(z_matrix, na.rm = TRUE)
  max_coords <- which(z_matrix == max_value, arr.ind = TRUE)
  max_x <- x[max_coords[2]]  # Corresponding x value
  max_y <- y[max_coords[1]]  # Corresponding y value
  
  # Create surface plot and remove the color scale (legend)
  fig <- plot_ly() %>% 
    add_surface(x = ~x, y = ~y, z = z_matrix, 
                colorscale = colorscale, 
                name = plot_name,
                showscale = FALSE) %>%  # Remove the legend bar
    layout(
      scene = list(
        aspectmode = "manual",  # Use manual aspect ratio
        aspectratio = list(x = 1, y = 1, z = 0.7),  # Adjust the aspect ratio
        xaxis = list(title = "Actual Treatment Rate", range = c(0.1, 0.9), dtick = 0.05),
        yaxis = list(title = "Actual Control Rate", range = c(0.1, 0.9), dtick = 0.05),
        zaxis = list(title = z_title, range = z_axis_range, dtick = z_dtick, tickformat = ".3f"),  # Use dynamic dtick and format to 3 decimals
        camera = list(
          eye = list(x = 1.5, y = 1.5, z = 1.5)
        )
      ),
      annotations = list(
        list(
          x = max_x, 
          y = max_y, 
          z = max_value, 
          text = paste("Max:", round(max_value, 3)),
          showarrow = TRUE,
          arrowhead = 4,
          ax = 0,
          ay = -40,
          font = list(color = "black")
        )
      ),
      title = list(
        text = paste0(plot_name, " Interim\nAssumed Control and Treatment Rate: 0.7 vs ", treatment_rate),
        x = 0.5,
        y = 0.95,
        xanchor = 'center',
        yanchor = 'bottom'
      ),
      margin = list(b=100),  # Increase margins
      width = 680,  # Set width
      height = 680  # Set height
    )
  
  return(fig)
}

# Function to generate surface plots for a given treatment rate
data_fig <- function(treatment_rate) {
  
  # Filter the data
  filtered_data <- combined_data %>% 
    filter(`Treatment Rate Assumption` == treatment_rate, `Type` != "Logrank")
  
  # Ensure the required columns are numeric
  filtered_data$`New N vs Old N Ratio` <- as.numeric(filtered_data$`New N vs Old N Ratio`)
  filtered_data$`Significant at Final` <- as.numeric(filtered_data$`Significant at Final`)
  filtered_data$`Probability of Increasing Sample Size` <- as.numeric(filtered_data$`Probability of Increasing Sample Size`)
  
  # Create grids for surface plot
  x <- sort(unique(filtered_data$`Actual Treatment Rate`))
  y <- sort(unique(filtered_data$`Actual Control Rate`))
  
  # Create z-matrices for blinded, unblinded, and the difference for each variable
  z_blinded_N <- matrix(filtered_data$`New N vs Old N Ratio`[filtered_data$`Type` == "Blinded"], 
                        nrow = length(y), ncol = length(x), byrow = TRUE)
  z_unblinded_N <- matrix(filtered_data$`New N vs Old N Ratio`[filtered_data$`Type` == "Unblinded"], 
                          nrow = length(y), ncol = length(x), byrow = TRUE)
  z_diff_N <- z_blinded_N - z_unblinded_N
  
  z_blinded_final <- matrix(filtered_data$`Significant at Final`[filtered_data$`Type` == "Blinded"], 
                            nrow = length(y), ncol = length(x), byrow = TRUE)
  z_unblinded_final <- matrix(filtered_data$`Significant at Final`[filtered_data$`Type` == "Unblinded"], 
                              nrow = length(y), ncol = length(x), byrow = TRUE)
  z_diff_final <- z_blinded_final - z_unblinded_final
  
  z_blinded_prob <- matrix(filtered_data$`Probability of Increasing Sample Size`[filtered_data$`Type` == "Blinded"], 
                           nrow = length(y), ncol = length(x), byrow = TRUE)
  z_unblinded_prob <- matrix(filtered_data$`Probability of Increasing Sample Size`[filtered_data$`Type` == "Unblinded"], 
                             nrow = length(y), ncol = length(x), byrow = TRUE)
  z_diff_prob <- z_blinded_prob - z_unblinded_prob
  
  # Generate plots for each variable
  fig_blinded_N <- create_surface_plot(x, y, z_blinded_N, "Blinded", 'Blues', "New N vs Old N Ratio", treatment_rate)
  fig_unblinded_N <- create_surface_plot(x, y, z_unblinded_N, "Unblinded", 'Reds', "New N vs Old N Ratio", treatment_rate)
  fig_diff_N <- create_surface_plot(x, y, z_diff_N, "Difference", 'Greens', "Difference of New N vs Old N Ratio", treatment_rate)
  
  fig_blinded_final <- create_surface_plot(x, y, z_blinded_final, "Blinded", 'Blues', "Significant at Final", treatment_rate)
  fig_unblinded_final <- create_surface_plot(x, y, z_unblinded_final, "Unblinded", 'Reds', "Significant at Final", treatment_rate)
  fig_diff_final <- create_surface_plot(x, y, z_diff_final, "Difference", 'Greens', "Difference of Significant at Final", treatment_rate)
  
  fig_blinded_prob <- create_surface_plot(x, y, z_blinded_prob, "Blinded", 'Blues', "Probability of Increasing Sample Size", treatment_rate)
  fig_unblinded_prob <- create_surface_plot(x, y, z_unblinded_prob, "Unblinded", 'Reds', "Probability of Increasing Sample Size", treatment_rate)
  fig_diff_prob <- create_surface_plot(x, y, z_diff_prob, "Difference", 'Greens', "Difference of Probability of Increasing Sample Size", treatment_rate)
  
  # Return all nine figures separately
  list(
    fig_blinded_N, fig_unblinded_N, fig_diff_N,
    fig_blinded_final, fig_unblinded_final, fig_diff_final,
    fig_blinded_prob, fig_unblinded_prob, fig_diff_prob
  )
}

# Generate plots for specific treatment rates
fig_list_0_2 <- data_fig(treatment_rate = 0.2)
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
fig_list_0_3 <- data_fig(treatment_rate = 0.3)
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
fig_list_0_4 <- data_fig(treatment_rate = 0.4)
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
fig_list_0_5 <- data_fig(treatment_rate = 0.5)
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
fig_list_0_6 <- data_fig(treatment_rate = 0.6)
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
# In the R Markdown file, output the figures explicitly for New N vs Old N Ratio, Significant at Final, and Probability of Increasing Sample Size

# Example for treatment rate 0.2
fig_list_0_2[[1]]  # Blinded plot for New N vs Old N Ratio for treatment rate 0.2
fig_list_0_2[[2]]  # Unblinded plot for New N vs Old N Ratio for treatment rate 0.2
fig_list_0_2[[3]]  # Difference plot for New N vs Old N Ratio for treatment rate 0.2
fig_list_0_2[[4]]  # Blinded plot for Significant at Final for treatment rate 0.2
fig_list_0_2[[5]]  # Unblinded plot for Significant at Final for treatment rate 0.2
fig_list_0_2[[6]]  # Difference plot for Significant at Final for treatment rate 0.2
fig_list_0_2[[7]]  # Blinded plot for Probability of Increasing Sample Size for treatment rate 0.2
fig_list_0_2[[8]]  # Unblinded plot for Probability of Increasing Sample Size for treatment rate 0.2
fig_list_0_2[[9]]  # Difference plot for Probability of Increasing Sample Size for treatment rate 0.2 
# # Repeat for treatment rates
# fig_list_0_3[[1]] 
# fig_list_0_3[[2]]  
# fig_list_0_3[[3]]  
# fig_list_0_3[[4]] 
# fig_list_0_3[[5]]  
# fig_list_0_3[[6]]  
# fig_list_0_3[[7]] 
# fig_list_0_3[[8]]  
# fig_list_0_3[[9]]  
# 
# fig_list_0_4[[1]] 
# fig_list_0_4[[2]]  
# fig_list_0_4[[3]]  
# fig_list_0_4[[4]] 
# fig_list_0_4[[5]]  
# fig_list_0_4[[6]]  
# fig_list_0_4[[7]] 
# fig_list_0_4[[8]]  
# fig_list_0_4[[9]]  
# 
# fig_list_0_5[[1]]
# fig_list_0_5[[2]]
# fig_list_0_5[[3]]
# fig_list_0_5[[4]]
# fig_list_0_5[[5]]
# fig_list_0_5[[6]]
# fig_list_0_5[[7]]
# fig_list_0_5[[8]]
# fig_list_0_5[[9]]
# 
# fig_list_0_6[[1]]
# fig_list_0_6[[2]]
# fig_list_0_6[[3]]
# fig_list_0_6[[4]]
# fig_list_0_6[[5]]
# fig_list_0_6[[6]]
# fig_list_0_6[[7]]
# fig_list_0_6[[8]]
# fig_list_0_6[[9]]