dat <- read.csv("EmCAB_Biomarker_Data.csv")

library(survival)

(196+187)/2886


####################
dat$CVorMI = NA
for (i in 1:length(dat$CVDeath5yr)){
  dat$CVorMI[i] <- 0
  if (dat$CVDeath5yr[i] == 1){
    dat$CVorMI[i] <- 1
  }
  else if (dat$MI5yr[i] == 1){
    dat$CVorMI[i] <- 1
  }
}
dat$DeathorMI <- 0
for (i in 1:length(dat$DeathorMI)){
  dat$DeathorMI[i] <- 0
  if (dat$death5yr[i] == 1){
    dat$DeathorMI[i] <- 1
  }
  else if (dat$MI5yr[i] == 1){
    dat$DeathorMI[i] <- 1
  }
}


dat$SurvObj3 <- with(dat, Surv(dat$timetoMI5yr, DeathorMI==1))

###################################
# Correcting zeros
###################################
min(dat$hsTroponin_pgmL[which(dat$hsTroponin_pgmL > 0)])
for (i in 1:nrow(dat)){
  if (dat$hsTroponin_pgmL[i] == 0){
    dat$hsTroponin_pgmL[i] <- 0.3*(1/sqrt(2))
  }
}

min(dat$HSP70_ngmL[which(dat$HSP70_ngmL > 0)])
for (i in 1:nrow(dat)){
  if (dat$HSP70_ngmL[i] == 0){
    dat$HSP70_ngmL[i] <- 1*(1/sqrt(2))
  }
}

min(dat$FDP_ugmL[which(dat$FDP_ugmL > 0)])
for (i in 1:nrow(dat)){
  if (dat$FDP_ugmL[i] == 0){
    dat$FDP_ugmL[i] <- 0.005*(1/sqrt(2))
  }
}

min(dat$suPAR_ngmL[which(dat$suPAR_ngmL > 0)])
for (i in 1:nrow(dat)){
  if (dat$suPAR_ngmL[i] == 0){
    dat$suPAR_ngmL[i] <- 0.78*(1/sqrt(2))
  }
}

min(dat$CRP_mgL[which(dat$CRP_mgL > 0)])
for (i in 1:nrow(dat)){
  if (dat$CRP_mgL[i] == 0){
    dat$CRP_mgL[i] <- 0.05*(1/sqrt(2))
  }
}

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


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

dat$hsTroponin_pgmL.log <- log(dat$hsTroponin_pgmL)
dat$HSP70_ngmL.log <- log(dat$HSP70_ngmL)
dat$suPAR_ngmL.log <-log(dat$suPAR_ngmL)
dat$FDP_ugmL.log <- log(dat$FDP_ugmL)
dat$CRP_mgL.log <-log(dat$CRP_mgL)


################
smp_size <- floor(0.5*(nrow(dat)))
set.seed(1836)
train_ind  <- sample(seq_len(nrow(dat)),size=smp_size)


training <- dat[train_ind,]
testing <- dat[-train_ind,]


dat<- training


length(dat$HSP70_ngmL[which(dat$HSP70_ngmL<=7.072e-1)])/length(dat$HSP70_ngmL)
########## Demographics 
# Age

length(dat$uniqueid)
mean(dat$age) 
sd(dat$age)

# Male sex
length(dat$male[which(dat$male == 1)])
(length(dat$male[which(dat$male == 1)])/ length(dat$male)) * 100

# White 
length(dat$white[which(dat$white == 1)])
(length(dat$white[which(dat$white == 1)])/ length(dat$white)) * 100

# Black
length(dat$black[which(dat$black == 1)])
(length(dat$black[which(dat$black == 1)])/ length(dat$black)) * 100

# Blood Pressure (NOT REPORTED)

# BMI
mean(dat$BMI[which(is.na(dat$BMI) == FALSE)])
sd(dat$BMI[which(is.na(dat$BMI) == FALSE)])

######## Follow-Up Events

# Myocardial Infarction
length(dat$MI[which(dat$MI == 1)])
(length(dat$MI[which(dat$MI == 1)])/ length(dat$MI)) * 100

# All-cause Death
length(dat$death[which(dat$death == 1)])
(length(dat$death[which(dat$death == 1)])/ length(dat$death)) * 100

# CV Death
length(dat$CVdeath[which(dat$CVdeath == 1)])
(length(dat$CVdeath[which(dat$CVdeath == 1)])/ length(dat$CVdeath)) * 100




c("timetoMI5yr", "DeathorMI", "BRS", "age.Cat" , "BMI.Cat", "eGFR.Cat", "Gensini.Cat", "male", 
                          "black", "PastSmoking", "CurrentSmoking", "HTN", "DM", "HC", "MIhx",
                          "HFhx", "EverStat", "EverARBACE", "BRS.CENTERED")

length(dat$CADhx[which(dat$HC == 1)])
(length(dat$CADhx[which(dat$HC == 1)])/ length(dat$HC)) * 100

length(testing$CADhx[which(testing$HC == 1)])
(length(testing$CADhx[which(testing$HC == 1)])/ length(testing$HC)) * 100


table(testing$eGFR.Cat)/(1084+358)
table(dat$Gensini.Cat)/(853+501)

summary(testing$age)

sum(dat$DeathorMI)
mean(dat$DeathorMI)

sum(testing$DeathorMI)
mean(testing$DeathorMI)
############# Biomarkers

# HS Trop.
quantile(dat$hsTroponin_pgmL)
# HS 70
quantile(dat$HSP70_ngmL)
# FDP 
quantile(dat$FDP_ugmL)
# suPAR
quantile(dat$suPAR_ngmL)
# CRP 
quantile(dat$CRP_mgL)

############## Histories

# Coronary Artery Disease
length(dat$CADhx[which(dat$CADhx == 1)])
(length(dat$CADhx[which(dat$CADhx == 1)])/ length(dat$CADhx)) * 100

# Hist of MI
length(dat$MIhx[which(dat$MIhx == 1)])
(length(dat$MIhx[which(dat$MIhx == 1)])/ length(dat$MIhx)) * 100

# Hist of PCI
length(dat$PCIhx[which(dat$PCIhx == 1)])
(length(dat$PCIhx[which(dat$PCIhx == 1)])/ length(dat$PCIhx)) * 100

# History of CABG
length(dat$CABGhx[which(dat$CABGhx == 1)])
(length(dat$CABGhx[which(dat$CABGhx == 1)])/ length(dat$CABGhx)) * 100

# diabetes 
length(dat$DM[which(dat$DM == 1)])
(length(dat$DM[which(dat$DM == 1)])/ length(dat$DM)) * 100

# Hupertension
length(dat$HTN[which(dat$HTN == 1)])
(length(dat$HTN[which(dat$HTN == 1)])/ length(dat$HTN)) * 100

# Smoking
length(dat$CurrentSmoking[which(dat$CurrentSmoking == 1)])
(length(dat$CurrentSmoking[which(dat$CurrentSmoking == 1)])/ length(dat$CurrentSmoking)) * 100

length(dat$PastSmoking[which(dat$PastSmoking == 1)])
(length(dat$PastSmoking[which(dat$PastSmoking == 1)])/ length(dat$PastSmoking)) * 100


####################### Medications

# ARB ACE-inh
length(dat$EverARBACE[which(dat$EverARBACE == 1)])
(length(dat$EverARBACE[which(dat$EverARBACE == 1)])/ length(dat$EverARBACE)) * 100

# Aspirin
length(dat$EverAsp[which(dat$EverAsp == 1)])
(length(dat$EverAsp[which(dat$EverAsp == 1)])/ length(dat$EverAsp)) * 100

# Plavix
length(dat$EverPlavix[which(dat$EverPlavix == 1)])
(length(dat$EverPlavix[which(dat$EverPlavix == 1)])/ length(dat$EverPlavix)) * 100

# Statins
length(dat$EverStat[which(dat$EverStat == 1)])
(length(dat$EverStat[which(dat$EverStat == 1)])/ length(dat$EverStat)) * 100

# Beta Blockers
length(dat$EverBB[which(dat$EverBB == 1)])
(length(dat$EverBB[which(dat$EverBB == 1)])/ length(dat$EverBB)) * 100

##################### Angiographic Findings

quantile(dat$Gensini_OM1[which(!is.na(dat$Gensini_OM1))])
table(dat$CAD50Vessels[which(!is.na(dat$CAD50Vessels))])
220/length((dat$CAD50Vessels[which(!is.na(dat$CAD50Vessels))]))

#*****************************************************************************
length(testing$uniqueid)
mean(testing$age) 
sd(testing$age)

# Male sex
length(testing$male[which(testing$male == 1)])
(length(testing$male[which(testing$male == 1)])/ length(testing$male)) * 100

# White 
length(testing$white[which(testing$white == 1)])
(length(testing$white[which(testing$white == 1)])/ length(testing$white)) * 100

# Black
length(testing$black[which(testing$black == 1)])
(length(testing$black[which(testing$black == 1)])/ length(testing$black)) * 100

# Blood Pressure (NOT REPORTED)

# BMI
mean(testing$BMI[which(is.na(testing$BMI) == FALSE)])
sd(testing$BMI[which(is.na(testing$BMI) == FALSE)])

######## Follow-Up Events

# Myocardial Infarction
length(testing$MI[which(testing$MI == 1)])
(length(testing$MI[which(testing$MI == 1)])/ length(testing$MI)) * 100

# All-cause Death
length(testing$death[which(testing$death == 1)])
(length(testing$death[which(testing$death == 1)])/ length(testing$death)) * 100

# CV Death
length(testing$CVdeath[which(testing$CVdeath == 1)])
(length(testing$CVdeath[which(testing$CVdeath == 1)])/ length(testing$CVdeath)) * 100

############# Biomarkers

# HS Trop.
quantile(testing$hsTroponin_pgmL)
# HS 70
quantile(testing$HSP70_ngmL)
# FDP 
quantile(testing$FDP_ugmL)
# suPAR
quantile(testing$suPAR_ngmL)
# CRP 
quantile(testing$CRP_mgL)

############## Histories

# Coronary Artery Disease
length(testing$CADhx[which(testing$CADhx == 1)])
(length(testing$CADhx[which(testing$CADhx == 1)])/ length(testing$CADhx)) * 100

# Hist of MI
length(testing$MIhx[which(testing$MIhx == 1)])
(length(testing$MIhx[which(testing$MIhx == 1)])/ length(testing$MIhx)) * 100

