#replace this directory with the one of the unzipped files
setwd("")



#install the "WWR" package from CRAN
install.packages("WWR")

library(WWR)


#install package "WR" from GitHub
install.packages("devtools")
library(devtools)
install_github("lmaowsic/Win-ratio/WR")
# or install it manually from "./R fun/WR_0.1.0.tar.gz"

library(WR)





# the non-ischemic patients data
data=non_ischemic

colnames(data)[4:16]=c(
  "Training vs Usual","Age (year)","Male vs Female","Black vs White", 
  "Other vs White", "BMI","LVEF","Hypertension","COPD","Diabetes",
  "ACE Inhibitor","Beta Blocker", "Smoker"
)



data[1:10,]

id=data[,1]
time=data[,2]
status=data[,3]
Z=as.matrix(data[,4:ncol(data)])
trt=Z[,1]

#total number of patients
n=length(unique(id))

#number of patients in each arm
n1=length(unique(id[trt==1&status<2]))
n0=length(unique(id[trt==0&status<2]))





## create dataset for TFE##############
## by taking the first record from each 

#assuming data have been sorted by ID and time
data.TFE = data[!duplicated(data[,"ID"]),]
data.TFE[1:100,]
#######################




######################################################
#   Two-sample analysis                              #
######################################################

### Cox Ph on TFE ####
summary(coxph(Surv(data.TFE[,"time"], data.TFE[,"status"]>0) ~ data.TFE[,"Training vs Usual"]))
#HR=0.804 (95% CI: 0.643, 1.006)

### WR ####

### First use the WWR package ###
library(WWR)

#unique id list
uid=unique(id)

#create the input variables for the function wwratio
y1=rep(NA,n)
y2=rep(NA,n)
d1=rep(NA,n)
d2=rep(NA,n)
z=rep(NA,n)


for (i in 1:n){
  tmp=data[id==uid[i],]
  
  z[i]=tmp[1,"Training vs Usual"]
  
  
  #survival time
  y2.ind=(tmp[,"status"]<2)
  y2[i]=tmp[y2.ind,"time"]
  d2[i]=tmp[y2.ind,"status"]
  
  #non-fatal event time
  if (2 %in% tmp[,"status"]){
    y1[i]=min(tmp[tmp[,"status"]==2,"time"])
    d1[i]=1
  }else{
    y1[i]=y2[i]
    d1[i]=0
  }
}

#unweighted win ratio
wtest.unweighted=winratio(y1=y1,y2=y2,d1=d1,d2=d2,z=z)
summary(wtest.unweighted)
# WR=1.2349 (95% CI: 0.9742, 1.5652)

### Then use the WR package with Training vs Usual as the only covariate###
obj=pwreg(time=time,status=status,Z=Z[,1],ID=id)
obj
# WR=1.2333 (95% CI:0.9719, 1.564854)





######################################################
#   Regression analysis                              #
######################################################
#### Cox PH on TFE#######
summary(coxph(Surv(time, status>0) 
  ~ `Training vs Usual`+`Age (year)`+`Male vs Female`+
 `Black vs White`+`Other vs White`+BMI+ LVEF
   +Hypertension+COPD+Diabetes+`ACE Inhibitor`+`Beta Blocker`+Smoker,
    data=data.TFE))

#### Proportional mean model for recurrent CE#######
CompoML(id=id,time=time,status=status,Z=Z,c(1,1))#unweighted


### (Priority-adjusted) proportional win-fractions regression model
obj=pwreg(time=time,status=status,Z=Z,ID=id)
obj

#Conduct a 2-df test on the effect of race
beta =matrix(obj$beta[4:5])
#extract estimated covariance matrix for (\beta_4,\beta_5)
Sigma =obj$Var[4:5,4:5]
#compute chisq statistic in quadratic form
chistats =t(beta) %*% solve(Sigma) %*% beta  

#compare the Wald statistic with the reference
# distribution of chisq(2) to obtain the p-value
1 - pchisq(chistats, df=2)



#plot the standardized score processes
#to check the proportionality assumption
score.obj =score.proc(obj)

par(mfrow = c(4,4))
for(i in c(1:13)){
  plot(score.obj, k = i, lwd=2)
}