TADPOLE_MCI_to_AD_v2.Rmd

---
title: "TADPOLE"
author: "José Tamez-Peña"
date: "October 20, 2017"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
```{r, echo = FALSE}
knitr::opts_chunk$set(collapse = TRUE, warning = FALSE, message = FALSE,comment = "#>")
```

```{r Libraries, echo = FALSE}

library("epiR")
library("FRESA.CAD")
library(network)
library(GGally)
library("e1071")

TADPOLE_BASIC <- read.delim("AllSubjectsWithBasicInfo.txt")
TADPOLE_CROSSMRI <- read.delim("CompleMRIObs.txt",na.strings = c("NA","","#DIV/0!"))
TADPOLE_predictors <- read.delim("PredictorsList.txt")
rownames(TADPOLE_BASIC) <- as.character(TADPOLE_BASIC$IDTIME)
rownames(TADPOLE_CROSSMRI) <- as.character(TADPOLE_CROSSMRI$ID_TP)


```


```{r imputescoreData}
print(sum(is.na(TADPOLE_CROSSMRI)))
toimputeNA <- c("RID","APOE4","CDRSB","ADAS11","ADAS13","MMSE","RAVLT_immediate","RAVLT_learning","RAVLT_forgetting","RAVLT_perc_forgetting","FAQ")
basictoImpute <- TADPOLE_BASIC[,toimputeNA]
rownames(basictoImpute) <- TADPOLE_BASIC$IDTIME
ids <- unique(basictoImpute$RID)
print(sum(is.na(basictoImpute)))
for (id in ids)
{
  setid <- subset(basictoImpute,RID==id)
  setrows <- nrow(setid)
#  print(id)
  if (setrows>1)
  {
    for (var in toimputeNA[-1])
    {
      if (is.na(setid[1,var]))
      {
        setid[1,var] <- setid[2,var];
      }
      for (n in 2:setrows)
      {
        if (n<setrows)
        {
         if (is.na(setid[n,var]))
         {
           mv <- mean(c(setid[n-1,var],setid[n+1,var]),na.rm=TRUE)
           if (length(mv)>0)
           {
            setid[n,var] <- mv
           }
         }
        }
        else
        {
         if (is.na(setid[n,var]))
         {
            setid[n,var] <- setid[n-1,var]
         }
        }
      }
    }
  }
  basictoImpute[rownames(setid),] <- setid
#  print(setid)
}
print(sum(is.na(basictoImpute)))
rownames(TADPOLE_CROSSMRI) <- TADPOLE_CROSSMRI$ID_TP

TADPOLE_CROSSMRI[,toimputeNA[-1]] <- basictoImpute[rownames(TADPOLE_CROSSMRI),toimputeNA[-1]]

print(subset(TADPOLE_CROSSMRI,RID==5160)[,toimputeNA[-1]])



toimputeNA <- c("RID","Ventricle_ICV","ICV")
basictoImpute <- TADPOLE_CROSSMRI[,toimputeNA]
rownames(basictoImpute) <- TADPOLE_CROSSMRI$ID_TP
ids <- unique(basictoImpute$RID)
print(sum(is.na(basictoImpute)))
for (id in ids)
{
  setid <- subset(basictoImpute,RID==id)
  setrows <- nrow(setid)
#  print(id)
  if (setrows>1)
  {
    for (var in toimputeNA[-1])
    {
      if (is.na(setid[1,var]))
      {
        setid[1,var] <- setid[2,var];
      }
      for (n in 2:setrows)
      {
        if (n<setrows)
        {
         if (is.na(setid[n,var]))
         {
           mv <- mean(c(setid[n-1,var],setid[n+1,var]),na.rm=TRUE)
           if (length(mv)>0)
           {
            setid[n,var] <- mv
           }
         }
        }
        else
        {
         if (is.na(setid[n,var]))
         {
            setid[n,var] <- setid[n-1,var]
         }
        }
      }
    }
  }
  basictoImpute[rownames(setid),] <- setid
#  print(setid)
}
print(sum(is.na(basictoImpute)))
TADPOLE_CROSSMRI[,toimputeNA[-1]] <- basictoImpute[rownames(TADPOLE_CROSSMRI),toimputeNA[-1]]
print(sum(is.na(TADPOLE_CROSSMRI)))

```


```{r subsets}
months <- c(0,3,6,12,18,24,30,36,42,48,54,60,66,72,78,84,90,96,102,108,114,120)


TRAIN_TADPOLEMRI <- subset(TADPOLE_CROSSMRI,TADPOLE_CROSSMRI$D1==1)
TRAIN_TADPOLEBASIC <- subset(TADPOLE_BASIC,TADPOLE_BASIC$D1==1)
TEST_TADPOLEMRI <- subset(TADPOLE_CROSSMRI,TADPOLE_CROSSMRI$D2==1)
TEST_TADPOLEBASIC <- subset(TADPOLE_BASIC,TADPOLE_BASIC$D2==1)



```



