# 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/CP_b_u_50.xlsx"

# Read the data from the Excel file
blinded_data <- read_excel(file_path, sheet = "Blinded")
unblinded_data <- read_excel(file_path, sheet = "Unblinded")

# Add a column to indicate the source of the data
blinded_data$`Interim Type` <- "Blinded"
unblinded_data$`Interim Type` <- "Unblinded"

# Combine the data
combined_data <- bind_rows(blinded_data, unblinded_data)

# Filter the data
filtered_data <- combined_data[combined_data$`Treatment Rate Assumption` == 0.3, ]


filtered_data$`New N vs Old N Ratio` <- as.numeric(filtered_data$`New N vs Old N Ratio`)

# Create a grid for surface plot
x <- sort(unique(filtered_data$`Actual Treatment Rate`))
y <- sort(unique(filtered_data$`Actual Control Rate`))


# Create a matrix for the z-values (Significant at Final)
z_blinded <- matrix(NA, nrow = length(y), ncol = length(x))
z_unblinded <- matrix(NA, nrow = length(y), ncol = length(x))

for(i in seq_along(x)){
  for(j in seq_along(y)){
    z_blinded[j, i] <- filtered_data$`New N vs Old N Ratio`[filtered_data$`Interim Type` == "Blinded" & 
                                                              filtered_data$`Actual Treatment Rate` == x[i] & 
                                                              filtered_data$`Actual Control Rate` == y[j]]
    z_unblinded[j, i] <- filtered_data$`New N vs Old N Ratio`[filtered_data$`Interim Type` == "Unblinded" & 
                                                                filtered_data$`Actual Treatment Rate` == x[i] & 
                                                                filtered_data$`Actual Control Rate` == y[j]]
  }
}

# Convert matrices to numeric
z_blinded <- as.numeric(z_blinded)
z_unblinded <- as.numeric(z_unblinded)

z_diff <- z_blinded - z_unblinded

# Create the surface plot
fig1 <- plot_ly() %>%
  add_surface(x = ~x, y = ~y, z = matrix(z_blinded, nrow = length(y), ncol = length(x)), 
              colorscale = 'Blues', 
              name = 'Blinded') %>%
  layout(
    scene = list(
      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 = "New N vs Old N Ratio", range = c(1, 1.5), dtick = 0.05),
      camera = list(
        eye = list(x = 1.5, y = 1.5, z = 1.5)
      )
    ),
    title = list(
      text = "Blinded Interim\nNew and Original Sample Size Ratio vs Actual Rates\n(Assumption: 70% vs 30%; Interim at 50%)",
      x = 0.5,
      y = 0.95,
      xanchor = 'center',
      yanchor = 'top'
    ),
    margin = list(t = 10, r = 50, b = 50, l = 50), 
    width = 1200,
    height = 800
  )
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
# Save the plot
htmlwidgets::saveWidget(fig1, "C:/Users/szhuo/Downloads/PhD/final_CP/final_plots/blinded_70V30_N_t50_surface.html")

fig1
# Create the surface plot
fig2 <- plot_ly() %>%
  add_surface(x = ~x, y = ~y, z = matrix(z_unblinded, nrow = length(y), ncol = length(x)), 
              colorscale = 'Reds', 
              name = 'Unblinded') %>%
  layout(
    scene = list(
      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 = "New N vs Old N Ratio", range = c(1, 1.5), dtick = 0.05),
      camera = list(
        eye = list(x = 1.5, y = 1.5, z = 1.5)
      )
    ),
    title = list(
      text = "Unblinded Interim\nNew and Original Sample Size Ratio vs Actual Rates\n(Assumption: 70% vs 30%; Interim at 50%)",
      x = 0.5,
      y = 0.95,
      xanchor = 'center',
      yanchor = 'top'
    ),
    margin = list(t = 10, r = 50, b = 50, l = 50), 
    width = 1200,
    height = 800
  )
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
# Save the plot
htmlwidgets::saveWidget(fig2, "C:/Users/szhuo/Downloads/PhD/final_CP/final_plots/unblinded_70V30_N_t50_surface.html")

fig2
# Create the surface plot
fig3 <- plot_ly() %>%
  add_surface(x = ~x, y = ~y, z = matrix(z_diff, nrow = length(y), ncol = length(x)), 
              colorscale = 'Greens', 
              name = 'Difference') %>%
  layout(
    scene = list(
      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 = "Dfference of New N vs Old N Ratio", range = c(-0.5, 0.5), dtick = 0.05),
      camera = list(
        eye = list(x = 1.5, y = 1.5, z = 1.5)
      )
    ),
    title = list(
      text = "Difference of Blinded - Unblinded Interim\nNew and Original Sample Size Ratio vs Actual Rates\n(Assumption: 70% vs 30%; Interim at 50%)",
      x = 0.5,
      y = 0.95,
      xanchor = 'center',
      yanchor = 'top'
    ),
    margin = list(t = 10, r = 50, b = 50, l = 50), 
    width = 1200,
    height = 800
  )
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
# Save the plot
htmlwidgets::saveWidget(fig3, "C:/Users/szhuo/Downloads/PhD/final_CP/final_plots/diff_70V30_N_t50_surface.html")

fig3