Power Simulations.Rmd

---
title: "Thesis Power Simulations"
author: "Stephen Ogbodo"
date: "2024-03-19"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


```{r}
library(tidyverse) #The epidemiologist's paradise
library(scales)    #To beautify axis labels
```

# **My data**

## *Groups and outcome*
### Treatment (A=1): Aromatase inhibitor (AI)
### Comparator (A=0): Tamoxifen
### Outcome for objective 1: Diabetes (Dia)
### Outcome for objective 2: MACE

## *Hypotheses*
### Null hypothesis: HR = 1
### Alternative hypothesis: HR < 1 (treatment AI doing better than Tamoxifen)

## *Sample size*
```{r}
# Number of patients in each group
AI <- 85000
Tam <- 87500

# No linkage to HES and ONS required for objective 1. Only CPRD needed.

# Account for 5% censoring
AI_cens <- AI*0.95
Tam_cens <- Tam*0.95

#Further 5% lost to propensity trimming
AI_trimmed <- AI_cens*0.95
Tam_trimmed <- Tam_cens*0.95

# True power is based on the size of the smaller group.
# Calculate total size as double of AI group size
n_final <- AI_trimmed*2

# Print the final sample size
print(n_final)
```


# **Objective 1: Diabetes risk**

## *Simulation for HR 1.25*
```{r}
# Simulation 
possible.ns <- seq(from=80000, to=180000, by=10000) # The sample sizes we'll be considering (selected around final sample size of 153,425 women calculated above)
powers <- rep(NA, length(possible.ns))              # Empty object to collect simulation estimates
alpha <- 0.05                                       # Standard significance level
sims <- 500                                         # Number of simulations to conduct for each N

# Write function to run simulation on
# Diabetes, 5Y FUP, HR 0.75, rate of event in control 18.3 per 1000 person-year (Lipscombe 2013)
set.seed(12345)
for (j in 1:length(possible.ns)){
  N <- possible.ns[j]                              # Pick the jth value for N
  significant.experiments <- rep(NA, sims)         # Empty object to count significant experiments
  for (i in 1:sims){
    Y0 <- rbinom(n=N, 1, 0.0183*4)                 # control potential outcome rate * mean years of follow-up (5 minus 1 year to conservatively account for depletion of susceptible individuals over time.)
    tau <- 1.25                                    # Hypothesized treatment effect
    Y1 <- Y0*tau                                   # treatment potential outcome
    Z.sim <- rbinom(n=N, size=1, prob=.5)          # Do a random assignment
    Y.sim <- Y1*Z.sim + Y0*(1-Z.sim)               # Reveal outcomes according to assignment
    fit.sim <- lm(Y.sim ~ Z.sim)                   # Do analysis (Simple regression)
    p.value <- summary(fit.sim)$coefficients[2,4]  # Extract p-values
    significant.experiments[i] <- (p.value <= alpha) #Determine significance according to p <= 0.05
  }
  powers[j] <- mean(significant.experiments)       # Store average success rate (power) for each N
}
pow_1 <- data.frame(possible.ns, powers)
pow_1$fup <- as.factor(c(2))
pow_1$effect <- as.factor(c(1.25))
pow_1
```