```{r Adjusting for age and ICV}


TRAIN_TADPOLEMRI$cICV <- TRAIN_TADPOLEMRI$ICV^(1/3)
TEST_TADPOLEMRI$cICV <- TEST_TADPOLEMRI$ICV^(1/3)

rownames(TADPOLE_predictors) <- as.character(TADPOLE_predictors[,1])
prednames <- as.character(TADPOLE_predictors[-c(1:17),1])

TRAIN_TADPOLEMRI <- TRAIN_TADPOLEMRI[complete.cases(TRAIN_TADPOLEMRI[,prednames]),]
ControlNormal <- subset(TRAIN_TADPOLEMRI,BXBL=="CN" & LastDX2=="NL" & LastVisit>2 & Month==0)
sum(ControlNormal$PTGENDER==1)
sum(ControlNormal$PTGENDER==2)
hist(ControlNormal$MeanT)
hist(subset(ControlNormal,PTGENDER==1)$MeanT)
hist(subset(ControlNormal,PTGENDER==2)$MeanT)

trainTadploe.adj <- featureAdjustment(TADPOLE_predictors[prednames,], baseModel="1+AGE+cICV",data=TRAIN_TADPOLEMRI,referenceframe=ControlNormal,strata="PTGENDER", type = "LM", pvalue = 0.001)

testTadploe.adj <- featureAdjustment(TADPOLE_predictors[prednames,], baseModel="1+AGE+cICV",data=TEST_TADPOLEMRI,referenceframe=ControlNormal,strata="PTGENDER", type = "LM", pvalue = 0.001)

ControlNormal.adj <- featureAdjustment(TADPOLE_predictors[prednames,], baseModel="1+AGE+cICV",data=ControlNormal,referenceframe=ControlNormal,strata="PTGENDER", type = "LM", pvalue = 0.001)


hist(subset(ControlNormal.adj,PTGENDER==1)$MeanT,breaks = 15)
hist(subset(ControlNormal.adj,PTGENDER==2)$MeanT,breaks = 15)
hist(subset(ControlNormal.adj,PTGENDER==1)$MeanSAD,breaks = 15)
hist(subset(ControlNormal.adj,PTGENDER==2)$MeanSAD,breaks = 15)


#testTadploe.norm <- rankInverseNormalDataFrame(TADPOLE_predictors[prednames,], testTadploe.adj, ControlNormal.adj,strata="PTGENDER")

trainTadploe.norm <- rankInverseNormalDataFrame(TADPOLE_predictors[prednames,], trainTadploe.adj, ControlNormal.adj,strata="PTGENDER")
testTadploe.norm <- rankInverseNormalDataFrame(TADPOLE_predictors[prednames,], testTadploe.adj, ControlNormal.adj,strata="PTGENDER")

mean(testTadploe.norm$MeanT)
mean(trainTadploe.norm$MeanT)
mean(testTadploe.norm$MeanSAD)
mean(trainTadploe.norm$MeanSAD)

plot(testTadploe.norm$MeanT~testTadploe.adj$MeanT)
plot(trainTadploe.norm$MeanT~trainTadploe.adj$MeanT)
plot(testTadploe.norm$MeanSAD~testTadploe.adj$MeanSAD)
plot(trainTadploe.norm$MeanSAD~trainTadploe.adj$MeanSAD)

```


```{r spliting by visit}

VISIT_TRAINCROSSMRI <- list()
VISIT_TRAINBASIC <- list()
VISIT_TESTCROSSMRI <- list()
VISIT_TESTBASIC <- list()
i = 1;
for (j in months )
{
  VISIT_TRAINCROSSMRI[[i]] <- subset(trainTadploe.norm,Month==j)
  VISIT_TRAINBASIC[[i]] <- subset(TRAIN_TADPOLEBASIC,Month2==j)
  VISIT_TESTCROSSMRI[[i]] <- subset(testTadploe.norm,Month==j)
  VISIT_TESTBASIC[[i]] <- subset(TEST_TADPOLEBASIC,Month2==j)
  i = i + 1
}
sampledcolumns <- c("RID",as.character(TADPOLE_predictors[,1]));

testLastTimePointWithNL <- subset(testTadploe.norm,LastCompleteObs==1 & LastDX2=="NL")[,sampledcolumns]
testLastTimePointWithMCI <- subset(testTadploe.norm,LastCompleteObs==1 & LastDX2=="MCI")[,sampledcolumns]
testLastTimePointWithAD <- subset(testTadploe.norm,LastCompleteObs==1 & LastDX2=="Dementia")[,sampledcolumns]

testLastTimePointWithMCI["55_41030",]

testWithNL <- subset(testTadploe.norm,LastDX2=="NL")[,sampledcolumns]
testWithMCI <- subset(testTadploe.norm,LastDX2=="MCI")[,sampledcolumns]
testWithAD <- subset(testTadploe.norm,LastDX2=="Dementia")[,sampledcolumns]

#checkIDs <- c(as.character(testLastTimePointWithNL$RID),as.character(testLastTimePointWithMCI$RID),as.character(testLastTimePointWithAD$RID))
#write.csv(checkIDs,file="TestIDsWithCO.csv")


```

