ftvdat <- read.csv("6657_FTV_20170817_noPHI_NoSite.csv")
preddat <- read.csv("I-SPY 1 level2a Patient Clinical and Outcome Data_predictors.csv")
outcdat <- read.csv("I-SPY 1 level2a Patient Clinical and Outcome Data_outcomes.csv")


## Plotting parameter values
cexval <- 1.9
cex.axisval <- cexval
cex.labval <- cexval
pchval <- 19
youngCol <- "blue"
oldCol <- "red"
plotWidth <- 720
plotHeight <- 720
plotBG <- "transparent"
marvals <- c(5,5,2,2)
lwdval <- 3
##########################

names(ftvdat)[1] <- "SUBJECTID"

dat <- merge(ftvdat,preddat)
dat <- merge(dat,outcdat)

names(dat) <- c("SubjectID","DataExtractDt","examT0","examT1","examT2","examT3","daysT0_T1","daysT0_T2","daysT0_T3","pixelVolT0","pixelVolT1","pixelVolT2","pixelVolT3","pixelVolTfinal","pixelVolPctChgT0_T1","pixelVolPctChgT0_T2","pixelVolPctChgT0_T3","ftvPeT0","ftvPeT1","ftvPeT2","ftvPeT3","ftvPeTfinal","ftvPePctChgT0_T1","ftvPePctChgT0_T2","ftvPePctChgT0_T3","age","race","ERpos","PgRpos","HRpos","Her2MostPos","HR_HER2cat","HR_HER2status","bilatCa","laterality","MRI_LD_T0","MRI_LD_T1","MRI_LD_T2","MRI_LD_Tfinal","sstat","survDtD2tx.","RFS","rfsInd","pCR","RCBclass")

write.csv(dat,file="ispy1.csv")

logit <- function(x,b0,b1){
  explinpred <- exp(b0 + b1*x)
  explinpred/(1+explinpred)
}

glm0 <- glm(pCR ~ MRI_LD_T0, family = "binomial", data=dat, na.action="na.exclude")
summary(glm0)
confint(glm0)

library(jtools)
summ(glm0,confint=TRUE,digits=3)

glm0_probs <- predict(glm0,type="response") 