## *Simulation for HR 1.20*
```{r}
# Simulation 
possible.ns <- seq(from=80000, to=180000, by=10000) # The sample sizes we'll be considering
powers <- rep(NA, length(possible.ns))              # Empty object to collect simulation estimates
alpha <- 0.05                                       # Standard significance level
sims <- 500                                         # Number of simulations to conduct for each N

# Write function to run simulation on
# Diabetes, 5Y FUP, HR 0.80, rate of event in control 18.3 per 1000 py (Lipscombe 2013)
set.seed(12345)
for (j in 1:length(possible.ns)){
  N <- possible.ns[j]                              # Pick the jth value for N
  significant.experiments <- rep(NA, sims)         # Empty object to count significant experiments
  for (i in 1:sims){
    Y0 <- rbinom(n=N, 1, 0.0183*4)                 # control potential outcome rate * mean years of follow-up (5 minus 1 year to conservatively account for depletion of susceptible individuals over time.)
    tau <- 1.20                                  # Hypothesize treatment effect
    Y1 <- Y0*tau                                   # treatment potential outcome
    Z.sim <- rbinom(n=N, size=1, prob=.5)          # Do a random assignment
    Y.sim <- Y1*Z.sim + Y0*(1-Z.sim)               # Reveal outcomes according to assignment
    fit.sim <- lm(Y.sim ~ Z.sim)                   # Do analysis (Simple regression)
    p.value <- summary(fit.sim)$coefficients[2,4]  # Extract p-values
    significant.experiments[i] <- (p.value <= alpha) # Determine significance according to p <= 0.05
  }
  powers[j] <- mean(significant.experiments)       # store average success rate (power) for each N
}
pow_2 <- data.frame(possible.ns, powers)
pow_2$fup <- as.factor(c(2))
pow_2$effect <- as.factor(c(1.20))
pow_2
```

## *Simulation for HR 1.15*
```{r}
# Simulation 
possible.ns <- seq(from=80000, to=180000, by=10000) # The sample sizes we'll be considering
powers <- rep(NA, length(possible.ns))              # Empty object to collect simulation estimates
alpha <- 0.05                                       # Standard significance level
sims <- 500                                         # Number of simulations to conduct for each N

# Write function to run simulation on
# Diabetes, 5Y FUP, HR 0.85, rate of event in control 18.3 per 1000 py (Lipscombe 2013)
set.seed(12345)
for (j in 1:length(possible.ns)){
  N <- possible.ns[j]                              # Pick the jth value for N
  significant.experiments <- rep(NA, sims)         # Empty object to count significant experiments
  for (i in 1:sims){
    Y0 <- rbinom(n=N, 1, 0.0183*4)                 # control potential outcome rate * mean years of follow-up (5 minus 1 year to conservatively account for depletion of susceptible individuals over time.)
    tau <- 1.15                                    # Hypothesize treatment effect
    Y1 <- Y0*tau                                   # treatment potential outcome
    Z.sim <- rbinom(n=N, size=1, prob=.5)          # Do a random assignment
    Y.sim <- Y1*Z.sim + Y0*(1-Z.sim)               # Reveal outcomes according to assignment
    fit.sim <- lm(Y.sim ~ Z.sim)                   # Do analysis (Simple regression)
    p.value <- summary(fit.sim)$coefficients[2,4]  # Extract p-values
    significant.experiments[i] <- (p.value <= alpha) # Determine significance according to p <= 0.05
  }
  powers[j] <- mean(significant.experiments)       # store average success rate (power) for each N
}
pow_3 <- data.frame(possible.ns, powers)
pow_3$fup <- as.factor(c(2))
pow_3$effect <- as.factor(c(1.15))
pow_3
```

## *Simulation for HR 1.10*
```{r}
# Simulation 
possible.ns <- seq(from=80000, to=180000, by=10000)   # The sample sizes we'll be considering
powers <- rep(NA, length(possible.ns))                # Empty object to collect simulation estimates
alpha <- 0.05                                         # Standard significance level
sims <- 500                                           # Number of simulations to conduct for each N

# Write function to run simulation on
# Diabetes, 5Y FUP, HR 0.90, rate of event in control 18.3 per 1000 py (Lipscombe 2013)
set.seed(12345)
for (j in 1:length(possible.ns)){
  N <- possible.ns[j]                                # Pick the jth value for N
  significant.experiments <- rep(NA, sims)           # Empty object to count significant experiments
  for (i in 1:sims){
    Y0 <- rbinom(n=N, 1, 0.0183*4)                   # control potential outcome rate * mean years of follow-up (5 minus 1 year to conservatively account for depletion of susceptible individuals over time.)
    tau <- 1.10                                      # Hypothesize treatment effect
    Y1 <- Y0*tau                                     # treatment potential outcome 
    Z.sim <- rbinom(n=N, size=1, prob=.5)            # Do a random assignment
    Y.sim <- Y1*Z.sim + Y0*(1-Z.sim)                 # Reveal outcomes according to assignment
    fit.sim <- lm(Y.sim ~ Z.sim)                     # Do analysis (Simple regression)
    p.value <- summary(fit.sim)$coefficients[2,4]    # Extract p-values
    significant.experiments[i] <- (p.value <= alpha) # Determine significance according to p <= 0.05
  }
  powers[j] <- mean(significant.experiments)         # store average success rate (power) for each N
}
pow_4 <- data.frame(possible.ns, powers)
pow_4$fup <- as.factor(c(2))
pow_4$effect <- as.factor(c(1.10))
pow_4
```

