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

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

####################
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$CVorMI[i] <- 0
  if (dat$death5yr[i] == 1){
    dat$DeathorMI[i] <- 1
  }
  else if (dat$MI5yr[i] == 1){
    dat$DeathorMI[i] <- 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==0)])/length(dat$HSP70_ngmL)
########## Demographics 
# Age
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

############# 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, na.rm=TRUE)
quantile(dat$Gensini_OM0, na.rm=TRUE)

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


###########################
###########################
########################### Distributions
###########################

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


# Hs Troponin

# Male vs Female
Male_hsTrop <- dat$hsTroponin_pgmL[which(dat$male == 1)]
Female_hsTrop <- dat$hsTroponin_pgmL[which(dat$male != 1)]
test_1 <- wilcox.test(Male_hsTrop, Female_hsTrop)
test_1$p.value

# Black VS Non-Black
Black_hsTrop <- dat$hsTroponin_pgmL[which(dat$black == 1)]
NonBlack_hsTrop <- dat$hsTroponin_pgmL[which(dat$black != 1)]
test_2 <- wilcox.test(Black_hsTrop, NonBlack_hsTrop)
test_2$p.value

# DM vs No DM
DM_hsTrop <- dat$hsTroponin_pgmL[which(dat$DM == 1)]
NoDM_hsTrop <- dat$hsTroponin_pgmL[which(dat$DM != 1)]
test_3 <- wilcox.test(DM_hsTrop, NoDM_hsTrop)
test_3$p.value

# Ckd vs No CKD
CKD_hsTrop <- dat$hsTroponin_pgmL[which(dat$ckd == 1)]
NoCKD_hsTrop <- dat$hsTroponin_pgmL[which(dat$ckd != 1)]
test_4 <- wilcox.test(CKD_hsTrop, NoCKD_hsTrop)
test_4$p.value

# MIhx vs No MIHx
MI_hsTrop <- dat$hsTroponin_pgmL[which(dat$MIhx == 1)]
NoMI_hsTrop <- dat$hsTroponin_pgmL[which(dat$MIhx != 1)]
test_5 <- wilcox.test(MI_hsTrop, NoMI_hsTrop)
test_5$p.value

# All hs Trop
boxplot(Male_hsTrop, 
        Female_hsTrop,
        Black_hsTrop, 
        NonBlack_hsTrop,
        DM_hsTrop, 
        NoDM_hsTrop,
        CKD_hsTrop, 
        NoCKD_hsTrop,
        MI_hsTrop, 
        NoMI_hsTrop,
        names = c("Male", "Female", "Black", "Non-Black", 
                  "DM", "No DM", "CKD", "No CKD", 
                  "MIhx", "No MIhx"),
        col = c(2,2,3,3,4,4,5,5,6,6),
        outline = FALSE,
        main=("HS Troponin Distributions"),
        ylim=c(0,40),
        ylab=("pg/mL")
        )
text(1.5,22,signif(test_1$p.value, digits= 3))
text(3.5,22,signif(test_2$p.value, digits= 3))
text(5.5,22,signif(test_3$p.value, digits= 3))
text(7.5,22,signif(test_4$p.value, digits= 3))
text(9.5,22,signif(test_5$p.value, digits= 3))


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

# HSP 70
quantile(dat$HSP70_ngmL, 0.80)


# Male vs Female
Male_hs70 <- dat$HSP70_ngmL[which(dat$male == 1)]
Female_hs70 <- dat$HSP70_ngmL[which(dat$male != 1)]
test_1 <- wilcox.test(Male_hs70, Female_hs70)
test_1$p.value

# Black VS Non-Black
Black_hs70 <- dat$HSP70_ngmL[which(dat$black == 1)]
NonBlack_hs70 <- dat$HSP70_ngmL[which(dat$black != 1)]
test_2 <- wilcox.test(Black_hs70, NonBlack_hs70)
test_2$p.value

# DM vs No DM
DM_hs70 <- dat$HSP70_ngmL[which(dat$DM == 1)]
NoDM_hs70 <- dat$HSP70_ngmL[which(dat$DM != 1)]
test_3 <- wilcox.test(DM_hs70, NoDM_hs70)
test_3$p.value

# Ckd vs No CKD
CKD_hs70 <- dat$HSP70_ngmL[which(dat$ckd == 1)]
NoCKD_hs70 <- dat$HSP70_ngmL[which(dat$ckd != 1)]
test_4 <- wilcox.test(CKD_hs70, NoCKD_hs70)
test_4$p.value

# MIhx vs No MIHx
MI_hs70 <- dat$HSP70_ngmL[which(dat$MIhx == 1)]
NoMI_hs70 <- dat$HSP70_ngmL[which(dat$MIhx != 1)]
test_5 <- wilcox.test(MI_hs70, NoMI_hs70)
test_5$p.value


# 
#########################
# FDP

# Male vs Female
Male_FDP <- dat$FDP_ugmL[which(dat$male == 1)]
Female_FDP <- dat$FDP_ugmL[which(dat$male != 1)]
test_1 <- wilcox.test(Male_FDP, Female_FDP)
test_1$p.value

# Black VS Non-Black
Black_FDP<- dat$FDP_ugmL[which(dat$black == 1)]
NonBlack_FDP <- dat$FDP_ugmL[which(dat$black != 1)]
test_2 <- wilcox.test(Black_FDP, NonBlack_FDP)
test_2$p.value

# DM vs No DM
DM_FDP <- dat$FDP_ugmL[which(dat$DM == 1)]
NoDM_FDP <- dat$FDP_ugmL[which(dat$DM != 1)]
test_3 <- wilcox.test(DM_FDP, NoDM_FDP)
test_3$p.value

# Ckd vs No CKD
CKD_FDP<- dat$FDP_ugmL[which(dat$ckd == 1)]
NoCKD_FDP <- dat$FDP_ugmL[which(dat$ckd != 1)]
test_4 <- wilcox.test(CKD_FDP, NoCKD_FDP)
test_4$p.value

# MIhx vs No MIHx
MI_FDP <- dat$FDP_ugmL[which(dat$MIhx == 1)]
NoMI_FDP <- dat$FDP_ugmL[which(dat$MIhx != 1)]
test_5 <- wilcox.test(MI_FDP, NoMI_FDP)
test_5$p.value

# All FDP
boxplot(Male_FDP, 
        Female_FDP,
        Black_FDP, 
        NonBlack_FDP,
        DM_FDP, 
        NoDM_FDP,
        CKD_FDP, 
        NoCKD_FDP,
        MI_FDP, 
        NoMI_FDP,
        names = c("Male", "Female", "Black", "Non-Black", 
                  "DM", "No DM", "CKD", "No CKD", 
                  "MIhx", "No MIhx"),
        col = c(2,2,3,3,4,4,5,5,6,6),
        outline = FALSE,
        main=("FDP Distributions"),
        ylim=c(0,2),
        ylab=("ug/mL")
)
text(1.5,1.2,signif(test_1$p.value, digits= 3))
text(3.5,1.2,signif(test_2$p.value, digits= 3))
text(5.5,1.2,signif(test_3$p.value, digits= 3))
text(7.5,1.2,signif(test_4$p.value, digits= 3))
text(9.5,1.2,signif(test_5$p.value, digits= 3))


###################
##################
quantile(dat$suPAR_ngmL)

# Male vs Female
Male_suPAR <- dat$suPAR_ngmL[which(dat$male == 1)]
Female_suPAR <- dat$suPAR_ngmL[which(dat$male != 1)]
test_1 <- wilcox.test(Male_suPAR, Female_suPAR)
test_1$p.value

# Black VS Non-Black
Black_suPAR<- dat$suPAR_ngmL[which(dat$black == 1)]
NonBlack_suPAR <- dat$suPAR_ngmL[which(dat$black != 1)]
test_2 <- wilcox.test(Black_suPAR, NonBlack_suPAR)
test_2$p.value

# DM vs No DM
DM_suPAR<- dat$suPAR_ngmL[which(dat$DM == 1)]
NoDM_suPAR <- dat$suPAR_ngmL[which(dat$DM != 1)]
test_3 <- wilcox.test(DM_suPAR, NoDM_suPAR)
test_3$p.value

# Ckd vs No CKD
CKD_suPAR <- dat$suPAR_ngmL[which(dat$ckd == 1)]
NoCKD_suPAR <- dat$suPAR_ngmL[which(dat$ckd != 1)]
test_4 <- wilcox.test(CKD_suPAR, NoCKD_suPAR)
test_4$p.value

# MIhx vs No MIHx
MI_suPAR <- dat$suPAR_ngmL[which(dat$MIhx == 1)]
NoMI_suPAR <- dat$suPAR_ngmL[which(dat$MIhx != 1)]
test_5 <- wilcox.test(MI_suPAR, NoMI_suPAR)
test_5$p.value

# All hs Trop
boxplot(Male_suPAR, 
        Female_suPAR,
        Black_suPAR, 
        NonBlack_suPAR,
        DM_suPAR, 
        NoDM_suPAR,
        CKD_suPAR, 
        NoCKD_suPAR,
        MI_suPAR, 
        NoMI_suPAR,
        names = c("Male", "Female", "Black", "Non-Black", 
                  "DM", "No DM", "CKD", "No CKD", 
                  "MIhx", "No MIhx"),
        col = c(2,2,3,3,4,4,5,5,6,6),
        outline = FALSE,
        main=("suPAR Distributions"),
        ylim=c(0,10),
        ylab=("ng/mL")
)
text(1.5,8,signif(test_1$p.value, digits= 3))
text(3.5,8,signif(test_2$p.value, digits= 3))
text(5.5,8,signif(test_3$p.value, digits= 3))
text(7.5,8,signif(test_4$p.value, digits= 3))
text(9.5,8,signif(test_5$p.value, digits= 3))