# Hist of PCI
length(testing$PCIhx[which(testing$PCIhx == 1)])
(length(testing$PCIhx[which(testing$PCIhx == 1)])/ length(testing$PCIhx)) * 100

# History of CABG
length(testing$CABGhx[which(testing$CABGhx == 1)])
(length(testing$CABGhx[which(testing$CABGhx == 1)])/ length(testing$CABGhx)) * 100

# diabetes 
length(testing$DM[which(testing$DM == 1)])
(length(testing$DM[which(testing$DM == 1)])/ length(testing$DM)) * 100

# Hupertension
length(testing$HTN[which(testing$HTN == 1)])
(length(testing$HTN[which(testing$HTN == 1)])/ length(testing$HTN)) * 100

# Smoking
length(testing$CurrentSmoking[which(testing$CurrentSmoking == 1)])
(length(testing$CurrentSmoking[which(testing$CurrentSmoking == 1)])/ length(testing$CurrentSmoking)) * 100

length(testing$PastSmoking[which(testing$PastSmoking == 1)])
(length(testing$PastSmoking[which(testing$PastSmoking == 1)])/ length(testing$PastSmoking)) * 100


####################### Medications

# ARB ACE-inh
length(testing$EverARBACE[which(testing$EverARBACE == 1)])
(length(testing$EverARBACE[which(testing$EverARBACE == 1)])/ length(testing$EverARBACE)) * 100

# Aspirin
length(testing$EverAsp[which(testing$EverAsp == 1)])
(length(testing$EverAsp[which(testing$EverAsp == 1)])/ length(testing$EverAsp)) * 100

# Plavix
length(testing$EverPlavix[which(testing$EverPlavix == 1)])
(length(testing$EverPlavix[which(testing$EverPlavix == 1)])/ length(testing$EverPlavix)) * 100

# Statins
length(testing$EverStat[which(testing$EverStat == 1)])
(length(testing$EverStat[which(testing$EverStat == 1)])/ length(testing$EverStat)) * 100

# Beta Blockers
length(testing$EverBB[which(testing$EverBB == 1)])
(length(testing$EverBB[which(testing$EverBB == 1)])/ length(testing$EverBB)) * 100

##################### Angiographic Findings

quantile(testing$Gensini_OM1[which(!is.na(testing$Gensini_OM1))])
table(testing$CAD50Vessels[which(!is.na(testing$CAD50Vessels))])
214/length((testing$CAD50Vessels[which(!is.na(testing$CAD50Vessels))]))



#*****************************************************************************
pairs(~ rank(hsTroponin_pgmL) + rank(HSP70_ngmL) + rank(FDP_ugmL) + rank(suPAR_ngmL) + rank(CRP_mgL), data=testing, pch=16)
text(1,1,"0.1056")


testing$hsTroponin_pgmL



cors <- cor(dat[,2:6],method = "spearman")


# Survival Analysis 



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


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


library(survival)
library(KMsurv)
library(xtable)
library(splines)
library("pROC")
library("aod")

head(dat)

#########################
dat.f <- dat[which(testing$male==0),]
dat.m <- dat[which(testing$male==1),]


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

dat.f <- dat[which(dat$male==0),]
dat.m <- dat[which(dat$male==1),]
##########################




dat$A3 <- NA
for (i in 1:length(dat$uniqueid)){
  if (dat$hsTroponin_pgmL[i] <= 1.6){
    dat$A3[i] = 1
  }
  else if (dat$hsTroponin_pgmL[i] <= 2.1){
    dat$A3[i] = 2
  }
  else if (dat$hsTroponin_pgmL[i] <= 2.6){
    dat$A3[i] = 3
  }
  else if(dat$hsTroponin_pgmL[i] <= 6.9){
    dat$A3[i] = 4
  }
  else if(dat$hsTroponin_pgmL[i] > 6.9){
    dat$A3[i] = 5
  }
}



length(dat$uniqueid[which(dat$B3==2)])

dat$B3 <- NA
for (i in 1:length(dat$uniqueid)){
  if (dat$HSP70_ngmL.log[i] <= 0.04116454){
    dat$B3[i] = 1
  }
  else if (dat$HSP70_ngmL.log[i] <= 5.017524){
    dat$B3[i] = 2
  }
  else if(dat$HSP70_ngmL.log[i] > 5.017524){
    dat$B3[i] = 3
  }
}

table(dat$B3)


dat$D3 <- NA
for (i in 1:length(dat$uniqueid)){
  if (dat$FDP_ugmL[i] <= 0.27){
    dat$D3[i] = 1
  }
  else if (dat$FDP_ugmL[i] <= 0.48){
    dat$D3[i] = 2
  }
  else if (dat$FDP_ugmL[i] <= 0.88){
    dat$D3[i] = 3
  }
  else if(dat$FDP_ugmL[i] > 0.88){
    dat$D3[i] = 4
  }
}




dat$E3 <- NA
for (i in 1:length(dat$uniqueid)){
  if (dat$CRP_mgL[i] <= 1.6){
    dat$E3[i] = 1
  }
  else if (dat$CRP_mgL[i] <= 5.8){
    dat$E3[i] = 2
  }
  else if (dat$CRP_mgL[i] <= 16.9){
    dat$E3[i] = 3
  }
  else if(dat$CRP_mgL[i] > 16.9){
    dat$E3[i] = 4
  }
}




dat.f$C3 <- NA
for (i in 1:length(dat.f$uniqueid)){
  if (dat.f$suPAR_ngmL[i] <= 3.6){
    dat.f$C3[i] = 1
  }
  else if(dat.f$suPAR_ngmL[i] > 3.6){
    dat.f$C3[i] = 2
  }
}

table(dat.f$C3)




dat.m$C3 <- NA
for (i in 1:length(dat.m$uniqueid)){
  if (dat.m$suPAR_ngmL[i]  <= 3.2){
    dat.m$C3[i] = 1
  }
  else if(dat.m$suPAR_ngmL[i] > 3.2){
    dat.m$C3[i] = 2
  }
}

dat$C3 = NA
dat$C3[which(testing$male==0 & dat$suPAR_ngmL <= 3.6)] <- 1
dat$C3[which(testing$male==0 & dat$suPAR_ngmL > 3.6)] <- 2

dat$C3[which(testing$male !=0 & dat$suPAR_ngmL <= 3.2)] <- 1
dat$C3[which(testing$male !=0 & dat$suPAR_ngmL > 3.2)] <- 2









table(dat$C3)

dat.m$SurvObj3 <- with(dat.m, Surv(dat.m$timetoMI5yr, DeathorMI==1))

KM_C3m <- survfit(SurvObj3~C3, data=dat.m, conf.type="log")
plot(KM_C3m, col=c(1:4), main=("CV Death or MI by suPAR (men)"), ylim=c(0.5, 1))
legend("bottomleft", c("Low","Medium","High","Very High"), col=c(1:4), lty=1, bty="n")

dat.f$SurvObj3 <- with(dat.f, Surv(dat.f$timetoMI5yr, DeathorMI==1))

KM_C3f <- survfit(SurvObj3~C3, data=dat.f, conf.type="log")
plot(KM_C3f, col=c(1:4), main=("CV Death or MI by suPAR (women)"), ylim=c(0.5, 1))
legend("bottomleft", c("Low","Medium","High","Very High"), col=c(1:4), lty=1, bty="n")


table(dat$A3)
table(dat$B3)
table(dat$D3)
table(dat$E3) 
table(dat.f$C3)
table(dat.m$C3)


###################################
###################################
testing$SurvObj3 <- with(testing, Surv(testing$timetoMI5yr, DeathorMI==1))

testing.f <- dat[which(testing$male==0),]
testing.m <- dat[which(testing$male==1),]

testing.f$SurvObj3 <- with(testing.f, Surv(testing.f$timetoMI5yr, DeathorMI==1))
testing.m$SurvObj3 <- with(testing.m, Surv(testing.m$timetoMI5yr, DeathorMI==1))



testing$A3 <- NA
for (i in 1:length(testing$uniqueid)){
  if (testing$hsTroponin_pgmL[i] <= 1.6){
    testing$A3[i] = 1
  }
  else if (testing$hsTroponin_pgmL[i] <= 2.1){
    testing$A3[i] = 2
  }
  else if (testing$hsTroponin_pgmL[i] <= 2.6){
    testing$A3[i] = 3
  }
  else if(testing$hsTroponin_pgmL[i] <= 6.9){
    testing$A3[i] = 4
  }
  else if(testing$hsTroponin_pgmL[i] > 6.9){
    testing$A3[i] = 5
  }
}



length(testing$uniqueid[which(testing$B3==2)])

testing$B3 <- NA
for (i in 1:length(testing$uniqueid)){
  if (testing$HSP70_ngmL.log[i] <= 0.04116454){
    testing$B3[i] = 1
  }
  else if (testing$HSP70_ngmL.log[i] <= 5.017524){
    testing$B3[i] = 2
  }
  else if(testing$HSP70_ngmL.log[i] > 5.017524){
    testing$B3[i] = 3
  }
}




testing$D3 <- NA
for (i in 1:length(testing$uniqueid)){
  if (testing$FDP_ugmL[i] <= 0.27){
    testing$D3[i] = 1
  }
  else if (testing$FDP_ugmL[i] <= 0.48){
    testing$D3[i] = 2
  }
  else if (testing$FDP_ugmL[i] <= 0.88){
    testing$D3[i] = 3
  }
  else if(testing$FDP_ugmL[i] > 0.88){
    testing$D3[i] = 4
  }
}




testing$E3 <- NA
for (i in 1:length(testing$uniqueid)){
  if (testing$CRP_mgL[i] <= 1.6){
    testing$E3[i] = 1
  }
  else if (testing$CRP_mgL[i] <= 5.8){
    testing$E3[i] = 2
  }
  else if (testing$CRP_mgL[i] <= 16.9){
    testing$E3[i] = 3
  }
  else if(testing$CRP_mgL[i] > 16.9){
    testing$E3[i] = 4
  }
}



testing.f$C3 <- NA
for (i in 1:length(testing.f$uniqueid)){
  if (testing.f$suPAR_ngmL[i] <= 3.6){
    testing.f$C3[i] = 3
  }
  else if(testing.f$suPAR_ngmL[i] > 3.6){
    testing.f$C3[i] = 4
  }
}