## *Simulation for HR 1.08*
```{r}
# Simulation 
possible.ns <- seq(from=80000, to=180000, by=10000)     # The sample sizes we'll be considering
powers <- rep(NA, length(possible.ns))                  # Empty object to collect simulation estimates
alpha <- 0.05                                           # Standard significance level
sims <- 500                                             # Number of simulations to conduct for each N

# Write function to run simulation on
# Diabetes, 5Y FUP, HR 0.93, rate of event in control 18.3 per 1000 py (Lipscombe 2013)
set.seed(12345)
for (j in 1:length(possible.ns)){
  N <- possible.ns[j]                              # Pick the jth value for N
  significant.experiments <- rep(NA, sims)         # Empty object to count significant experiments
  for (i in 1:sims){
    Y0 <- rbinom(n=N, 1, 0.0183*4)                 # control potential outcome rate * mean years of follow-up (5 minus 1 year to conservatively account for depletion of susceptible individuals over time.)
    tau <- 1.08                                   # Hypothesize treatment effect
    Y1 <- Y0*tau                                   # treatment potential outcome
    Z.sim <- rbinom(n=N, size=1, prob=.5)          # Do a random assignment
    Y.sim <- Y1*Z.sim + Y0*(1-Z.sim)               # Reveal outcomes according to assignment
    fit.sim <- lm(Y.sim ~ Z.sim)                   # Do analysis (Simple regression)
    p.value <- summary(fit.sim)$coefficients[2,4]  # Extract p-values
    significant.experiments[i] <- (p.value <= alpha) # Determine significance according to p <= 0.05
  }
  powers[j] <- mean(significant.experiments)       # store average success rate (power) for each N
}
pow_5 <- data.frame(possible.ns, powers)
pow_5$fup <- as.factor(c(2))
pow_5$effect <- as.factor(c(1.08))
pow_5
```

## *Simulation for HR 1.05*
```{r}
# Simulation 
possible.ns <- seq(from=80000, to=180000, by=10000)     # The sample sizes we'll be considering
powers <- rep(NA, length(possible.ns))                  # Empty object to collect simulation estimates
alpha <- 0.05                                           # Standard significance level
sims <- 500                                             # Number of simulations to conduct for each N

# Write function to run simulation on
# Diabetes, 5Y FUP, HR 0.95, rate of event in control 18.3 per 1000 py (Lipscombe 2013)
set.seed(12345)
for (j in 1:length(possible.ns)){
  N <- possible.ns[j]                              # Pick the jth value for N
  significant.experiments <- rep(NA, sims)         # Empty object to count significant experiments
  for (i in 1:sims){
    Y0 <- rbinom(n=N, 1, 0.0183*4)                 # control potential outcome rate * mean years of follow-up (minus 1 year to conservatively account for depletion of susceptible individuals over time.)
    tau <- 1.05                                    # Hypothesize treatment effect
    Y1 <- Y0*tau                                   # treatment potential outcome
    Z.sim <- rbinom(n=N, size=1, prob=.5)          # Do a random assignment
    Y.sim <- Y1*Z.sim + Y0*(1-Z.sim)               # Reveal outcomes according to assignment
    fit.sim <- lm(Y.sim ~ Z.sim)                   # Do analysis (Simple regression)
    p.value <- summary(fit.sim)$coefficients[2,4]  # Extract p-values
    significant.experiments[i] <- (p.value <= alpha) # Determine significance according to p <= 0.05
  }
  powers[j] <- mean(significant.experiments)       # store average success rate (power) for each N
}
pow_6 <- data.frame(possible.ns, powers)
pow_6$fup <- as.factor(c(2))
pow_6$effect <- as.factor(c(1.05))
pow_6
```

