TADPOLE_Time_to_AD.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}
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]]


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]

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;
j=1
for (j in 1:length(VISIT_TRAINCROSSMRI) )
{
  tmpss <-  subset(VISIT_TRAINCROSSMRI[[j]],MCI_AD==1 & BDX==2)
  YearMCItoAD <- c(YearMCItoAD,as.numeric(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")

MCIADUniRankFeaturesRaw["YearsFromEvent",]

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

CASESNLtoAD <- NULL
YearNLtoAD <- NULL
for (j in 1:length(VISIT_TRAINCROSSMRI) )
{
  tmpss <-  subset(VISIT_TRAINCROSSMRI[[j]],NL_AD==1 & BDX == 1)
  YearNLtoAD <- c(YearNLtoAD,(tmpss$YearNLtoAD))
  CASESNLtoAD <- rbind(CASESNLtoAD,tmpss[,sampledcolumns])
}

CasesTimetoAD <- cbind(YearNLtoAD,CASESNLtoAD)
colnames(CasesTimetoAD) <- c("theTimeToAD",colnames(CasesTimetoAD)[-1])
colnames(CASESMCItoAD) <- colnames(CasesTimetoAD)
CasesTimetoAD <- rbind(CasesTimetoAD,CASESMCItoAD)
CasesTimetoAD <- CasesTimetoAD[complete.cases(CasesTimetoAD),]
#CasesTimetoAD$theTimeToAD <- log(CasesTimetoAD$theTimeToAD/2)
CasesTimetoAD$theTimeToAD <- sqrt(CasesTimetoAD$theTimeToAD)


ToADTIme <- univariateRankVariables(variableList = TADPOLE_predictors,
	                                            formula = "theTimeToAD ~ 1",
	                                            Outcome = "theTimeToAD",
	                                            data = CasesTimetoAD, 
	                                            categorizationType = "Raw", 
	                                            type = "LM", 
	                                            rankingTest = "Ztest",
	                                            description = "Description",
                                               uniType="Regression")
```


```{r Modeling time to AD}

numberofRepp=20
numberOfTrials=5
allformulas <- character();
logitcor <- numeric();
svmcor <- numeric();
svmMAD <- numeric();
logitMAD <- numeric();
tsamples <- numberofRepp*numberOfTrials

testWithMCI$theTimeToAD <- numeric(nrow(testWithMCI))

for ( repp in 1:numberofRepp)
{

  ids <- unique(as.character(CasesTimetoAD$RID))
  sampleIDS_TimeToAD <- ids[sample(length(ids),9*length(ids)/10)]
  save(sampleIDS_TimeToAD,file=paste("sampleIDS_TimeToAD",numberOfTrials,repp,".RDATA",sep="_"))
  
  trainTimeToAD <- subset(CasesTimetoAD,RID %in% sampleIDS_TimeToAD)
  testTimeToAD <- subset(CasesTimetoAD,!(RID %in% sampleIDS_TimeToAD))
  
  print(length(unique(as.character(trainTimeToAD$RID))))
  print(length(unique(as.character(testTimeToAD$RID))))
  
  ids <- unique(as.character(trainTimeToAD$RID))
  sys_times <- NULL;
  TimeToAD <- list();
  n=1
  for (n in 1:numberOfTrials)
  {
    singlecaseTimeToAD <- NULL;
    for (i in ids)
    {
      case1 <- subset(trainTimeToAD,RID==i)
      caserows <- nrow(case1)
      if (caserows>1)
      {
        singlecaseTimeToAD <- rbind(singlecaseTimeToAD,case1[sample(caserows, 1),])
      }
      else
      {
        singlecaseTimeToAD <- rbind(singlecaseTimeToAD,case1)
      }
    }
    singlecaseTimeToAD$RID <- NULL;
    singlecaseTimeToAD$BDX <- NULL;
    print(nrow(singlecaseTimeToAD))
    sys_times <- rbind(sys_times,system.time(TimeToAD[[n]] <- FRESA.Model(theTimeToAD ~ 1,singlecaseTimeToAD)))
  }
  save(sys_times,file=paste("sys_times_TimeToAD3",numberOfTrials,repp,".RDATA",sep="_"))
  save(TimeToAD,file=paste("TimeToAD3",numberOfTrials,repp,".RDATA",sep="_"))
  
  #load(file=paste("TimeToAD",numberOfTrials,repp,".RDATA",sep="_"))
  #load(file=paste("sampleIDS_TimeToAD",numberOfTrials,repp,".RDATA",sep="_"))
  #trainTimeToAD <- subset(CasesTimetoAD,RID %in% sampleIDS_TimeToAD)
  #testTimeToAD <- subset(CasesTimetoAD,!(RID %in% sampleIDS_TimeToAD))
  #trainTimeToAD <- trainTimeToAD[complete.cases(trainTimeToAD),]
  #trainTimeToAD$theTimeToAD <- log(trainTimeToAD$theTimeToAD+1)
  #testTimeToAD <- testTimeToAD[complete.cases(testTimeToAD),]
  #testTimeToAD$theTimeToAD <- log(testTimeToAD$theTimeToAD+1)
  
  baggedformula <- character();
  baggedformula2 <- character();
  for (n in 1:numberOfTrials)
  {
    if (length(TimeToAD[[n]]$bagging$bagged.model$coefficients)>1)
    {
      baggedformula <- append(baggedformula,TimeToAD[[n]]$bagging$formula)
    }
    if (length(TimeToAD[[n]]$BSWiMS.model$coefficients)>1)
    {
      baggedformula2 <- append(baggedformula2,TimeToAD[[n]]$BSWiMS.model$formula)
    }
  }
  
  mp <- medianPredict(as.character(baggedformula2),trainTimeToAD,testTimeToAD,type="SVM")
  plot(testTimeToAD$theTimeToAD~mp$medianPredict)
  ct <- cor.test(testTimeToAD$theTimeToAD,mp$medianPredict)
  svmcor <- append(svmcor,ct$estimate)
  res <- testTimeToAD$theTimeToAD^2-mp$medianPredict^2
  res <- res[testTimeToAD$theTimeToAD<1]
  svmMAD <- append(svmMAD,mean(abs(res)))

  mp <- medianPredict(as.character(baggedformula2),trainTimeToAD,testTimeToAD,predictType ="linear",type="LM")
  plot(testTimeToAD$theTimeToAD~mp$medianPredict)
  ct <- cor.test(testTimeToAD$theTimeToAD,mp$medianPredict)
  logitcor <- append(logitcor,ct$estimate)
  res <- testTimeToAD$theTimeToAD^2-mp$medianPredict^2
  res <- res[testTimeToAD$theTimeToAD<1]
  logitMAD <- append(logitMAD,mean(abs(res)))
  print(svmMAD)
  print(logitMAD)

  mp <- medianPredict(as.character(baggedformula2),CasesTimetoAD,testWithMCI,type="SVM")
  as.character(baggedformula2)

#  mp$predictions[,-1] <- 2*exp(mp$predictions[,-1])
  mp$predictions[,-1] <- mp$predictions[,-1]^2
  cdft <- NULL
  x <- (1:60)/12
  for (i in 1:nrow(mp$predictions))
  {
    nd <- mp$predictions[i,-1]
    if (sna <- sum(is.na(nd))>0)
    {
      if (sna<numberOfTrials/2)
      {
        nd <- numeric(numberOfTrials) 
      }
      else
      {
        nd <- nd[!is.na(nd)];
      }
    }
    fn <- ecdf(as.numeric(nd))
    cdft <- rbind(cdft,fn(x))
  }
  rownames(cdft) <- rownames(mp$predictions)
  write.csv(cdft,file=paste(repp,"SVMTimeToADTimeCDF.csv",sep="_"))
  write.csv(mp$predictions,file=paste(repp,"SVMTimeToADTimePrediction.csv",sep="_"))
  
  mpl <- medianPredict(as.character(baggedformula2),CasesTimetoAD,testWithMCI,predictType ="linear",type="LM")
#  mpl$predictions[,-1] <- 2*exp(mpl$predictions[,-1])
  mpl$predictions[,-1] <- mpl$predictions[,-1]^2
  cdft <- NULL
  x <- (1:60)/12
  mpre <- mpl$predictions
  for (i in 1:nrow(mpl$predictions))
  {
    nd <- mpl$predictions[i,-1]
    if (sna <- sum(is.na(nd))>0)
    {
      if (sna<numberOfTrials/2)
      {
        nd <- numeric(numberOfTrials) 
      }
      else
      {
        nd <- nd[!is.na(nd)];
      }
    }
    fn <- ecdf(as.numeric(nd))
    cdft <- rbind(cdft,fn(x))
  }
  rownames(cdft) <- rownames(mpl$predictions)
  write.csv(cdft,file=paste(repp,"LMTimeToADTimeCDF.csv",sep="_"))
  write.csv(mpl$predictions,file=paste(repp,"LMTimeToADTimePrediction.csv",sep="_"))
  allformulas <- append(allformulas,as.character(baggedformula2))
}

mpl <- medianPredict(allformulas,CasesTimetoAD,testWithMCI,predictType ="linear",type="LM")
#mpl$predictions[,-1] <- 2*exp(mpl$predictions[,-1])
mpl$predictions[,-1] <- mpl$predictions[,-1]^2
cdft <- NULL
x <- (1:60)/12
sd(mpl$predictions[30,-1])
nelem <- nrow(mpl$predictions[,-1])*ncol(mpl$predictions[,-1])

mpl$predictions[,-1] <-  mpl$predictions[,-1]+as.data.frame(matrix(rnorm(nelem,0,mean(logitMAD)),nrow=nrow(mpl$predictions)));
sd(mpl$predictions[30,-1])

for (i in 1:nrow(mpl$predictions))
{
  nd <- mpl$predictions[i,-1]
  if (sna <- sum(is.na(nd))>0)
  {
    if (sna<numberOfTrials/2)
    {
      nd <- numeric(numberOfTrials) 
    }
    else
    {
      nd <- nd[!is.na(nd)];
    }
  }
  fn <- ecdf(as.numeric(nd))
  cdft <- rbind(cdft,fn(x))
}
rownames(cdft) <- rownames(mpl$predictions)
write.csv(cdft,file=paste(repp,"FinalADTimeCDF.csv",sep="_"))
write.csv(mpl$predictions,file=paste(repp,"FINALTimeToADTimePrediction.csv",sep="_"))
write.csv(allformulas,file=paste(repp,"AllFormulasFinalTimetoAD.csv",sep="_"))

MCIandNL <- rbind(testLastTimePointWithMCI,testLastTimePointWithNL)
MCIandNL$theTimeToAD <- numeric(nrow(MCIandNL))

mpl <- medianPredict(allformulas,CasesTimetoAD,MCIandNL,predictType ="linear",type="LM")
#mpl$predictions[,-1] <- 2*exp(mpl$predictions[,-1])
mpl$predictions[,-1] <- mpl$predictions[,-1]^2
cdft <- NULL
x <- (1:60)/12
sd(mpl$predictions[10,-1])
nelem <- nrow(mpl$predictions[,-1])*ncol(mpl$predictions[,-1])

mpl$predictions[,-1] <-  mpl$predictions[,-1]+as.data.frame(matrix(rnorm(nelem,0,mean(logitMAD)),nrow=nrow(mpl$predictions)));
sd(mpl$predictions[10,-1])
sum(is.na(mpl$predictions[,-1]))
i=nrow(mpl$predictions)
for (i in 1:nrow(mpl$predictions))
{
  nd <- mpl$predictions[i,-1]
  if (sna <- sum(is.na(nd))>0)
  {
    if (sna>tsamples/2)
    {
      nd <- numeric(tsamples) 
    }
    else
    {
      nd <- nd[!is.na(nd)];
    }
  }
  fn <- ecdf(as.numeric(nd))
  cdft <- rbind(cdft,fn(x))
}
rownames(cdft) <- rownames(mpl$predictions)
write.csv(cdft,file=paste(repp,"FinalLastObservationADTimeCDF.csv",sep="_"))
write.csv(mpl$predictions,file=paste(repp,"FINALLastObservationTimeToADTimePrediction.csv",sep="_"))



mpl <- medianPredict(allformulas,CasesTimetoAD,testWithMCI,type="SVM")
#mpl$predictions[,-1] <- 2*exp(mpl$predictions[,-1])
mpl$predictions[,-1] <- mpl$predictions[,-1]^2
cdft <- NULL
x <- (1:60)/12


mpl$predictions[,-1] <-  mpl$predictions[,-1]+as.data.frame(matrix(rnorm(nelem,0,mean(svmMAD)),nrow=nrow(mpl$predictions)));

for (i in 1:nrow(mpl$predictions))
{
  nd <- mpl$predictions[i,-1]
  if (sna <- sum(is.na(nd))>0)
  {
    if (sna<numberOfTrials/2)
    {
      nd <- numeric(numberOfTrials) 
    }
    else
    {
      nd <- nd[!is.na(nd)];
    }
  }
  fn <- ecdf(as.numeric(nd))
  cdft <- rbind(cdft,fn(x))
}
rownames(cdft) <- rownames(mpl$predictions)
write.csv(cdft,file=paste(repp,"FinalADTimeCDF_SVM.csv",sep="_"))
write.csv(mpl$predictions,file=paste(repp,"FINALTimeToADTimePrediction_SVM.csv",sep="_"))


barplot(cbind(svmcor,logitcor),beside=TRUE)
barplot(cbind(svmMAD,logitMAD),beside=TRUE)


```


```{r Modeling short time to AD}

numberofRepp=10
numberOfTrials=5
allformulas <- character();
logitcor <- numeric();
svmcor <- numeric();

ShortCasesTimeToAD <- subset(CasesTimetoAD,theTimeToAD<1.5)

testWithMCI$theTimeToAD <- numeric(nrow(testWithMCI))
repp=1
for ( repp in 1:numberofRepp)
{

  ids <- unique(as.character(ShortCasesTimeToAD$RID))
  sampleIDS_ShortTimeToAD <- ids[sample(length(ids),7*length(ids)/8)]
  save(sampleIDS_ShortTimeToAD,file=paste("sampleIDS_ShortTimeToAD",numberOfTrials,repp,".RDATA",sep="_"))
  
  trainShortTimeToAD <- subset(ShortCasesTimeToAD,RID %in% sampleIDS_ShortTimeToAD)
  testShortTimeToAD <- subset(ShortCasesTimeToAD,!(RID %in% sampleIDS_ShortTimeToAD))
  
  print(length(unique(as.character(trainShortTimeToAD$RID))))
  print(length(unique(as.character(testShortTimeToAD$RID))))
  
  ids <- unique(as.character(trainShortTimeToAD$RID))
  sys_times <- NULL;
  ShortTimeToAD <- list();
  n=1
  for (n in 1:numberOfTrials)
  {
    singlecaseShortTimeToAD <- NULL;
    for (i in ids)
    {
      case1 <- subset(trainShortTimeToAD,RID==i)
      caserows <- nrow(case1)
      if (caserows>1)
      {
        singlecaseShortTimeToAD <- rbind(singlecaseShortTimeToAD,case1[sample(caserows, 1),])
      }
      else
      {
        singlecaseShortTimeToAD <- rbind(singlecaseShortTimeToAD,case1)
      }
    }
    singlecaseShortTimeToAD$RID <- NULL;
    singlecaseShortTimeToAD$BDX <- NULL;
    print(nrow(singlecaseShortTimeToAD))
    sys_times <- rbind(sys_times,system.time(ShortTimeToAD[[n]] <- FRESA.Model(theTimeToAD ~ 1,singlecaseShortTimeToAD)))
  }
  save(sys_times,file=paste("sys_times_ShortTimeToAD3",numberOfTrials,repp,".RDATA",sep="_"))
  save(ShortTimeToAD,file=paste("ShortTimeToAD",numberOfTrials,repp,".RDATA",sep="_"))
  
  #load(file=paste("ShortTimeToAD",numberOfTrials,repp,".RDATA",sep="_"))
  #load(file=paste("sampleIDS_ShortTimeToAD",numberOfTrials,repp,".RDATA",sep="_"))
  #trainShortTimeToAD <- subset(ShortCasesTimeToAD,RID %in% sampleIDS_ShortTimeToAD)
  #testShortTimeToAD <- subset(ShortCasesTimeToAD,!(RID %in% sampleIDS_ShortTimeToAD))
  #trainShortTimeToAD <- trainShortTimeToAD[complete.cases(trainShortTimeToAD),]
  #trainShortTimeToAD$theTimeToAD <- log(trainShortTimeToAD$theTimeToAD+1)
  #testShortTimeToAD <- testShortTimeToAD[complete.cases(testShortTimeToAD),]
  #testShortTimeToAD$theTimeToAD <- log(testShortTimeToAD$theTimeToAD+1)
  
  baggedformula <- character();
  baggedformula2 <- character();
  for (n in 1:numberOfTrials)
  {
    if (length(ShortTimeToAD[[n]]$bagging$bagged.model$coefficients)>1)
    {
      baggedformula <- append(baggedformula,ShortTimeToAD[[n]]$bagging$formula)
    }
    if (length(ShortTimeToAD[[n]]$BSWiMS.model$coefficients)>1)
    {
      baggedformula2 <- append(baggedformula2,ShortTimeToAD[[n]]$BSWiMS.model$formula)
    }
  }
  
  mp <- medianPredict(as.character(baggedformula2),trainShortTimeToAD,testShortTimeToAD,type="SVM")
  plot(testShortTimeToAD$theTimeToAD~mp$medianPredict)
  ct <- cor.test(testShortTimeToAD$theTimeToAD,mp$medianPredict)
  svmcor <- append(svmcor,ct$estimate)
  
  mp <- medianPredict(as.character(baggedformula2),trainShortTimeToAD,testShortTimeToAD,predictType ="linear",type="LM")
  plot(testShortTimeToAD$theTimeToAD~mp$medianPredict)
  ct <- cor.test(testShortTimeToAD$theTimeToAD,mp$medianPredict)
  logitcor <- append(logitcor,ct$estimate)

  mp <- medianPredict(as.character(baggedformula2),ShortCasesTimeToAD,testWithMCI,type="SVM")
#  mp$predictions[,-1] <- 2*exp(mp$predictions[,-1])
  mp$predictions[,-1] <- mp$predictions[,-1]^2
  cdft <- NULL
  x <- (1:60)/12
  for (i in 1:nrow(mp$predictions))
  {
    nd <- mp$predictions[i,-1]
    if (sna <- sum(is.na(nd))>0)
    {
      if (sna<numberOfTrials/2)
      {
        nd <- numeric(numberOfTrials) 
      }
      else
      {
        nd <- nd[!is.na(nd)];
      }
    }
    fn <- ecdf(as.numeric(nd))
    cdft <- rbind(cdft,fn(x))
  }
  rownames(cdft) <- rownames(mp$predictions)
  write.csv(cdft,file=paste(repp,"SVMShortTimeToADTimeCDF.csv",sep="_"))
  write.csv(mp$predictions,file=paste(repp,"SVMShortTimeToADTimePrediction.csv",sep="_"))
  
  mpl <- medianPredict(as.character(baggedformula2),ShortCasesTimeToAD,testWithMCI,predictType ="linear",type="LM")
#  mpl$predictions[,-1] <- 2*exp(mpl$predictions[,-1])
  mpl$predictions[,-1] <- mpl$predictions[,-1]^2
  cdft <- NULL
  x <- (1:60)/12
  mpre <- mpl$predictions
  for (i in 1:nrow(mpl$predictions))
  {
    nd <- mpl$predictions[i,-1]
    if (sna <- sum(is.na(nd))>0)
    {
      if (sna<numberOfTrials/2)
      {
        nd <- numeric(numberOfTrials) 
      }
      else
      {
        nd <- nd[!is.na(nd)];
      }
    }
    fn <- ecdf(as.numeric(nd))
    cdft <- rbind(cdft,fn(x))
  }
  rownames(cdft) <- rownames(mpl$predictions)
  write.csv(cdft,file=paste(repp,"LMShortTimeToADTimeCDF.csv",sep="_"))
  write.csv(mpl$predictions,file=paste(repp,"LMShortTimeToADTimePrediction.csv",sep="_"))
  allformulas <- append(allformulas,as.character(baggedformula2))
}

mpl <- medianPredict(allformulas,ShortCasesTimeToAD,testWithMCI,predictType ="linear",type="LM")
#mpl$predictions[,-1] <- 2*exp(mpl$predictions[,-1])
mpl$predictions[,-1] <- mpl$predictions[,-1]^2
cdft <- NULL
x <- (1:60)/12
for (i in 1:nrow(mpl$predictions))
{
  nd <- mpl$predictions[i,-1]
  if (sna <- sum(is.na(nd))>0)
  {
    if (sna<numberOfTrials/2)
    {
      nd <- numeric(numberOfTrials) 
    }
    else
    {
      nd <- nd[!is.na(nd)];
    }
  }
  fn <- ecdf(as.numeric(nd))
  cdft <- rbind(cdft,fn(x))
}
rownames(cdft) <- rownames(mpl$predictions)
write.csv(cdft,file=paste(repp,"FinalShortADTimeCDF.csv",sep="_"))
write.csv(mpl$predictions,file=paste(repp,"FINALShortTimeToADTimePrediction.csv",sep="_"))

barplot(svmcor)
barplot(logitcor)

```

```{r Modeling long time to AD}

numberofRepp=10
numberOfTrials=5
allformulas <- character();
logitcor <- numeric();
svmcor <- numeric();

longCasesTimeToAD <- subset(CasesTimetoAD,theTimeToAD>=1.5)

testWithMCI$theTimeToAD <- numeric(nrow(testWithMCI))
repp=1
for ( repp in 1:numberofRepp)
{

  ids <- unique(as.character(longCasesTimeToAD$RID))
  sampleIDS_longTimeToAD <- ids[sample(length(ids),7*length(ids)/8)]
  save(sampleIDS_longTimeToAD,file=paste("sampleIDS_longTimeToAD",numberOfTrials,repp,".RDATA",sep="_"))
  
  trainlongTimeToAD <- subset(longCasesTimeToAD,RID %in% sampleIDS_longTimeToAD)
  testlongTimeToAD <- subset(longCasesTimeToAD,!(RID %in% sampleIDS_longTimeToAD))
  
  print(length(unique(as.character(trainlongTimeToAD$RID))))
  print(length(unique(as.character(testlongTimeToAD$RID))))
  
  ids <- unique(as.character(trainlongTimeToAD$RID))
  sys_times <- NULL;
  longTimeToAD <- list();
  n=1
  for (n in 1:numberOfTrials)
  {
    singlecaselongTimeToAD <- NULL;
    for (i in ids)
    {
      case1 <- subset(trainlongTimeToAD,RID==i)
      caserows <- nrow(case1)
      if (caserows>1)
      {
        singlecaselongTimeToAD <- rbind(singlecaselongTimeToAD,case1[sample(caserows, 1),])
      }
      else
      {
        singlecaselongTimeToAD <- rbind(singlecaselongTimeToAD,case1)
      }
    }
    singlecaselongTimeToAD$RID <- NULL;
    singlecaselongTimeToAD$BDX <- NULL;
    print(nrow(singlecaselongTimeToAD))
    sys_times <- rbind(sys_times,system.time(longTimeToAD[[n]] <- FRESA.Model(theTimeToAD ~ 1,singlecaselongTimeToAD)))
  }
  save(sys_times,file=paste("sys_times_longTimeToAD3",numberOfTrials,repp,".RDATA",sep="_"))
  save(longTimeToAD,file=paste("longTimeToAD",numberOfTrials,repp,".RDATA",sep="_"))
  
  #load(file=paste("longTimeToAD",numberOfTrials,repp,".RDATA",sep="_"))
  #load(file=paste("sampleIDS_longTimeToAD",numberOfTrials,repp,".RDATA",sep="_"))
  #trainlongTimeToAD <- subset(longCasesTimeToAD,RID %in% sampleIDS_longTimeToAD)
  #testlongTimeToAD <- subset(longCasesTimeToAD,!(RID %in% sampleIDS_longTimeToAD))
  #trainlongTimeToAD <- trainlongTimeToAD[complete.cases(trainlongTimeToAD),]
  #trainlongTimeToAD$theTimeToAD <- log(trainlongTimeToAD$theTimeToAD+1)
  #testlongTimeToAD <- testlongTimeToAD[complete.cases(testlongTimeToAD),]
  #testlongTimeToAD$theTimeToAD <- log(testlongTimeToAD$theTimeToAD+1)
  
  baggedformula <- character();
  baggedformula2 <- character();
  for (n in 1:numberOfTrials)
  {
    if (length(longTimeToAD[[n]]$bagging$bagged.model$coefficients)>1)
    {
      baggedformula <- append(baggedformula,longTimeToAD[[n]]$bagging$formula)
    }
    if (length(longTimeToAD[[n]]$BSWiMS.model$coefficients)>1)
    {
      baggedformula2 <- append(baggedformula2,longTimeToAD[[n]]$BSWiMS.model$formula)
    }
  }
  
  mp <- medianPredict(as.character(baggedformula2),trainlongTimeToAD,testlongTimeToAD,type="SVM")
  plot(testlongTimeToAD$theTimeToAD~mp$medianPredict)
  ct <- cor.test(testlongTimeToAD$theTimeToAD,mp$medianPredict)
  svmcor <- append(svmcor,ct$estimate)
  
  mp <- medianPredict(as.character(baggedformula2),trainlongTimeToAD,testlongTimeToAD,predictType ="linear",type="LM")
  plot(testlongTimeToAD$theTimeToAD~mp$medianPredict)
  ct <- cor.test(testlongTimeToAD$theTimeToAD,mp$medianPredict)
  logitcor <- append(logitcor,ct$estimate)

  mp <- medianPredict(as.character(baggedformula2),longCasesTimeToAD,testWithMCI,type="SVM")
#  mp$predictions[,-1] <- 2*exp(mp$predictions[,-1])
  mp$predictions[,-1] <- mp$predictions[,-1]^2
  cdft <- NULL
  x <- (1:60)/12
  for (i in 1:nrow(mp$predictions))
  {
    nd <- mp$predictions[i,-1]
    if (sna <- sum(is.na(nd))>0)
    {
      if (sna<numberOfTrials/2)
      {
        nd <- numeric(numberOfTrials) 
      }
      else
      {
        nd <- nd[!is.na(nd)];
      }
    }
    fn <- ecdf(as.numeric(nd))
    cdft <- rbind(cdft,fn(x))
  }
  rownames(cdft) <- rownames(mp$predictions)
  write.csv(cdft,file=paste(repp,"SVMlongTimeToADTimeCDF.csv",sep="_"))
  write.csv(mp$predictions,file=paste(repp,"SVMlongTimeToADTimePrediction.csv",sep="_"))
  
  mpl <- medianPredict(as.character(baggedformula2),longCasesTimeToAD,testWithMCI,predictType ="linear",type="LM")
#  mpl$predictions[,-1] <- 2*exp(mpl$predictions[,-1])
  mpl$predictions[,-1] <- mpl$predictions[,-1]^2
  cdft <- NULL
  x <- (1:60)/12
  mpre <- mpl$predictions
  for (i in 1:nrow(mpl$predictions))
  {
    nd <- mpl$predictions[i,-1]
    if (sna <- sum(is.na(nd))>0)
    {
      if (sna<numberOfTrials/2)
      {
        nd <- numeric(numberOfTrials) 
      }
      else
      {
        nd <- nd[!is.na(nd)];
      }
    }
    fn <- ecdf(as.numeric(nd))
    cdft <- rbind(cdft,fn(x))
  }
  rownames(cdft) <- rownames(mpl$predictions)
  write.csv(cdft,file=paste(repp,"LMlongTimeToADTimeCDF.csv",sep="_"))
  write.csv(mpl$predictions,file=paste(repp,"LMlongTimeToADTimePrediction.csv",sep="_"))
  allformulas <- append(allformulas,as.character(baggedformula2))
}

mpl <- medianPredict(allformulas,longCasesTimeToAD,testWithMCI,predictType ="linear",type="LM")
#mpl$predictions[,-1] <- 2*exp(mpl$predictions[,-1])
mpl$predictions[,-1] <- mpl$predictions[,-1]^2
cdft <- NULL
x <- (1:60)/12
for (i in 1:nrow(mpl$predictions))
{
  nd <- mpl$predictions[i,-1]
  if (sna <- sum(is.na(nd))>0)
  {
    if (sna<numberOfTrials/2)
    {
      nd <- numeric(numberOfTrials) 
    }
    else
    {
      nd <- nd[!is.na(nd)];
    }
  }
  fn <- ecdf(as.numeric(nd))
  cdft <- rbind(cdft,fn(x))
}
rownames(cdft) <- rownames(mpl$predictions)
write.csv(cdft,file=paste(repp,"FinalongADTimeCDF.csv",sep="_"))
write.csv(mpl$predictions,file=paste(repp,"FINALlongTimeToADTimePrediction.csv",sep="_"))

barplot(svmcor)
barplot(logitcor)

```