---
title: Impact of COVID-19 Pandemic on Oncology Clinical Trials - Simulation Results for Time Off-Treatment (Section 3.4) and Missed Visits (Section 3.3) on Efficacy Assessment
author:   | 
          |
          | Jiabu Ye, Binbing Yu, Zsolt Szijgyarto and Helen Mann
          | for the COVID-19 Oncology Biometrics Working Group
          | AstraZeneca Oncology Biometrics
date: "July 10, 2021"
output: slidy_presentation
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE,fig.align = "center", fig.width = 12)
```

<!-- # Example simulation ouput for assessing impact of 2-missed-visit rules (missing data) -->

<!-- - We consider 5000 simulations for each of the two settings. Each simulation takes about 1 min. We have used internal Scientific Computing Cluster (SCP) to generate the simulation results -->

<!-- - To demonstrate reproducibility, we output the simulation results for /output/MissingVisits/MissingVisitSimSetting2/MissingVisitSimout1.csv (See Line 229-235 in the R markdown file COVID19 Supplementary Tables and Figures v2.Rmd) -->

```{r,warning=FALSE, message=FALSE, echo=FALSE}
library(coin)
library(kableExtra)
library(flexsurvcure)
library(xlsx)
library(survival)
#library(survminer)
library(ggplot2)
library(doBy)
library(gsDesign)
library(PWEALL)
#library(LogrankPower)

#library(bindata)
setwd("C:/Users/klrd547/AZCollaboration/COVID-19 Stats Workstream - Documents/Simmulations/Simulation scripts/JDS Code")

# Simulate survival data according to current hypothetical design
sim.data<-function(sample.size,medS,HR,rr,study.cutoff) {
  
   # sample.size<-530
   # medS<-26
   # HR<-0.69
   # target.event<-233
   # rr<-1 
   # study.cutoff<-45
   
   accrual.duration<-29 
   accrual.weight<-2    # quadratic accrual
   lambda0 = log(2)/medS
   lambda1 = lambda0*HR
   entry.time <- runif(sample.size, 0, 1)
   entry.time <- (((accrual.duration)^accrual.weight)*entry.time)^(1/accrual.weight)
  
   n0<-round(sample.size/(rr+1))
   n1<-sample.size-n0
   event.time<-c(rexp(n0,lambda0),rexp(n1,lambda1))
   trt<-c(rep(0,n0),rep(1,n1))
   caltime<-entry.time+event.time
   Time<-pmin(study.cutoff-entry.time,event.time)
   Event<-as.numeric(event.time<=(study.cutoff-entry.time))  
   SUBJID<-1:sample.size
   data.frame(SUBJID,trt,event.time,entry.time,caltime,Time,Event) # 0 control 1 treatment
}
# number of consecutive NA
num.consecutive.NA<-function(x) {
  r <- rle(x==1)
  if (sum(x==1)==0) num.NA<-0 else num.NA<-max(r$lengths[r$values])
  num.NA
}

# Introducing missing visits to the simulated survival data
simulate.missing.visit<-function(sdata,study.cutoff,pandemic.parameter,fix.target.duration) {
  # study.cutoff<-45
  # pandemic.parameter<-pandemic.assumption[14,]
  pandemic.parameter<-as.numeric(pandemic.parameter)
  Pan.before.cutoff<-pandemic.parameter[1] # Time (month) of pandemic before study cutoff
  Pan.Duration<-pandemic.parameter[2]      # Duration of pandemic (month)
  prob.missing<-pandemic.parameter[3]      # probability of missing a visit
  missing.type<-pandemic.parameter[4]      # 1=Missing all visits during pandemic, 2=missing with higher prob 
  visit.schedule<-c(c(8,8+12*(1:9),116+16*(1:3))/(52/12),36+6*(1:6)) # schedule in months
  max.schedule.time<-max(visit.schedule)
  
  #print(pandemic.parameter)
  sdata$num.missing<-0
  #sdata$num.any.missing<-0
  sdata$adjusted.event.time<-sdata$Time
  sdata$max.visit.time<-sdata$Time
  sdata$event.censored.by.missing.2visits<-0
  sdata$length.non.missing.visit.time<-0
  
  subject.missing<-rbinom(nrow(sdata),1,prob=prob.missing)
  for (i in 1:nrow(sdata)) {
    visit.times<-visit.schedule[visit.schedule<=sdata$Time[i]]
    if (fix.target.duration==1) time.pan.start<-59-Pan.before.cutoff-sdata$entry.time[i]
    if (fix.target.duration!=1) time.pan.start<-study.cutoff-Pan.before.cutoff-sdata$entry.time[i]
    time.pan.end<-time.pan.start+Pan.Duration
    within.pandemic<-((visit.times<=time.pan.end) & (visit.times>=time.pan.start))
    if (missing.type==1) missing.visit<-rep(subject.missing[i],length(visit.times))
    if (missing.type==2) missing.visit<-rbinom(length(visit.times),1,prob=prob.missing)
    missing.within.pandemic<-(missing.visit & within.pandemic)
    #num.missing<-sum(missing.within.pandemic)
    #sdata$num.any.missing[i]<-sum(missing.within.pandemic)
    num.missing<-num.consecutive.NA(missing.within.pandemic)
    non.missing.visit.times<-visit.times[!missing.within.pandemic] 
    sdata$length.non.missing.visit.time[i]<-length(non.missing.visit.times)
 
    if (length(non.missing.visit.times)>0) sdata$max.visit.time[i]<-max(non.missing.visit.times)
    sdata$num.missing[i]<-num.missing
    
    sdata$event.censored.by.missing.2visits[i]<-((sdata$Event[i]==1) & (num.missing>=2))
    if (num.missing>=2 & sdata$Event[i]==1) {
      max.missing.time.within.pandemic<-min(visit.times[missing.within.pandemic])
      sdata$adjusted.event.time[i]<-min(visit.times[visit.times>max.missing.time.within.pandemic])
      #sdata$max.visit.time.before.missing[i]<-max(visit.times[visit.times<sdata$adjusted.event.time[i]]
    }
  }
  sdata
}

  
# fit cox model to the survival data after introduce possible missing visits  
fit.cox.model<-function(analysis.data,alpha) {
  HR.null<-1
  cox.fit<-coxph(Surv(Time,Event)~trt,ties="efron",data=analysis.data)
  cox.est<-summary(cox.fit)$coef[1:5]
  HRest<-cox.est[2]
  #CIupper<-exp(cox.est[1]+qnorm(1-alpha/2)*cox.est[3])
  #RejH0<-(CIupper<hr.null)      
  RejH0<-(cox.est[5]<(alpha/2))
  maturity<-summary(cox.fit)$nevent/sample.size
  crHR<-exp(log(HR.null)-qnorm(1-alpha/2)*cox.est[3]) 
  c(HRest,cox.est[5],maturity,crHR)
}