## *(Optional) Add some noise to perfect powers to separate them out for visualization (limit absolute overlap)*
```{r}
pow_2n<-
  pow_2%>%
  mutate(power_n = case_when(possible.ns==80000~powers-0.01, possible.ns==90000~powers-0.009, possible.ns==100000~powers-0.008, possible.ns==110000~powers-0.007, possible.ns==120000~powers-0.006, possible.ns==130000~powers-0.005, possible.ns==140000~powers-0.004, possible.ns==150000~powers-0.003, possible.ns==160000~powers-0.002, possible.ns==170000~powers-0.001, .default = powers))

pow_3n<-
  pow_3%>%
  mutate(power_n = case_when(possible.ns==80000~powers-0.015, possible.ns==90000~powers-0.0135, possible.ns==100000~powers-0.012, possible.ns==110000~powers-0.0105, possible.ns==120000~powers-0.009, possible.ns==130000~powers-0.0075, possible.ns==140000~powers-0.006, possible.ns==150000~powers-0.0045, possible.ns==160000~powers-0.003, possible.ns==170000~powers-0.0015, .default = powers))

pow_4n<-
  pow_4%>%
  mutate(power_n = case_when(possible.ns==80000~powers-0.02, possible.ns==90000~powers-0.018, possible.ns==100000~powers-0.016, possible.ns==110000~powers-0.014, possible.ns==120000~powers-0.012, possible.ns==130000~powers-0.010, possible.ns==140000~powers-0.008, possible.ns==150000~powers-0.006, possible.ns==160000~powers-0.004, possible.ns==170000~powers-0.002, .default = powers))

```


## *Graphical representation*
```{r}
# Plot power curves
ggplot() + 
  geom_line(data = pow_1, aes(x = possible.ns, y = powers, color = "a) HR=1.25"), linewidth=1)+
  geom_smooth(data = pow_2n, aes(x = possible.ns, y = power_n, color = "b) HR=1.20"), 
              se = FALSE, method = "gam", formula = y ~ s(log(x), k = 4)) +
  geom_smooth(data = pow_3n, aes(x = possible.ns, y = power_n, color = "c) HR=1.15"), 
              se = FALSE, method = "gam", formula = y ~ s(log(x), k = 4)) +
  geom_smooth(data = pow_4n, aes(x = possible.ns, y = power_n, color = "d) HR=1.10"), 
              se = FALSE, method = "gam", formula = y ~ s(log(x), k = 4)) +
  geom_smooth(data = pow_5, aes(x = possible.ns, y = powers, color = "e) HR=1.08"), 
              se = FALSE, method = "gam", formula = y ~ s(log(x), k = 4)) +
  geom_smooth(data = pow_6, aes(x = possible.ns, y = powers, color = "f)  HR=1.05"), 
              se = FALSE, method = "gam", formula = y ~ s(log(x), k = 4)) +
  labs(x = "Study Sample Size", y = "Power", title = "Power for Objective II (5-year follow-up)") + 
  theme_classic() +  
  theme(plot.title = element_text(hjust = 0.5),  
        panel.grid = element_blank()) +
  scale_x_continuous(labels = scales::comma) +
  scale_y_continuous(labels=scales::percent) +
  geom_hline(yintercept = 0.8, linetype = "dashed", color = "black") +
  geom_vline(xintercept = 153425, linetype = "dashed", color = "black") +
  scale_color_manual(name = "Effect Size", 
                     values = c("a) HR=1.25" = "blue", "b) HR=1.20" = "purple", "c) HR=1.15" = "green", "d) HR=1.10" = "brown", "e) HR=1.08" = "darkgreen" , "f)  HR=1.05" = "darkgrey"))

ggsave("fig1.jpg", units="in", width=10, height=7, dpi=300)
```