#################
#################
quantile(dat$CRP_mgL)

# Male vs Female
Male_CRP <- dat$CRP_mgL[which(dat$male == 1)]
Female_CRP <- dat$CRP_mgL[which(dat$male != 1)]
test_1 <- wilcox.test(Male_CRP, Female_CRP)
test_1$p.value

# Black VS Non-Black
Black_CRP<- dat$CRP_mgL[which(dat$black == 1)]
NonBlack_CRP <- dat$CRP_mgL[which(dat$black != 1)]
test_2 <- wilcox.test(Black_CRP, NonBlack_CRP)
test_2$p.value

# DM vs No DM
DM_CRP<- dat$CRP_mgL[which(dat$DM == 1)]
NoDM_CRP <- dat$CRP_mgL[which(dat$DM != 1)]
test_3 <- wilcox.test(DM_CRP, NoDM_CRP)
test_3$p.value

# Ckd vs No CKD
CKD_CRP <- dat$CRP_mgL[which(dat$ckd == 1)]
NoCKD_CRP <- dat$CRP_mgL[which(dat$ckd != 1)]
test_4 <- wilcox.test(CKD_CRP, NoCKD_CRP)
test_4$p.value

# MIhx vs No MIHx
MI_CRP <- dat$CRP_mgL[which(dat$MIhx == 1)]
NoMI_CRP <- dat$CRP_mgL[which(dat$MIhx != 1)]
test_5 <- wilcox.test(MI_CRP, NoMI_CRP)
test_5$p.value

# All hs Trop
boxplot(Male_CRP, 
        Female_CRP,
          Black_CRP, 
        NonBlack_CRP,
        DM_CRP, 
        NoDM_CRP,
        CKD_CRP, 
        NoCKD_CRP,
        MI_CRP, 
        NoMI_CRP,
        names = c("Male", "Female", "Black", "Non-Black", 
                  "DM", "No DM", "CKD", "No CKD", 
                  "MIhx", "No MIhx"),
        col = c(2,2,3,3,4,4,5,5,6,6),
        outline = FALSE,
        main=("CRP Distributions"),
        ylim=c(0,20),
        ylab=("mg/L")
)
text(1.5,10,signif(test_1$p.value, digits= 3))
text(3.5,10,signif(test_2$p.value, digits= 3))
text(5.5,10,signif(test_3$p.value, digits= 3))
text(7.5,10,signif(test_4$p.value, digits= 3))
text(9.5,10,signif(test_5$p.value, digits= 3))

#################
#################
# HSP Proportions of zero
length(dat$HSP70_ngmL[which(dat$male == 1 & dat$HSP70_ngmL != 0)]) / 
  length(dat$HSP70_ngmL[which(dat$male == 1)])
length(dat$HSP70_ngmL[which(dat$male != 1 & dat$HSP70_ngmL != 0)]) / 
  length(dat$HSP70_ngmL[which(dat$male != 1)])
# Males have 21.1 % Nonzero, Females have 20.3 % Nonzero

length(dat$HSP70_ngmL[which(dat$black == 1 & dat$HSP70_ngmL != 0)]) / 
  length(dat$HSP70_ngmL[which(dat$black == 1)])
length(dat$HSP70_ngmL[which(dat$black != 1 & dat$HSP70_ngmL != 0)]) / 
  length(dat$HSP70_ngmL[which(dat$black != 1)])
# Black people have 20.3 % Nonzero, Non-Black have 20.9 % Nonzero

length(dat$HSP70_ngmL[which(dat$DM == 1 & dat$HSP70_ngmL != 0)]) / 
  length(dat$HSP70_ngmL[which(dat$DM == 1)])
length(dat$HSP70_ngmL[which(dat$DM != 1 & dat$HSP70_ngmL != 0)]) / 
  length(dat$HSP70_ngmL[which(dat$DM != 1)])
# Those with DM have 24.2 % Nonzero, Non-DM have 19.2 % Nonzero

length(dat$HSP70_ngmL[which(dat$ckd == 1 & dat$HSP70_ngmL != 0)]) / 
  length(dat$HSP70_ngmL[which(dat$ckd == 1)])
length(dat$HSP70_ngmL[which(dat$ckd != 1 & dat$HSP70_ngmL != 0)]) / 
  length(dat$HSP70_ngmL[which(dat$ckd != 1)])
# CKD people have 26.5 % Nonzero, Non-ckd have 18.9 % Nonzero

length(dat$HSP70_ngmL[which(dat$MIhx == 1 & dat$HSP70_ngmL != 0)]) / 
  length(dat$HSP70_ngmL[which(dat$MIhx == 1)])
length(dat$HSP70_ngmL[which(dat$MIhx != 1 & dat$HSP70_ngmL != 0)]) / 
  length(dat$HSP70_ngmL[which(dat$MIhx != 1)])
# MIhx people have 22.1 % Nonzero, Non-MI have 20.4 % Nonzero

# Male vs Female
Male_HSP70 <- dat$HSP70_ngmL[which(dat$male == 1 & dat$HSP70_ngmL != 0)]
Female_HSP70 <- dat$HSP70_ngmL[which(dat$male != 1 & dat$HSP70_ngmL != 0)]
test_1 <- wilcox.test(Male_HSP70, Female_HSP70)
test_1$p.value

# Black VS Non-Black
Black_HSP70<- dat$HSP70_ngmL[which(dat$black == 1 & dat$HSP70_ngmL != 0)]
NonBlack_HSP70 <- dat$HSP70_ngmL[which(dat$black != 1 & dat$HSP70_ngmL != 0)]
test_2 <- wilcox.test(Black_HSP70, NonBlack_HSP70)
test_2$p.value

# DM vs No DM
DM_HSP70 <- dat$HSP70_ngmL[which(dat$DM == 1 & dat$HSP70_ngmL != 0)]
NoDM_HSP70 <- dat$HSP70_ngmL[which(dat$DM != 1 & dat$HSP70_ngmL != 0)]
test_3 <- wilcox.test(DM_HSP70, NoDM_HSP70)
test_3$p.value

# Ckd vs No CKD
CKD_HSP70 <- dat$HSP70_ngmL[which(dat$ckd == 1 & dat$HSP70_ngmL != 0)]
NoCKD_HSP70 <- dat$HSP70_ngmL[which(dat$ckd != 1 & dat$HSP70_ngmL != 0)]
test_4 <- wilcox.test(CKD_HSP70, NoCKD_HSP70)
test_4$p.value

# MIhx vs No MIHx
MI_HSP70 <- dat$HSP70_ngmL[which(dat$MIhx == 1 & dat$HSP70_ngmL != 0)]
NoMI_HSP70 <- dat$HSP70_ngmL[which(dat$MIhx != 1 & dat$HSP70_ngmL != 0)]
test_5 <- wilcox.test(MI_HSP70, NoMI_HSP70)
test_5$p.value

# All hs Trop
boxplot(Male_HSP70, 
        Female_HSP70,
        Black_HSP70, 
        NonBlack_HSP70,
        DM_HSP70, 
        NoDM_HSP70,
        CKD_HSP70, 
        NoCKD_HSP70,
        MI_HSP70, 
        NoMI_HSP70,
        names = c("Male", "Female", "Black", "Non-Black", 
                  "DM", "No DM", "CKD", "No CKD", 
                  "MIhx", "No MIhx"),
        col = c(2,2,3,3,4,4,5,5,6,6),
        outline = FALSE,
        main=("HSP-70 Distributions"),
        ylim=c(0,1000),
        ylab=("ng/mL")
)
text(1.5,600,signif(test_1$p.value, digits= 3))
text(3.5,600,signif(test_2$p.value, digits= 3))
text(5.5,600,signif(test_3$p.value, digits= 3))
text(7.5,600,signif(test_4$p.value, digits= 3))
text(9.5,600,signif(test_5$p.value, digits= 3))



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


# Survival Analysis 

library(survival)

# Outcome 1, All-Cause Death
length(dat$timetodeath5yr)
length(dat$death)

dat$SurvObj1 <- with(dat, Surv(timetodeath5yr, death5yr==1))
km1_ckd <- survfit(SurvObj1 ~ ckd, data = dat, conf.type = "log-log")
plot(km1_ckd, col = c(1,2), xlim=c(0,2000), main="All-cause Death", ylab=("S(x)"))
legend("bottomleft", c("CKD", "No CKD"), col = c(1,2), lty=1)

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$SurvObj2 <- with(dat, Surv(dat$timetoMI5yr, dat$CVorMI==1))
km2_ckd <- survfit(SurvObj2 ~ ckd, data = dat, conf.type = "log-log")
plot(km2_ckd, col = c(1,2), xlim=c(0,2000), main="CV death or MI", ylab=("S(x)"))
legend("bottomleft", c("CKD", "No CKD"), col = c(1,2), lty=1)

dat$DeathorMI <- 0
for (i in 1:length(dat$DeathorMI)){
  dat$CVorMI[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))
km3_ckd <- survfit(SurvObj3 ~ ckd, data = dat, conf.type = "log-log")
plot(km3_ckd, col = c(1,2), xlim=c(0,2000), main="All-cause death or MI", ylab=("S(x)"))
legend("bottomleft", c("CKD", "No CKD"), col = c(1,2), lty=1)


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