testing.m$C3 <- NA
for (i in 1:length(testing.m$uniqueid)){
  if (testing.m$suPAR_ngmL[i]  <= 3.2){
    testing.m$C3[i] = 1
  }
  else if(testing.m$suPAR_ngmL[i] > 3.2){
    testing.m$C3[i] = 2
  }
}


for (i in 1:length(testing$uniqueid)){
  if (is.na(testing$male[i])) {
    testing$male[i] <- 0
  }
}

testing$C3 = NA
for (i in 1:length(dat$uniqueid)){
  if (testing$male[i] == 1 & testing$suPAR_ngmL[i] <= 3.2){
    testing$C3[i] <- 1
  }
  else if (testing$male[i] == 1 & testing$suPAR_ngmL[i] > 3.2){
    testing$C3[i] <- 2
  }
  else if (testing$male[i] == 0  & testing$suPAR_ngmL[i] <= 3.6){
    testing$C3[i] <- 1
  }
  else if (testing$male[i] == 0 & testing$suPAR_ngmL[i] > 3.6){
    testing$C3[i] <- 2
  }
}




testing$C3[which(testing$male==0 & testing$suPAR_ngmL <= 3.6)] <- 1
testing$C3[which(testing$male==0 & testing$suPAR_ngmL > 3.6)] <- 2

testing$C3[which(testing$male !=0 & testing$suPAR_ngmL <= 3.2)] <- 1
testing$C3[which(testing$male !=0 & testing$suPAR_ngmL > 3.2)] <- 2




table(testing$A3)
table(testing$B3)
table(testing$D3)
table(testing$E3)
table(testing.m$C3)
table(testing.f$C3)



KM_A3 <- survfit(SurvObj3~A3, data=testing, conf.type="log")
plot(KM_A3, col=c(1:5), main=("Death or MI by HS-Trop Levels"), ylim=c(0.5, 1))
legend("bottomleft", c("Low","Medium","High","Very High","Extremely High"), col=c(1:5), lty=1, bty="n")


KM_B3 <- survfit(SurvObj3~B3, data=testing, conf.type="log")
plot(KM_B3, col=c(1:4), main=("Death or MI by HSP-70"), ylim=c(0.5, 1))
legend("bottomleft", c("Low","Medium","High"), col=c(1:3), lty=1, bty="n")

KM_D3 <- survfit(SurvObj3~D3, data=testing, conf.type="log")
plot(KM_D3, col=c(1:4), main=("Death or MI by FDP"), ylim=c(0.5, 1))
legend("bottomleft", c("Low","Medium","High","Very High"), col=c(1:4), lty=1, bty="n")

KM_E3 <- survfit(SurvObj3~E3, data=testing, conf.type="log")
plot(KM_E3, col=c(1:4), main=("Death or MI by CRP"), ylim=c(0.5, 1))
legend("bottomleft", c("Low","Medium","High","Very High"), col=c(1:4), lty=1, bty="n")

KM_C3 <- survfit(SurvObj3~C3, data=testing, conf.type="log")
plot(KM_C3, col=c(1:2), main=("Death or MI by SuPAR"), ylim=c(0.5, 1))
legend("bottomleft", c("Low","High"), col=c(1:2), lty=1, bty="n")





table(testing$A3)
table(testing$B3)
table(testing$D3)
table(testing$E3)
table(testing.f$C3)
table(testing.m$C3)

C3t = c(testing.m$C3, testing.f$C3)


#************************* REGRESSION

#Age  can try categories (e.g., <60, 60-70, >70, etc.)
#BMI  can try to create categories (e.g., normal, overweight, obese, etc.)
#eGFR  eGFR<60 is defined as chronic kidney disease
#Gensini_OM0 ?an define severe CAD (binary) using Gensini score > 20
#male
#race (black vs. nonblack)
#Smoking (current, past, never) ?check if past can be combined with never smokers
#HTN
#DM
#HC
#MIHx
#HFhx
#everstat
#everARBACE
table(dat$BMI.Cat)

dat$age.Cat <- NA
for (i in 1:length(dat$uniqueid)){
  if (dat$age[i]  < 60){
    dat$age.Cat[i] = 1
  }
  else if(dat$age[i] <= 70){
    dat$age.Cat[i] = 2
  }
  else if(dat$age[i] > 70){
    dat$age.Cat[i] = 3
  }
}

dat$BMI.Cat <- NA
for (i in 1:length(dat$uniqueid)){
  if (dat$BMI_Impute[i]  < 18.5){
    dat$BMI.Cat[i] = 1
  }
  else if(dat$BMI_Impute[i] < 25){
    dat$BMI.Cat[i] = 2
  }
  else if(dat$BMI_Impute[i] < 30){
    dat$BMI.Cat[i] = 3
  }
  else if(dat$BMI_Impute[i] >= 30){
    dat$BMI.Cat[i] = 4
  }
}

dat$eGFR.Cat <- NA
for (i in 1:length(dat$uniqueid)){
  if (is.na(dat$eGFR_max120[i])) {
    dat$eGFR.Cat[i] = NA
  }
  else if (dat$eGFR_max120[i]  < 60){
    dat$eGFR.Cat[i] = 1
  }
  else if(dat$eGFR_max120[i] >= 60){
    dat$eGFR.Cat[i] = 2
  }
}

dat$Gensini.Cat <- NA
for (i in 1:length(dat$uniqueid)){
  if (is.na(dat$Gensini_OM0[i])) {
    dat$Gensini.Cat[i] = NA
  }
  else if (dat$Gensini_OM0[i]  < 20){
    dat$Gensini.Cat[i] = 1
  }
  else if(dat$Gensini_OM0[i] >= 20){
    dat$Gensini.Cat[i] = 2
  }
}

dat$Smoke.Cat <- NA
for (i in 1:length(dat$uniqueid)){
  if (dat$PastSmoking[i]  == 1){
    dat$Smoke.Cat[i] = 1
  }
  else if(dat$CurrentSmoking[i] == 1){
    dat$Smoke.Cat[i] = 2
  }
  else {dat$Smoke.Cat[i] = 0}
}

shapiro.test(dat$hsTroponin_pgmL.log)
shapiro.test(dat$HSP70_ngmL.log)
shapiro.test(dat$FDP_ugmL.log)
shapiro.test(dat$CRP_mgL.log)
shapiro.test(dat$suPAR_ngmL.log)

citation("stats")

#



KM <- coxph(SurvObj3 ~ hsTroponin_pgmL.log + HSP70_ngmL.log + FDP_ugmL.log + CRP_mgL.log + suPAR_ngmL.log
            , data=dat)
summary(KM)

KM <- coxph(SurvObj3 ~ age + BMI_Impute + eGFR.Cat + Gensini.Cat + male + 
              black + factor(Smoke.Cat) + HTN + DM + HC + MIhx + HFhx + EverStat+ 
              EverARBACE + hsTroponin_pgmL.log + HSP70_ngmL.log + FDP_ugmL.log + CRP_mgL.log + suPAR_ngmL.log
            , data=dat)
summary(KM)



dat$FDP_log_time = dat$FDP_ugmL.log * dat$timetoMI5yr

KM <- coxph(SurvObj3 ~ hsTroponin_pgmL.log + HSP70_ngmL.log + FDP_ugmL.log+ FDP_log_time + CRP_mgL.log + suPAR_ngmL.log
            , data=dat)
summary(KM)

KM <- coxph(SurvObj3 ~ age + BMI_Impute + eGFR.Cat + Gensini.Cat + male + 
              black + factor(Smoke.Cat) + HTN + DM + HC + MIhx + HFhx + EverStat+ 
              EverARBACE + hsTroponin_pgmL.log + HSP70_ngmL.log + FDP_ugmL.log + FDP_log_time + CRP_mgL.log + suPAR_ngmL.log
            , data=dat)
summary(KM)
basehaz(KM)

coxsnellres= dat$DeathorMI - (resid(KM,type="martingale"))

summary(dat$eGFR.Cat)

dat.complete = dat[which(is.na(dat$age + dat$BMI_Impute + dat$eGFR.Cat + dat$Gensini.Cat + dat$male + 
                           dat$black + dat$Smoke.Cat + dat$HTN + dat$DM + dat$HC + dat$MIhx + 
                             dat$HFhx + dat$EverStat+ 
                           dat$EverARBACE + dat$hsTroponin_pgmL.log + 
                             dat$HSP70_ngmL.log + dat$FDP_ugmL.log + dat$CRP_mgL.log + 
                             dat$suPAR_ngmL.log)==FALSE),]

coxsnellres= dat.complete$DeathorMI - (resid(KM,type="martingale"))

plot(coxsnellres)

FIT4=survfit(Surv(coxsnellres,dat.complete$DeathorMI)~1)
Htilde=cumsum(FIT4$n.event/FIT4$n.risk)
plot(FIT4$time,Htilde,type='s',col='blue',
     ylab = ("Estimated Cumulative Hazard"), 
     xlab = ("Cox-Snell Residuals")
     #,xlim=c(0,.5)
     )
abline(0,1,col='red',lty=2)


#


fit2=coxph(Surv(timetoMI5yr,DeathorMI)~age + BMI_Impute + eGFR.Cat + Gensini.Cat + male + 
             black + factor(Smoke.Cat) + HTN + DM + HC + MIhx + HFhx + EverStat+ 
             EverARBACE + hsTroponin_pgmL.log + HSP70_ngmL.log + FDP_ugmL.log  + CRP_mgL.log + suPAR_ngmL.log,
           data=dat,method='breslow')

resid(fit2,type='martingale')
plot(dat.complete$FDP_ugmL.log, resid(fit2), ylab=("Martingale Residuals"), xlab=("log(FDP (ug/mL)"), pch=16)
lines(lowess(dat.complete$FDP_ugmL.log, resid(fit2)),col='red')

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

# LETS SEE IF IT WORKS BETTER WITH CATEGORIES

KM <- coxph(SurvObj3 ~ age.Cat + BMI.Cat + eGFR.Cat + Gensini.Cat + male + 
              black + factor(Smoke.Cat) + HTN + DM + HC + MIhx + HFhx + EverStat+ 
              EverARBACE + hsTroponin_pgmL.log + HSP70_ngmL.log + FDP_ugmL.log + CRP_mgL.log + suPAR_ngmL.log
            , data=dat)
summary(KM)

coxsnellres= dat.complete$DeathorMI - (resid(KM,type="martingale"))

plot(coxsnellres)