```{r univariate MCI to AD}
CASESMCItoAD <- NULL
YearMCItoAD <- NULL
controlYearMCItoAD <- NULL;
CONTROLMCItoAD <- NULL;
TestSetMCItoAD <- NULL;

for (j in 1:length(VISIT_TRAINCROSSMRI) )
{
  tmpss <-  subset(VISIT_TRAINCROSSMRI[[j]],MCI_AD==1 & BDX==2)
  YearMCItoAD <- c(YearMCItoAD,(tmpss$YearMCItoAD))
  CASESMCItoAD <- rbind(CASESMCItoAD,tmpss[,sampledcolumns])

  tmpss <- subset(VISIT_TRAINCROSSMRI[[j]],BXBL=="MCI" & LastDX2=="MCI" & LastVisit>4.0)
  controlYearMCItoAD <- c(controlYearMCItoAD,(tmpss$YearMCItoAD))
  CONTROLMCItoAD <- rbind(CONTROLMCItoAD,tmpss[,sampledcolumns])
}

CASESMCItoAD <- cbind(YearMCItoAD,CASESMCItoAD)
CASESMCItoAD <- CASESMCItoAD[complete.cases(CASESMCItoAD),]
CONTROLMCItoAD <- cbind(controlYearMCItoAD,CONTROLMCItoAD)
CONTROLMCItoAD <- CONTROLMCItoAD[complete.cases(CONTROLMCItoAD),]
colnames(CONTROLMCItoAD) <- colnames(CASESMCItoAD)

MCIADUniRankFeaturesRaw <- univariateRankVariables(variableList = TADPOLE_predictors,
	                                            formula = "YearMCItoAD ~ 1",
	                                            Outcome = "YearMCItoAD",
	                                            data = CASESMCItoAD, 
	                                            categorizationType = "Raw", 
	                                            type = "LM", 
	                                            rankingTest = "Ztest",
	                                            description = "Description",
                                               uniType="Regression")


mean(CASESMCItoAD$YearMCItoAD)
sum(CASESMCItoAD$YearMCItoAD<1)/nrow(CASESMCItoAD)
sum(CASESMCItoAD$YearMCItoAD<2)/nrow(CASESMCItoAD)



```

