#colAUC {caTools}
library(caTools)
library(ROCR)

setwd("I://Spinale//Diagnostic Code")

LVH4 <- read.table("LVH_4markers.csv", sep = ",", header=T)


CV <- function(data,K){
	N <- nrow(data)
	M <- ncol(data)
	
	modulus <- N%%K 

	true <-predict<-list()
	id <- sample(N,N,replace=FALSE)
	cut <- floor(N/K)

	data.test <- array(0,dim=c(cut,M,(K-1)))
	data.train <- array(0,dim=c((N-cut),M,(K-1)))

	data <- as.matrix(data)
	
	# FIRST K-1 CUTS OF THE DATA
	for(k in 1:(K-1)) {
	  data.test[,,k] <- data[id[((k-1)*cut+1):(cut*k)],]
	  data.train[,,k] <- data[-id[((k-1)*cut+1):(cut*k)],]	
 	  colnames(data.test[,,k]) <- colnames(data.train[,,k]) <- colnames(data)
	}
	
	#Kth PORTION OF THE DATA
	test.K <- data[id[(cut*(K-1)+1):N],] 
	train.K <- data[-id[(cut*(K-1)+1):N],]
	
	dimnames(data.test)=list(a=1:cut,b=colnames(data),c=1:(K-1))
	dimnames(data.train)=list(a=1:(N-cut),b=colnames(data),c=1:(K-1))

	beta <- pred <- truth <- AUC <- 0
	x.val <- y.val <- list()
	color <- c("red","green","black","blue","yellow")

	for(k in 1:(K-1)) {
	Xtrain <- data.train[,,k]
	Xtest <- data.test[,,k]
	logistic.fit <- glm(Xtrain[,1] ~ Xtrain[,2:M], family = binomial(logit))

	x.matrix <- cbind(Xtest[,2:M])
	beta <- logistic.fit$coefficients
	pred <- antilogit(beta[1] + x.matrix %*% beta[2:M])
	truth <- Xtest[,1]

	# Computes a simple ROC curve from predicted values
	predict <- prediction(pred, truth)
	perf <- performance(predict, measure = "tpr", x.measure = "fpr")
	AUC[k] <- performance(predict,'auc')@y.values[[1]]
	x.val[[k]] <- perf@x.values
	y.val[[k]] <- perf@y.values
	}

	logistic.fit <- glm(train.K[,1] ~ train.K[,2:M], family = binomial(logit))
	x.matrix <- cbind(test.K[,2:M])
	beta <- logistic.fit$coefficients
	pred <- antilogit(beta[1] + x.matrix %*% beta[2:M])
	truth <- test.K[,1]
	
	predict <- prediction(pred, truth)
	perf <- performance(predict, measure = "tpr", x.measure = "fpr")
	AUC[K] <- performance(predict,'auc')@y.values[[1]]
	x.val[[K]] <- perf@x.values
	y.val[[K]] <- perf@y.values
	
	plot(x.val[[1]][[1]],y.val[[1]][[1]],type="l",col=color[1],xlab=
		"1-Specificity", ylab="Sensitivity",cex=1.7,cex.axis=1.7)
	for(k in 2:K){
	points(x.val[[k]][[1]],y.val[[k]][[1]],type="l",col=color[k])
	}

	return(c(AUC=mean(AUC), quant1 = quantile(AUC,.05), quant2 = quantile(AUC,.95)))
}
antilogit <- function(u){ 
		exp(u)/(1+exp(u))
}


par(mfrow=c(3,3))

Nsims <- 50
AUC.LVH <- AUC.LVH.BNP <- AUC.CHF <- AUC.CHF.BNP<- matrix(c(0),Nsims,3)
#par(mfrow=c(2,2))#CHANGE FOR PANEL EFFECT FOR SIMULATION
for(i in 1:Nsims){
	
	AUC.LVH[i,] <- CV(LVH4,5)
	title("LVH 4 Biomarker Subset")
}

CV(LVH4,5)

#FUNCTION CALL
#CV(data,K)