FIT4=survfit(Surv(coxsnellres,dat.complete$DeathorMI)~1)
Htilde=cumsum(FIT4$n.event/FIT4$n.risk)
plot(FIT4$time,Htilde,type='s',col='blue',
     ylab = ("Estimated Cumulative Hazard"), 
     xlab = ("Cox-Snell Residuals")
     #,xlim=c(0,.5)
)
abline(0,1,col='red',lty=2)

#

fit2=coxph(Surv(timetoMI5yr,DeathorMI)~age.Cat + BMI.Cat + eGFR.Cat + Gensini.Cat + male + 
             black + factor(Smoke.Cat) + HTN + DM + HC + MIhx + HFhx + EverStat+ 
             EverARBACE + hsTroponin_pgmL.log + HSP70_ngmL.log + FDP_ugmL.log + CRP_mgL.log + suPAR_ngmL.log,
           data=dat,method='breslow')

resid(fit2,type='martingale')
plot(dat.complete$FDP_ugmL.log, resid(fit2))
lines(lowess(dat.complete$FDP_ugmL.log, resid(fit2)),col='red')


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

# HOW ABOUT SPLINES?

dat$AGESp1 <- (dat$age - 60)*as.numeric(dat$age >= 60)
dat$AGESp2 <- (dat$age - 70)*as.numeric(dat$age >= 70)

KM <- coxph(SurvObj3 ~ age + AGESp1 + AGESp2 + BMI.Cat + eGFR.Cat + Gensini.Cat + male + 
              black + factor(Smoke.Cat) + HTN + DM + HC + MIhx + HFhx + EverStat+ 
              EverARBACE + hsTroponin_pgmL.log + HSP70_ngmL.log + FDP_ugmL.log + CRP_mgL.log + suPAR_ngmL.log
            , data=dat)
sum <- summary(KM)



coxsnellres= dat.complete$DeathorMI - (resid(KM,type="martingale"))


FIT4=survfit(Surv(coxsnellres,dat.complete$DeathorMI)~1)
Htilde=cumsum(FIT4$n.event/FIT4$n.risk)
plot(FIT4$time,Htilde,type='s',col='blue',
     ylab = ("Estimated Cumulative Hazard"), 
     xlab = ("Cox-Snell Residuals")
     ,xlim=c(0,2.0)
)
abline(0,1,col='red',lty=2)


fit2=coxph(Surv(timetoMI5yr,DeathorMI)~ age + AGESp1+ AGESp2 + BMI.Cat + eGFR.Cat + Gensini.Cat + male + 
             black + factor(Smoke.Cat) + HTN + DM + HC + MIhx + HFhx + EverStat+ 
             EverARBACE + hsTroponin_pgmL.log + HSP70_ngmL.log + FDP_ugmL.log + CRP_mgL.log + suPAR_ngmL.log,
           data=dat,method='breslow')

resid(fit2,type='martingale')
plot(dat.complete$FDP_ugmL.log, resid(fit2))
lines(lowess(dat.complete$FDP_ugmL.log, resid(fit2)),col='red')

fit2


resid.schoenfeld <- residuals(KM, type = "schoenfeld")

plot(resid.schoenfeld[,1] ~ dat.complete$timetoMI5yr[which(dat.complete$DeathorMI == 1)])

par(mfrow=c(2,3))
plot(resid.schoenfeld[,18] ~ dat.complete$timetoMI5yr[which(dat.complete$DeathorMI == 1)], xlab = ("Time to Event"), ylab=("Partial Residual"), main = ("HS-Troponin"), pch=16)
plot(resid.schoenfeld[,19] ~ dat.complete$timetoMI5yr[which(dat.complete$DeathorMI == 1)], xlab = ("Time to Event"), ylab=("Partial Residual"), main = ("HSP 70"), pch=16)
plot(resid.schoenfeld[,20] ~ dat.complete$timetoMI5yr[which(dat.complete$DeathorMI == 1)], xlab = ("Time to Event"), ylab=("Partial Residual"), main = ("FDP"), pch=16)
plot(resid.schoenfeld[,21] ~ dat.complete$timetoMI5yr[which(dat.complete$DeathorMI == 1)], xlab = ("Time to Event"), ylab=("Partial Residual"), main = ("CRP"), pch=16)
plot(resid.schoenfeld[,22] ~ dat.complete$timetoMI5yr[which(dat.complete$DeathorMI == 1)], xlab = ("Time to Event"), ylab=("Partial Residual"), main = ("SuPAR"), pch=16)


summary(dat$HSP70_ngmL)

dat$HSP70_ngmL[which(dat$HSP70_ngmL<=7.072e-01)]

##########

par(mfrow=c(2,3))  

KM_A3 <- survfit(SurvObj3~A3, data=dat, conf.type="log")
plot(KM_A3, col=c(1:4,6), main=("Death or MI by HS-Trop Levels"), ylim=c(0.5, 1), ylab = ("Survival"), xlab = ("time to event (days)"))
legend("bottomleft", c("<1.6 (pg/mL)","1.6-2.1","2.1-2.6","2.6-6.9","≥6.9"), 
       col=c(1:4,6), lty=1, bty="n")

KM_B3 <- survfit(SurvObj3~B3, data=dat, conf.type="log")
plot(KM_B3, col=c(1:4), main=("Death or MI by HSP-70 Levels"), ylim=c(0.5, 1), ylab = ("Survival"), xlab = ("time to event (days)"))
legend("bottomleft", c("<1.0 (ng/mL)","1.0-151.0","≥151.0"), col=c(1:3), lty=1, bty="n")

KM_D3 <- survfit(SurvObj3~D3, data=dat, conf.type="log")
plot(KM_D3, col=c(1:4), main=("Death or MI by FDP"), ylim=c(0.5, 1), ylab = ("Survival"), xlab = ("time to event (days)"))
legend("bottomleft", c("<0.27 (ug/mL)","0.27-0.48","0.48-0.88", "≥0.88"), col=c(1:4), lty=1, bty="n")

KM_E3 <- survfit(SurvObj3~E3, data=dat, conf.type="log")
plot(KM_E3, col=c(1:4), main=("Death or MI by CRP"), ylim=c(0.5, 1), ylab = ("Survival"), xlab = ("time to event (days)"))
legend("bottomleft", c("<1.6 (ug/mL)","1.6-5.8","5.8-16.9", "≥16.9"), col=c(1:4), lty=1, bty="n")

KM_C3 <- survfit(SurvObj3~C3, data=dat, conf.type="log")
plot(KM_C3, col=c(1:2), main=("Death or MI by suPAR"), ylim=c(0.5, 1), ylab = ("Survival"), xlab = ("time to event (days)"))
legend("bottomleft", c("<3.2 (M), 3.6 (F)","≥3.2(M), 3.6(F) (ng/mL)"), col=c(1:2), lty=1, bty="n")










#******************************#
# Making the testing stuff too

testing$age.Cat <- NA
for (i in 1:length(testing$uniqueid)){
  if (testing$age[i]  < 60){
    testing$age.Cat[i] = 1
  }
  else if(testing$age[i] <= 70){
    testing$age.Cat[i] = 2
  }
  else if(testing$age[i] > 70){
    testing$age.Cat[i] = 3
  }
}



testing$BMI.Cat <- NA
for (i in 1:length(testing$uniqueid)){
  if (testing$BMI_Impute[i]  < 18.5){
    testing$BMI.Cat[i] = 1
  }
  else if(testing$BMI_Impute[i] < 25){
    testing$BMI.Cat[i] = 2
  }
  else if(testing$BMI_Impute[i] < 30){
    testing$BMI.Cat[i] = 3
  }
  else if(testing$BMI_Impute[i] >= 30){
    testing$BMI.Cat[i] = 4
  }
}

testing$eGFR.Cat <- NA
for (i in 1:length(testing$uniqueid)){
  if (is.na(testing$eGFR_max120[i])) {
    testing$eGFR.Cat[i] = NA
  }
  else if (testing$eGFR_max120[i]  < 60){
    testing$eGFR.Cat[i] = 1
  }
  else if(testing$eGFR_max120[i] >= 60){
    testing$eGFR.Cat[i] = 2
  }
}

testing$Gensini.Cat <- NA
for (i in 1:length(testing$uniqueid)){
  if (is.na(testing$Gensini_OM0[i])) {
    testing$Gensini.Cat[i] = NA
  }
  else if (testing$Gensini_OM0[i]  < 20){
    testing$Gensini.Cat[i] = 1
  }
  else if(testing$Gensini_OM0[i] >= 20){
    testing$Gensini.Cat[i] = 2
  }
}

testing$Smoke.Cat <- NA
for (i in 1:length(testing$uniqueid)){
  if (testing$PastSmoking[i]  == 1){
    testing$Smoke.Cat[i] = 1
  }
  else if(testing$CurrentSmoking[i] == 1){
    testing$Smoke.Cat[i] = 2
  }
  else {testing$Smoke.Cat[i] = 0}
}


KM <- coxph(SurvObj3 ~ age + BMI_Impute + eGFR.Cat + Gensini.Cat + male + 
              black + factor(Smoke.Cat) + HTN + DM + HC + MIhx + HFhx + EverStat+ 
              EverARBACE + hsTroponin_pgmL.log + HSP70_ngmL.log + FDP_ugmL.log + CRP_mgL.log + suPAR_ngmL.log
            , data=testing)


sumKM <- summary(KM)
sumKM
RISK.NoBRS <- exp(KM$linear.predictors)





#**************#
KMf <- coxph(SurvObj3 ~ factor(A3) + factor(B3) + factor(C3) + factor(D3) + factor(E3) 
              , data=testing)


coefs <- KMf$coefficients
coefs