library(survival)
findcutnum = function (factor,outcome,datatype,nmin=20,segment=100)
{
  if (missing(factor)) 
    stop("The argument factor is missing")
  if (missing(outcome)) 
    stop("The argument outcome is missing")
  if (missing(datatype)) 
    stop("The argument datatype is missing")
  if (missing(datatype)) 
    stop("The argument datatype is missing")
  if (datatype=="survival"){
    if(!is.matrix(outcome)){
      stop("The outcome must be matrix")
    }
    if(dim(outcome)[2]!=2){
      stop("The outcome's column dimensions must be two ")
    }
  }
  if (datatype=="logistic"){
    if (!is.vector(outcome))
      stop("The argument outcome must be  vector")
  }
  if (datatype == "survival"){
    delta<-data.frame(outcome)
    colnames(delta)=c("event","time")
    userdata=na.omit(data.frame(delta$event,delta$time,factor))  #Collating of survival data
    colnames(userdata)=c("event","time","factor")
    index=order(userdata$factor)
    userdata=userdata[index,]
    n=dim(userdata)[1] #number of subject
    
    range=range(userdata$factor) #range of continous predictor
    cutunit=diff(range)/segment
    p=1
    coxfit = coxph(Surv(userdata$time,userdata$event) ~ userdata$factor)
    originAIC=(-2*coxfit$loglik[2])+2*p
    #originAICc=AIC+(2*p*(p+1)/(n-p-1))
    
    
    cut1AIC=c()
    start1=userdata$factor[nmin] #confirmed group1 has 20 person and where cut1 to start cutting
    end1=userdata$factor[n-nmin]#cut1 end
    while(start1<end1)
    {	   
      group1=sum(userdata$factor<=start1) #nummber of group1
      factor_status=c(rep(0,group1),rep(1,n-group1))
      coxfit = coxph(Surv(userdata$time,userdata$event) ~ factor(factor_status))
      AIC=(-2*coxfit$loglik[2])+2*p
      AICc=AIC+(2*p*(p+1)/(n-p-1))
      cut1AIC=c(cut1AIC,AIC)
      start1=start1+cutunit
    }
    #-------------------------------------cut2------------------------------------------------------------#
    p=2
    cut2AIC=c()
    start1=userdata$factor[nmin] #confirmed group1 has 20 person and where cut1 to start cutting
    end1=userdata$factor[n-nmin*2]#cut1 end  
    
    start2=userdata$factor[which(userdata$factor>start1)[nmin]] #confirmed group2 has 20 person and where cut2 to start cutting
    end2=userdata$factor[n-nmin] #cut2 end
    
    while(start1<end1)
    {
      group1=sum(userdata$factor<=start1) #nummber of group1
      while(start2<end2)
      {
        group2=sum(userdata$factor<=start2)-group1
        factor_status=c(rep(0,group1),rep(1,group2),rep(2,n-group1-group2))
        coxfit = coxph(Surv(userdata$time,userdata$event) ~ factor(factor_status))
        AIC=(-2*coxfit$loglik[2])+2*p
        AICc=AIC+(2*p*(p+1)/(n-p-1))
        cut2AIC=c(cut2AIC,AIC) 	   				
        start2=start2+cutunit
      }
      start1 = start1+cutunit
      start2 = userdata$factor[sum(userdata$factor<=start1)+nmin]
    }
    #-------------------------------------cut3------------------------------------------------------------#
    p=3
    cut3AIC=c()
    start1=userdata$factor[nmin] #confirmed group1 has 20 person and where cut1 to start cutting
    end1=userdata$factor[n-nmin*3]#cut1 end  
    
    start2=userdata$factor[which(userdata$factor>start1)[nmin]] #confirmed group2 has 20 person and where cut2 to start cutting
    end2=userdata$factor[n-nmin*2] #cut2 end
    
    start3=userdata$factor[which(userdata$factor>start2)[nmin]] #confirmed group2 has 20 person and where cut2 to start cutting
    end3=userdata$factor[n-nmin]
    
    while(start1<end1)
    {
      group1=sum(userdata$factor<=start1) #nummber of group1
      while(start2<end2)
      {
        group2=sum(userdata$factor<=start2)-group1
        while(start3<end3)
        {
          group3=sum(userdata$factor<=start3)-group1-group2
          factor_status=c(rep(0,group1),rep(1,group2),rep(2,group3),rep(3,n-group1-group2-group3))
          coxfit=coxph(Surv(userdata$time,userdata$event) ~ factor(factor_status))
          AIC=(-2*coxfit$loglik[2])+2*p
          AICc=AIC+(2*p*(p+1)/(n-p-1))
          start3=start3+cutunit
          cut3AIC=c(cut3AIC,AIC)
        }
        start2 = start2+cutunit
        start3 = userdata$factor[sum(userdata$factor<=start2)+nmin]
      }
      start1 = start1+cutunit
      start2 = userdata$factor[sum(userdata$factor<=start1)+nmin]
    }
    #-----------------------------------------cut4---------------------------------------------#
    p=4
    cut4AIC=c()
    start1=userdata$factor[nmin] #confirmed group1 has 20 person and where cut1 to start cutting
    end1=userdata$factor[n-nmin*4]#cut1 end  
    
    start2=userdata$factor[which(userdata$factor>start1)[nmin]] #confirmed group2 has 20 person and where cut2 to start cutting
    end2=userdata$factor[n-nmin*3] #cut2 end
    
    start3=userdata$factor[which(userdata$factor>start2)[nmin]] #confirmed group2 has 20 person and where cut2 to start cutting
    end3=userdata$factor[n-nmin*2]
    
    start4=userdata$factor[which(userdata$factor>start3)[nmin]] #confirmed group2 has 20 person and where cut2 to start cutting
    end4=userdata$factor[n-nmin]
    
    while(start1<end1)
    {
      group1=sum(userdata$factor<=start1) #nummber of group1
      while(start2<end2)
      {
        group2=sum(userdata$factor<=start2)-group1
        while(start3<end3)
        { 
          group2=sum(userdata$factor<=start2)-group1
          while(start4<end4)
          { 
            group4=sum(userdata$factor<=start4)-group1-group2-group3
            factor_status=c(rep(0,group1),rep(1,group2),rep(2,group3),rep(3,group4),rep(4,n-group1-group2-group3-group4))
            coxfit=coxph(Surv(userdata$time,userdata$event) ~ factor(factor_status))
            AIC=(-2*coxfit$loglik[2])+2*p
            AICc=AIC+(2*p*(p+1)/(n-p-1))
            start4=start4+cutunit
            cut4AIC=c(cut4AIC,AIC)
          }
          start3 = start3+cutunit
          start4 = userdata$factor[sum(userdata$factor<=start3)+nmin]
        }
        start2 = start2+cutunit
        start3 = userdata$factor[sum(userdata$factor<=start2)+nmin]
      }
      start1 = start1+cutunit
      start2 = userdata$factor[sum(userdata$factor<=start1)+nmin]
    }
    
  }#type=sur
  
  
  #-------------------------------------------------------------------------------------------------------------------#
  if (datatype == "logistic"){
    userdata=na.omit(data.frame(outcome,factor))
    index=order(userdata$factor)
    userdata=userdata[index,]
    
    n=dim(userdata)[1] #number of subject
    
    range=range(userdata$factor) #range of continous predictor
    cutunit=diff(range)/segment
    logist_red <- glm(userdata$outcome ~userdata$factor, family=binomial())
    originAIC=logist_red$aic
    
    #----------------------------------------cut 1 ---------------------------------------------------------#
    cut1AIC=c()
    start1=userdata$factor[nmin] #confirmed group1 has 20 person and where cut1 to start cutting
    end1=userdata$factor[n-nmin]#cut1 end  
    
    while(start1<end1)
    {
      group1=sum(userdata$factor<=start1) #nummber of group1
      factor_status=c(rep(0,group1),rep(1,n-group1))
      logist_fit =glm(userdata$outcome ~factor(factor_status), family=binomial())
      AIC=logist_fit$aic
      start1=start1+cutunit
      cut1AIC=c(cut1AIC,AIC) 
    }
    #----------------------------------------cut 2 ---------------------------------------------------------#
    cut2AIC=c()
    start1=userdata$factor[nmin] #confirmed group1 has 20 person and where cut1 to start cutting
    end1=userdata$factor[n-nmin*2]#cut1 end  
    
    start2=userdata$factor[which(userdata$factor>start1)[nmin]] #confirmed group2 has 20 person and where cut2 to start cutting
    end2=userdata$factor[n-nmin] #cut2 end
    while(start1<end1)
    {
      group1=sum(userdata$factor<=start1) #nummber of group1
      while(start2<end2)
      {
        group2=sum(userdata$factor<=start2)-group1
        factor_status=c(rep(0,group1),rep(1,group2),rep(2,n-group1-group2))
        logist_fit =glm(userdata$outcome ~factor(factor_status),family=binomial())
        start2=start2+cutunit
        AIC=logist_fit$aic
        cut2AIC=c(cut2AIC,AIC) 
      }
      
      start1 = start1+cutunit
      start2 = userdata$factor[sum(userdata$factor<=start1)+nmin]
    }
    #----------------------------------------cut 3 ---------------------------------------------------------#
    cut3AIC=c()
    start1=userdata$factor[nmin] #confirmed group1 has 20 person and where cut1 to start cutting
    end1=userdata$factor[n-nmin*3]#cut1 end  
    
    start2=userdata$factor[which(userdata$factor>start1)[nmin]] #confirmed group2 has 20 person and where cut2 to start cutting
    end2=userdata$factor[n-nmin*2] #cut2 end
    
    start3=userdata$factor[which(userdata$factor>start2)[nmin]] #confirmed group2 has 20 person and where cut2 to start cutting
    end3=userdata$factor[n-nmin]
    while(start1<end1)
    {
      group1=sum(userdata$factor<=start1) #nummber of group1
      
      while(start2<end2)
      {
        group2=sum(userdata$factor<=start2)-group1
        
        while(start3<end3)
        {
          group3=sum(userdata$factor<=start3)-group1-group2
          factor_status=c(rep(0,group1),rep(1,group2),rep(2,group3),rep(3,n-group1-group2-group3))
          logist_fit =glm(userdata$outcome ~factor(factor_status), family=binomial())
          AIC=logist_fit$aic
          cut3AIC=c(cut3AIC,AIC) 
          start3=start3+cutunit
          
        }
        
        start2 = start2+cutunit
        start3 = userdata$factor[sum(userdata$factor<=start2)+nmin]
        
      }
      
      start1 = start1+cutunit
      start2 = userdata$factor[sum(userdata$factor<=start1)+nmin]
    }
    #----------------------------------------cut 4 ---------------------------------------------------------#
    cut4AIC=c()
    start1=userdata$factor[nmin] #confirmed group1 has 20 person and where cut1 to start cutting
    end1=userdata$factor[n-nmin*4]#cut1 end  
    
    start2=userdata$factor[which(userdata$factor>start1)[nmin]] #confirmed group2 has 20 person and where cut2 to start cutting
    end2=userdata$factor[n-nmin*3] #cut2 end
    
    start3=userdata$factor[which(userdata$factor>start2)[nmin]] #confirmed group2 has 20 person and where cut2 to start cutting
    end3=userdata$factor[n-nmin*2]
    
    start4=userdata$factor[which(userdata$factor>start3)[nmin]] 
    end4=userdata$factor[n-nmin]
    while(start1<end1)
    {
      group1=sum(userdata$factor<=start1) #nummber of group1
      
      while(start2<end2)
      {
        group2=sum(userdata$factor<=start2)-group1
        
        while(start3<end3)
        {
          group3=sum(userdata$factor<=start3)-group1-group2
          while(start4<end4)
          {
            group4=sum(userdata$factor<=start4)-group1-group2-group3
            factor_status=c(rep(0,group1),rep(1,group2),rep(2,group3),rep(3,group4),rep(4,n-group1-group2-group3-group4))
            logist_fit =glm(userdata$outcome ~factor(factor_status), family=binomial())
            AIC=logist_fit$aic
            cut4AIC=c(cut4AIC,AIC) 
            start4=start4+cutunit
            
          }
          start3 = start3+cutunit
          start4 = userdata$factor[sum(userdata$factor<=start3)+nmin]
          
        }
        
        start2 = start2+cutunit
        start3 = userdata$factor[sum(userdata$factor<=start2)+nmin]
        
      }
      
      start1 = start1+cutunit
      start2 = userdata$factor[sum(userdata$factor<=start1)+nmin]
    }
    
  } #log end
  AICmin=which.min(c(min(originAIC),min(cut1AIC),min(cut2AIC),min(cut3AIC),min(cut4AIC)))
  cat(paste("AIC have a minimum value ","When cutnum = ",AICmin-1,"\n",sep=""))
  return(list(min=c(min(originAIC),min(cut1AIC),min(cut2AIC),min(cut3AIC),min(cut4AIC))))
}