Simulation.design.file<-"Simulation design parameters.xlsx"
design.parameter<-read.xlsx(Simulation.design.file,sheetIndex = 1,rowIndex=1:11,colIndex = 1:8)
design.parameter<-design.parameter[design.parameter$Number.of.patients==530,]
design.null.parameter<-design.parameter
design.null.parameter$Target.HR<-1
design.parameter<-rbind(design.parameter,design.null.parameter)
pandemic.assumption.v1<-expand.grid(pan.before.cutoff=c(6,12,24),pan.duration=c(3,6,9),prob.missing=c(0.25,0.5),missing.type=1:2)
pandemic.assumption.v1<-pandemic.assumption.v1[!(pandemic.assumption.v1$pan.before.cutoff<=12 & pandemic.assumption.v1$pan.duration==9),]

#pandemic.assumption<-expand.grid(pan.before.cutoff=c(6,20,24),pan.duration=c(3,6,9),prob.missing=c(0.25,0.5),missing.type=1:2)
pandemic.assumption<-expand.grid(pan.before.cutoff=c(20,24),pan.duration=c(3,6,9),prob.missing=c(0.25,0.5),missing.type=1:2)

n.sim<-5000
sample.size<-530
missing.data.design.parameter<-design.parameter[1:2,]
for (Fup in c(12,14,seq(18,28,2),32,34,36)) {
  added.DOC.for.PFS<-design.parameter[1,]
  added.DOC.for.PFS$Analysis.Stage<-""
  added.DOC.for.PFS$HR.cv<-NA
  added.DOC.for.PFS$Number.events<-NA
  added.DOC.for.PFS$Followup.Time<-Fup+29
  missing.data.design.parameter<-rbind(missing.data.design.parameter,added.DOC.for.PFS)
}
missing.data.design.parameter<-orderBy(~Followup.Time,data=missing.data.design.parameter)
missing.data.design.parameter<-rbind(missing.data.design.parameter,design.null.parameter[1:2,])

sim.power.with.missing.data<-function(sample.size,medS,HR,rr,study.cutoff,alpha, seed){
  # sample.size<-530
  # medS<-26
  # HR<-0.69
  # rr<-1
  # study.cutoff<-45
  # alpha<-0.05
  
  set.seed(seed)
  sdata<-sim.data(sample.size,medS,HR,rr,study.cutoff)
  est.out<-matrix(NA,nrow=1+nrow(pandemic.assumption),ncol=9)
  colnames(est.out)<-c("Pandemic.assumption","n.missing","n.2missing","n.event","n.event.2missing",'HR.M2','pvalue.M2','HR.MR','pvalue.MR')
  fit.FD<-fit.cox.model(sdata,alpha)[1:2]
  est.out[1,]<-c(0,0,0,sum(sdata$Event),0,fit.FD,fit.FD)
  
  for (i in 1:nrow(pandemic.assumption)) {
    adata<-simulate.missing.visit(sdata,study.cutoff,pandemic.assumption[i,],1)
    n.missing<-as.numeric(table(adata$num.missing)[1])
    n.2missing<-sum(adata$num.missing>=2)
    n.event<-sum(adata$Event)
    n.event.censored.by.2missing<-sum(adata$event.censored.by.missing.2visits)
    missing.pattern<-c(n.missing,n.2missing,n.event,n.event.censored.by.2missing)
    
    # if # missing visits>=2, use the missing-2 visit rule & record last visit time as censoring time
    bdata<-adata
    index<-(bdata$num.missing>=2 & bdata$length.non.missing.visit.time>1)
    if (sum(index)>0) {
      bdata[index,]$Time<-bdata[index,]$max.visit.time
      bdata[index,]$Event<-0
    }
    fit.missing.2visit.rule<-fit.cox.model(bdata,alpha)[1:2]
    
    # not using missing-2-visit rule, events are observed, but interval censored b/c missing visits
    bdata<-adata
    index<-(bdata$num.missing>=2 & bdata$event.censored.by.missing.2visits==1)
    if (sum(index)>0) bdata[index,]$Time<-bdata[index,]$adjusted.event.time
    fit.with.missing<-fit.cox.model(bdata,alpha)[1:2]
    est.out[i+1,]<-c(i, missing.pattern,fit.missing.2visit.rule,fit.with.missing)
  }
  est.out
}