# **Objective 2: MACE outcome among women with comorbid diabetes at baseline**

## *Sample size*
```{r}
# Number of patients in each group
AI <- 15500
Tam <- 17000

# Participants that are linkable to HES and ONS (75%)
AI_linked <- AI*0.75
Tam_linked <- Tam*0.75

# Account for 5% censoring
AI_cens <- AI_linked*0.95
Tam_cens <- Tam_linked*0.95

#Further 5% lost to trimming
AI_trimmed <- AI_cens*0.95
Tam_trimmed <- Tam_cens*0.95

# True power is based on the size of the smaller group.
# Calculate total size as double of AI group size
n_final <- AI_trimmed*2

# Print the data frame
print(n_final)
```


## *Simulation for HR 1.25*
```{r}
# Simulation 
possible.ns <- seq(from=15000, to=25000, by=1000) # The sample sizes we'll be considering (selected around final sample size of 20,983 women calculated above)
powers <- rep(NA, length(possible.ns))           # Empty object to collect simulation estimates
alpha <- 0.05                                    # Standard significance level
sims <- 500                                      # Number of simulations to conduct for each N

# Write function to run simulation on
# MACE, 5Y FUP, HR 0.75, rate of event in control 60 per 1000 person-year (Angoulvant et al, 2021)
set.seed(12345)
for (j in 1:length(possible.ns)){
  N <- possible.ns[j]                              # Pick the jth value for N
  significant.experiments <- rep(NA, sims)         # Empty object to count significant experiments
  for (i in 1:sims){
    Y0 <- rbinom(n=N, 1, 0.06*4)                   # control potential outcome rate * mean years of follow-up (5 minus 1 year to conservatively account for depletion of susceptible individuals over time.)
    tau <- 1.25                                    # Hypothesize treatment effect
    Y1 <- Y0*tau                                   # treatment potential outcome
    Z.sim <- rbinom(n=N, size=1, prob=.5)          # Do a random assignment
    Y.sim <- Y1*Z.sim + Y0*(1-Z.sim)               # Reveal outcomes according to assignment
    fit.sim <- lm(Y.sim ~ Z.sim)                   # Do analysis (Simple regression)
    p.value <- summary(fit.sim)$coefficients[2,4]  # Extract p-values
    significant.experiments[i] <- (p.value <= alpha) # Determine significance according to p <= 0.05
  }
  powers[j] <- mean(significant.experiments)       # store average success rate (power) for each N
}
MACE1 <- data.frame(possible.ns, powers)
MACE1$fup <- as.factor(c(2))
MACE1$effect <- as.factor(c(1.25))
MACE1
```

## *Simulation for HR 1.20*
```{r}
# Simulation 
possible.ns <- seq(from=15000, to=25000, by=1000) # The sample sizes we'll be considering
powers <- rep(NA, length(possible.ns))            # Empty object to collect simulation estimates
alpha <- 0.05                                     # Standard significance level
sims <- 500                                       # Number of simulations to conduct for each N

# Write function to run simulation on
# MACE, 5Y FUP, HR 0.80, rate of event in control 60 per 1000 py (Angoulvant et al, 2021)
set.seed(12345)
for (j in 1:length(possible.ns)){
  N <- possible.ns[j]                              # Pick the jth value for N
  significant.experiments <- rep(NA, sims)         # Empty object to count significant experiments
  for (i in 1:sims){
    Y0 <- rbinom(n=N, 1, 0.06*4)                   # control potential outcome rate * mean years of follow-up (5 minus 1 year to conservatively account for depletion of susceptible individuals over time.)
    tau <- 1.20                                    # Hypothesize treatment effect
    Y1 <- Y0*tau                                   # treatment potential outcome
    Z.sim <- rbinom(n=N, size=1, prob=.5)          # Do a random assignment
    Y.sim <- Y1*Z.sim + Y0*(1-Z.sim)               # Reveal outcomes according to assignment
    fit.sim <- lm(Y.sim ~ Z.sim)                   # Do analysis (Simple regression)
    p.value <- summary(fit.sim)$coefficients[2,4]  # Extract p-values
    significant.experiments[i] <- (p.value <= alpha) # Determine significance according to p <= 0.05
  }
  powers[j] <- mean(significant.experiments)       # store average success rate (power) for each N
}
MACE2 <- data.frame(possible.ns, powers)
MACE2$fup <- as.factor(c(2))
MACE2$effect <- as.factor(c(1.20))
MACE2
```