##################
##################
##################
##################
##################
#The package was developed by Meng-Ke Hsieh, Wen-Yi Chang, and Chung Chang
#A step-by-step manual for using our package can be found here: http://www.math.nsysu.edu.tw/~cchang/cutoff/manual.doc

#findcut is a function that finds the optimal location of cutoff points. #The output is a list of locations of # This function is implemented in R; please save this file to findcut.r
#This function can deal with both survival outcome and binary one. 
#Although there are many criteria in this function, for survival outcome, we suggest using criteria of "likelihood ratio test" and "logrank test", and for binary one, using criteria of "AUC" and "Likelihood ratio test".

# input arguments for findcut: (Let N be the number of individuals) :findcut(factor,outcome,cutnum,datatype,nmin,segment)
#"factor" (Nx1 vector) is the continuous risk factor that needs to find the optimal location of cutoffs.
# cutnum (a postive integer) is the number of cutoffs;
# datatype (a string) indicates the type of data (for survival outcome, put "survival"; for binary one, put "binary"); 
# nmin (a positive integer) is the minimum number of individuals in each group and the default number is 20; segment is the total number of pieces and the default number is 100. 


#for survival outcome: outcome is a survival object (event, time), where 
# event (Nx1 vector) is the event indicator (1: event occurs, 0:right censoring) and  
#time (Nx1 vector) is the observed time (i.e. minimum of event time and censoring time). 
#Example for survival outcome: BMI vs overall survival (OS): findcut(BMI,(event,os),2,"survival",30,100) 
#Here event is the event indicator, os is the overall survial time;  2 cutoffs, minimum number of individuals in each group is 30; 100 is the number of pieces.

# for binary outcome: outcome is a binary variable (1: yes, 0: no)
# Example for binary outcome: fraction of tumor invasion (a risk factor) vs. lymphavascular space invasion (LVSI), a binary outcome. 
# findcut(invasion,LVSI,2,"logistic",30,100) 
# Here invasion is the risk factor that represents fraction of tumor invasion and LVSI is the binary outcome; number of cuts is 2; minimum number of individuals in each group is 30; 100 is the number of pieces.