BRS <- 0
for (i in 1:length(testing$uniqueid)){
  BRS[i] <- 0
  if (testing$A3[i] == 2) {
    BRS[i] <- BRS[i] + KMf$coefficients[1]
  }
  if (testing$A3[i] == 3) {
    BRS[i] <- BRS[i] + KMf$coefficients[2]
  }
  if (testing$A3[i] == 4) {
    BRS[i] <- BRS[i] + KMf$coefficients[3]
  }
  if (testing$A3[i] == 5) {
    BRS[i] <- BRS[i] + KMf$coefficients[4]
  }
  
  ##########################################
  
  if (testing$B3[i] == 2) {
    BRS[i] <- BRS[i] + KMf$coefficients[5]
  }
  if (testing$B3[i] == 3) {
    BRS[i] <- BRS[i] + KMf$coefficients[6]
  }
  
  ##########################################
  
  if (testing$D3[i] == 2) {
    BRS[i] <- BRS[i] + KMf$coefficients[7]
  }
  if (testing$D3[i] == 3) {
    BRS[i] <- BRS[i] + KMf$coefficients[8]
  }
  if (testing$D3[i] == 4) {
    BRS[i] <- BRS[i] + KMf$coefficients[9]
  }
  
  #########################################
  
  if (testing$E3[i] == 2) {
    BRS[i] <- BRS[i] + KMf$coefficients[10]
  }
  if (testing$E3[i] == 3) {
    BRS[i] <- BRS[i] + KMf$coefficients[11]
  }
  if (testing$E3[i] == 4) {
    BRS[i] <- BRS[i] + KMf$coefficients[12]
  }
}

hist(BRS)

KMscore <- coxph(testing$SurvObj3 ~ BRS)
sum <- summary(KMscore)

sum

BRS.CENTERED <- BRS-mean(BRS)
testing$BRS.CENTERED = BRS.CENTERED
######################################################
testing$age.Cat


list.of.packages <- c("foreign", "survC1", "survIDINRI")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
library(foreign)
library(survC1)
library(survIDINRI)

model1 <- c(BRS)
mydata = cbind(testing$timetoMI5yr, testing$DeathorMI, model1)

mod2dat <- cbind(testing$age.Cat , testing$BMI.Cat, testing$eGFR.Cat, testing$Gensini.Cat, testing$male, 
                   testing$black, testing$PastSmoking, testing$CurrentSmoking, testing$HTN, testing$DM, testing$HC, testing$MIhx,
                   testing$HFhx, testing$EverStat, testing$EverARBACE, testing$BRS.CENTERED)

head(mod2dat)

mydata = cbind(mydata, mod2dat)
head(mydata)

mydata[,1:2]

head(mydata)
Cfun = Est.Cval(mydata[,1:3], 1826, nofit=FALSE)
Cfun$Dhat 


mydatacomp = mydata[complete.cases(mydata),]

xout <- IDI.INF(indata = mydatacomp[,1:2],
        covs0  = mydatacomp[,4:18],
        covs1  = mydatacomp[,3:18],  
        t0   = 1826,                                             
        npert    = 200) 

IDI.INF.OUT(xout)
IDI.INF.GRAPH(xout)


C.12 = Inf.Cval.Delta(mydata = mydatacomp[,1:2],
                      covs0  = mydatacomp[,4:18],
                      covs1  = mydatacomp[,3:18],  
                      tau    = 1826,                                             
                      itr    = 20) 
C.12


################################################
mydatacomp <- as.data.frame(mydatacomp)
head(mydatacomp)
colnames(mydatacomp) <- c("timetoMI5yr", "DeathorMI", "BRS", "age.Cat" , "BMI.Cat", "eGFR.Cat", "Gensini.Cat", "male", 
                 "black", "PastSmoking", "CurrentSmoking", "HTN", "DM", "HC", "MIhx",
                 "HFhx", "EverStat", "EverARBACE", "BRS.CENTERED")

head(mydatacomp)



mydatacomp$SurvObj3 <- with(mydatacomp, Surv(mydatacomp$timetoMI5yr, mydatacomp$DeathorMI==1))

mydatacomp$Smoke.Cat <- NA
for (i in 1:length(mydatacomp$BRS)){
  if (mydatacomp$PastSmoking[i]  == 1){
    mydatacomp$Smoke.Cat[i] = 1
  }
  else if(mydatacomp$CurrentSmoking[i] == 1){
    mydatacomp$Smoke.Cat[i] = 2
  }
  else {mydatacomp$Smoke.Cat[i] = 0}
}


head(mydatacomp)
mydatacomp$age.Cat <- mydatacomp$age.Cat - 1
mydatacomp$BMI.Cat <- mydatacomp$BMI.Cat - 1
mydatacomp$Gensini.Cat <- mydatacomp$Gensini.Cat - 1
mydatacomp$eGFR.Cat <- mydatacomp$eGFR.Cat - 1


KM0 <- coxph(SurvObj3 ~ factor(age.Cat) + factor(BMI.Cat) + factor(eGFR.Cat) + factor(Gensini.Cat) + male + 
              black + factor(Smoke.Cat) + HTN + DM + HC + MIhx + HFhx + EverStat+ 
              EverARBACE
            , data=mydatacomp)

sumKM0 <- summary(KM0)
KM0_coeffs <- sumKM0$coefficients


RISK.NoBRS <- exp(KM$linear.predictors)

head(mydatacomp)
head(testing$SurvObj3)


KMfull <- coxph(SurvObj3 ~ factor(age.Cat) + factor(BMI.Cat) + factor(eGFR.Cat) + factor(Gensini.Cat) + male + 
               black + factor(Smoke.Cat) + HTN + DM + HC + MIhx + HFhx + EverStat+ 
               EverARBACE  + BRS
             , data=testing)

coeffs <- KMfull

KM1 <- coxph(SurvObj3 ~ factor(age.Cat) + factor(BMI.Cat) + factor(eGFR.Cat) + factor(Gensini.Cat) + male + 
              black + factor(Smoke.Cat) + HTN + DM + HC + MIhx + HFhx + EverStat+ 
              EverARBACE  + BRS.CENTERED
            , data=mydatacomp)

KM1

sumKM1 <- summary(KM1)
KM1_coeffs <- sumKM1$coefficients






sumKM1 <- summary(KM1)
cumhaz <- basehaz(KM1)
cumhaz <- cumhaz$hazard
BASELINESURVIVAL <- exp(-1*cumhaz)
BS5yr1 <- BASELINESURVIVAL[length(BASELINESURVIVAL)]

sumKM0 <- summary(KM0)
cumhaz <- basehaz(KM0)
cumhaz <- cumhaz$hazard
BASELINESURVIVAL <- exp(-1*cumhaz)
BS5yr0 <- BASELINESURVIVAL[length(BASELINESURVIVAL)]


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

# BASELINE SURVIVAL FINDINGS


head(mydatacomp)
KM0
BetaKM0 = rep(0,1345)

for (i in 1:1345){
  if (mydatacomp$age.Cat[i] == 2){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[1]
  }
  if (mydatacomp$age.Cat[i] == 3){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[2]
  }
  if (mydatacomp$BMI.Cat[i] == 2){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[3]
  }
  if (mydatacomp$BMI.Cat[i] == 3){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[4]
  }
  if (mydatacomp$age.Cat[i] == 4){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[5]
  }
  if (mydatacomp$eGFR.Cat[i] == 2){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[6]
  }
  if (mydatacomp$Gensini.Cat[i] == 2){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[7]
  }
  if (mydatacomp$male[i] == 1){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[8]
  }
  if (mydatacomp$black[i] == 1){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[9]
  }
  if (mydatacomp$Smoke.Cat[i] == 1){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[10]
  }
  if (mydatacomp$Smoke.Cat[i] == 2){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[11]
  }
  if (mydatacomp$HTN[i] == 1){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[12]
  }
  if (mydatacomp$DM[i] == 1){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[13]
  }
  if (mydatacomp$HC[i] == 1){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[14]
  }
  if (mydatacomp$MIhx[i] == 1){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[15]
  }
  if (mydatacomp$HFhx[i] == 1){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[16]
  }
  if (mydatacomp$EverStat[i] == 1){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[17]
  }
  if (mydatacomp$EverARBACE[i] == 1){
    BetaKM0[i] = BetaKM0[i]+KM0_coeffs[18]
  }
}

BetaKM0


Est.Risk0 <- rep(0,1345)
for (i in 1:1345){
  Est.Risk0[i] = 1-BS5yr0^(exp(BetaKM0[i]))
} 

mean(Est.Risk0)

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

head(mydatacomp)
KM1
BetaKM1 = rep(0,1345)

for (i in 1:1345){
  if (mydatacomp$age.Cat[i] == 1){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[1]
  }
  if (mydatacomp$age.Cat[i] == 2){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[2]
  }
  if (mydatacomp$BMI.Cat[i] == 1){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[3]
  }
  if (mydatacomp$BMI.Cat[i] == 2){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[4]
  }
  if (mydatacomp$age.Cat[i] == 3){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[5]
  }
  if (mydatacomp$eGFR.Cat[i] == 1){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[6]
  }
  if (mydatacomp$Gensini.Cat[i] == 1){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[7]
  }
  if (mydatacomp$male[i] == 1){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[8]
  }
  if (mydatacomp$black[i] == 1){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[9]
  }
  if (mydatacomp$Smoke.Cat[i] == 1){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[10]
  }
  if (mydatacomp$Smoke.Cat[i] == 2){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[11]
  }
  if (mydatacomp$HTN[i] == 1){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[12]
  }
  if (mydatacomp$DM[i] == 1){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[13]
  }
  if (mydatacomp$HC[i] == 1){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[14]
  }
  if (mydatacomp$MIhx[i] == 1){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[15]
  }
  if (mydatacomp$HFhx[i] == 1){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[16]
  }
  if (mydatacomp$EverStat[i] == 1){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[17]
  }
  if (mydatacomp$EverARBACE[i] == 1){
    BetaKM1[i] = BetaKM1[i]+KM1_coeffs[18]
  }
  
  BetaKM1[i] = BetaKM1[i] + (mydatacomp$BRS.CENTERED[i]*KM1_coeffs[19])
}


Est.Risk1 <- rep(0,1345)
for (i in 1:1345){
  Est.Risk1[i] = 1-BS5yr1^(exp(BetaKM1[i]))
} 

mean(Est.Risk1)

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

# MAKE THE 10 GROUPS


mydatacomp <- cbind(mydatacomp,Est.Risk0)
mydatacomp <- cbind(mydatacomp,Est.Risk1)



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

# OLD MODEL

mydataORDER0 <- mydatacomp[order(mydatacomp$Est.Risk0),]
head(mydataORDER0)

DEATHS.by.group0 = rep(0,10)
MEAN.EST.RISK.BY.GROUP0 = rep(0,10)

for (i in 1:10){
  DEATHS.by.group0[i] <- sum(mydataORDER0$DeathorMI[which(mydataORDER0$GROUPS10==i)]) / table(mydataORDER0$GROUPS10)[i]
  MEAN.EST.RISK.BY.GROUP0[i]  <- mean(mydataORDER0$Est.Risk0[which(mydataORDER0$GROUPS10==i)])
  
}