## *Simulation for HR 1.15*
```{r}
# Simulation 
possible.ns <- seq(from=15000, to=25000, by=1000) # The sample sizes we'll be considering
powers <- rep(NA, length(possible.ns))            # Empty object to collect simulation estimates
alpha <- 0.05                                     # Standard significance level
sims <- 500                                       # Number of simulations to conduct for each N

# Write function to run simulation on
# MACE, 5Y FUP, HR 0.85, rate of event in control 60 per 1000 py (Angoulvant et al, 2021)
set.seed(12345)
for (j in 1:length(possible.ns)){
  N <- possible.ns[j]                              # Pick the jth value for N
  significant.experiments <- rep(NA, sims)         # Empty object to count significant experiments
  for (i in 1:sims){
    Y0 <- rbinom(n=N, 1, 0.06*4)                   # control potential outcome rate * mean years of follow-up (5 minus 1 year to conservatively account for depletion of susceptible individuals over time.)
    tau <- 1.15                                    # Hypothesize treatment effect
    Y1 <- Y0*tau                                   # treatment potential outcome
    Z.sim <- rbinom(n=N, size=1, prob=.5)          # Do a random assignment
    Y.sim <- Y1*Z.sim + Y0*(1-Z.sim)               # Reveal outcomes according to assignment
    fit.sim <- lm(Y.sim ~ Z.sim)                   # Do analysis (Simple regression)
    p.value <- summary(fit.sim)$coefficients[2,4]  # Extract p-values
    significant.experiments[i] <- (p.value <= alpha) # Determine significance according to p <= 0.05
  }
  powers[j] <- mean(significant.experiments)       # store average success rate (power) for each N
}
MACE3 <- data.frame(possible.ns, powers)
MACE3$fup <- as.factor(c(2))
MACE3$effect <- as.factor(c(1.15))
MACE3
```

## *Simulation for HR 1.10*
```{r}
# Simulation 
possible.ns <- seq(from=15000, to=25000, by=1000)  # The sample sizes we'll be considering
powers <- rep(NA, length(possible.ns))             # Empty object to collect simulation estimates
alpha <- 0.05                                      # Standard significance level
sims <- 500                                        # Number of simulations to conduct for each N

# Write function to run simulation on
# MACE, 5Y FUP, HR 0.90, rate of event in control 90 per 1000 py (Angoulvant et al, 2021)
set.seed(12345)
for (j in 1:length(possible.ns)){
  N <- possible.ns[j]                              # Pick the jth value for N
  significant.experiments <- rep(NA, sims)         # Empty object to count significant experiments
  for (i in 1:sims){
    Y0 <- rbinom(n=N, 1, 0.06*4)                   # control potential outcome rate * mean years of follow-up (5 minus 1 year to conservatively account for depletion of susceptible individuals over time.)
    tau <- 1.10                                   # Hypothesize treatment effect
    Y1 <- Y0*tau                                   # treatment potential outcome
    Z.sim <- rbinom(n=N, size=1, prob=.5)          # Do a random assignment
    Y.sim <- Y1*Z.sim + Y0*(1-Z.sim)               # Reveal outcomes according to assignment
    fit.sim <- lm(Y.sim ~ Z.sim)                   # Do analysis (Simple regression)
    p.value <- summary(fit.sim)$coefficients[2,4]  # Extract p-values
    significant.experiments[i] <- (p.value <= alpha) # Determine significance according to p <= 0.05
  }
  powers[j] <- mean(significant.experiments)       # store average success rate (power) for each N
}
MACE4 <- data.frame(possible.ns, powers)
MACE4$fup <- as.factor(c(3))
MACE4$effect <- as.factor(c(1.10))
MACE4
```