#===================================================================================================
# Repeated simulation (N=n.sim) to estimate the p-value, HR
#===================================================================================================

run.missing.data.simulation<-function(n.sim, missing.data.design.parameter,output.csv) {
  for (iS in 1:n.sim) {
    simout<-NULL
    for (k in 1:nrow(missing.data.design.parameter)) {
      sample.size<-missing.data.design.parameter$Number.of.patients[k]
      medS<-missing.data.design.parameter$SOC.median[k]
      HR<-missing.data.design.parameter$Target.HR[k]
      study.cutoff<-missing.data.design.parameter$Followup.Time[k]
      alpha<-0.05
      rr<-1
      seed<-iS
      simout<-rbind(simout,cbind(k,sim.power.with.missing.data(sample.size,medS,HR,rr,study.cutoff,alpha, seed)))
    }
    if (output.csv==TRUE) write.csv(simout,paste("output/MissingData/MissingVisitSimout",iS,".csv",sep=""),row.names = FALSE)
    #if (output.csv==FALSE) print(head(simout))
  }
  simout
}

#====================================================================================================
# Because each repeatation takes about 1 min. The simulation was performed on AZ Scientific Computing 
# Cluster. For the purpose of reproducibility. We demonstrate the results are reproducible using the 
# 1st simulation (iS=1) for current pandemic assumption (Set 2. pandemic.assumption)
# The output is the /output/MissingVisits/MissingVisitSimSetting2/MissingVisitSimout1.csv
#======================================================================================================
simout.iS1<-run.missing.data.simulation(1,missing.data.design.parameter,FALSE)

# aPFS<-gsDesign(k=2,  alpha=0.05/2,  timing=233/311,  sfu="OF", test.type=2)$upper$bound 
# aOS<-gsDesign(k=3,  alpha=0.05/2,  timing=c(175/300,242/300),  sfu="OF", test.type=2)$upper$bound 
# alpha.PFS<-2*(1-pnorm(aPFS))
# alpha.OS<-2*(1-pnorm(aOS))
# alpha.list<-c(alpha.PFS,alpha.OS) # alphal level from R package gsDesign
#alpha.list<-c(0.019221,0.05,0.006677728,0.0251497315,0.05) # alpha level from EAST (Cytel)
alpha.list<-c(0.019221,0.05)
calc.power<-function(simdat) {
  alpha.level<-alpha.list[2-(simdat$k %in% c(3,14))]
  simdat$Rej.M2<-(simdat$pvalue.M2<alpha.level)
  simdat$Rej.MR<-(simdat$pvalue.MR<alpha.level)
  simdat$prob.2missing<-(simdat$n.2missing/530)
  simdat
}
#simdat<-read.csv(paste("output/MissingData/MissingVisitSimout",1,".csv",sep=""))
# Functions to read and output simulation results
extract.simout<-function(output.version,pandemic.assumption) {
  #output.version<-1
  if (output.version==1) output.dir<-"Output/MissingVisits/MissingVisitSimSetting1/MissingVisitSimout";
  if (output.version==2) output.dir<-"Output/MissingVisits/MissingVisitSimSetting2/MissingVisitSimout";
  
  simout<-calc.power(read.csv(paste(output.dir,1,".csv",sep="")))
  for (iS in 2:n.sim) {
    simout<-simout+calc.power(read.csv(paste(output.dir,iS,".csv",sep="")))
  }
  simout<-simout/n.sim
  simout<-merge(simout,data.frame(Pandemic.assumption=1:nrow(pandemic.assumption),pandemic.assumption),by="Pandemic.assumption",all.x=TRUE)
  simout<-merge(simout,data.frame(k=1:nrow(missing.data.design.parameter),missing.data.design.parameter),by="k")
  simout<-orderBy(~k+Pandemic.assumption,data=simout)
  simout$bias.M2<-(simout$HR.M2/simout$Target.HR-1)*100
  simout$bias.MR<-(simout$HR.MR/simout$Target.HR-1)*100
  simout$Power.M2<-simout$Rej.M2*100
  simout$Power.MR<-simout$Rej.MR*100
  
  # For the convenience of output only to select some rows
  simout[simout$Pandemic.assumption==0,]$pan.before.cutoff<-0
  simout[simout$Pandemic.assumption==0,]$pan.duration<-0
  simout[simout$Pandemic.assumption==0,]$prob.missing<-0 
  simout[simout$Pandemic.assumption==0,]$missing.type<-0
  #simout<-simout[!(simout$pan.before.cutoff==20 & simout$pan.duration==9),]

  simout
}

# Read the simulation output and combine them into one csv file
# simout.v1<-extract.simout(1,pandemic.assumption.v1)
# simout.v2<-extract.simout(2,pandemic.assumption)
# 
# simout.m6<-simout.v1[simout.v1$pan.before.cutoff==6,]
# simout.m24<-simout.v2[simout.v2$pan.before.cutoff %in% c(0,24),]
# simout<-rbind(simout.m6,simout.m24)
# write.csv(simout,"Missing data simulation output.csv",row.names=FALSE)
simout<-read.csv("Missing data simulation output.csv")