findcut= function(factor=NULL,outcome=NULL,cutnum=NA,datatype=c("survival","logistic"),nmin=20,segment=100)
{   
  if (missing(factor)) 
    stop("The argument factor is missing")
  if (missing(outcome)) 
    stop("The argument outcome is missing")
  if (missing(datatype)) 
    stop("The argument datatype is missing")
  if (datatype=="survival"){
    if(!is.matrix(outcome)){
      stop("The outcome must be matrix")
    }
    if(dim(outcome)[2]!=2){
      stop("The outcome's column dimensions must be two ")
    }
  }
  if (datatype=="logistic"){
    if (!is.vector(outcome))
      stop("The argument outcome must be vector")
  }
  #if (cutnum>4)
  # stop("The argument cutnum must be less than or equal to 4")  
  z=qnorm(1-(1-0.95)/2)#caculate 95% confidence interval(cut1 need)
  if (datatype == "survival"){ 
    delta<-data.frame(outcome)
    colnames(delta)=c("event","time")
    userdata=na.omit(data.frame(delta$event,delta$time,factor))  #Collating of survival data
    colnames(userdata)=c("event","time","factor")
    index=order(userdata$factor)
    userdata=userdata[index,]
    n=dim(userdata)[1] #number of subject
    range=range(userdata$factor) #range of continous predictor
    cutunit=diff(range)/segment
    k=1 #count for result
    start <- c() 
    end <- c()   
    group <- rep(0,cutnum)#dummy variable
    ###############################
    if(cutnum==1){
      start[1]=userdata$factor[nmin]
      end[1]=userdata$factor[n-nmin*cutnum]
      result=matrix(NA,ncol=11,nrow=100000)
      colnames(result)=c("Cut1","Log-rank test","Gehan-Wilcoxon test",
                         "PH_assumption1","HR1","HRlow","HRup","P1","Likelihood ratio test","Waldtest","Scoretest")
      
    }else{
      for(i in 2:cutnum){ 
        start[1]=userdata$factor[nmin]
        end[1]=userdata$factor[n-nmin*cutnum]
        start[i]=userdata$factor[which(userdata$factor>start[i-1])[nmin]] #confirmed group2 has 20 person and where cut2 to start cutting
        end[i]=userdata$factor[n-nmin*(cutnum-i+1)] 
        result=matrix(NA,ncol=4*cutnum+5,nrow=100000)
        colnames(result)=c(paste("Cut",1:cutnum,sep=""),"Log-rank test","Gehan-Wilcoxon test",
                           paste("PH_assumption",1:cutnum,sep=""),paste("HR",1:cutnum,sep=""),paste("P",1:cutnum,sep=""),"Likelihood ratio test","Waldtest","Scoretest")}
    }
    process <- function(userdata,start,end,i,cutnum,group,n,cutunit,nmin,result,k){
      if(i==cutnum)
      {
        while(start[i]<end[i])
        {
          if((i-1)==0){group[i]=sum(userdata$factor<=start[i])} else {group[i]=sum(userdata$factor<=start[i])-sum(group[1:(i-1)])}
          j <- 0:cutnum
          factor_status=c(rep(j[1:(tail(j+1,1)-1)],group[1:(tail(j+1,1)-1)]),rep(tail(j,1),n-sum(group[1:tail(j,1)]))) 
          coxfit=coxph(Surv(userdata$time,userdata$event) ~ factor(factor_status))   #˘Xt1˘F?nĄLO-ĄPAI?O??
          model= summary(coxfit)
          coef = model$coefficients
          z=qnorm(1-(1-0.95)/2)
          lrtest_fit= survdiff(Surv(userdata$time,userdata$event)~ factor(factor_status)) #log-rank test
          wilcoxon_fit=survdiff(Surv(userdata$time,userdata$event)~factor(factor_status),rho=1)  #wilcoxon test
          if(i==1){
            result[k,]=c(start[1],1-pchisq(lrtest_fit$chisq,1), 1-pchisq(wilcoxon_fit$chisq,1)
                         ,cox.zph(coxfit)$table[,3],coef[,2],exp(coef[1] - z * coef[3]),exp(coef[1] + z * coef[3]),coef[5]
                         ,model[["logtest"]]["pvalue"],model[["waldtest"]]["pvalue"],model[["sctest"]]["pvalue"])
          }else{
            result[k,]=c(start[1:cutnum],1-pchisq(lrtest_fit$chisq,cutnum), 1-pchisq(wilcoxon_fit$chisq,cutnum)
                         ,cox.zph(coxfit)$table[,3][1:cutnum],coef[,2],coef[,5],model[["logtest"]]["pvalue"],model[["waldtest"]]["pvalue"],model[["sctest"]]["pvalue"])
          }
          start[i]=start[i]+cutunit
          k=k+1
        }
        if(cutnum==1){
          return(result)
        }else{
          return(list(start=start,group=group,result=result,k=k))
        }
      }
      else
      {
        while(start[i]<end[i])     
        {
          if(i==1){group[i]=sum(userdata$factor<=start[i])} else {group[i]=sum(userdata$factor<=start[i])-sum(group[1:(i-1)])}
          processin=process(userdata,start,end,i+1,cutnum,group,n,cutunit,nmin,result,k)
          group=processin$group
          start=processin$start
          result=processin$result
          k=processin$k;
          start[i] = start[i]+cutunit
          start[i+1] = userdata$factor[sum(userdata$factor<=start[i])+nmin]
        }
        if(i!=1)
        {
          return(list(start=start,group=group,result=result,k=k))
        }
        else
        {
          return(result)
        }
      }
    }      
    allcut=na.omit(data.frame(process(userdata,start,end,i=1,cutnum,group,n,cutunit,nmin,result,k)))
    if(cutnum==1){
      par(mfrow=c(1,2))  
      plotdata=data.frame(cut=allcut[,"Cut1"],HR=allcut[,"HR1"],HRlow=allcut[,"HRlow"],HRup=allcut[,"HRup"],p=allcut[,"Log.rank.test"])
      plot(plotdata$cut,plotdata$HR,type="l",lwd=2,xlab="Cutoff point",ylab="Hazard ratio",main="Cutoff point and Hazard Ratio",ylim=c(min(allcut[,"HRlow"]),max(allcut[,"HR1"])))
      lines(plotdata$cut,plotdata$HRlow,lty=2)
      lines(plotdata$cut,plotdata$HRup,lty=2)
      plot(plotdata$cut,plotdata$p,type="l",lwd=2,xlab="Cutoff point",ylab="p-value",main="Cutoff point and p-value")
      phokin=which(allcut[,"PH_assumption1"]>0.05)
    }else{
      ph <- paste("PH_assumption",1:cutnum,sep="") 
      phok<-matrix(data=NA,ncol=10000,nrow=4)
      for(i in 1:cutnum)
      {
        maxph <- length(which(allcut[,ph[i]]>0.05))
        for(x in 1:maxph){
          phok[i,x]=which(allcut[,ph[i]]>0.05)[x]
        }
      }
      phokin <- as.vector(phok[1:cutnum,])
      for(i in 1:(cutnum-1)){
        phokin=na.exclude(intersect(intersect(phok[i,],phokin),phok[i+1,]))
      }
    }
    ################################################################################
    bestcut <- matrix(NA,ncol=cutnum,nrow=1)
    colnames(bestcut)=c(paste("Cut",1:cutnum,sep=""))
    cutpoint <- c(allcut[which.min(allcut[,"Log.rank.test"]),][1:cutnum],allcut[which.min(allcut[,"Gehan.Wilcoxon.test"]),][1:cutnum],
                  allcut[phokin[which.min(allcut[phokin,"Likelihood.ratio.test"])],][1:cutnum],allcut[phokin[which.min(allcut[phokin,"Waldtest"])],][1:cutnum],
                  allcut[phokin[which.min(allcut[phokin,"Scoretest"])],][1:cutnum])
    cutmatrix <- matrix(as.numeric(cutpoint),ncol=cutnum,byrow=T)
    count=c()
    for(i in 1:5){
      count[i]=paste(cutmatrix[i,],collapse = "")
    }
    tablenumber <- data.frame(table(count))
    bestcut[1,] <- cutmatrix[which(count==tablenumber[which.max(tablenumber$Freq),1])[1],]
    if(cutnum==1){
      output=list(allcut=allcut,
                  logranktest=allcut[which.min(allcut[,"Log.rank.test"]),], 
                  wilcoxon=allcut[which.min(allcut[,"Gehan.Wilcoxon.test"]),],       
                  logtest=allcut[phokin[which.min(allcut[phokin,"Likelihood.ratio.test"])],],   
                  waldtest=allcut[phokin[which.min(allcut[phokin,"Waldtest"])],],               
                  scoretest=allcut[phokin[which.min(allcut[phokin,"Scoretest"])],],             
                  HRabsmax=allcut[phokin[which.max(abs(allcut[phokin,"HR1"]-1))],])
      #,bestcut=bestcut)                 
    } 
    else {
      output=list(allcut=allcut,
                  logranktest=allcut[which.min(allcut[,"Log.rank.test"]),],
                  wilcoxon=allcut[which.min(allcut[,"Gehan.Wilcoxon.test"]),],       
                  logtest_likelihood.ratio.test=allcut[phokin[which.min(allcut[phokin,"Likelihood.ratio.test"])],],
                  waldtest=allcut[phokin[which.min(allcut[phokin,"Waldtest"])],],
                  scoretest=allcut[phokin[which.min(allcut[phokin,"Scoretest"])],])
      #,bestcut=bestcut)
    }
    return(output)
  }#datatype=survival end
  #----------------------------------------logit regression---------------------------------------------------------#
  if (datatype == "logistic"){
    userdata=na.omit(data.frame(outcome,factor))
    index=order(userdata$factor)
    userdata=userdata[index,]
    n=dim(userdata)[1] #number of subject
    range=range(userdata$factor) #range of continous predictor
    cutunit=diff(range)/segment
    logist_red <- glm(userdata$outcome ~1, family=binomial())# no other variable so the natural prediction is the mean of userdata$outcome ,reduce model
    k=1 #count for result
    start <- c() 
    end <- c()   
    group <- rep(0,cutnum)  #dummy variable
    #----------------------------------------cut---------------------------------------------------------#
    if(cutnum==1){
      start[1]=userdata$factor[nmin]
      end[1]=userdata$factor[n-nmin*cutnum]
      result=matrix(NA,ncol=13,nrow=100000)
      colnames(result)=c("Cut1","OR1","OR_up","OR_low","P1","Likelihood ratio test","Waldtest","Scoretest","Specificity","Sensitivity","Youden",
                         "Fisher's Exact Test","AUC")
    }else{
      for(i in 2:cutnum){ 
        start[1]=userdata$factor[nmin]
        end[1]=userdata$factor[n-nmin*cutnum]
        start[i]=userdata$factor[which(userdata$factor>start[i-1])[nmin]] #confirmed group2 has 20 person and where cut2 to start cutting
        end[i]=userdata$factor[n-nmin*(cutnum-i+1)] 
        result=matrix(NA,ncol=3*cutnum+11,nrow=100000)
        colnames(result)=c(paste("Cut",1:cutnum,sep=""),paste("OR",1:cutnum,sep=""),paste("P",1:cutnum,sep=""),"Likelihood ratio test","Waldtest","Scoretest","Specificity","Sensitivity","Lower","Higher","Youden","euclidean","manhattan","AUC")
      }
    }
    process <- function(userdata,start,end,i,cutnum,group,n,cutunit,nmin,result,k){
      if(i==cutnum)
      { 
        while(start[i]<end[i])
        {
          if(i==1){group[i]=sum(userdata$factor<=start[i])} else {group[i]=sum(userdata$factor<=start[i])-sum(group[1:(i-1)])}
          j <- 0:cutnum
          factor_status=c(rep(j[1:(tail(j+1,1)-1)],group[1:cutnum]),rep(tail(j,1),n-sum(group[1:tail(j,1)])))  
          logist_fit =glm(userdata$outcome ~factor(factor_status), family=binomial())             
          model= summary(logist_fit)
          coef = model$coefficients
          table<- table(userdata$outcome, factor_status)  
          Fisher.test=fisher.test(table)$p.value     
          roc <- roc(userdata$outcome,factor_status) 
          spe=roc$specificities*100 
          sen=roc$sensitivities*100 
          euc <- sqrt((rep(100,cutnum+2)-spe)^2+(rep(100,cutnum+2)-sen)^2)   
          manhan <- abs(rep(100,cutnum+2)-spe)+abs(sen-rep(100,cutnum+2))  
          LRT=anova(logist_fit,logist_red,test="Chisq")                      
          Scoretest=anova(logist_fit,logist_red,test="Rao")                  
          z=qnorm(1-(1-0.95)/2)
          if(i==1){
            Waldtest=wald.test(b = coef(logist_fit), Sigma = vcov(logist_fit), Terms = 2)
            result[k,]=c(start[1],exp(coef[2]),exp((coef[2]+z*coef[4])),exp((coef[2]-z*coef[4])),coef[8],LRT$P[2],Waldtest$result$chi2[3],Scoretest$P[2],
                         spe[2],sen[2],sen[2]+spe[2]-1,Fisher.test,as.numeric(roc$auc))
          }else{
            Waldtest=wald.test(b = coef(logist_fit), Sigma = vcov(logist_fit), Terms = 2:(cutnum+1))
            maxroc=which.max(sen+spe)    
            mineuc=which.min(euc)        
            minman=which.min(manhan)  
            result[k,]=c(start[1:cutnum],exp(coef[2:(cutnum+1)]),coef[(3*cutnum+5):(4*cutnum+4)],LRT$P[2], Waldtest$result$chi2[3],Scoretest$P[2]
                         ,spe[maxroc],sen[maxroc],sen[maxroc]-100*z*sqrt((sen[maxroc]/100)*(1-sen[maxroc]/100)/(sum(na.exclude(outcome)==1)))
                         ,sen[maxroc]+100*z*sqrt((sen[maxroc]/100)*(1-sen[maxroc]/100)/(sum(na.exclude(outcome)==1)))
                         ,sen[maxroc]+spe[maxroc]-1,sqrt((sen[mineuc]-100)^2 + (spe[mineuc]-100)^2),abs(sen[minman]-100)+ abs(100-spe[minman]),as.numeric(roc$auc))
            
          }
          start[i]=start[i]+cutunit
          k=k+1
        }
        if(cutnum==1){
          return(result)
        }else{
          return(list(start=start,group=group,result=result,k=k))
        }
      }
      else
      {
        while(start[i]<end[i])
        {
          if(i==1){group[i]=sum(userdata$factor<=start[i])} else {group[i]=sum(userdata$factor<=start[i])-sum(group[1:(i-1)])}
          processin=process(userdata,start,end,i+1,cutnum,group,n,cutunit,nmin,result,k)
          group=processin$group
          start=processin$start
          result=processin$result
          k=processin$k;
          start[i] = start[i]+cutunit
          start[i+1] = userdata$factor[sum(userdata$factor<=start[i])+nmin]
        }
        if(i!=1)
        {
          return(list(start=start,group=group,result=result,k=k))
        }
        else
        {
          return(result)
        }
      }
    }      
    
    allcut=na.omit(data.frame(process(userdata,start,end,i=1,cutnum,group,n,cutunit,nmin,result,k)))
    if(cutnum==1){
      par(mfrow=c(1,2))
      plotdata=data.frame(cut=allcut[,"Cut1"],OR=allcut[,"OR1"],ORup=allcut[,"OR_up"],ORlow=allcut[,"OR_low"],p=allcut[,"Likelihood.ratio.test"])
      plot(plotdata$cut,plotdata$OR,type="l",lwd=2,xlab="Cutoff point",ylab="Odds ratio",main="Cutoff point and Odds Ratio", ylim=c(min(allcut[,"OR_low"]),max(allcut[,"OR1"])))
      lines(plotdata$cut,plotdata$ORlow,lty=2)
      lines(plotdata$cut,plotdata$ORup,lty=2)
      plot(plotdata$cut,plotdata$p,type="l",lwd=2,xlab="Cutoff point",ylab="p-value",main="Cutoff point and p-value")
      allcut$euclidean=sqrt((100-allcut$Sensitivity)^2 + (100-allcut$Specificity)^2)
      allcut$manhattan=100-allcut$Sensitivity+ 100-allcut$Specificity
    }
    #-----------------------------------------------------------------------------------------------#
    bestcut <- matrix(NA,ncol=cutnum,nrow=1) 
    colnames(bestcut)=c(paste("Cut",1:cutnum,sep=""))
    if(cutnum==1){
      cutpoint <- c(allcut[which.min(allcut[,"Likelihood.ratio.test"]),][1:cutnum],allcut[which.min(allcut[,"Waldtest"]),][1:cutnum],
                    allcut[which.min(allcut[,"Scoretest"]),][1:cutnum],allcut[which.min(allcut$Fisher.s.Exact.Test),][1:cutnum],allcut[which.max(allcut[,"Youden"]),][1:cutnum],
                    allcut[which.min(allcut[,"euclidean"]),][1:cutnum],allcut[which.max(allcut[,"manhattan"]),][1:cutnum],allcut[which.max(allcut[,"AUC"]),][1:cutnum])
    }
    else{
      cutpoint <- c(allcut[which.min(allcut[,"Likelihood.ratio.test"]),][1:cutnum],allcut[which.min(allcut[,"Waldtest"]),][1:cutnum],
                    allcut[which.min(allcut[,"Scoretest"]),][1:cutnum],allcut[which.max(allcut[,"Youden"]),][1:cutnum],
                    allcut[which.min(allcut[,"euclidean"]),][1:cutnum],allcut[which.max(allcut[,"manhattan"]),][1:cutnum],allcut[which.max(allcut[,"AUC"]),][1:cutnum])
    }
    cutmatrix <- matrix(as.numeric(cutpoint),ncol=cutnum,byrow=T)
    count=c()
    if(cutnum==1){
      for(i in 1:8){
        count[i]=paste(cutmatrix[i,],collapse = "")
      }
    }
    else{
      for(i in 1:7){
        count[i]=paste(cutmatrix[i,],collapse = "")
      }
    }
    tablenumber <- data.frame(table(count))
    bestcut[1,] <- cutmatrix[which(count==tablenumber[which.max(tablenumber$Freq),1])[1],]
    if(cutnum==1){
      output=list(allcut=allcut,
                  Likelihood.ratio.test=allcut[which.min(allcut[,"Likelihood.ratio.test"]),], 
                  waldtest=allcut[which.min(allcut[,"Waldtest"]),],                           
                  scoretest=allcut[which.min(allcut[,"Scoretest"]),],                       
                  ORabsmax=allcut[which.max(abs(allcut$OR1)),],   
                  youden=allcut[which.max(allcut$Youden),],       
                  Fisher.test=allcut[which.min(allcut$Fisher.s.Exact.Test),],  
                  euclidean=allcut[which.min(allcut$euclidean),],    
                  manhattan=allcut[which.min(allcut$manhattan),],    
                  AUC=allcut[which.max(allcut[,"AUC"]),])
      #,bestcut=bestcut)
    }
    else 
    {
      output=list(allcut=allcut,
                  Likelihood.ratio.test=allcut[which.min(allcut[,"Likelihood.ratio.test"]),],
                  waldtest=allcut[which.min(allcut[,"Waldtest"]),],
                  scoretest=allcut[which.min(allcut[,"Scoretest"]),],
                  youden=allcut[which.max(allcut[,"Youden"]),],
                  euclidean=allcut[which.min(allcut[,"euclidean"]),],
                  manhattan=allcut[which.min(allcut[,"manhattan"]),],
                  AUC=allcut[which.max(allcut[,"AUC"]),])
      #,bestcut=bestcut)
    }
  }
  return(output)
}

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


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