## *Simulation for HR 1.08*
```{r}
# Simulation 
possible.ns <- seq(from=15000, to=25000, by=1000)  # The sample sizes we'll be considering
powers <- rep(NA, length(possible.ns))             # Empty object to collect simulation estimates
alpha <- 0.05                                      # Standard significance level
sims <- 500                                        # Number of simulations to conduct for each N

# Write function to run simulation on
# MACE, 5Y FUP, HR 0.93, rate of event in control 60 per 1000 person-years (Angoulvant et al, 2021)
set.seed(12345)
for (j in 1:length(possible.ns)){
  N <- possible.ns[j]                              # Pick the jth value for N
  significant.experiments <- rep(NA, sims)         # Empty object to count significant experiments
  for (i in 1:sims){
    Y0 <- rbinom(n=N, 1, 0.06*4)                   # control potential outcome rate * mean years of follow-up (5 minus 1 year to conservatively account for depletion of susceptible individuals over time.)
    tau <- 1.08                                   # Hypothesize treatment effect
    Y1 <- Y0*tau                                   # treatment potential outcome
    Z.sim <- rbinom(n=N, size=1, prob=.5)          # Do a random assignment
    Y.sim <- Y1*Z.sim + Y0*(1-Z.sim)               # Reveal outcomes according to assignment
    fit.sim <- lm(Y.sim ~ Z.sim)                   # Do analysis (Simple regression)
    p.value <- summary(fit.sim)$coefficients[2,4]  # Extract p-values
    significant.experiments[i] <- (p.value <= alpha) # Determine significance according to p <= 0.05
  }
  powers[j] <- mean(significant.experiments)       # store average success rate (power) for each N
}
MACE5 <- data.frame(possible.ns, powers)
MACE5$fup <- as.factor(c(3))
MACE5$effect <- as.factor(c(1.08))
MACE5
```


## *Simulation for HR 1.05*
```{r}
# Simulation 
possible.ns <- seq(from=15000, to=25000, by=1000) # The sample sizes we'll be considering
powers <- rep(NA, length(possible.ns))            # Empty object to collect simulation estimates
alpha <- 0.05                                     # Standard significance level
sims <- 500                                       # Number of simulations to conduct for each N

# Write function to run simulation on
# MACE, 5Y FUP, HR 0.95, rate of event in control 60 per 1000 py (Angoulvant et al, 2021)
set.seed(12345)
for (j in 1:length(possible.ns)){
  N <- possible.ns[j]                              # Pick the jth value for N
  significant.experiments <- rep(NA, sims)         # Empty object to count significant experiments
  for (i in 1:sims){
    Y0 <- rbinom(n=N, 1, 0.06*4)                   # control potential outcome rate * mean years of follow-up (5 minus 1 year to conservatively account for depletion of susceptible individuals over time.)
    tau <- 1.05                                    # Hypothesize treatment effect
    Y1 <- Y0*tau                                   # treatment potential outcome
    Z.sim <- rbinom(n=N, size=1, prob=.5)          # Do a random assignment
    Y.sim <- Y1*Z.sim + Y0*(1-Z.sim)               # Reveal outcomes according to assignment
    fit.sim <- lm(Y.sim ~ Z.sim)                   # Do analysis (Simple regression)
    p.value <- summary(fit.sim)$coefficients[2,4]  # Extract p-values
    significant.experiments[i] <- (p.value <= alpha) # Determine significance according to p <= 0.05
  }
  powers[j] <- mean(significant.experiments)       # store average success rate (power) for each N
}
MACE6 <- data.frame(possible.ns, powers)
MACE6$fup <- as.factor(c(2))
MACE6$effect <- as.factor(c(1.05))
MACE6
```