print.missing.simout<-function(PFS.stage,missing.type) {
  #PFS.stage<-1
  #missing.type<-1
  if (PFS.stage==1) {
    k.list<-c(3,14)
    n.row<-7
  }  
  if (PFS.stage==2) {
    k.list<-c(10,15)
    n.row<-11
  }  
  print(kable(missing.data.design.parameter[k.list[1],],row.names=FALSE,caption="Design characteristics")%>%kable_styling())
  merge.by.var<-names(pandemic.assumption)[1:3]
  simtab1<-simout[simout$k==k.list[1] & simout$missing.type!=3-missing.type,c(merge.by.var,'prob.2missing','bias.M2','Power.M2','bias.MR','Power.MR')]
  simtab2<-simout[simout$k==k.list[2] & simout$missing.type!=3-missing.type,c(merge.by.var,'prob.2missing','bias.M2','Power.M2','bias.MR','Power.MR')]
  simtab<-merge(simtab1,simtab2,by=merge.by.var)
  simtab<-simtab[1:n.row,]
  names(simtab)[1:3]<-c("Time COVID to DCO","Duration","P(Miss)")
  names(simtab)[4:8]<-c("P(2 missing)","% Bias","Power","%Bias","Power")
  names(simtab)[9:13]<-c("P(2 missing)","% Bias","Alpha","%Bias","Alpha")
  print(kable(simtab,digits=2,align=rep("c",ncol(simtab)),row.names=FALSE)%>%add_header_above(c("Use missing-2-visit rule?"=3," "=1,"Yes"=2,"No"=2," "=1,"Yes"=2,"No"=2))%>%add_header_above(c(" "=3,"H1: HR=0.69"=5,"H0: HR=1.0"=5)))
}
```

# Simulated Clinical trial and Questions

- A Phase III randomized controlled trial of sample size n with r:1 randomization ratio. Primary endpoint is time to event (OS or PFS). Assume the median survival time is m years for the control with an hypothesised treatment (statistically measured by a hazard ratio  and associated p value).

- Design characteristics considered in the simulation

```{r,warning=FALSE, message=FALSE,results='asis'}
kable(data.frame(Analysis=c(1,2),missing.data.design.parameter[c(3,10),]),row.names=FALSE)%>%kable_styling()
```

- Question. Effect of missing visits during pandemic on interim and final analyses

    - Section 3.3: Effect of missing-2 visit rule if the missing probability during pandemic is higher. If we apply the missing 2 visits censoring rule on interuppted tumour assessments, how will the power and hypothesised HR change for PFS/EFS/DFS.

    - Section 3.2: The time off-treatment/missing on-study procedure may be associated with OS/PFS. We set the multiplicative factor of HR due to COVID-19 is 1.0 or 1.1 for the control group and 1.4 (1/0.72) for the treatment group. This creates differential COVID-19 impact on the treatment and control groups. 

# Simulation for Section 3.3. Impact of missing visits on HR, Type 1 error and Power

The missingness indicator for a clinical visit follows a Bernoulli distribution during the COVID-19 pandemic (a,b)

- Missing type: 1=missing by patient, 2=missing by visit

- Probability of missing: 0.25 (moderate), 0.5 (high)

- Time that pandemic occurs before DOC: 6 and 24 months before the FA PFS time (59m)

- Pandemic duration: 3, 6 and 9 months

# Supplementary Table 5. Impact of 2-missed visits censoring rule on bias of HR estimates and power/Type 1 error of the interim PFS  (patient-level missing)

- The COVID start time (months after FPR) = total study time (59m) - Time from COVID to DCO (data cutoff) 
- Columns under H1: HR=0.69 are the power reported in Table 5.
- Columns under H0: HR=1.0 are not reported

```{r,warning=FALSE, message=FALSE,results='asis'}
print.missing.simout(1,1)
```

<!-- # Table Not Presented. 5Impact of 2-missed visits censoring rule on bias of HR estimates and power/Type 1 error of the interim PFS  (visit-level missing) -->
<!-- ```{r,warning=FALSE, message=FALSE,results='asis'} -->
<!-- print.missing.simout(1,2) -->
<!-- ``` -->

# Supplementary Tabels 6 and 7. Impact of 2-missed visits censoring rule on bias of HR estimates and power/Type 1 error of the final PFS  (patient-level missing)

-The COVID start time (months after FPR) = total study time (59m) - Time from COVID to DCO (data cutoff) 

- Columns under H1: HR=0.69 are the % bias and power reported in Table 6.

- Columns under H0: HR=1.0 are the % bias and Type 1 error reported in Table 7

```{r,warning=FALSE, message=FALSE,results='asis'}
print.missing.simout(2,1)
```

<!-- # Impact of missing visits on HR and power/Type 1 error for Analysis 2 - Final PFS (visit-level missing) -->
<!-- ```{r,warning=FALSE, message=FALSE,results='asis'} -->
<!-- print.missing.simout(2,2) -->
<!-- ``` -->

# Figure 2A. Plot of \# events, \# of having 2 missing and \# events being impacted by missing 2 visit rule by total study duration (COVID-19 pandemics starts 24 months before the target DCO (59m) and has duration of 9 months, patient-leve missing probably 50\% during pandemic)

```{r,warning=FALSE, message=FALSE,results='asis'}
plot.trend.of.missing<-function(sample.size,medS,HR,rr,pandemic.parameter) {
 
  #study.cutoff<-59
  #sample.size<-10000
  #medS<-26
  #HR<-0.69
  #rr<-1
  #pandemic.parameter<-c(24,9,0.5,1)
  study.cutoff.range<-seq(35,65,2)
  missing.trend<-matrix(NA,nrow=length(study.cutoff.range),ncol=4)
  for (i in 1:length(study.cutoff.range)) {
    #study.cutoff<-45
    set.seed(123)
    sdata<-sim.data(sample.size,medS,HR,rr,study.cutoff.range[i])
    sim.missing.data<-simulate.missing.visit(sdata,study.cutoff.range[i],pandemic.parameter,1)
    sim.missing.data$num.2missing<-(sim.missing.data$num.missing>1)
    missing.trend[i,]<-c(study.cutoff.range[i],as.numeric(colMeans(sim.missing.data)[c(7,13,11)]*530))
  }
  plot(study.cutoff.range,missing.trend[,2],type="o",lty=1,lwd=2,xlab="Study duration (month)",ylab="Number of patients",ylim=c(0,375),xlim=c(35,65),pch=0)
  lines(study.cutoff.range,missing.trend[,3],type="o",lty=2,lwd=2,pch=1)
  lines(study.cutoff.range,missing.trend[,4],type="o",lty=2,lwd=2,pch=2)
  abline(v=c(59-24,59-24+pandemic.parameter[2],59),lty=3)
  text(x=37,y=250,"COVID Starts")
  text(x=45,y=250,"COVID Ends")
  text(x=58,y=250,"Target DCO")
  legend("topleft",lty=1:3,c("# PFS events","# patients having 2 missing","# patients with events impacted by 2 missing visits"),bty="n",lwd=2,pch=0:2)
  #print(missing.trend[nrow(missing.trend)-3,])
  missing.trend
}
Event.Trend<-plot.trend.of.missing(5000,26,0.69,1,c(24,9,0.5,1))
```

# Figure 2B. Plot of patients being impacted: a sensitivity analysis by changing patient-level missing probability from 0.5 to 0.25
```{r,warning=FALSE, message=FALSE,results='asis'}
TmpTrend<-plot.trend.of.missing(5000,26,0.69,1,c(24,9,0.25,1))
Event.Trend<-cbind(Event.Trend,TmpTrend[,4])
```

# Figure 2C. Plot of patients being impacted: a sensitivity analysis by changing Pandemic duration 9m to 3m
```{r,warning=FALSE, message=FALSE,results='asis'}
TmpTrend<-plot.trend.of.missing(5000,26,0.69,1,c(24,3,0.5,1))
Event.Trend<-cbind(Event.Trend,TmpTrend[,4])
```

# Figure 2D. Plot of patients being impacted: a sensitivity analysis by changing patient-level missing to visit-level missing
```{r,warning=FALSE, message=FALSE,results='asis'}
TmpTrend<-plot.trend.of.missing(5000,26,0.69,1,c(24,9,0.5,2))
Event.Trend<-cbind(Event.Trend,TmpTrend[,4])
```

<!-- # Plot of the power for final PFS analysis (alpha=0.05) by study duration after applying the missing 2-visit rules -->

<!-- - The power of logrank test is $\Phi(|\log(\mbox{HR})|\times \sqrt{Dp(1-p)}-Z_{1-\alpha/2})$ where $p=0.5$ for 1:1 randomization, $D$ is number of events and $HR=0.69$ for the PFS analysis. -->

<!-- - Basecase: COVID-19 starts 24 month before final analysis and lasts for 9 month, patient-level missing probably 0.5 during pandemic -->

<!-- ```{r,warning=FALSE, message=FALSE,results='asis'} -->
<!-- event.column.index<-c(2,4:7) -->
<!-- power.trend.matrix<-matrix(NA,nrow=nrow(Event.Trend),ncol=5) -->
<!-- for (i in 1:5) { -->
<!--   n.event<-Event.Trend[,2]-(i>1)*Event.Trend[,event.column.index[i]] -->
<!--   #print(n.event) -->
<!--   for (j in 1:nrow(power.trend.matrix)) { -->
<!--     #alpha<-alpha.list[1+(Event.Trend[j,1]>45)] -->
<!--     alpha<-0.05 -->
<!--     power.trend.matrix[j,i]<-pnorm(abs(log(0.69))*sqrt(0.5*0.5*n.event[j])-qnorm(1-alpha/2)) -->
<!--   }   -->
<!-- } -->
<!-- power.trend.matrix<-cbind(Event.Trend[,1],power.trend.matrix) -->
<!-- par(mar=c(4,4,1,1)) -->
<!-- plot(power.trend.matrix[,1],power.trend.matrix[,2],type="o",lty=1,lwd=2,xlab="Study duration (month)",ylab="Power",ylim=c(0.6,1),xlim=c(35,65)) -->
<!-- color.list<-c("black","red","blue","purple","orange") -->
<!-- for (i in 3:6) lines(power.trend.matrix[,1],power.trend.matrix[,i],lty=i-1,lwd=2,pch=1,col=color.list[i-1]) -->
<!-- abline(v=c(59-24,44,45,59),lty=3) -->
<!-- text(x=37,y=0.9,"COVID starts",cex=1.5) -->
<!-- text(x=44,y=0.9,"COVID ends",cex=1.5) -->
<!-- text(x=46,y=0.95,"Interim (45m)",cex=1.5) -->
<!-- text(x=59,y=0.95,"Target DCO (59m)",cex=1.5) -->
<!-- legend.list<-c("Full data","Base case","Pandemic duration 3 months","patient-level missing prob. 0.25","Visit-level missing prob. 0.5") -->
<!-- legend("bottomright",lty=1:5,col=color.list,legend.list,bty="n",lwd=2,pch=c(1,rep(NA,4))) -->
<!-- ``` -->


# Simulation for Section 3.2. Impact of excess/reduced hazard due to time off-treatment during COVID-19 on HR and power

- Both treatment and control groups have the same excess/reduced hazard during COVID19 pandemic

- Treatment group has more excess hazard than the control during COVID19, i.e., a mulitiplicative factor of HR due to pandemic is 1.0 and 1.1 among the control and a multiplicative factor of HR is 1.1 or 1/0.72 in the treatment group, respectively. For example, when the multiplicative factor is 1.0 in the control and 1/0.72 in the treatment group, this means the patients in the treatment group has the same mortality rate as those in the control group during the COVID-19 period.

- We assume that 10\% or 30% of the patients in the control and treatment groups will be impacted by COVID-19

- Measure of impact (see results in the tables hereafter):

     - Delta: Number of additional number of events in the treatment group vs. control due to COVID-19 if 100\% patients are impacted.
     - Bias of the HR estimate
     - Power/Type 1 error of the log-rank test

```{r,warning=FALSE, message=FALSE, echo=FALSE}
# Simulate survival data according to current hypothetical design
sim.excess.hazard.data<-function(sample.size,medS,HR,rr,study.cutoff,pandemic.EH.parameter) {
  
   # sample.size<-530
   # medS<-26
   # HR<-0.69
   # rr<-1
   # study.cutoff<-45
   # pandemic.EH.parameter<-c(24,6,0.95,1.05)
   # print(pandemic.EH.parameter)
   pandemic.EH.parameter<-as.numeric(pandemic.EH.parameter)
   pan.before.cutoff<-pandemic.EH.parameter[1]
   pan.duration<-pandemic.EH.parameter[2]
   EHR0<-pandemic.EH.parameter[3]
   EHR1<-pandemic.EH.parameter[4]

   accrual.duration<-29 
   accrual.weight<-2    # quadratic accrual
   lambda0 = log(2)/medS
   lambda1 = lambda0*HR
   entry.time <- runif(sample.size)
   entry.time <- (((accrual.duration)^accrual.weight)*entry.time)^(1/accrual.weight)
  
   n0<-round(sample.size/(rr+1))
   n1<-sample.size-n0
   ExpD0<-sum(pexp(study.cutoff-entry.time,rate=lambda0)[1:n0])
   ExpD1<-sum(pexp(study.cutoff-entry.time,rate=lambda1)[n0+(1:n1)])
   
   t0.Exp<-rexp(n0,lambda0)
   t1.Exp<-rexp(n1,lambda1)
   if (pan.before.cutoff>0) {
    IMP<-runif(n0+n1) 
    t0<-rep(NA,n0)
    t1<-rep(NA,n1)
    d0<-rep(NA,n0)
    d1<-rep(NA,n1)
    #t.int<-cbind(pmax(0,study.cutoff-entry.time-pan.before.cutoff),pmax(0,study.cutoff-entry.time-pan.before.cutoff+pan.duration))
    t.int<-cbind(pmax(0,59-entry.time-pan.before.cutoff),pmax(0,59-entry.time-pan.before.cutoff+pan.duration))
    for (i in 1:n0) {
      lambda.pwe<-lambda0*c(1,EHR0,1)
      tchange<-c(0,t.int[i,1],t.int[i,2])
      t0[i]<-rpwe(nr=1,rate=lambda.pwe,tchange=tchange)$r
      d0[i]<-pwe(t=study.cutoff-entry.time[i],rate=lambda.pwe,tchange=tchange)$dist
    }  
    for (i in 1:n1) {
      lambda.pwe<-lambda1*c(1,EHR1,1)
      tchange<-c(0,t.int[i+n0,1],t.int[i+n0,2])
      t1[i]<-rpwe(nr=1,rate=lambda.pwe,tchange=tchange)$r
      d1[i]<-pwe(t=study.cutoff-entry.time[i+n0],rate=lambda.pwe,tchange=tchange)$dist
    }  
    event.time<-(IMP<PERCENT.IMPACTED)*c(t0,t1)+(IMP>=PERCENT.IMPACTED)*c(t0.Exp,t1.Exp)
    ActD0<-(1-PERCENT.IMPACTED)*ExpD0+PERCENT.IMPACTED*sum(d0)
    ActD1<-(1-PERCENT.IMPACTED)*ExpD1+PERCENT.IMPACTED*sum(d1)
    num.death<-c(ExpD0,ExpD1,ActD0,ActD1)
   } else {
    event.time<-c(t0.Exp,t1.Exp)
    num.death<-c(ExpD0,ExpD1,ExpD0,ExpD1)
   }   
   trt<-c(rep(0,n0),rep(1,n1))
   Time<-pmin(study.cutoff-entry.time,event.time)
   Event<-as.numeric(event.time<=(study.cutoff-entry.time))  
   SUBJID<-1:sample.size
   sdata<-data.frame(SUBJID,trt,Time,Event) # 0 control 1 treatment
   c(num.death,fit.cox.model(sdata,alpha))
}