```{r predicting MCI to AD conversion}

numberofRepp = 20
numberOfTrials = 5
case_classMCItoAD <- CASESMCItoAD
case_classMCItoAD$Event <- 1
control_classMCItoAD <- CONTROLMCItoAD
control_classMCItoAD$Event <- 0
fulltrain <-rbind(case_classMCItoAD,control_classMCItoAD)
fulltrain <- fulltrain[complete.cases(fulltrain),] 
testLastTimePointWithMCI$Event <- rep(0,nrow(testLastTimePointWithMCI))
testWithMCI$Event <- rep(0,nrow(testWithMCI))
logitauc <- NULL;
svmauc <- NULL;
allformulas <- character();

for (repp in 1:numberofRepp)
{
  
  ids <- unique(as.character(case_classMCItoAD$RID))
  sampleCasesIDS_MCItoAD <- ids[sample(length(ids),4*length(ids)/5)]
  save(sampleCasesIDS_MCItoAD,file=paste("sampleCasesIDS_MCItoAD",numberOfTrials,repp,".RDATA",sep="_"));
  ids <- unique(as.character(control_classMCItoAD$RID))
  sampleControlsIDS_MCItoAD <- ids[sample(length(ids),4*length(ids)/5)]
  save(sampleControlsIDS_MCItoAD,file=paste("sampleControlsIDS_MCItoAD",numberOfTrials,repp,".RDATA",sep="_"));
  
  
  trainDataMCItoAD <- rbind(subset(case_classMCItoAD,RID %in% sampleCasesIDS_MCItoAD),
                            subset(control_classMCItoAD,RID %in% sampleControlsIDS_MCItoAD))
  testDataMCItoAD <- rbind(subset(case_classMCItoAD,!(RID %in% sampleCasesIDS_MCItoAD)),
                            subset(control_classMCItoAD,!(RID %in% sampleControlsIDS_MCItoAD)))
  
  trainDataMCItoAD$YearMCItoAD <- NULL
  trainDataMCItoAD$BDX <- NULL
  trainDataMCItoAD <- trainDataMCItoAD[complete.cases(trainDataMCItoAD),]
  testDataMCItoAD <- testDataMCItoAD[complete.cases(testDataMCItoAD),]
  
  
  ids <- unique(as.character(trainDataMCItoAD$RID))
  preMCItoAD  <- list();
  #load(file="preMCItoAD.RDATA")
  n=1
  sys_times <- NULL;
  for (n in 1:numberOfTrials)
  {
    singlecaseMCItoAD <- NULL;
    for (i in ids)
    {
      case1 <- subset(trainDataMCItoAD,RID==i)
      caserows <- nrow(case1)
      if (caserows>1)
      {
        singlecaseMCItoAD <- rbind(singlecaseMCItoAD,case1[sample(caserows, 1),])
      }
      else
      {
        singlecaseMCItoAD <- rbind(singlecaseMCItoAD,case1)
      }
    }
    singlecaseMCItoAD$RID <- NULL;
    singlecaseMCItoAD$BDX <- NULL;
    print(sum(singlecaseMCItoAD$Event))
    print(sum(singlecaseMCItoAD$Event==0))
    sys_times <- rbind(sys_times,system.time(preMCItoAD[[n]] <- FRESA.Model(Event ~ 1,singlecaseMCItoAD)))
  }
  save(sys_times,file=paste("sys_times_preMCItoAD",numberOfTrials,repp,".RDATA",sep="_"))
  save(preMCItoAD,file=paste("preMCItoAD",numberOfTrials,repp,".RDATA",sep="_"))
  
#  load(file=paste("preMCItoAD2",numberOfTrials,repp,".RDATA",sep="_"))
#  load(file=paste("sampleCasesIDS_MCItoAD2",numberOfTrials,repp,".RDATA",sep="_"))
#  load(file=paste("sampleControlsIDS_MCItoAD2",numberOfTrials,repp,".RDATA",sep="_"))
#  trainDataMCItoAD <- rbind(subset(case_classMCItoAD,RID %in% sampleCasesIDS_MCItoAD),
#                            subset(control_classMCItoAD,RID %in% sampleControlsIDS_MCItoAD))
#  testDataMCItoAD <- rbind(subset(case_classMCItoAD,!(RID %in% sampleCasesIDS_MCItoAD)),
#                            subset(control_classMCItoAD,!(RID %in% sampleControlsIDS_MCItoAD)))
#  trainDataMCItoAD <- trainDataMCItoAD[complete.cases(trainDataMCItoAD),]
#  testDataMCItoAD <- testDataMCItoAD[complete.cases(testDataMCItoAD),]
  
#  trainDataMCItoAD$YearMCItoAD <- NULL
#  trainDataMCItoAD$BDX <- NULL
  
  baggedformula <- character();
  baggedformula2 <- character();
  for (n in 1:numberOfTrials)
  {
    if (length(preMCItoAD[[n]]$bagging$bagged.model$coefficients)>1)
    {
      baggedformula <- append(baggedformula,preMCItoAD[[n]]$bagging$formula)
    }
    if (length(preMCItoAD[[n]]$BSWiMS.model$coefficients)>1)
    {
      baggedformula2 <- append(baggedformula2,preMCItoAD[[n]]$BSWiMS.model$formula)
    }
  }
  
  mp2 <- medianPredict(as.character(baggedformula2),trainDataMCItoAD,testDataMCItoAD,predictType ="prob",type="LOGIT")
  pm2 <- plotModels.ROC(mp2$predictions,main="Bagging",cex=0.90)
  logitauc <- append(logitauc,pm2$ensemble.auc)
  
  mpS2 <- medianPredict(as.character(baggedformula2),trainDataMCItoAD,testDataMCItoAD,type="SVM")
  mpS2$predictions[-1] <- mpS2$predictions[-1] - 0.5
  pmS2 <- plotModels.ROC(mpS2$predictions,main="SVM",cex=0.90)
  svmauc <- append(svmauc,pmS2$ensemble.auc)


  predictionMCIToAD <- medianPredict(as.character(baggedformula2),fulltrain,testLastTimePointWithMCI,predictType ="prob",type="SVM")
  write.csv(predictionMCIToAD$medianPredict,file=paste(repp,"ADFromMCIPrediction.csv",sep="_"))
  
  predictionAllMCIToAD <- medianPredict(as.character(baggedformula2),fulltrain,testWithMCI,predictType ="prob",type="SVM")
  
  write.csv(predictionAllMCIToAD$medianPredict,file=paste(repp,"ALLADFromMCIPrediction.csv",sep="_"))
  write.csv(predictionAllMCIToAD$predictions,file=paste(repp,"ALLPredADFromMCIPrediction.csv",sep="_"))

  predictionMCIToAD <- medianPredict(as.character(baggedformula2),fulltrain,testLastTimePointWithMCI,predictType ="prob",type="LOGIT")
  write.csv(predictionMCIToAD$medianPredict,file=paste(repp,"LOGITADFromMCIPrediction.csv",sep="_"))
  
  logpredictionAllMCIToAD <- medianPredict(as.character(baggedformula2),fulltrain,testWithMCI,predictType ="prob",type="LOGIT")
  
  write.csv(logpredictionAllMCIToAD$medianPredict,file=paste(repp,"LogitALLADFromMCIPrediction.csv",sep="_"))
  write.csv(logpredictionAllMCIToAD$predictions,file=paste(repp,"LogitALLPredADFromMCIPrediction.csv",sep="_"))
  allformulas <- append(allformulas,as.character(baggedformula2))
  
}
barplot(svmauc)
barplot(logitauc)
pluscts <- str_count(allformulas,"\\+")
write.csv(allformulas,file=paste(repp,"AllFormulasFinalMCItoADPred.csv",sep="_"))

logpredictionAllMCIToAD <- medianPredict(allformulas,fulltrain,testWithMCI,predictType ="prob",type="LOGIT")
  
write.csv(logpredictionAllMCIToAD$medianPredict,file=paste(repp,"FINAL_LogitALLADFromMCIPrediction.csv",sep="_"))
write.csv(logpredictionAllMCIToAD$predictions,file=paste(repp,"FINAL_LogitALLPredADFromMCIPrediction.csv",sep="_"))

predictionMCIToAD <- medianPredict(allformulas,fulltrain,testLastTimePointWithMCI,predictType ="prob",type="LOGIT")
write.csv(predictionMCIToAD$medianPredict,file=paste(repp,"FINAL_LOGITADFromMCIPrediction.csv",sep="_"))

logpredictionAllMCIToAD <- medianPredict(allformulas,fulltrain,testWithMCI,type="SVM")
  
write.csv(logpredictionAllMCIToAD$medianPredict,file=paste(repp,"FINAL_SVMALLADFromMCIPrediction.csv",sep="_"))
write.csv(logpredictionAllMCIToAD$predictions,file=paste(repp,"FINAL_SVMALLPredADFromMCIPrediction.csv",sep="_"))

predictionMCIToAD <- medianPredict(allformulas,fulltrain,testLastTimePointWithMCI,type="SVM")
write.csv(predictionMCIToAD$medianPredict,file=paste(repp,"FINAL_SVMADFromMCIPrediction.csv",sep="_"))

```