DEATHS.by.group0

plot(DEATHS.by.group0)

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

# GOOD MODEL



mydataORDER1 <- mydatacomp[order(mydatacomp$Est.Risk1),]
head(mydataORDER1)

mydataORDER1$GROUPS10 = 0

for (i in 1:135){
  mydataORDER1$GROUPS10[i] <- 1
}
for (i in 136:270){
  mydataORDER1$GROUPS10[i] <- 2
}
for (i in 271:405){
  mydataORDER1$GROUPS10[i] <- 3
}
for (i in 406:540){
  mydataORDER1$GROUPS10[i] <- 4
}
for (i in 541:675){
  mydataORDER1$GROUPS10[i] <- 5
}
for (i in 676:809){
  mydataORDER1$GROUPS10[i] <- 6
}
for (i in 810:943){
  mydataORDER1$GROUPS10[i] <- 7
}
for (i in 944:1077){
  mydataORDER1$GROUPS10[i] <- 8
}
for (i in 1078:1211){
  mydataORDER1$GROUPS10[i] <- 9
}
for (i in 1212:1345){
  mydataORDER1$GROUPS10[i] <- 10
}


#######

######





decilegroups = findInterval(Est.Risk1, quantile(Est.Risk1, probs=(0:9)/10))
ESTIMATED_RISK1 <- as.vector(tapply(Est.Risk1, decilegroups, mean))

mydatacomp <- cbind(mydatacomp, decilegroups)

KM_EST_RISK = rep(0,10)
for (i in 1:10){
  KMEST <- survfit(mydatacomp$SurvObj3[which(mydatacomp$decilegroups == i)]~1)
  sumEST <- summary(KMEST)
  KM_EST_RISK[i] <- 1-min(sumEST$surv)
}
KM_EST_RISK



mean(ESTIMATED_RISK1)
mean(KM_EST_RISK)
BOTH <- rbind(ESTIMATED_RISK1,KM_EST_RISK)

par(mfrow=c(1,1))
barplot(BOTH,beside=T, xlab=("Risk Group"),legend=c("Mean Est Risk", "KM Observed Risk"),args.legend = list(x = "topleft"),
        names.arg=c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"),
        ylab=("Proportion of Individuals with Event / Estimated Risk"))
abline(h=0)

mean(mydatacomp[,2])




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

# NRI?
head(mydatacomp)

mydataCENSORED <- mydatacomp[which((mydatacomp$DeathorMI==0 & mydatacomp$timetoMI5yr <1826)==FALSE),]
length(mydataCENSORED$timetoMI5yr)
##################




RISK.BRS <- exp(KM$linear.predictors)


plot(RISK.NoBRS, RISK.BRS)

head(testing)

Cases <- which(mydatacomp$DeathorMI==1)
NonCases <- which(mydatacomp$DeathorMI==0)

RISK.NoBRS[Cases]

plot(RISK.NoBRS[NonCases], RISK.BRS[NonCases], pch=16,col=1, xlim=c(0,10), ylim=c(0,10), ylab=("Risk, Model with BRS"),  xlab=("Risk, Model with no BRS"))
points(RISK.NoBRS[Cases], RISK.BRS[Cases], pch=16, col=3, xlim=c(0,10), ylim=c(0,10))

legend("topleft", c("Cases","Non-cases"), col=c(3,1), pch=16, bty="n")
abline(v=c(2,4))
abline(h=c(2,4))

RISK.Cases <- as.data.frame(cbind(RISK.NoBRS[Cases], RISK.BRS[Cases]))
RISK.NonCases <- as.data.frame(cbind(RISK.NoBRS[NonCases], RISK.BRS[NonCases]))

# NonCases

# Low to Low
length(RISK.NonCases$V1[which(RISK.NonCases$V1<2 & RISK.NonCases$V2<2)])/length(RISK.NonCases$V1)
# High to High
length(RISK.NonCases$V1[which(RISK.NonCases$V1>4 & RISK.NonCases$V2>4)])/length(RISK.NonCases$V1)
# Int to Int
length(RISK.NonCases$V1[which(RISK.NonCases$V1 >2 & RISK.NonCases$V1 <4 & RISK.NonCases$V2 >2 & RISK.NonCases$V2 <4)])/length(RISK.NonCases$V1)
# Int to Low
length(RISK.NonCases$V1[which(RISK.NonCases$V1 >2 & RISK.NonCases$V1 <4 & RISK.NonCases$V2 <2 )])/length(RISK.NonCases$V1)
# High to Low
length(RISK.NonCases$V1[which(RISK.NonCases$V1 >4 & RISK.NonCases$V2 <2 )])/length(RISK.NonCases$V1)
# High to Int
length(RISK.NonCases$V1[which(RISK.NonCases$V1 >4 & RISK.NonCases$V2 >2 & RISK.NonCases$V2 <4)])/length(RISK.NonCases$V1)
# Low to Int
length(RISK.NonCases$V1[which(RISK.NonCases$V1 <2 & RISK.NonCases$V2 >2 & RISK.NonCases$V2 <4)])/length(RISK.NonCases$V1)
# Low to High
length(RISK.NonCases$V1[which(RISK.NonCases$V1 <2 & RISK.NonCases$V2 >4)])/length(RISK.NonCases$V1)
# Int to High
length(RISK.NonCases$V1[which(RISK.NonCases$V1 >2 & RISK.NonCases$V1 <4 & RISK.NonCases$V2 >4)])/length(RISK.NonCases$V1)

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

plot(Est.Risk0, Est.Risk1)

head(dat.complete)






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

####################################### The iterations

######################################################
set.seed(122)

ss <- sample(1:5,size=length(testing$uniqueid),replace=TRUE,prob=c(1/5,1/5,1/5,1/5,1/5))
testing$GROUP = ss

testing$TEST.GROUP.1 <- NA
for (i in 1:length(testing$GROUP)){
  if (testing$GROUP[i] %in% c(2:5)) {
    testing$TEST.GROUP.1[i] <- 1
  }
  else{
    testing$TEST.GROUP.1[i] <- 2
  }
}


KM1 <- coxph(SurvObj3 ~ factor(A3) + factor(B3) + factor(C3) + factor(D3) + factor(E3) 
             , data=testing[which(testing$TEST.GROUP.1==1),])

table(testing$A3[which(testing$TEST.GROUP.1==1)])
table(testing$B3[which(testing$TEST.GROUP.1==1)])
table(testing$C3[which(testing$TEST.GROUP.1==1)])
table(testing$D3[which(testing$TEST.GROUP.1==1)])
table(testing$E3[which(testing$TEST.GROUP.1==1)])



KM1

KM1$coefficients



BRS1 <- 0
for (i in 1:length(testing$uniqueid)){
  BRS1[i] <- 0
   if (testing$A3[i] == 2) {
    BRS1[i] <- BRS1[i] + KM1$coefficients[1]
    }
    if (testing$A3[i] == 3) {
      BRS1[i] <- BRS1[i] + KM1$coefficients[2]
    }
    if (testing$A3[i] == 4) {
      BRS1[i] <- BRS1[i] + KM1$coefficients[3]
    }
    if (testing$A3[i] == 5) {
      BRS1[i] <- BRS1[i]  + KM1$coefficients[4]
    }
    
    ##########################################
    
    if (testing$B3[i] == 2) {
      BRS1[i] <- BRS1[i] + KM1$coefficients[5]
    }
    if (testing$B3[i] == 3) {
      BRS1[i] <- BRS1[i] + KM1$coefficients[6]
    }
    
    ##########################################
  if (testing$C3[i] == 2) {
    BRS1[i] <- BRS1[i] + KM1$coefficients[7]
  }
  
  
  
    if (testing$D3[i] == 2) {
      BRS1[i] <- BRS1[i] + KM1$coefficients[8]
    }
    if (testing$D3[i] == 3) {
      BRS1[i] <- BRS1[i] + KM1$coefficients[9]
    }
    if (testing$D3[i] == 4) {
      BRS1[i] <- BRS1[i] + KM1$coefficients[10]
    }
    
    #########################################
    
    if (testing$E3[i] == 2) {
      BRS1[i] <- BRS1[i] + KM1$coefficients[11]
    }
    if (testing$E3[i] == 3) {
      BRS1[i] <- BRS1[i] + KM1$coefficients[12]
    }
    if (testing$E3[i] == 4) {
      BRS1[i] <- BRS1[i] + KM1$coefficients[13]
  }
}

hist(BRS1)
testing <- cbind(testing,BRS1)
mydata = cbind(testing$timetoMI5yr[which(testing$TEST.GROUP.1==2)], 
               testing$DeathorMI[which(testing$TEST.GROUP.1==2)], 
               testing$BRS1[which(testing$TEST.GROUP.1==2)]
               )

mod2dat <- cbind(testing$age[which(testing$TEST.GROUP.1==2)] , 
                 testing$BMI_Impute[which(testing$TEST.GROUP.1==2)], 
                 testing$eGFR.Cat[which(testing$TEST.GROUP.1==2)], 
                 testing$Gensini.Cat[which(testing$TEST.GROUP.1==2)], 
                 testing$male[which(testing$TEST.GROUP.1==2)], 
                 testing$black[which(testing$TEST.GROUP.1==2)], 
                 testing$PastSmoking[which(testing$TEST.GROUP.1==2)], 
                 testing$CurrentSmoking[which(testing$TEST.GROUP.1==2)], 
                 testing$HTN[which(testing$TEST.GROUP.1==2)], 
                 testing$DM[which(testing$TEST.GROUP.1==2)], 
                 testing$HC[which(testing$TEST.GROUP.1==2)], 
                 testing$MIhx[which(testing$TEST.GROUP.1==2)],
                 testing$HFhx[which(testing$TEST.GROUP.1==2)], 
                 testing$EverStat[which(testing$TEST.GROUP.1==2)], 
                 testing$EverARBACE[which(testing$TEST.GROUP.1==2)])

head(mod2dat)

mydata = cbind(mydata, mod2dat)
head(mydata)

mydata[,1:2]

head(mydata)
Cfun = Est.Cval(mydata[,1:3], 1826, nofit=FALSE)
Cfun$Dhat 


mydata = mydata[complete.cases(mydata),]

C.Test1 = Inf.Cval.Delta(mydata = mydata[,1:2],
                      covs0  = mydata[,4:18],
                      covs1  = mydata[,3:18],  
                      tau    = 1826,                                             
                      itr    = 100) 
