## R function #############

meandurationAUC=function(PFS, STATUS_PFS, ETIME, STATUS_TTR, tau){
  ETIME[STATUS_TTR==0]=Inf
  
  y1=pmin(PFS, ETIME, tau)
  delta1=1*(STATUS_TTR+STATUS_PFS>0)
  delta1[y1==tau]=1
  
  y2=pmin(PFS, tau)
  delta2=STATUS_PFS
  delta2[y2==tau]=1
  
  fit1=survfit(Surv(y1, delta1)~1)
  n1=length(fit1$time)
  auc1=sum(c(1, fit1$surv[-n1])*(fit1$time-c(0, fit1$time[-n1])))
  
  
  fit2=survfit(Surv(y2, delta2)~1)
  n2=length(fit2$time)
  auc2=sum(c(1, fit2$surv[-n2])*(fit2$time-c(0, fit2$time[-n2])))
  
  muhat=auc2-auc1
  
  N=length(y1)
  prisk1=pauc1=rep(0, N)
  for(i in 1:N)
  {k1=sum(fit1$time<=y1[i])
  if(k1>1)
    pauc1[i]=sum(c(1, fit1$surv[1:(k1-1)])*(fit1$time[1:k1]-c(0, fit1$time[1:(k1-1)])))
  if(k1==1)
    pauc1[i]=fit1$time[1]
  prisk1[i]=mean(y1>=y1[i])
  }
  pauc1=auc1-pauc1
  
  xi1=delta1*pauc1/prisk1
  
  tau1=rep(0, N)
  for(i in 1:N)
  {tau1[i]=xi1[i]-mean(xi1*(y1[i]>=y1)/prisk1)}
  
  prisk2=pauc2=rep(0, N)
  for(i in 1:N)
  {k2=sum(fit2$time<=y2[i])
  if(k2>1)
    pauc2[i]=sum(c(1, fit2$surv[1:(k2-1)])*(fit2$time[1:k2]-c(0, fit2$time[1:(k2-1)])))
  if(k2==1)
    pauc2[i]=fit2$time[1]
  prisk2[i]=mean(y2>=y2[i])
  }
  pauc2=auc2-pauc2
  
  xi2=delta2*pauc2/prisk2
  
  tau2=rep(0, N)
  for(i in 1:N)
  {tau2[i]=xi2[i]-mean(xi2*(y2[i]>=y2)/prisk2)}
  
  sigma=sqrt(var(tau1-tau2)/N)
  
  return(list(meandur.est=muhat, meandur.se=sigma))
}

##### INPUT

# ETIME: time to response (month)
# STATUS_TTR: status for response (=1 for responders; =0 for non-responders)
# PFS: progression-free survival time (month)
# STATUS_PFS: 1=event; 0=censored
# tau:  truncation time for defining the time-window of interest

##### Output

# meandur.est:  the mean response-duration during the time window [0, tau]
# meandur.es :  the standard error of meandur.est

#### Example #####

library(survival)
data=read.table(data.csv)
meandurationAUC(data$PFS, data$STATUS_PFS, data$ETIME, data$STATUS_TTR, 24)

##################