pandemic.excess.hazard.parameter1<-expand.grid(pan.before.cutoff=c(6,-9,24),pan.duration=c(3,6,9),EHR0=c(1.0),EHR1=c(1.1,1/0.72))
pandemic.excess.hazard.parameter2<-expand.grid(pan.before.cutoff=c(6,-9,24),pan.duration=c(3,6,9),EHR0=c(1.1),EHR1=c(1.1,1/0.72))
pandemic.excess.hazard.parameter<-rbind(pandemic.excess.hazard.parameter1,pandemic.excess.hazard.parameter2)
pandemic.excess.hazard.parameter<-pandemic.excess.hazard.parameter[!(pandemic.excess.hazard.parameter$pan.before.cutoff<=12 & pandemic.excess.hazard.parameter$pan.duration==9),]

run.excess.hazard.simulation<-function(n.sim,design.parameter,output.csv) {
  for (iS in 1:n.sim) {
    simout<-NULL
    for (k in 1:nrow(design.parameter)) {
      sample.size<-design.parameter$Number.of.patients[k]
      medS<-design.parameter$SOC.median[k]
      HR<-design.parameter$Target.HR[k]
      study.cutoff<-design.parameter$Followup.Time[k]
      alpha<-0.05
      rr<-1
      for (iP in 0:nrow(pandemic.excess.hazard.parameter)) {
        if (iP==0) sfit<-sim.excess.hazard.data(sample.size,medS,HR,rr,study.cutoff,c(0,0,0,0))
        if (iP>0)  sfit<-sim.excess.hazard.data(sample.size,medS,HR,rr,study.cutoff,pandemic.excess.hazard.parameter[iP,])
        simout<-rbind(simout,c(k,iP,sfit))
      }
    }
    colnames(simout)<-c("iDesign","iPandemic","ExpD0","ExpD1","ActD0","ActD1",'HREst','RejH0', 'Maturity','CriticalHR')
    if (output.csv==TRUE) {
      if (PERCENT.IMPACTED==0.1) write.csv(simout,paste("output/TimeOffTreatment/Impact10Percent/ExcessHazardSimout",iS,".csv",sep=""),row.names = FALSE)
      if (PERCENT.IMPACTED==0.3) write.csv(simout,paste("output/TimeOffTreatment/Impact30Percent/ExcessHazardSimout",iS,".csv",sep=""),row.names = FALSE)
      
    }
    if (output.csv==FALSE) {
      #print(head(simout))
    }
  }
  simout
}