## * (Optional) Add some noise to perfect powers to separate them out for visualization (limit absolute overlap)*
```{r}
library(tidyverse)
MACE2n<-
  MACE2%>%
  mutate(power_n = case_when(possible.ns==15000~powers-0.01, possible.ns==16000~powers-0.009, possible.ns==17000~powers-0.008, possible.ns==18000~powers-0.007, possible.ns==19000~powers-0.006, possible.ns==20000~powers-0.005, possible.ns==21000~powers-0.004, possible.ns==22000~powers-0.003, possible.ns==23000~powers-0.002, possible.ns==24000~powers-0.001, .default = powers))

MACE3n<-
  MACE3%>%
  mutate(power_n = case_when(possible.ns==15000~powers-0.015, possible.ns==16000~powers-0.0135, possible.ns==17000~powers-0.012, possible.ns==18000~powers-0.0105, possible.ns==19000~powers-0.009, possible.ns==20000~powers-0.0075, possible.ns==21000~powers-0.006, possible.ns==22000~powers-0.0045, possible.ns==23000~powers-0.003, possible.ns==24000~powers-0.0015, .default = powers))

MACE4n<-
  MACE4%>%
  mutate(power_n = case_when(possible.ns==15000~powers-0.02, possible.ns==16000~powers-0.018, possible.ns==17000~powers-0.016, possible.ns==18000~powers-0.014, possible.ns==19000~powers-0.012, possible.ns==20000~powers-0.010, possible.ns==21000~powers-0.008, possible.ns==22000~powers-0.006, possible.ns==23000~powers-0.004, possible.ns==24000~powers-0.002, .default = powers))

```


## *Graphical representation*
```{r}
# Plot power curves
ggplot() + 
  geom_line(data = MACE1, aes(x = possible.ns, y = powers, color = "a) HR=1.25"), linewidth = 1) +
  geom_smooth(data = MACE2n, aes(x = possible.ns, y = power_n, color = "b) HR=1.20"), 
              se = FALSE, method = "gam", formula = y ~ s(log(x), k = 4)) +
  geom_smooth(data = MACE3n, aes(x = possible.ns, y = power_n, color = "c) HR=1.15"), 
              se = FALSE, method = "gam", formula = y ~ s(log(x), k = 4)) +
  geom_smooth(data = MACE4, aes(x = possible.ns, y = powers, color = "d) HR=1.10"), 
              se = FALSE, method = "gam", formula = y ~ s(log(x), k = 4)) +
  geom_smooth(data = MACE5, aes(x = possible.ns, y = powers, color = "e) HR=1.08"), 
              se = FALSE, method = "gam", formula = y ~ s(log(x), k = 4)) +
  geom_smooth(data = MACE6, aes(x = possible.ns, y = powers, color = "f) HR=1.05"), 
              se = FALSE, method = "gam", formula = y ~ s(log(x), k = 4)) +
  labs(x = "Study Sample Size", y = "Power", title = "Power for Objective III (5-year follow-up)") + 
  theme_classic() +  
  theme(plot.title = element_text(hjust = 0.5),  
        panel.grid = element_blank()) +
  scale_x_continuous(labels = scales::comma) +
  scale_y_continuous(labels=scales::percent) +
  geom_hline(yintercept = 0.8, linetype = "dashed", color = "black") +
  geom_vline(xintercept = 20983, linetype = "dashed", color = "black") +
  scale_color_manual(name = "Effect Size", 
                     values = c("a) HR=1.25" = "blue", "b) HR=1.20" = "purple", "c) HR=1.15" = "green", "d) HR=1.10" = "brown", "e) HR=1.08" = "darkgreen" , "f) HR=1.05" = "darkgrey"))

ggsave("fig2.jpg", units="in", width=10, height=7, dpi=300)
```