head(dat)

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

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







##########################
# OUTCOME 1
###
a1 <- findcutnum(dat$hsTroponin_pgmL.log, cbind(dat$death5yr,
            dat$timetodeath5yr),"survival",50,50)
#4
b1 <- findcutnum(dat$HSP70_ngmL.log, cbind(dat$death5yr,
            dat$timetodeath5yr),"survival",100,100)
#2
d1 <- findcutnum(dat$FDP_ugmL.log, cbind(dat$death5yr,
             dat$timetodeath5yr),"survival",100,100)
#3
e1 <- findcutnum(dat$CRP_mgL.log, cbind(dat$death5yr,
              dat$timetodeath5yr),"survival",100,50)
# 3
c1m <- findcutnum(dat.m$suPAR_ngmL.log, cbind(dat.m$death5yr,
          dat.m$timetodeath5yr),"survival",50,100)
#4
c1f <- findcutnum(dat.f$suPAR_ngmL.log, cbind(dat.f$death5yr,
          dat.f$timetodeath5yr),"survival",50,100)
#3




a1.cuts <- findcut(dat$hsTroponin_pgmL.log, cbind(dat$death5yr,
          dat$timetodeath5yr), cutnum = 3, "survival",50,50)
b1.cuts <- findcut(dat$HSP70_ngmL.log, cbind(dat$death5yr,
          dat$timetodeath5yr), cutnum = 2,"survival",100,100)
d1.cuts <- findcut(dat$FDP_ugmL.log, cbind(dat$death5yr,
          dat$timetodeath5yr), cutnum = 3,"survival",100,100)
c1m.cuts <- findcut(dat.m$suPAR_ngmL.log, cbind(dat.m$death5yr,
                  dat.m$timetodeath5yr),cutnum=4,"survival",50,100)
c1f.cuts <- findcut(dat.f$suPAR_ngmL.log, cbind(dat.f$death5yr,
                   dat.f$timetodeath5yr),cutnum=3,"survival",50,100)
e1.cuts <- findcut(dat$CRP_mgL.log, cbind(dat$death5yr,
                   dat$timetodeath5yr),cutnum=3,"survival",100,50)
exp(0.6931472)
exp(0.8329091)
exp(1.44096)
quantile(dat$FDP_ugmL)
##########################




##########################
# OUTCOME 2
###
a2 <- findcutnum(dat$hsTroponin_pgmL.log, cbind(dat$CVorMI,
              dat$timetoMI5yr),"survival",100,100)
#4
b2 <- findcutnum(dat$HSP70_ngmL.log, cbind(dat$CVorMI,
             dat$timetoMI5yr),"survival",100,100)
#2
d2 <- findcutnum(dat$FDP_ugmL.log, cbind(dat$CVorMI,
              dat$timetoMI5yr),"survival",100,100)
#3
e2 <- findcutnum(dat$CRP_mgL.log, cbind(dat$CVorMI,
                 dat$timetoMI5yr),"survival",100,50)
# 3
c2m <- findcutnum(dat.m$suPAR_ngmL.log, cbind(dat.m$CVorMI,
                 dat.m$timetoMI5yr),"survival",50,100)
#4
c2f <- findcutnum(dat.f$suPAR_ngmL.log, cbind(dat.f$CVorMI,
                 dat.f$timetoMI5yr),"survival",50,100)