PERCENT.IMPACTED<-0.1
#EH.simout<-run.excess.hazard.simulation(1,pandemic.excess.hazard.parameter,FALSE)

# Functions to read and output simulation results
#scp.directory<-"//samba.scp.astrazeneca.net/klrd547/Oncology/Projects/COVID19 Simulation/output"
alpha.list<-c(0.019221,0.05,0.006677728,0.0251497315,0.05) # alpha level from EAST (Cytel)
mod<-function(n, m) { n - m * floor(n/m) }
comp.power<-function(simdat) {
  alpha.level<-alpha.list[1+mod(simdat$iDesign-1,5)]
  simdat$RejH0<-(simdat$RejH0<alpha.level)
  simdat
}

# dir.output.10pct.impacted<-"output/TimeOffTreatment/Impact10Percent/ExcessHazardSimout";
# dir.output.30pct.impacted<-"output/TimeOffTreatment/Impact30Percent/ExcessHazardSimout";
# 
# simout30<-comp.power(read.csv(paste(dir.output.30pct.impacted,1,".csv",sep="")))
# for (iS in 2:n.sim) {
# simout30<-simout30+comp.power(read.csv(paste(dir.output.30pct.impacted,iS,".csv",sep="")))
# }
# simout30<-simout30/n.sim
# 
# simout10<-comp.power(read.csv(paste(dir.output.10pct.impacted,1,".csv",sep="")))
# for (iS in 2:n.sim) {
# simout10<-simout10+comp.power(read.csv(paste(dir.output.10pct.impacted,iS,".csv",sep="")))
# }
# simout10<-simout10/n.sim
# 
# # Only use one set of 5000 simulation under null hypothesis (iPandemic==0)
# simout30<-rbind(simout10[simout10$iPandemic==0,],simout30[simout30$iPandemic!=0,])
# 
# write.csv(simout10,"Excess hazard simulation output 10percent.csv",row.names=FALSE)
# write.csv(simout30,"Excess hazard simulation output 30percent.csv",row.names=FALSE)
simout10<-read.csv("Excess hazard simulation output 10percent.csv")
simout30<-read.csv("Excess hazard simulation output 30percent.csv")