C.Test1

xout <- IDI.INF(indata = mydata[,1:2],
                covs0  = mydata[,4:18],
                covs1  = mydata[,3:18],  
                t0   = 1826,                                             
                npert    = 200) 

IDI.INF.OUT(xout)
#######################################################################

# Testing group 2


testing$TEST.GROUP.2 <- NA
for (i in 1:length(testing$GROUP)){
  if (testing$GROUP[i] %in% c(1,3,4,5)) {
    testing$TEST.GROUP.2[i] <- 1
  }
  else{
    testing$TEST.GROUP.2[i] <- 2
  }
}


KM1 <- coxph(SurvObj3 ~ factor(A3) + factor(B3) + factor(C3) + factor(D3) + factor(E3) 
             , data=testing[which(testing$TEST.GROUP.2==1),])

table(testing$A3[which(testing$TEST.GROUP.1==1)])
table(testing$B3[which(testing$TEST.GROUP.1==1)])
table(testing$C3[which(testing$TEST.GROUP.1==1)])
table(testing$D3[which(testing$TEST.GROUP.1==1)])
table(testing$E3[which(testing$TEST.GROUP.1==1)])



KM1

KM1$coefficients



BRS2 <- 0
for (i in 1:length(testing$uniqueid)){
  BRS2[i] <- 0
  if (testing$A3[i] == 2) {
    BRS2[i] <- BRS2[i] + KM1$coefficients[1]
  }
  if (testing$A3[i] == 3) {
    BRS2[i] <- BRS2[i] + KM1$coefficients[2]
  }
  if (testing$A3[i] == 4) {
    BRS2[i] <- BRS2[i] + KM1$coefficients[3]
  }
  if (testing$A3[i] == 5) {
    BRS2[i] <- BRS2[i]  + KM1$coefficients[4]
  }
  
  ##########################################
  
  if (testing$B3[i] == 2) {
    BRS2[i] <- BRS2[i] + KM1$coefficients[5]
  }
  if (testing$B3[i] == 3) {
    BRS2[i] <- BRS2[i] + KM1$coefficients[6]
  }
  
  ##########################################
  if (testing$C3[i] == 2) {
    BRS2[i] <- BRS2[i] + KM1$coefficients[7]
  }
  
  
  
  if (testing$D3[i] == 2) {
    BRS2[i] <- BRS2[i] + KM1$coefficients[8]
  }
  if (testing$D3[i] == 3) {
    BRS2[i] <- BRS2[i] + KM1$coefficients[9]
  }
  if (testing$D3[i] == 4) {
    BRS2[i] <- BRS2[i] + KM1$coefficients[10]
  }
  
  #########################################
  
  if (testing$E3[i] == 2) {
    BRS2[i] <- BRS2[i] + KM1$coefficients[11]
  }
  if (testing$E3[i] == 3) {
    BRS2[i] <- BRS2[i] + KM1$coefficients[12]
  }
  if (testing$E3[i] == 4) {
    BRS2[i] <- BRS2[i] + KM1$coefficients[13]
  }
}

hist(BRS2)
testing <- cbind(testing,BRS2)
mydata = cbind(testing$timetoMI5yr[which(testing$TEST.GROUP.2==2)], 
               testing$DeathorMI[which(testing$TEST.GROUP.2==2)], 
               testing$BRS2[which(testing$TEST.GROUP.2==2)]
)

mod2dat <- cbind(testing$age[which(testing$TEST.GROUP.2==2)] , 
                 testing$BMI_Impute[which(testing$TEST.GROUP.2==2)], 
                 testing$eGFR.Cat[which(testing$TEST.GROUP.2==2)], 
                 testing$Gensini.Cat[which(testing$TEST.GROUP.2==2)], 
                 testing$male[which(testing$TEST.GROUP.2==2)], 
                 testing$black[which(testing$TEST.GROUP.2==2)], 
                 testing$PastSmoking[which(testing$TEST.GROUP.2==2)], 
                 testing$CurrentSmoking[which(testing$TEST.GROUP.2==2)], 
                 testing$HTN[which(testing$TEST.GROUP.2==2)], 
                 testing$DM[which(testing$TEST.GROUP.2==2)], 
                 testing$HC[which(testing$TEST.GROUP.2==2)], 
                 testing$MIhx[which(testing$TEST.GROUP.2==2)],
                 testing$HFhx[which(testing$TEST.GROUP.2==2)], 
                 testing$EverStat[which(testing$TEST.GROUP.2==2)], 
                 testing$EverARBACE[which(testing$TEST.GROUP.2==2)])

head(mod2dat)

mydata = cbind(mydata, mod2dat)
head(mydata)

mydata[,1:2]


head(mydata)
Cfun = Est.Cval(mydata[,1:3], 1826, nofit=FALSE)
Cfun$Dhat 

mydata = mydata[complete.cases(mydata),]

C.Test2 = Inf.Cval.Delta(mydata = mydata[,1:2],
                         covs0  = mydata[,4:18],
                         covs1  = mydata[,3:18],  
                         tau    = 1826,                                             
                         itr    = 100) 
C.Test2

xout <- IDI.INF(indata = mydata[,1:2],
                covs0  = mydata[,4:18],
                covs1  = mydata[,3:18],  
                t0   = 1826,                                             
                npert    = 200) 

IDI.INF.OUT(xout)
#######################################################################

# Testing group 3


testing$TEST.GROUP.3 <- NA
for (i in 1:length(testing$GROUP)){
  if (testing$GROUP[i] %in% c(1,2,4,5)) {
    testing$TEST.GROUP.3[i] <- 1
  }
  else{
    testing$TEST.GROUP.3[i] <- 2
  }
}


KM1 <- coxph(SurvObj3 ~ factor(A3) + factor(B3) + factor(C3) + factor(D3) + factor(E3) 
             , data=testing[which(testing$TEST.GROUP.3==1),])

table(testing$A3[which(testing$TEST.GROUP.1==1)])
table(testing$B3[which(testing$TEST.GROUP.1==1)])
table(testing$C3[which(testing$TEST.GROUP.1==1)])
table(testing$D3[which(testing$TEST.GROUP.1==1)])
table(testing$E3[which(testing$TEST.GROUP.1==1)])



KM1

KM1$coefficients



BRS3 <- 0
for (i in 1:length(testing$uniqueid)){
  BRS3[i] <- 0
  if (testing$A3[i] == 2) {
    BRS3[i] <- BRS3[i] + KM1$coefficients[1]
  }
  if (testing$A3[i] == 3) {
    BRS3[i] <- BRS3[i] + KM1$coefficients[2]
  }
  if (testing$A3[i] == 4) {
    BRS3[i] <- BRS3[i] + KM1$coefficients[3]
  }
  if (testing$A3[i] == 5) {
    BRS3[i] <- BRS3[i]  + KM1$coefficients[4]
  }
  
  ##########################################
  
  if (testing$B3[i] == 2) {
    BRS3[i] <- BRS3[i] + KM1$coefficients[5]
  }
  if (testing$B3[i] == 3) {
    BRS3[i] <- BRS3[i] + KM1$coefficients[6]
  }
  
  ##########################################
  if (testing$C3[i] == 2) {
    BRS3[i] <- BRS3[i] + KM1$coefficients[7]
  }
  
  
  
  if (testing$D3[i] == 2) {
    BRS3[i] <- BRS3[i] + KM1$coefficients[8]
  }
  if (testing$D3[i] == 3) {
    BRS3[i] <- BRS3[i] + KM1$coefficients[9]
  }
  if (testing$D3[i] == 4) {
    BRS3[i] <- BRS3[i] + KM1$coefficients[10]
  }
  
  #########################################
  
  if (testing$E3[i] == 2) {
    BRS3[i] <- BRS3[i] + KM1$coefficients[11]
  }
  if (testing$E3[i] == 3) {
    BRS3[i] <- BRS3[i] + KM1$coefficients[12]
  }
  if (testing$E3[i] == 4) {
    BRS3[i] <- BRS3[i] + KM1$coefficients[13]
  }
}

hist(BRS3)
testing <- cbind(testing,BRS3)
mydata = cbind(testing$timetoMI5yr[which(testing$TEST.GROUP.3==2)], 
               testing$DeathorMI[which(testing$TEST.GROUP.3==2)], 
               testing$BRS3[which(testing$TEST.GROUP.3==2)]
)

mod2dat <- cbind(testing$age[which(testing$TEST.GROUP.3==2)] , 
                 testing$BMI_Impute[which(testing$TEST.GROUP.3==2)], 
                 testing$eGFR.Cat[which(testing$TEST.GROUP.3==2)], 
                 testing$Gensini.Cat[which(testing$TEST.GROUP.3==2)], 
                 testing$male[which(testing$TEST.GROUP.3==2)], 
                 testing$black[which(testing$TEST.GROUP.3==2)], 
                 testing$PastSmoking[which(testing$TEST.GROUP.3==2)], 
                 testing$CurrentSmoking[which(testing$TEST.GROUP.3==2)], 
                 testing$HTN[which(testing$TEST.GROUP.3==2)], 
                 testing$DM[which(testing$TEST.GROUP.3==2)], 
                 testing$HC[which(testing$TEST.GROUP.3==2)], 
                 testing$MIhx[which(testing$TEST.GROUP.3==2)],
                 testing$HFhx[which(testing$TEST.GROUP.3==2)], 
                 testing$EverStat[which(testing$TEST.GROUP.3==2)], 
                 testing$EverARBACE[which(testing$TEST.GROUP.3==2)])

head(mod2dat)

mydata = cbind(mydata, mod2dat)
head(mydata)

mydata[,1:2]


head(mydata)
Cfun = Est.Cval(mydata[,1:3], 1826, nofit=FALSE)
Cfun$Dhat 

mydata = mydata[complete.cases(mydata),]

C.Test3 = Inf.Cval.Delta(mydata = mydata[,1:2],
                         covs0  = mydata[,4:18],
                         covs1  = mydata[,3:18],  
                         tau    = 1826,                                             
                         itr    = 100) 
C.Test3


xout <- IDI.INF(indata = mydata[,1:2],
                covs0  = mydata[,4:18],
                covs1  = mydata[,3:18],  
                t0   = 1826,                                             
                npert    = 200) 

IDI.INF.OUT(xout)
#######################################################################

# Testing group 4