```{r predicting short term MCI to AD conversion}

numberofRepp = 10
numberOfTrials = 5
case_classShtMCItoAD <- subset(CASESMCItoAD,YearMCItoAD<1.5)
case_classShtMCItoAD$Event <- 1
control_classShtMCItoAD <- CONTROLMCItoAD
control_classShtMCItoAD$Event <- 0
fulltrain <-rbind(case_classShtMCItoAD,control_classShtMCItoAD)
fulltrain <- fulltrain[complete.cases(fulltrain),] 
testLastTimePointWithMCI$Event <- rep(0,nrow(testLastTimePointWithMCI))
testWithMCI$Event <- rep(0,nrow(testWithMCI))
logitauc <- NULL;
svmauc <- NULL;
allformulas <- character();
repp <- 1
for (repp in 1:numberofRepp)
{
  
  ids <- unique(as.character(case_classShtMCItoAD$RID))
  sampleCasesIDS_ShtMCItoAD <- ids[sample(length(ids),3*length(ids)/4)]
  save(sampleCasesIDS_ShtMCItoAD,file=paste("short_sampleCasesIDS_ShtMCItoAD",numberOfTrials,repp,".RDATA",sep="_"));
  ids <- unique(as.character(control_classShtMCItoAD$RID))
  sampleControlsIDS_ShtMCItoAD <- ids[sample(length(ids),3*length(ids)/4)]
  save(sampleControlsIDS_ShtMCItoAD,file=paste("short_sampleControlsIDS_ShtMCItoAD",numberOfTrials,repp,".RDATA",sep="_"));
  
  
  trainDataShtMCItoAD <- rbind(subset(case_classShtMCItoAD,RID %in% sampleCasesIDS_ShtMCItoAD),
                            subset(control_classShtMCItoAD,RID %in% sampleControlsIDS_ShtMCItoAD))
  testDataShtMCItoAD <- rbind(subset(case_classShtMCItoAD,!(RID %in% sampleCasesIDS_ShtMCItoAD)),
                            subset(control_classShtMCItoAD,!(RID %in% sampleControlsIDS_ShtMCItoAD)))
  
  trainDataShtMCItoAD$YearMCItoAD <- NULL
  trainDataShtMCItoAD$BDX <- NULL
  trainDataShtMCItoAD <- trainDataShtMCItoAD[complete.cases(trainDataShtMCItoAD),]
  testDataShtMCItoAD <- testDataShtMCItoAD[complete.cases(testDataShtMCItoAD),]
  
  
  ids <- unique(as.character(trainDataShtMCItoAD$RID))
  preShtMCItoAD  <- list();
  #load(file="preShtMCItoAD.RDATA")
  n=1
  sys_times <- NULL;
  for (n in 1:numberOfTrials)
  {
    singlecaseShtMCItoAD <- NULL;
    for (i in ids)
    {
      case1 <- subset(trainDataShtMCItoAD,RID==i)
      caserows <- nrow(case1)
      if (caserows>1)
      {
        singlecaseShtMCItoAD <- rbind(singlecaseShtMCItoAD,case1[sample(caserows, 1),])
      }
      else
      {
        singlecaseShtMCItoAD <- rbind(singlecaseShtMCItoAD,case1)
      }
    }
    singlecaseShtMCItoAD$RID <- NULL;
    singlecaseShtMCItoAD$BDX <- NULL;
    print(sum(singlecaseShtMCItoAD$Event))
    print(sum(singlecaseShtMCItoAD$Event==0))
    sys_times <- rbind(sys_times,system.time(preShtMCItoAD[[n]] <- FRESA.Model(Event ~ 1,singlecaseShtMCItoAD)))
  }
  save(sys_times,file=paste("short_sys_times_preShtMCItoAD",numberOfTrials,repp,".RDATA",sep="_"))
  save(preShtMCItoAD,file=paste("short_preShtMCItoAD",numberOfTrials,repp,".RDATA",sep="_"))
  
#  load(file=paste("preShtMCItoAD2",numberOfTrials,repp,".RDATA",sep="_"))
#  load(file=paste("sampleCasesIDS_ShtMCItoAD2",numberOfTrials,repp,".RDATA",sep="_"))
#  load(file=paste("sampleControlsIDS_ShtMCItoAD2",numberOfTrials,repp,".RDATA",sep="_"))
#  trainDataShtMCItoAD <- rbind(subset(case_classShtMCItoAD,RID %in% sampleCasesIDS_ShtMCItoAD),
#                            subset(control_classShtMCItoAD,RID %in% sampleControlsIDS_ShtMCItoAD))
#  testDataShtMCItoAD <- rbind(subset(case_classShtMCItoAD,!(RID %in% sampleCasesIDS_ShtMCItoAD)),
#                            subset(control_classShtMCItoAD,!(RID %in% sampleControlsIDS_ShtMCItoAD)))
#  trainDataShtMCItoAD <- trainDataShtMCItoAD[complete.cases(trainDataShtMCItoAD),]
#  testDataShtMCItoAD <- testDataShtMCItoAD[complete.cases(testDataShtMCItoAD),]
  
#  trainDataShtMCItoAD$YearMCItoAD <- NULL
#  trainDataShtMCItoAD$BDX <- NULL
  
  baggedformula <- character();
  baggedformula2 <- character();
  for (n in 1:numberOfTrials)
  {
    if (length(preShtMCItoAD[[n]]$bagging$bagged.model$coefficients)>1)
    {
      baggedformula <- append(baggedformula,preShtMCItoAD[[n]]$bagging$formula)
    }
    if (length(preShtMCItoAD[[n]]$BSWiMS.model$coefficients)>1)
    {
      baggedformula2 <- append(baggedformula2,preShtMCItoAD[[n]]$BSWiMS.model$formula)
    }
  }
  
  mp2 <- medianPredict(as.character(baggedformula2),trainDataShtMCItoAD,testDataShtMCItoAD,predictType ="prob",type="LOGIT")
  pm2 <- plotModels.ROC(mp2$predictions,main="Bagging",cex=0.90)
  logitauc <- append(logitauc,pm2$ensemble.auc)
  
  mpS2 <- medianPredict(as.character(baggedformula2),trainDataShtMCItoAD,testDataShtMCItoAD,type="SVM")
  mpS2$predictions[-1] <- mpS2$predictions[-1] - 0.5
  pmS2 <- plotModels.ROC(mpS2$predictions,main="SVM",cex=0.90)
  svmauc <- append(svmauc,pmS2$ensemble.auc)

  predictionShtMCItoAD <- medianPredict(as.character(baggedformula2),fulltrain,testLastTimePointWithMCI,predictType ="prob",type="SVM")
  write.csv(predictionShtMCItoAD$medianPredict,file=paste(repp,"short_ADFromMCIPrediction.csv",sep="_"))
  
  predictionAllShtMCItoAD <- medianPredict(as.character(baggedformula2),fulltrain,testWithMCI,predictType ="prob",type="SVM")
  
  write.csv(predictionAllShtMCItoAD$medianPredict,file=paste(repp,"short_ALLADFromMCIPrediction.csv",sep="_"))
  write.csv(predictionAllShtMCItoAD$predictions,file=paste(repp,"short_ALLPredADFromMCIPrediction.csv",sep="_"))

  predictionShtMCItoAD <- medianPredict(as.character(baggedformula2),fulltrain,testLastTimePointWithMCI,predictType ="prob",type="LOGIT")
  write.csv(predictionShtMCItoAD$medianPredict,file=paste(repp,"short_LOGITADFromMCIPrediction.csv",sep="_"))
  
  logpredictionAllShtMCItoAD <- medianPredict(as.character(baggedformula2),fulltrain,testWithMCI,predictType ="prob",type="LOGIT")
  
  write.csv(logpredictionAllShtMCItoAD$medianPredict,file=paste(repp,"short_LogitALLADFromMCIPrediction.csv",sep="_"))
  write.csv(logpredictionAllShtMCItoAD$predictions,file=paste(repp,"short_LogitALLPredADFromMCIPrediction.csv",sep="_"))
  allformulas <- append(allformulas,as.character(baggedformula2))
  
}
barplot(svmauc)
barplot(logitauc)

logpredictionAllShtMCItoAD <- medianPredict(allformulas,fulltrain,testWithMCI,predictType ="prob",type="LOGIT")
  
write.csv(logpredictionAllShtMCItoAD$medianPredict,file=paste(repp,"short_FINAL_LogitALLADFromMCIPrediction.csv",sep="_"))
write.csv(logpredictionAllShtMCItoAD$predictions,file=paste(repp,"short_FINAL_LogitALLPredADFromMCIPrediction.csv",sep="_"))

predictionShtMCItoAD <- medianPredict(allformulas,fulltrain,testLastTimePointWithMCI,predictType ="prob",type="LOGIT")
write.csv(predictionShtMCItoAD$medianPredict,file=paste(repp,"short_FINAL_LOGITADFromMCIPrediction.csv",sep="_"))


```