merge.pandemic.excess.hazard.parameter<-function(simout,prop.impacted) {
  simout$k<-simout$iDesign
  simout$Sim<-NULL
  simout$iDesign<-NULL
  simout<-merge(simout,data.frame(iPandemic=1:nrow(pandemic.excess.hazard.parameter),pandemic.excess.hazard.parameter),by="iPandemic",all.x=TRUE)
  
  simout<-merge(simout,data.frame(k=1:nrow(design.parameter),design.parameter),by="k")
  simout<-orderBy(~k+iPandemic,data=simout)
  simout$bias<-(simout$HREst/simout$Target.HR-1)*100
  simout$Power<-simout$RejH0*100
  
  # For the convenience of output only to select some rows
  simout[simout$iPandemic==0,]$pan.before.cutoff<-0
  simout[simout$iPandemic==0,]$pan.duration<-0
  simout[simout$iPandemic==0,]$EHR0<-1
  simout[simout$iPandemic==0,]$EHR1<-1
  simout<-simout[simout$pan.before.cutoff!=12,]
  #simout$Diff.Dth<-((simout$ActD1-simout$ExpD1)-(simout$ActD0-simout$ExpD0))*10
  simout$Diff.Dth<-((simout$ActD1-simout$ExpD1)-(simout$ActD0-simout$ExpD0))/prop.impacted
  #simout[simout$pan.before.cutoff==6,]$pan.before.cutoff<-20
  simout
}