testing$TEST.GROUP.4 <- NA
for (i in 1:length(testing$GROUP)){
  if (testing$GROUP[i] %in% c(1,2,3,5)) {
    testing$TEST.GROUP.4[i] <- 1
  }
  else{
    testing$TEST.GROUP.4[i] <- 2
  }
}


KM1 <- coxph(SurvObj3 ~ factor(A3) + factor(B3) + factor(C3) + factor(D3) + factor(E3) 
             , data=testing[which(testing$TEST.GROUP.4==1),])

table(testing$A3[which(testing$TEST.GROUP.1==1)])
table(testing$B3[which(testing$TEST.GROUP.1==1)])
table(testing$C3[which(testing$TEST.GROUP.1==1)])
table(testing$D3[which(testing$TEST.GROUP.1==1)])
table(testing$E3[which(testing$TEST.GROUP.1==1)])



KM1

KM1$coefficients



BRS4 <- 0
for (i in 1:length(testing$uniqueid)){
  BRS4[i] <- 0
  if (testing$A3[i] == 2) {
    BRS4[i] <- BRS4[i] + KM1$coefficients[1]
  }
  if (testing$A3[i] == 3) {
    BRS4[i] <- BRS4[i] + KM1$coefficients[2]
  }
  if (testing$A3[i] == 4) {
    BRS4[i] <- BRS4[i] + KM1$coefficients[3]
  }
  if (testing$A3[i] == 5) {
    BRS4[i] <- BRS4[i]  + KM1$coefficients[4]
  }
  
  ##########################################
  
  if (testing$B3[i] == 2) {
    BRS4[i] <- BRS4[i] + KM1$coefficients[5]
  }
  if (testing$B3[i] == 3) {
    BRS4[i] <- BRS4[i] + KM1$coefficients[6]
  }
  
  ##########################################
  if (testing$C3[i] == 2) {
    BRS4[i] <- BRS4[i] + KM1$coefficients[7]
  }
  
  
  
  if (testing$D3[i] == 2) {
    BRS4[i] <- BRS4[i] + KM1$coefficients[8]
  }
  if (testing$D3[i] == 3) {
    BRS4[i] <- BRS4[i] + KM1$coefficients[9]
  }
  if (testing$D3[i] == 4) {
    BRS4[i] <- BRS4[i] + KM1$coefficients[10]
  }
  
  #########################################
  
  if (testing$E3[i] == 2) {
    BRS4[i] <- BRS4[i] + KM1$coefficients[11]
  }
  if (testing$E3[i] == 3) {
    BRS4[i] <- BRS4[i] + KM1$coefficients[12]
  }
  if (testing$E3[i] == 4) {
    BRS4[i] <- BRS4[i] + KM1$coefficients[13]
  }
}

hist(BRS4)
testing <- cbind(testing,BRS4)
mydata = cbind(testing$timetoMI5yr[which(testing$TEST.GROUP.4==2)], 
               testing$DeathorMI[which(testing$TEST.GROUP.4==2)], 
               testing$BRS4[which(testing$TEST.GROUP.4==2)]
)

mod2dat <- cbind(testing$age[which(testing$TEST.GROUP.4==2)] , 
                 testing$BMI_Impute[which(testing$TEST.GROUP.4==2)], 
                 testing$eGFR.Cat[which(testing$TEST.GROUP.4==2)], 
                 testing$Gensini.Cat[which(testing$TEST.GROUP.4==2)], 
                 testing$male[which(testing$TEST.GROUP.4==2)], 
                 testing$black[which(testing$TEST.GROUP.4==2)], 
                 testing$PastSmoking[which(testing$TEST.GROUP.4==2)], 
                 testing$CurrentSmoking[which(testing$TEST.GROUP.4==2)], 
                 testing$HTN[which(testing$TEST.GROUP.4==2)], 
                 testing$DM[which(testing$TEST.GROUP.4==2)], 
                 testing$HC[which(testing$TEST.GROUP.4==2)], 
                 testing$MIhx[which(testing$TEST.GROUP.4==2)],
                 testing$HFhx[which(testing$TEST.GROUP.4==2)], 
                 testing$EverStat[which(testing$TEST.GROUP.4==2)], 
                 testing$EverARBACE[which(testing$TEST.GROUP.4==2)])

head(mod2dat)

mydata = cbind(mydata, mod2dat)
head(mydata)

mydata[,1:2]


head(mydata)
Cfun = Est.Cval(mydata[,1:3], 1826, nofit=FALSE)
Cfun$Dhat 

mydata = mydata[complete.cases(mydata),]

C.Test4 = Inf.Cval.Delta(mydata = mydata[,1:2],
                         covs0  = mydata[,4:18],
                         covs1  = mydata[,3:18],  
                         tau    = 1826,                                             
                         itr    = 100) 
C.Test4

xout <- IDI.INF(indata = mydata[,1:2],
                covs0  = mydata[,4:18],
                covs1  = mydata[,3:18],  
                t0   = 1826,                                             
                npert    = 200) 

IDI.INF.OUT(xout)
#######################################################################

# Testing group 5


testing$TEST.GROUP.5 <- NA
for (i in 1:length(testing$GROUP)){
  if (testing$GROUP[i] %in% c(1,2,3,4)) {
    testing$TEST.GROUP.5[i] <- 1
  }
  else{
    testing$TEST.GROUP.5[i] <- 2
  }
}


KM1 <- coxph(SurvObj3 ~ factor(A3) + factor(B3) + factor(C3) + factor(D3) + factor(E3) 
             , data=testing[which(testing$TEST.GROUP.5==1),])

table(testing$A3[which(testing$TEST.GROUP.1==1)])
table(testing$B3[which(testing$TEST.GROUP.1==1)])
table(testing$C3[which(testing$TEST.GROUP.1==1)])
table(testing$D3[which(testing$TEST.GROUP.1==1)])
table(testing$E3[which(testing$TEST.GROUP.1==1)])



KM1

KM1$coefficients



BRS5 <- 0
for (i in 1:length(testing$uniqueid)){
  BRS5[i] <- 0
  if (testing$A3[i] == 2) {
    BRS5[i] <- BRS5[i] + KM1$coefficients[1]
  }
  if (testing$A3[i] == 3) {
    BRS5[i] <- BRS5[i] + KM1$coefficients[2]
  }
  if (testing$A3[i] == 4) {
    BRS5[i] <- BRS5[i] + KM1$coefficients[3]
  }
  if (testing$A3[i] == 5) {
    BRS5[i] <- BRS5[i]  + KM1$coefficients[4]
  }
  
  ##########################################
  
  if (testing$B3[i] == 2) {
    BRS5[i] <- BRS5[i] + KM1$coefficients[5]
  }
  if (testing$B3[i] == 3) {
    BRS5[i] <- BRS5[i] + KM1$coefficients[6]
  }
  
  ##########################################
  if (testing$C3[i] == 2) {
    BRS5[i] <- BRS5[i] + KM1$coefficients[7]
  }
  
  
  
  if (testing$D3[i] == 2) {
    BRS5[i] <- BRS5[i] + KM1$coefficients[8]
  }
  if (testing$D3[i] == 3) {
    BRS5[i] <- BRS5[i] + KM1$coefficients[9]
  }
  if (testing$D3[i] == 4) {
    BRS5[i] <- BRS5[i] + KM1$coefficients[10]
  }
  
  #########################################
  
  if (testing$E3[i] == 2) {
    BRS5[i] <- BRS5[i] + KM1$coefficients[11]
  }
  if (testing$E3[i] == 3) {
    BRS5[i] <- BRS5[i] + KM1$coefficients[12]
  }
  if (testing$E3[i] == 4) {
    BRS5[i] <- BRS5[i] + KM1$coefficients[13]
  }
}

hist(BRS5)
testing <- cbind(testing,BRS5)
mydata = cbind(testing$timetoMI5yr[which(testing$TEST.GROUP.5==2)], 
               testing$DeathorMI[which(testing$TEST.GROUP.5==2)], 
               testing$BRS5[which(testing$TEST.GROUP.5==2)]
)

mod2dat <- cbind(testing$age[which(testing$TEST.GROUP.5==2)] , 
                 testing$BMI_Impute[which(testing$TEST.GROUP.5==2)], 
                 testing$eGFR.Cat[which(testing$TEST.GROUP.5==2)], 
                 testing$Gensini.Cat[which(testing$TEST.GROUP.5==2)], 
                 testing$male[which(testing$TEST.GROUP.5==2)], 
                 testing$black[which(testing$TEST.GROUP.5==2)], 
                 testing$PastSmoking[which(testing$TEST.GROUP.5==2)], 
                 testing$CurrentSmoking[which(testing$TEST.GROUP.5==2)], 
                 testing$HTN[which(testing$TEST.GROUP.5==2)], 
                 testing$DM[which(testing$TEST.GROUP.5==2)], 
                 testing$HC[which(testing$TEST.GROUP.5==2)], 
                 testing$MIhx[which(testing$TEST.GROUP.5==2)],
                 testing$HFhx[which(testing$TEST.GROUP.5==2)], 
                 testing$EverStat[which(testing$TEST.GROUP.5==2)], 
                 testing$EverARBACE[which(testing$TEST.GROUP.5==2)])

head(mod2dat)

mydata = cbind(mydata, mod2dat)
head(mydata)

mydata[,1:2]


head(mydata)
Cfun = Est.Cval(mydata[,1:3], 1826, nofit=FALSE)
Cfun$Dhat 

mydata = mydata[complete.cases(mydata),]

C.Test5 = Inf.Cval.Delta(mydata = mydata[,1:2],
                         covs0  = mydata[,4:18],
                         covs1  = mydata[,3:18],  
                         tau    = 1826,                                             
                         itr    = 20) 
C.Test5



xout <- IDI.INF(indata = mydata[,1:2],
                covs0  = mydata[,4:18],
                covs1  = mydata[,3:18],  
                t0   = 1826,                                             
                npert    = 200) 

IDI.INF.OUT(xout)
#####################################


head(testing)


MOD1_Cstats <- c(C.Test1[1], C.Test2[1], C.Test3[1], C.Test4[1], C.Test5[1])

MOD0_Cstats <- c(C.Test1[2], C.Test2[2], C.Test3[2], C.Test4[2], C.Test5[2])

t.test(MOD1_Cstats, MOD0_Cstats, conf.level=.95)

mean(MOD1_Cstats - MOD0_Cstats)

#