#3


#############################
a2.cuts <- findcut(dat$hsTroponin_pgmL.log, cbind(dat$CVorMI,
                                                dat$timetoMI5yr),cutnum=3, "survival",100,100)
b2.cuts <- findcut(dat$HSP70_ngmL.log, cbind(dat$CVorMI,
                                           dat$timetoMI5yr),cutnum = 2, "survival",100,100)

d2.cuts <- findcut(dat$FDP_ugmL.log, cbind(dat$CVorMI,
                                         dat$timetoMI5yr),cutnum = 3, "survival",100,100)

e2.cuts <- findcut(dat$CRP_mgL.log, cbind(dat$CVorMI,
                                        dat$timetoMI5yr),cutnum = 3, "survival",100,50)


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


##########################
# OUTCOME 3
###
a3 <- findcutnum(dat$hsTroponin_pgmL.log, cbind(dat$DeathorMI,
                                                dat$timetoMI5yr),"survival",100,100)
#4
b3 <- findcutnum(dat$HSP70_ngmL.log, cbind(dat$DeathorMI,
                                           dat$timetoMI5yr),"survival",75,100)
#2
d3 <- findcutnum(dat$FDP_ugmL.log, cbind(dat$DeathorMI,
                                         dat$timetoMI5yr),"survival",100,100)
#3
e3 <- findcutnum(dat$CRP_mgL.log, cbind(dat$DeathorMI,
                                        dat$timetoMI5yr),"survival",100,50)
# 3
c3m <- findcutnum(dat.m$suPAR_ngmL.log, cbind(dat.m$DeathorMI,
                                              dat.m$timetoMI5yr),"survival",50,100)
#4
c3f <- findcutnum(dat.f$suPAR_ngmL.log, cbind(dat.f$DeathorMI,
                                              dat.f$timetoMI5yr),"survival",50,100)
#3
c3t <- findcutnum(dat$suPAR_ngmL.log, cbind(dat$DeathorMI,
                                              dat$timetoMI5yr),"survival",100,100)

####################################
###
a3cuts <- findcut(dat$hsTroponin_pgmL.log, cbind(dat$DeathorMI,
                                                dat$timetoMI5yr),cutnum= 3,"survival",100,100)
#4
b3cuts <- findcut(dat$HSP70_ngmL.log, cbind(dat$DeathorMI,
                                           dat$timetoMI5yr),cutnum = 2,"survival",75,100)
#2
d3cuts <- findcut(dat$FDP_ugmL.log, cbind(dat$DeathorMI,
                                         dat$timetoMI5yr),cutnum= 3,"survival",100,100)
#3
e3cuts <- findcut(dat$CRP_mgL.log, cbind(dat$DeathorMI,
                                        dat$timetoMI5yr),cutnum=3,"survival",100,50)

c3m <- findcut(dat.m$suPAR_ngmL.log, cbind(dat.m$DeathorMI,
                                              dat.m$timetoMI5yr),cutnum=3,"survival",50,100)
c3f <- findcut(dat.f$suPAR_ngmL.log, cbind(dat.f$DeathorMI,
                                              dat.f$timetoMI5yr),cutnum=3,"survival",50,100)
#3



exp(c(0.5933, 0.7129498, 0.7975072, 1.159098))
exp(c(0.8372475, 1.071948, 1.669522))



dat.f$C3 <- NA
for (i in 1:length(dat.f$uniqueid)){
  if (dat.f$suPAR_ngmL.log[i] <= 0.7870739){
    dat.f$C3[i] = 1
  }
  else if (dat.f$suPAR_ngmL.log[i] <= 1.031232){
    dat.f$C3[i] = 2
  }
  else if (dat.f$suPAR_ngmL.log[i] <= 1.770012){
    dat.f$C3[i] = 3
  }
  else if(dat.f$suPAR_ngmL.log[i] > 1.770012){
    dat.f$C3[i] = 4
  }
}

exp(c(0.5306283,0.6471032,1.202025))


hist(dat.f$C3)






exp(c(0.78707, 1.031232, 1.7700))

hist(dat.f$suPAR_ngmL.log, breaks=20)
abline(v=exp(c(0.78707, 1.031232, 1.7700)), col=3)

hist(dat.f$suPAR_ngmL.log, breaks=20)
abline(v=exp(c(1.074377,1.455215)), col=3)


dat.f$suPAR_ngmL.log
quantile(dat.f$supa)
########################################
#######################################
######################################


# OUTCOME 3

par(mfrow=c(1,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
  }
}

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

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



testing$B3 <- NA
for (i in 1:length(testing$uniqueid)){
  if (testing$HSP70_ngmL[i] <= 1.0){
    testing$B3[i] = 1
  }
  else if (testing$HSP70_ngmL[i] <= 37.4){
    testing$B3[i] = 2
  }
  else if(dat$hsTroponin_pgmL[i] <= 144.9){
    testing$B3[i] = 3
  }
  else if(dat$hsTroponin_pgmL[i] > 144.9){
    testing$B3[i] = 4
  }
}


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", "Very High"), col=c(1:3), lty=1, bty="n")



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
  }
}


KM_D3 <- survfit(SurvObj3~D3, data=testing, conf.type="log")
plot(KM_D3, col=c(1:4), main=("CV 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")



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
  }
}