```{r predicting late term MCI to AD conversion}

numberofRepp = 10
numberOfTrials = 5
case_classLatMCItoAD <- subset(CASESMCItoAD,YearMCItoAD>=1.5)
case_classLatMCItoAD$Event <- 1
control_classLatMCItoAD <- CONTROLMCItoAD
control_classLatMCItoAD$Event <- 0
fulltrain <-rbind(case_classLatMCItoAD,control_classLatMCItoAD)
fulltrain <- fulltrain[complete.cases(fulltrain),] 
testLastTimePointWithMCI$Event <- rep(0,nrow(testLastTimePointWithMCI))
testWithMCI$Event <- rep(0,nrow(testWithMCI))
logitauc <- NULL;
svmauc <- NULL;
allformulas <- character();
repp <- 1
for (repp in 1:numberofRepp)
{
  
  ids <- unique(as.character(case_classLatMCItoAD$RID))
  sampleCasesIDS_LatMCItoAD <- ids[sample(length(ids),3*length(ids)/4)]
  save(sampleCasesIDS_LatMCItoAD,file=paste("late_sampleCasesIDS_LatMCItoAD",numberOfTrials,repp,".RDATA",sep="_"));
  ids <- unique(as.character(control_classLatMCItoAD$RID))
  sampleControlsIDS_LatMCItoAD <- ids[sample(length(ids),3*length(ids)/4)]
  save(sampleControlsIDS_LatMCItoAD,file=paste("late_sampleControlsIDS_LatMCItoAD",numberOfTrials,repp,".RDATA",sep="_"));
  
  
  trainDataLatMCItoAD <- rbind(subset(case_classLatMCItoAD,RID %in% sampleCasesIDS_LatMCItoAD),
                            subset(control_classLatMCItoAD,RID %in% sampleControlsIDS_LatMCItoAD))
  testDataLatMCItoAD <- rbind(subset(case_classLatMCItoAD,!(RID %in% sampleCasesIDS_LatMCItoAD)),
                            subset(control_classLatMCItoAD,!(RID %in% sampleControlsIDS_LatMCItoAD)))
  
  trainDataLatMCItoAD$YearMCItoAD <- NULL
  trainDataLatMCItoAD$BDX <- NULL
  trainDataLatMCItoAD <- trainDataLatMCItoAD[complete.cases(trainDataLatMCItoAD),]
  testDataLatMCItoAD <- testDataLatMCItoAD[complete.cases(testDataLatMCItoAD),]
  
  
  ids <- unique(as.character(trainDataLatMCItoAD$RID))
  preLatMCItoAD  <- list();
  #load(file="preLatMCItoAD.RDATA")
  n=1
  sys_times <- NULL;
  for (n in 1:numberOfTrials)
  {
    singlecaseLatMCItoAD <- NULL;
    for (i in ids)
    {
      case1 <- subset(trainDataLatMCItoAD,RID==i)
      caserows <- nrow(case1)
      if (caserows>1)
      {
        singlecaseLatMCItoAD <- rbind(singlecaseLatMCItoAD,case1[sample(caserows, 1),])
      }
      else
      {
        singlecaseLatMCItoAD <- rbind(singlecaseLatMCItoAD,case1)
      }
    }
    singlecaseLatMCItoAD$RID <- NULL;
    singlecaseLatMCItoAD$BDX <- NULL;
    print(sum(singlecaseLatMCItoAD$Event))
    print(sum(singlecaseLatMCItoAD$Event==0))
    sys_times <- rbind(sys_times,system.time(preLatMCItoAD[[n]] <- FRESA.Model(Event ~ 1,singlecaseLatMCItoAD)))
  }
  save(sys_times,file=paste("late_sys_times_preLatMCItoAD",numberOfTrials,repp,".RDATA",sep="_"))
  save(preLatMCItoAD,file=paste("late_preLatMCItoAD",numberOfTrials,repp,".RDATA",sep="_"))
  
#  load(file=paste("preLatMCItoAD2",numberOfTrials,repp,".RDATA",sep="_"))
#  load(file=paste("sampleCasesIDS_LatMCItoAD2",numberOfTrials,repp,".RDATA",sep="_"))
#  load(file=paste("sampleControlsIDS_LatMCItoAD2",numberOfTrials,repp,".RDATA",sep="_"))
#  trainDataLatMCItoAD <- rbind(subset(case_classLatMCItoAD,RID %in% sampleCasesIDS_LatMCItoAD),
#                            subset(control_classLatMCItoAD,RID %in% sampleControlsIDS_LatMCItoAD))
#  testDataLatMCItoAD <- rbind(subset(case_classLatMCItoAD,!(RID %in% sampleCasesIDS_LatMCItoAD)),
#                            subset(control_classLatMCItoAD,!(RID %in% sampleControlsIDS_LatMCItoAD)))
#  trainDataLatMCItoAD <- trainDataLatMCItoAD[complete.cases(trainDataLatMCItoAD),]
#  testDataLatMCItoAD <- testDataLatMCItoAD[complete.cases(testDataLatMCItoAD),]
  
#  trainDataLatMCItoAD$YearMCItoAD <- NULL
#  trainDataLatMCItoAD$BDX <- NULL
  
  baggedformula <- character();
  baggedformula2 <- character();
  for (n in 1:numberOfTrials)
  {
    if (length(preLatMCItoAD[[n]]$bagging$bagged.model$coefficients)>1)
    {
      baggedformula <- append(baggedformula,preLatMCItoAD[[n]]$bagging$formula)
    }
    if (length(preLatMCItoAD[[n]]$BSWiMS.model$coefficients)>1)
    {
      baggedformula2 <- append(baggedformula2,preLatMCItoAD[[n]]$BSWiMS.model$formula)
    }
  }
  
  mp2 <- medianPredict(as.character(baggedformula2),trainDataLatMCItoAD,testDataLatMCItoAD,predictType ="prob",type="LOGIT")
  pm2 <- plotModels.ROC(mp2$predictions,main="Bagging",cex=0.90)
  logitauc <- append(logitauc,pm2$ensemble.auc)
  
  mpS2 <- medianPredict(as.character(baggedformula2),trainDataLatMCItoAD,testDataLatMCItoAD,type="SVM")
  mpS2$predictions[-1] <- mpS2$predictions[-1] - 0.5
  pmS2 <- plotModels.ROC(mpS2$predictions,main="SVM",cex=0.90)
  svmauc <- append(svmauc,pmS2$ensemble.auc)

  predictionLatMCItoAD <- medianPredict(as.character(baggedformula2),fulltrain,testLastTimePointWithMCI,predictType ="prob",type="SVM")
  write.csv(predictionLatMCItoAD$medianPredict,file=paste(repp,"late_ADFromMCIPrediction.csv",sep="_"))
  
  predictionAllLatMCItoAD <- medianPredict(as.character(baggedformula2),fulltrain,testWithMCI,predictType ="prob",type="SVM")
  
  write.csv(predictionAllLatMCItoAD$medianPredict,file=paste(repp,"late_ALLADFromMCIPrediction.csv",sep="_"))
  write.csv(predictionAllLatMCItoAD$predictions,file=paste(repp,"late_ALLPredADFromMCIPrediction.csv",sep="_"))

  predictionLatMCItoAD <- medianPredict(as.character(baggedformula2),fulltrain,testLastTimePointWithMCI,predictType ="prob",type="LOGIT")
  write.csv(predictionLatMCItoAD$medianPredict,file=paste(repp,"late_LOGITADFromMCIPrediction.csv",sep="_"))
  
  logpredictionAllLatMCItoAD <- medianPredict(as.character(baggedformula2),fulltrain,testWithMCI,predictType ="prob",type="LOGIT")
  
  write.csv(logpredictionAllLatMCItoAD$medianPredict,file=paste(repp,"late_LogitALLADFromMCIPrediction.csv",sep="_"))
  write.csv(logpredictionAllLatMCItoAD$predictions,file=paste(repp,"late_LogitALLPredADFromMCIPrediction.csv",sep="_"))
  allformulas <- append(allformulas,as.character(baggedformula2))
  
}
barplot(svmauc)
barplot(logitauc)

logpredictionAllLatMCItoAD <- medianPredict(allformulas,fulltrain,testWithMCI,predictType ="prob",type="LOGIT")
  
write.csv(logpredictionAllLatMCItoAD$medianPredict,file=paste(repp,"late_FINAL_LogitALLADFromMCIPrediction.csv",sep="_"))
write.csv(logpredictionAllLatMCItoAD$predictions,file=paste(repp,"late_FINAL_LogitALLPredADFromMCIPrediction.csv",sep="_"))

predictionLatMCItoAD <- medianPredict(allformulas,fulltrain,testLastTimePointWithMCI,predictType ="prob",type="LOGIT")
write.csv(predictionLatMCItoAD$medianPredict,file=paste(repp,"late_FINAL_LOGITADFromMCIPrediction.csv",sep="_"))


```