png(filename="Lecture2plots/plot1.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(dat$MRI_LD_T0,dat$pCR,xlab="LD T0",ylab="Pr(pCR)",pch=19,cex=0.4, cex.axis=cex.axisval, cex.lab=cex.labval)
curve(logit(x,b0=as.numeric(coef(glm0)[1]),b1=as.numeric(coef(glm0)[2])),add=TRUE,lwd=lwdval,col="red")
dev.off()
glm1 <- glm(pCR ~ MRI_LD_T1, family = "binomial", data=dat, na.action="na.exclude")
summ(glm1,confint=TRUE,digits=3)

png(filename="Lecture2plots/plot2.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(dat$MRI_LD_T1,dat$pCR,xlab="LD T1",ylab="Pr(pCR)",pch=19,cex=0.4, cex.axis=cex.axisval, cex.lab=cex.labval)
curve(logit(x,b0=as.numeric(coef(glm1)[1]),b1=as.numeric(coef(glm1)[2])),add=TRUE,lwd=lwdval,col="red")
dev.off()

glmfinal <- glm(pCR ~ MRI_LD_Tfinal, family = "binomial", data=dat, na.action="na.exclude")
summ(glmfinal,confint=TRUE,digits=3)

png(filename="Lecture2plots/plot3.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(dat$MRI_LD_Tfinal,dat$pCR,xlab="LD Tfinal",ylab="Pr(pCR)",pch=19,cex=0.4, cex.axis=cex.axisval, cex.lab=cex.labval)
curve(logit(x,b0=as.numeric(coef(glmfinal)[1]),b1=as.numeric(coef(glmfinal)[2])),add=TRUE,lwd=lwdval,col="red")
dev.off()


png(filename="Lecture2plots/plot4.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(dat$MRI_LD_Tfinal,dat$pCR,xlab="LD Tfinal",ylab="Pr(pCR)",pch=19,cex=0.4, cex.axis=cex.axisval, cex.lab=cex.labval)
curve(logit(x,b0=as.numeric(coef(glmfinal)[1]),b1=as.numeric(coef(glmfinal)[2])),add=TRUE,lwd=lwdval,col="red")
abline(h=0.5,col="blue",lwd=lwdval)
dev.off()

glmfinal_probs <- predict(glmfinal,type="response") 
max(glmfinal_probs)


glmMult <- glm(pCR ~ MRI_LD_Tfinal * ftvPePctChgT0_T2, family = "binomial", data=dat, na.action="na.exclude")
summ(glmMult,confint=TRUE,digits=3)

glmMult_probs <- predict(glmMult,type="response") 
max(glmMult_probs,na.rm=TRUE)

################################################################
######## DISCRIMINANT ANALYSIS #################################
################################################################

set.seed(1)
sampSize <- 10^7 ## Size of each sample
mean1 <- -1.0 ## Mean distribution 1
mean2 <- 1.0 ## Mean distribution 2
sd1 <- sd2 <- 1.0 ## SD for each of sample 1 and 2
xlim1 <- -6  ## limits on the x-axis for the plots
xlim2 <- 6

samp1 <- rnorm(n=sampSize,mean=mean1,sd=sd1)
samp2 <- rnorm(n=sampSize,mean=mean2,sd=sd2)

# Histogram Colored (blue and red)
png(filename="Lecture2plots/plot4.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
curve(dnorm(x,mean=mean1,sd=sd1),from=-6,to=6,col=rgb(1,0,0,1),lty=0,xlab="X",ylab="Pr(X=x)",cex.axis=cex.axisval, cex.lab=cex.labval)
curve(dnorm(x,mean=mean2,sd=sd2),add=TRUE,col=rgb(0,0,1,1),lty=0)
hist(samp1, col=rgb(1,0,0,0.5),xlim=c(-6,6),prob=TRUE,add=TRUE)
hist(samp2, col=rgb(0,0,1,0.5), prob=TRUE, add=TRUE)
legend(x="topright",legend=c("cancer","not cancer"),lwd=lwdval,col=c(rgb(1,0,0,1),rgb(0,0,1,1)))
box()
dev.off()

png(filename="Lecture2plots/plot5.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
curve(dnorm(x,mean=mean1,sd=sd1),from=-6,to=6,col=rgb(1,0,0,1),lwd=lwdval,xlab="X",ylab="Pr(X=x|Y=k)",cex.axis=cex.axisval, cex.lab=cex.labval)
curve(dnorm(x,mean=mean2,sd=sd2),add=TRUE,col=rgb(0,0,1,1),lwd=lwdval)
hist(samp1, col=rgb(1,0,0,0.5),xlim=c(-6,6),prob=TRUE,add=TRUE)
hist(samp2, col=rgb(0,0,1,0.5), prob=TRUE, add=TRUE)
legend(x="topright",legend=c("cancer","not cancer"),lwd=lwdval,col=c(rgb(1,0,0,1),rgb(0,0,1,1)))
box()
dev.off()

rm(samp1,samp2) ## Free up the memory

png(filename="Lecture2plots/plot6.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
curve(dnorm(x,mean=mean1,sd=sd1),from=-6,to=6,col=rgb(1,0,0,1),lwd=lwdval,xlab="X",ylab="Pr(X=x|Y=k)",cex.axis=cex.axisval, cex.lab=cex.labval)
curve(dnorm(x,mean=mean2,sd=sd2),add=TRUE,col=rgb(0,0,1,1),lwd=lwdval)
abline(v=0,lwd=lwdval,lty=2)
legend(x="topright",legend=c("cancer","not cancer"),lwd=lwdval,col=c(rgb(1,0,0,1),rgb(0,0,1,1)))
box()
dev.off()


png(filename="Lecture2plots/plot7.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
curve(1.2*dnorm(x,mean=mean2,sd=sd2),from=-6,to=6,col=rgb(0,0,1,1),lwd=lwdval,xlab="X",ylab="",cex.axis=cex.axisval, cex.lab=cex.labval)
curve(0.8*dnorm(x,mean=mean1,sd=sd1),add=TRUE,col=rgb(1,0,0,1),lwd=lwdval)
abline(v=-0.2,lwd=lwdval,lty=2)
legend(x="topright",legend=c("cancer","not cancer"),lwd=lwdval,col=c(rgb(1,0,0,1),rgb(0,0,1,1)))
box()
dev.off()

############## Same simulated example with smaller dataset

set.seed(5)
sampSize2 <- 25

samp1 <- rnorm(n=sampSize2,mean=mean1,sd=sd1)
samp2 <- rnorm(n=sampSize2,mean=mean2,sd=sd2)

png(filename="Lecture2plots/plot8.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
curve(dnorm(x,mean=mean1,sd=sd1),ylim=c(0,0.5),from=-6,to=6,col=rgb(1,0,0,1),lwd=lwdval,xlab="X",ylab="Pr(X=x|Y=k)",cex.axis=cex.axisval, cex.lab=cex.labval)
curve(dnorm(x,mean=mean2,sd=sd2),add=TRUE,col=rgb(0,0,1,1),lwd=lwdval)
legend(x="topright",legend=c("cancer","not cancer"),lwd=lwdval,col=c(rgb(1,0,0,1),rgb(0,0,1,1)))
box()
dev.off()

png(filename="Lecture2plots/plot9.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
curve(dnorm(x,mean=mean1,sd=sd1),from=-6,to=6,col=rgb(1,0,0,1),lty=0,xlab="X",ylab="Pr(X=x|Y=k)",ylim=c(0,0.5),cex.axis=cex.axisval, cex.lab=cex.labval)
hist(samp1, col=rgb(1,0,0,0.5),breaks=10,xlim=c(-6,6),prob=TRUE,add=TRUE)
hist(samp2, col=rgb(0,0,1,0.5),breaks=10,prob=TRUE, add=TRUE)
legend(x="topright",legend=c("cancer","not cancer"),lwd=lwdval,col=c(rgb(1,0,0,1),rgb(0,0,1,1)))
box()
dev.off()

datLDAtrain <- data.frame(cancer=c(rep("cancer",sampSize2),rep("benign",sampSize2)),x=c(samp1,samp2))

head(datLDAtrain)
tail(datLDAtrain)

library(MASS)
ldafit <- lda(cancer ~ x, data=datLDAtrain)
ldafit
plot(ldafit)

png(filename="Lecture2plots/plot10.png", width=870, height=950, bg="transparent")
par(mar=marvals) ## setting margin width
plot(ldafit,cex.axis=cex.axisval, cex.lab=cex.labval)
dev.off()

png(filename="Lecture2plots/plot11.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
curve(dnorm(x,mean=mean1,sd=sd1),ylim=c(0,0.5),from=-6,to=6,col=rgb(1,0,0,1),lwd=lwdval,xlab="X",ylab="Pr(X=x|Y=k)",cex.axis=cex.axisval, cex.lab=cex.labval)
curve(dnorm(x,mean=mean2,sd=sd2),add=TRUE,col=rgb(0,0,1,1),lwd=lwdval)
legend(x="topright",legend=c("cancer","not cancer"),lwd=lwdval,col=c(rgb(1,0,0,1),rgb(0,0,1,1)))
abline(v=-1,col=rgb(1,0,0,1),lwd=lwdval,lty=2)
abline(v=1,col=rgb(0,0,1,1),lwd=lwdval,lty=2)
abline(v=0,lwd=lwdval,lty=2)
box()
dev.off()

png(filename="Lecture2plots/plot12.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
curve(dnorm(x,mean=mean1,sd=sd1),from=-6,to=6,col=rgb(1,0,0,1),lty=0,xlab="X",ylab="Pr(X=x|Y=k)",ylim=c(0,0.5),cex.axis=cex.axisval, cex.lab=cex.labval)
hist(samp1, col=rgb(1,0,0,0.5),breaks=10,xlim=c(-6,6),prob=TRUE,add=TRUE)
hist(samp2, col=rgb(0,0,1,0.5),breaks=10,prob=TRUE, add=TRUE)
abline(v=-1,col=rgb(1,0,0,1),lwd=lwdval,lty=2)
abline(v=1,col=rgb(0,0,1,1),lwd=lwdval,lty=2)
abline(v=0,lwd=lwdval,lty=2)
abline(v=ldafit$means[2],col=rgb(1,0,0,1),lwd=lwdval)
abline(v=ldafit$means[1],col=rgb(0,0,1,1),lwd=lwdval)
abline(v=mean(ldafit$means),lwd=lwdval)
legend(x="topright",legend=c("cancer","not cancer"),lwd=lwdval,col=c(rgb(1,0,0,1),rgb(0,0,1,1)))
box()
dev.off()

ldapred <- predict(ldafit)
names(ldapred)
head(ldapred$class)
head(ldapred$posterior)
head(ldapred$x)
ldaclass <- ldapred$class
table(ldaclass,datLDAtrain$cancer)
mean(ldaclass==datLDAtrain$cancer)
sum(ldapred$posterior[,1]>=.5)
sum(ldapred$posterior[,1]<.5)
ldapred$posterior[1:20,1]
ldaclass[1:20]
sum(ldapred$posterior[,1]>.9)


set.seed(6)
samp1b <- rnorm(n=sampSize2,mean=mean1,sd=sd1)
samp2b <- rnorm(n=sampSize2,mean=mean2,sd=sd2)

datLDAtest <- data.frame(cancer=c(rep("cancer",sampSize2),rep("benign",sampSize2)),x=c(samp1b,samp2b))

head(datLDAtest)
tail(datLDAtest)

ldapredTest <- predict(ldafit,newdata=datLDAtest)
names(ldapredTest)
head(ldapredTest$class)
head(ldapredTest$posterior)
head(ldapredTest$x)
ldaclassTest <- ldapredTest$class
table(ldaclassTest,datLDAtest$cancer)
mean(ldaclassTest==datLDAtest$cancer)
sum(ldapredTest$posterior[,1]>=.5)
sum(ldapredTest$posterior[,1]<.5)
ldapredTest$posterior[1:20,1]
ldaclassTest[1:20]
sum(ldapredTest$posterior[,1]>.9)

library(MASS)

ldafit2 <- lda(cancer ~ x, prior=c(0.8,0.2), data=datLDAtrain)
ldafit2

ldapred2 <- predict(ldafit2)
names(ldapred2)
head(ldapred2$class)
head(ldapred2$posterior)
head(ldapred2$x)
ldaclass2 <- ldapred2$class
table(ldaclass2,datLDAtrain$cancer)
mean(ldaclass2==datLDAtrain$cancer)
sum(ldapred2$posterior[,1]>=.5)
sum(ldapred2$posterior[,1]<.5)
ldapred2$posterior[1:20,1]
ldaclass2[1:20]
sum(ldapred2$posterior[,1]>.9)

############## Example on breast data with 2 predictors

## First normalize the predictors
dat$LDfz <- (dat$MRI_LD_Tfinal-mean(dat$MRI_LD_Tfinal,na.rm=TRUE))/sd(dat$MRI_LD_Tfinal,na.rm=TRUE)

dat$ftvChgT0_T2 <- (dat$ftvPePctChgT0_T2-mean(dat$ftvPePctChgT0_T2,na.rm=TRUE))/sd(dat$ftvPePctChgT0_T2,na.rm=TRUE)
                    
ldaMult <- lda(pCR ~ MRI_LD_Tfinal + ftvPePctChgT0_T2, data=dat)

ldaMult

png(filename="Lecture2plots/plot13.png", width=870, height=950, bg="transparent")
par(mar=marvals) ## setting margin width
plot(ldaMult,cex.axis=cex.axisval, cex.lab=cex.labval)
dev.off()