simout10<-merge.pandemic.excess.hazard.parameter(simout10,0.1)
simout30<-merge.pandemic.excess.hazard.parameter(simout30,0.3)

print.excess.hazard.simout<-function(simout10,simout30,k) {
  
  print(kable(design.parameter[k,],row.names=FALSE,caption="Design characteristics")%>%kable_styling())
  merge.by.var<-names(pandemic.excess.hazard.parameter)[1:2]
  
  select.simout<-function(simout,k) {
    extract.simtab<-function(EHR) {
      simtab0<-simout[simout$k==k & simout$EHR0==1 & simout$EHR1==1,c(merge.by.var,"Diff.Dth","bias","Power")]
      simtab1<-simout[simout$k==k & simout$EHR0==EHR[1] & simout$EHR1==EHR[2],c(merge.by.var,"Diff.Dth","bias","Power")]
      rbind(simtab0,simtab1)
    }
    EHR.list<-data.matrix(expand.grid(EHR0=c(1.0,1.1),EHR1=c(1.1,1/0.72)))
    simtab<-extract.simtab(EHR.list[3,]) # only use select HR=(1.0,1.4) and (1.1, 1.4)
    for (i in 4:4) simtab<-merge(simtab,extract.simtab(EHR.list[i,]),by=merge.by.var,all=TRUE)
    # Final PFS analysis
    if (k %in% c(2,7)) {
      simtab<-simtab[simtab$pan.before!=-9,]
    } 
    # Final OS analysis
    if (k %in% c(5,10)) {
      simtab[simtab$pan.before.cutoff!=0,]$pan.before.cutoff<-simtab[simtab$pan.before.cutoff!=0,]$pan.before.cutoff+15
    }
    simtab<-orderBy(~pan.before.cutoff+pan.duration,data=simtab)
  }

  simtab10<-select.simout(simout10,k)
  simtab30<-select.simout(simout30,k)
  simtab.comb<-merge(simtab10,simtab30,by=merge.by.var)
  simtab.comb<-orderBy(~pan.before.cutoff,data=simtab.comb)
  
  names(simtab.comb)[1:2]<-c("Time COVID-DCO","COVID duration")
  if (k<=5) names(simtab.comb)[3:14]<-rep(c("Delta","% Bias","Power"),4)
  if (k>5) names(simtab.comb)[3:14]<-rep(c("Delta","% Bias","Alpha"),4)
  print(kable(simtab.comb,digits=1,align=rep("c",ncol(simtab.comb)),row.names=FALSE)%>%add_header_above(c(" "=2,"(1.0,1.4)"=3,"(1.1,1.4)"=3,"(1.0,1.4)"=3,"(1.1,1.4)"=3))%>%add_header_above(c(" "=2,"10% Impacted"=6,"30% Impacted"=6))%>%add_header_above(c(" "=2,"Excess HR of patients during pandemic among control and treatment groups (HR0,HR1) "=12)))
}
```


# Supplementary Table 1. Impact of time off-treatment on PFS HR estimates and Type 1 error under the null hypothesis for PFS 

- In Table 1, The COVID start time (months after FPR) = total study time (59m) - Time from COVID to DCO (data cutoff) 

```{r,warning=FALSE, message=FALSE,results='asis'}
print.excess.hazard.simout(simout10,simout30,7)
```

# Supplementary Table 2. Impact of time off-treatment on PFS HR estimates and power under the alternative hypothesis for PFS .

- In Table 2, The COVID start time (months after FPR) = total study time (59m) - Time from COVID to DCO (data cutoff) 

```{r,warning=FALSE, message=FALSE,results='asis'}
print.excess.hazard.simout(simout10,simout30,2)
```

# Supplementary Table 3. Impact of time off-treatment on OS HR estimates and Type 1 error under the null hypothesis for OS 

- In Table 3, The COVID start time (months after FPR) = total study time (74m) - Time from COVID to DCO (data cutoff) 

```{r,warning=FALSE, message=FALSE,results='asis'}
print.excess.hazard.simout(simout10,simout30,10)
```

# Supplementary Table 4. Impact of time off-treatment on OS HR estimates and power under the alternative hypothesis for OS .

- In Table 4, The COVID start time (months after FPR) = total study time (74m) - Time from COVID to DCO (data cutoff) 

```{r,warning=FALSE, message=FALSE,results='asis'}
print.excess.hazard.simout(simout10,simout30,5)
```