KM_E3 <- survfit(SurvObj3~E3, data=testing, conf.type="log")
plot(KM_E3, col=c(1:4), main=("CV 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")



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




testing.f$C3 <- NA
for (i in 1:length(testing.f$uniqueid)){
  if (testing.f$FDP_ugmL[i] <= 1.7){
    testing.f$C3[i] = 1
  }
  else if (testing.f$FDP_ugmL[i] <= 1.9){
    testing.f$C3[i] = 2
  }
  else if (testing.f$FDP_ugmL[i] <= 3.3){
    testing.f$C3[i] = 3
  }
  else if(testing.f$FDP_ugmL[i] > 3.3){
    testing.f$C3[i] = 4
  }
}



KM_C3 <- survfit(SurvObj3~C3, data=testing.f, conf.type="log")
plot(KM_C3, col=c(1:4), main=("CV 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")



testing.m$C3 <- NA
for (i in 1:length(testing.m$uniqueid)){
  if (testing.m$CRP_mgL[i] <= 2.2){
    testing.m$C3[i] = 1
  }
  else if (testing.m$CRP_mgL[i] <= 2.8){
    testing.m$C3[i] = 2
  }
  else if (testing.m$CRP_mgL[i] <= 5.9){
    testing.m$C3[i] = 3
  }
  else if(testing.m$CRP_mgL[i] > 5.9){
    testing.m$C3[i] = 4
  }
}


KM_C3 <- survfit(SurvObj3~C3, data=testing.m, conf.type="log")
plot(KM_C3, col=c(1:4), main=("CV 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")



#####################################
###################################
# OUTCOME 1

par(mfrow=c(2,2))

dat$A1 <- NA
for (i in 1:length(dat$uniqueid)){
  if (dat$hsTroponin_pgmL[i] <= 2.0){
    dat$A1[i] = 1
  }
  else if (dat$hsTroponin_pgmL[i] <= 2.3){
    dat$A1[i] = 2
  }
  else if (dat$hsTroponin_pgmL[i] <= 4.2){
    dat$A1[i] = 3
  }
  else if(dat$hsTroponin_pgmL[i] > 4.2){
    dat$A1[i] = 4
  }
}

KM_A1 <- survfit(SurvObj1~A1, data=dat, conf.type="log")
plot(KM_A1, col=c(1:4), main=("All-Cause Death by HS Trop"), ylim=c(0.5, 1))
legend("bottomleft", c("Low","Medium","High","Very High"), col=c(1:4), lty=1, bty="n")



dat$B1 <- NA
for (i in 1:length(dat$uniqueid)){
  if (dat$HSP70_ngmL[i] <= 3.1){
    dat$B1[i] = 1
  }
  else if (dat$HSP70_ngmL[i] <= 151.1){
    dat$B1[i] = 2
  }
  else if(dat$hsTroponin_pgmL[i] > 151.1){
    dat$B1[i] = 3
  }
}


KM_B1 <- survfit(SurvObj1~B1, data=dat, conf.type="log")
plot(KM_B1, col=c(1:3), main=("All-Cause Death by HSP-70"), ylim=c(0.5, 1))
legend("bottomleft", c("Low","Medium","High"), col=c(1:3), lty=1, bty="n")



dat$C1 <- NA
for (i in 1:length(dat$uniqueid)){
  if (dat$FDP_ugmL[i] <= 1.17){
    dat$C1[i] = 1
  }
  else if (dat$FDP_ugmL[i] <= 1.22){
    dat$C1[i] = 2
  }
  else if (dat$FDP_ugmL[i] <= 1.51){
    dat$C1[i] = 3
  }
  else if(dat$FDP_ugmL[i] > 1.51){
    dat$C1[i] = 4
  }
}


KM_C1 <- survfit(SurvObj1~C1, data=dat, conf.type="log")
plot(KM_C1, col=c(1:4), main=("All-Cause Death by FDP"), ylim=c(0.5, 1))
legend("bottomleft", c("Low","Medium","High","Very High"), col=c(1:4), lty=1, bty="n")



dat$E1 <- NA
for (i in 1:length(dat$uniqueid)){
  if (dat$CRP_mgL[i] <= 1.4){
    dat$E1[i] = 1
  }
  else if (dat$CRP_mgL[i] <= 2.8){
    dat$E1[i] = 2
  }
  else if (dat$CRP_mgL[i] <= 14.1){
    dat$E1[i] = 3
  }
  else if(dat$CRP_mgL[i] > 14.1){
    dat$E1[i] = 4
  }
}


KM_E1 <- survfit(SurvObj1~E1, data=dat, conf.type="log")
plot(KM_E1, col=c(1:4), main=("All-Cause Death by CRP"), ylim=c(0.5, 1))
legend("bottomleft", c("Low","Medium","High","Very High"), col=c(1:4), lty=1, bty="n")


######################################
#######################################
########################################
# OUTCOME 2

par(mfrow=c(2,2))

dat$A2 <- NA
for (i in 1:length(dat$uniqueid)){
  if (dat$hsTroponin_pgmL[i] <= 2.3){
    dat$A2[i] = 1
  }
  else if (dat$hsTroponin_pgmL[i] <= 2.7){
    dat$A2[i] = 2
  }
  else if (dat$hsTroponin_pgmL[i] <= 4.3){
    dat$A2[i] = 3
  }
  else if(dat$hsTroponin_pgmL[i] > 4.3){
    dat$A2[i] = 4
  }
}

KM_A2 <- survfit(SurvObj2~A2, data=dat, conf.type="log")
plot(KM_A2, col=c(1:4), main=("CV Death or MI by HS Trop"), ylim=c(0.5, 1))
legend("bottomleft", c("Low","Medium","High","Very High"), col=c(1:4), lty=1, bty="n")



dat$B2 <- NA
for (i in 1:length(dat$uniqueid)){
  if (dat$HSP70_ngmL[i] <= 0){
    dat$B2[i] = 1
  }
  else if (dat$HSP70_ngmL[i] <= 384){
    dat$B2[i] = 2
  }
  else if(dat$hsTroponin_pgmL[i] > 384){
    dat$B2[i] = 3
  }
}


KM_B2 <- survfit(SurvObj2~B2, data=dat, conf.type="log")
plot(KM_B2, col=c(1:3), main=("CV 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")



dat$C2 <- NA
for (i in 1:length(dat$uniqueid)){
  if (dat$FDP_ugmL[i] <= 1.17){
    dat$C2[i] = 1
  }
  else if (dat$FDP_ugmL[i] <= 1.54){
    dat$C2[i] = 2
  }
  else if (dat$FDP_ugmL[i] <= 1.88){
    dat$C2[i] = 3
  }
  else if(dat$FDP_ugmL[i] > 1.88){
    dat$C2[i] = 4
  }
}


KM_C2 <- survfit(SurvObj2~C2, data=dat, conf.type="log")
plot(KM_C2, col=c(1:4), main=("CV 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")



dat$E2 <- NA
for (i in 1:length(dat$uniqueid)){
  if (dat$CRP_mgL[i] <= 2.46){
    dat$E2[i] = 1
  }
  else if (dat$CRP_mgL[i] <= 5.91){
    dat$E2[i] = 2
  }
  else if (dat$CRP_mgL[i] <= 20.04){
    dat$E2[i] = 3
  }
  else if(dat$CRP_mgL[i] > 20.04){
    dat$E2[i] = 4
  }
}


KM_E2 <- survfit(SurvObj2~E2, data=dat, conf.type="log")
plot(KM_E2, col=c(1:4), main=("CV 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")


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























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


system.time(findcut(dat$hsTroponin_pgmL.rank, cbind(dat$death5yr,
         dat$timetodeath5yr),cutnum=1,"survival",100,100))

system.time(findcut(dat$hsTroponin_pgmL.log, cbind(dat$death5yr,
      dat$timetodeath5yr),cutnum=1,"survival",100,100))

system.time(findcutnum(dat$hsTroponin_pgmL.log, cbind(dat$death5yr,
          dat$timetodeath5yr),"survival",100,100))

system.time(findcutnum(dat$hsTroponin_pgmL.log, cbind(dat$death5yr,
               dat$timetodeath5yr),"survival",100,500))

system.time(findcut(dat$hsTroponin_pgmL.log, cbind(dat$death5yr,
             dat$timetodeath5yr),cutnum = 4,"survival",100,100))

system.time(findcut(dat$hsTroponin_pgmL.log, cbind(dat$death5yr,
             dat$timetodeath5yr),cutnum = 3,"survival",100,100))



##################
##################
# Chop up data in order to do these functions
quantile(dat$hsTroponin_pgmL, c(0,.95))
quantile(dat$HSP70_ngmL, c(0.05,0.95))
quantile(dat$FDP_ugmL, c(0.05,0.95))
quantile(dat$suPAR_ngmL, c(0.05,0.95))
quantile(dat$CRP_mgL, c(0.05,0.95))
dat.90.hsTrop <- dat[which(dat$hsTroponin_pgmL <= 67.7 & dat$hsTroponin_pgmL >= 1.5),]
dat.90.hsp70 <- dat[which(dat$HSP70_ngmL <= 365 & dat$HSP70_ngmL >= 0),]
dat.90.FDP <- dat[which(dat$FDP_ugmL <= 5.01575 & dat$FDP_ugmL >= 0.196),]
dat.90.suPAR <- dat[which(dat$suPAR_ngmL <= 6.7025 & dat$suPAR_ngmL >= 1.74),]
dat.90.CRP <- dat[which(dat$CRP_mgL <= 23 & dat$CRP_mgL >= 0.4),]


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

a1 <- findcut(dat$hsTroponin_pgmL, cbind(dat$death5yr,
            dat$timetodeath5yr),cutnum=1,"survival",100,1000)
a1


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

a <- findcut(dat.90.hsTrop$hsTroponin_pgmL, cbind(dat.90.hsTrop$death5yr,
        dat.90.hsTrop$timetodeath5yr),cutnum=1,"survival",100,1000)
a
# 3.7522
b <- findcut(dat.90.hsp70$HSP70_ngmL, cbind(dat.90.hsp70$death5yr,
        dat.90.hsp70$timetodeath5yr),cutnum=1,"survival",100,1000)
b$wilcoxon
b$allcut[,1]
# 6.154
c <- findcut(dat.90.FDP$FDP_ugmL, cbind(dat.90.FDP$death5yr,
       dat.90.FDP$timetodeath5yr),cutnum=1,"survival",100,1000)
c$wilcoxon
c$allcut[,1]
# 0.4944
d <- findcut(dat.90.suPAR$suPAR_ngmL, cbind(dat.90.suPAR$death5yr,
       dat.90.suPAR$timetodeath5yr),cutnum=1,"survival",100,1000)
d$wilcoxon
d$allcut[,1]
# 2.70004
e <- findcut(dat.90.CRP$CRP_mgL, cbind(dat.90.CRP$death5yr,
       dat.90.CRP$timetodeath5yr),cutnum=1,"survival",100,1000)
e$wilcoxon
e$allcut[,1]
# 3.9126

#####################################
# Outcome 2
head(dat.90.hsTrop)
a <- findcut(dat.90.hsTrop$hsTroponin_pgmL, cbind(dat.90.hsTrop$CVorMI,dat.90.hsTrop$timetoMI5yr),
        cutnum=1,"survival",100,1000)
a$wilcoxon
# 4.017
b <- findcut(dat.90.hsp70$HSP70_ngmL, cbind(dat.90.hsp70$CVorMI,
      dat.90.hsp70$timetoMI5yr),cutnum=1,"survival",100,1000)
b$wilcoxon
b$allcut[,1]
# 1.086
c <- findcut(dat.90.FDP$FDP_ugmL, cbind(dat.90.FDP$CVorMI,
      dat.90.FDP$timetoMI5yr),cutnum=1,"survival",100,1000)
c$wilcoxon
c$allcut[,1]
# 0.8403
d <- findcut(dat.90.suPAR$suPAR_ngmL, cbind(dat.90.suPAR$CVorMI,
      dat.90.suPAR$timetoMI5yr),cutnum=1,"survival",100,1000)
d$wilcoxon
d$allcut[,1]
# 2.8433
e <- findcut(dat.90.CRP$CRP_mgL, cbind(dat.90.CRP$CVorMI,
      dat.90.CRP$timetoMI),cutnum=1,"survival",100,1000)
e$wilcoxon
e$allcut[,1]
# 5.2008

#####################################
# Outcome 3
head(dat.90.hsTrop)
a <- findcut(dat.90.hsTrop$hsTroponin_pgmL, cbind(dat.90.hsTrop$DeathorMI,dat.90.hsTrop$timetoMI5yr),
        cutnum=1,"survival",100,1000)
a$wilcoxon
# 3.6198
b <- findcut(dat.90.hsp70$HSP70_ngmL, cbind(dat.90.hsp70$DeathorMI,
       dat.90.hsp70$timetoMI5yr),cutnum=1,"survival",100,1000)
b$wilcoxon
b$allcut[,1]
# 1.086
c <- findcut(dat.90.FDP$FDP_ugmL, cbind(dat.90.FDP$DeathorMI,
        dat.90.FDP$timetoMI5yr),cutnum=1,"survival",100,1000)
c$wilcoxon
c$allcut[,1]
# 0.494416
d <- findcut(dat.90.suPAR$suPAR_ngmL, cbind(dat.90.suPAR$DeathorMI,
         dat.90.suPAR$timetoMI5yr),cutnum=1,"survival",100,1000)
d$wilcoxon
d$allcut[,1]
# 2.66052
e <- findcut(dat.90.CRP$CRP_mgL, cbind(dat.90.CRP$DeathorMI,
         dat.90.CRP$timetoMI),cutnum=1,"survival",100,1000)
e$wilcoxon
e$allcut[,1]
# 4.116