datamatrixModeling_binary.RMD

---
title: "Data Matrix Modeling"
output: pdf_document
params:
  csvPath:            datalocation/pyradiomicsout.csv
  target:             MutationalStatusp53
  positive_class:     mut
  inputs:            !r NULL 
  leaveOneOut:       !r TRUE
  test_csv:          !r NULL
  rescale:           !r FALSE
  removeCorrelated:  !r TRUE
  plot:              !r TRUE
  semisupervised:    !r FALSE
  kClusters:         !r as.numeric(9)
  genetic:           !r FALSE
  boruta:            !r FALSE
  univariate:        !r TRUE
  unipValThresh:     !r 0.05
---
## Version: BINARY CLASSIFIER

Update: Default behavior is now to take all non-target features as inpute

DOC: Usage - Rscript -e "rmarkdown::render( 'datamatrixModeling_binary.RMD', output_file = 'myfile.pdf',params = list( csvPath='datalocation/pyradiomics_modeling.csv', target='MutationalStatus', positive_class='mut', inputs=NULL , leaveOneOut=TRUE, test_csv='datalocation/pyradiomicsout_experimentalcohort.csv', rescale=FALSE, removeCorrelated=TRUE, plot=TRUE, semisupervised=FALSE, kClusters=as.numeric(9), genetic=FALSE, boruta=TRUE, univariate=TRUE, unipValThresh=0.05 ))"

DOC: Debugging - set params in R shell:  params = list( csvPath='datalocation/pyradiomics_modeling.csv', target='MutationalStatus', positive_class='mut', inputs=NULL , leaveOneOut=TRUE, test_csv='datalocation/pyradiomicsout_experimentalcohort.csv', rescale=FALSE, removeCorrelated=TRUE, plot=TRUE, semisupervised=FALSE, kClusters=as.numeric(9), genetic=FALSE, boruta=TRUE, univariate=TRUE, unipValThresh=0.05 ),  copy block by block

TODO: default syntax :set syntax=rmd

TODO: semi-supervised feature and genetic variable selection

TODO: Compare to null model for imbalanced datasets

BUG: variable correlation plots rendered off page.

BUG: preprocess crashes when trying too many image features

DOC: boxplot  https://www.r-bloggers.com/about-boxplot/ : the bottom and top of the box are always the 25th and 75th percentile (the lower and upper quartiles, respectively), and the band near the middle of the box is always the 50th percentile (the median). But the ends of the whiskers can represent several possible alternative values...”
In R’s default boxplot{graphics} code,

upper whisker = min(max(x), Q_3 + 1.5 * IQR) 
lower whisker = max(min(x), Q_1 – 1.5 * IQR)

where IQR = Q_3 – Q_1, the box length.
So the upper whisker is located at the *smaller* of the maximum x value and Q_3 + 1.5 IQR, 
whereas the lower whisker is located at the *larger* of the smallest x value and Q_1 – 1.5 IQR.


Compiled: `r format(Sys.time(), "%Y-%b-%d %H:%M:%S")` 

Target Variable: `r params$target`

Input File:      `r params$csvPath`

Target and inputs are column headings in csv file,
everything else is ignored


```{r Loading Packages + Data, echo=TRUE, warning=FALSE, message=TRUE}

params
options("width"=180)
# WIP reduce number of packages
libs <- c("caret", "kernlab", "knitr", "e1071", "magrittr", "rpart", "nnet", "parallel",
"randomForest", "xgboost", "Boruta", "leaps", "MASS",
"ranger", "cluster", "subselect", "corrplot", "gridExtra","pROC")

invisible(lapply(libs, require,character.only=T))
set.seed(25)

#Load data
datamatrix <- read.csv(params$csvPath)
allinputs <- colnames(datamatrix)
inputcols <- allinputs [grep("original_*",colnames(datamatrix))]
if(!is.null(params$inputs)){
  print(" manually select columns")
  inputcols <- allinputs [params$inputs]
} 

#error check
print(sprintf("target %s specified",params$target))
if(is.null(params$target)) stop("No target specified")

# Set model parameters
modelparams <- list(#tree = list(method  = "rpart",
                    #            #tuneGrid = data.frame(.cp = 0.01),
                    #            parms   = list(split="information")
                    #          ),
                    forest = list(method      = "rf",
                                  ntree       = 500,
                                 #tuneGrid    = data.frame(.mtry = mtry),
                                 #replace    = TRUE,
                                 #na.action  = randomForest::na.roughfix,
                                 importance  = FALSE,
                                 predict.all = FALSE
                                 ),
                     #na.action removed since "na.action" is used in caret
                    
                     #xgboost = list(method = "xgbLinear"),
 
                    # nnet = list(method = "nnet",
                    #             #tuneGrid=data.frame(.size = 10, .decay = 0),
                    #             #linout  = TRUE,
                    #             skip    = TRUE,
                    #             MaxNWts = 10000,
                    #             trace   = FALSE,
                    #             maxit   = 100),
                     
                     svm = list(method = "svmRadial"),
                     logit = list(method = "glm")
)

# Partition + Trim dataset for convenience
dsraw <- datamatrix[ c(params$target, inputcols ) ]

# TODO option for train/test/vallidate split
#set.seed(5)
#inTrain <- createDataPartition(as.factor(datamatrix[,params$target]), 1, p=0.7)
#  testing <- dsraw[-inTrain[[1]], ]
#  dsraw <- dsraw[inTrain[[1]], ]

```
## Pre-processing data: 
By default removes columns with zero variance and discards variables correlated >0.8
```{r Pre-Processing Data, echo=FALSE}
#  methods: zv removes zero-variance columns
#           corr removes highly correlated columns
#           center/scale recenters variables


#coerce to factor if needed
print(sprintf("target variable %s ",params$target) )  
if(!is.factor(dsraw[,params$target])){
  dsraw[params$target]<- as.factor(dsraw[,params$target])
  levels(dsraw[,params$target])<- c("control","case")
  params$positive_class <- "case"
}

#check if binary
if(dsraw[,params$target] %>% levels %>% length !=2){
    stop("Target must have two levels")
}   

#remove columns where all values are NA
# REMOVE ROWS WITH NA

narows <- sum(!complete.cases(dsraw))
nacols <- sum(colSums(is.na(dsraw)) == nrow(dsraw))
cat(sprintf("Removing %d rows for missing values\nRemoving %d NA columns",narows,nacols) )  

dsraw <- dsraw[complete.cases(dsraw),colSums(is.na(dsraw))<nrow(dsraw)]

#reduced input set
inputs <- setdiff(names(dsraw),params$target)

ppMethods <- "zv"
if(params$rescale)          {ppMethods <- c(ppMethods,"center","scale")}
if(params$removeCorrelated) {ppMethods <- c(ppMethods, "corr")}

#first remove zero variance (caret bug workaround)
pp <- preProcess(dsraw[inputs], method = "zv")
ds <- cbind(predict(pp, newdata = dsraw[inputs]),dsraw[params$target])
Filteredinputs <- setdiff(names(ds), params$target)

pp <- preProcess(predict(pp, newdata = ds[Filteredinputs]), method = ppMethods, cutoff=0.8)

#get reduced input set/dataset
ds <- cbind(predict(pp, newdata = ds[Filteredinputs]),ds[params$target])
Filteredinputs <- setdiff(names(ds), params$target)


# WIP: Clustering for semi-supervised analysis with kClusters
if(params$semisupervised){
cat(paste("Clustering into ", params$kClusters, " clusters... By unsupervised random forest\n\n"))
rfUL <- randomForest(x=ds[,Filteredinputs],
                     ntree = 500,
                     replace=FALSE)
ds <- cbind(ds, clusters.conv = pam(1-rfUL$proximity, k = params$kClusters, diss = TRUE, cluster.only = TRUE))
ds$clusters.conv <- as.factor(ds$clusters.conv) #change to factor not int
Filteredinputs <- setdiff(names(ds),params$target)
}

```

## Pre-Processing results: 
Started with `r length(inputs)` non-NA variables.
```{r,echo=FALSE}
pp
```
`r length(Filteredinputs)` remained after pre-processing

## Variable Selection: 

Default is Boruta and Wilcoxon test (P value cutoff 0.20/ `r length(Filteredinputs)`). Wilcoxon currently only tests numeric input variables

```{r Performing Variable Selection, echo=FALSE, warning=FALSE, fig.width=8, fig.height=8}

# Variable selection 
variableSelections <- list() #list(all=Filteredinputs)

#if <30 variables use all of them
if(Filteredinputs %>% length <= 30){
  variableSelections[["all"]] <- Filteredinputs

  correlations <- cor(Filter(is.numeric,ds[variableSelections$all]), use="pairwise")
   corrord <- order(correlations[1,])
   correlations <- correlations[corrord,corrord]
   corrplot(correlations,
   title = "Correlations for all variables",
   mar = c(1,2,2,0)  )
  }

#genetic
if(params$genetic){
gen <- genetic(cor(ds[,Filteredinputs]), 4) #manually selected 4 outputs
variableSelections$genetic <- names(ds[,Filteredinputs])[gen$bestsets]
print("Genetic variable selections:")
print(variableSelections$genetic)

correlations <- cor(Filter(is.numeric,ds[variableSelections$genetic]), use="pairwise")
   corrord <- order(correlations[1,])
   correlations <- correlations[corrord,corrord]
   corrplot(correlations,
   title = "Correlations for Genetic method",
   mar = c(1,2,2,0)  )

}

#univariate (currently only works for numeric variables)
if(params$univariate){

nums <- sapply(ds[Filteredinputs],is.numeric)
nums <- names(nums)[nums] #get names not T/F
pvals <- lapply(nums,
       function(var) {          
           formula    <- as.formula(paste(var, "~", params$target))
           test <- wilcox.test(formula, ds[, c(var, params$target)])
           test$p.value #could use 1-pchisq(test$statistic, df= test$parameter)
       })
names(pvals) <- nums

variableSelections$univariate <- setdiff(Filteredinputs, nums[pvals > params$unipValThresh]) #discard variables above threshold

if(variableSelections$univariate %>% length > 1 & params$plot){
correlations <- cor(Filter(is.numeric,ds[variableSelections$univariate]), use="pairwise")
#   corrord <- order(correlations[1,])
#   correlations <- correlations[corrord,corrord]
   corrplot(correlations,
   title = "Correlations for Univariate method",
   mar = c(1,2,2,0),
   order="hclust"  )
   

#par(mfrow=c(length(variableSelections$univariate),1))
for(iii in 1:length(variableSelections$univariate)){
boxplot( as.formula( paste(variableSelections$univariate[[iii]], "~", params$target) ),
         data=ds,
         ylab=paste(variableSelections$univariate[[iii]]),
         main=sprintf("Wilcoxon Rank Sum test P = %0.3e",pvals[variableSelections$univariate[[iii]]]),
         mar=c(12,1,0,0)
       )

 univariable= ds[variableSelections$univariate[[iii]]]
 cat(setdiff(levels(ds[,params$target]),params$positive_class), 'summary\n')
 print(summary(univariable[which(ds[,params$target] != params$positive_class ),]))
 cat('\n')
 cat(params$positive_class, 'summary\n')
 print(summary(univariable[which(ds[,params$target] == params$positive_class ),]))
 cat('\n')
}# end boxplot
} #end if >0 variables
} #end if params$univariate



#boruta
if(params$boruta){
bor <- Boruta(x=ds[,Filteredinputs], y=ds[,params$target], pValue=.35)
variableSelections$boruta <- names(ds[,Filteredinputs])[which(bor$finalDecision == "Confirmed")]

print("Finished Boruta variable selection")
print(bor)


if(length(variableSelections$boruta) > 1 & params$plot){

correlations <- cor(Filter(is.numeric,ds[variableSelections$boruta]), use="pairwise")
#   corrord <- order(correlations[1,])
#   correlations <- correlations[corrord,corrord]
   corrplot(correlations,
   title = "Correlations for Boruta method",
   mar = c(1,2,2,0),
   order = "hclust"  )
}


if(length(variableSelections$boruta) > 0 & params$plot){
#par(mfrow=c(length(variableSelections$boruta),1))
for(iii in 1:length(variableSelections$boruta)){
boxplot( as.formula( paste(variableSelections$boruta[[iii]], "~", params$target) ),
         data=ds,
         ylab=paste(variableSelections$boruta[[iii]]),
         mar=c(12,1,0,0),
         main=sprintf("Wilcoxon Rank Sum test P = %0.3e",pvals[variableSelections$boruta[[iii]]])
       )
 borutavariable= ds[variableSelections$boruta[[iii]]]
 cat(setdiff(levels(ds[,params$target]),params$positive_class), 'summary\n')
 print(summary(borutavariable[which(ds[,params$target] != params$positive_class ),]))
 cat('\n')
 cat(params$positive_class, 'summary\n')
 print(summary(borutavariable[which(ds[,params$target] == params$positive_class ),]))
 cat('\n')
}# end boxplots
} #end if >0 variables
} #end if params$boruta


```


```{r Univariate selection, echo=FALSE}

#print results and error handle if no significant variables
if(params$univariate){
  if(length(variableSelections$univariate) == 0){ 
    cat(sprintf("No P values < %.5f\n", 0.2/length(nums)))
    variableSelections$univariate <- NULL
    } 
   
    pvaltable <- data.frame(Variable = Filteredinputs[pvals<.2], 
                        P_value = unlist(pvals[pvals<.2] )
                        )
    pvaltable <- pvaltable[order(pvaltable$P_value, decreasing=F),]
    cat("Listing variables with WILCOXON RANK SUM test P value < 0.2")
    kable(pvaltable, format = "markdown", row.names=F)
    }

```




# Modeling using `r names(modelparams)`.
Use Leave-one-out cross validation: `r params$leaveOneOut`


```{r Modeling (may take some time), echo=FALSE, warning=FALSE, message=FALSE}

#WIP: consider using metric=roc for binary classification

#Modeling, dataparams are arguements to caret::train
modelformula <- as.formula(paste(params$target,"~."))
dataparams   <- list(form = modelformula,
#                      data = ds[,c(params$target,Filteredinputs)],
                      metric="Accuracy", #other option: AUC?
                      trControl=trainControl(allowParallel  = T,
                                             method = ifelse(params$leaveOneOut,"LOOCV", "repeatedcv"),
                                             classProbs=TRUE,
                                             #returnResamp = "final",
                                             #number = 10,
                                             #repeats= 5,
                                             verboseIter = F) # use method="none" to disable grid tuning for speed
                     )  
caretparams      <- lapply(modelparams,function(x) c(dataparams,x))

#initialize outputs
models    <- list()
acc <- list()

for(jjj in 1:length(variableSelections)){
modeldata  <- ds[,c(params$target,variableSelections[[jjj]])]

# RE seeding: Hawthorn et al, "The design and analysis of benchmark experiments" (2005)

for (iii in 1:length(modelparams)){
model_name <-paste(names(variableSelections)[jjj], names(modelparams)[iii],sep="_")
#print(paste("Training model", iii, "of", (length(modelparams)*length(variableSelections)), ":", model_name, sep=" "))
set.seed(3141) #seed before train to get same subsamples
invisible(  models[[model_name]] <- do.call(caret::train, c(caretparams[[iii]], list(data=modeldata)))  )

metric <- models[[model_name]]$metric

#get best acuracy manually for LOOCV
acc[[model_name]] <- max(models[[model_name]]$results[metric])

}


}


```


```{r Finding best model and plotting, echo = FALSE, fig.width=8, fig.height=10}
if(params$leaveOneOut){
  maxacc <- max(unlist(acc))
  maxmodels <- names(which(acc==maxacc))

  #plot
  par(mar=c(8,5,1,1))
  barplot(unlist(acc),las=2, ylim=c(0,1),
    ylab=paste("training set LOOCV",metric))
  
} else {
  rs <- resamples(models)
  print(summary(object = rs))
  
  acc        <- rs$values[,grepl(rs$metrics[1], names(rs$values))]
  names(acc) <- rs$models
  maxacc <- max(apply(acc,2,mean))
  maxmodels <- rs$models[apply(acc,2,mean)==maxacc]
  
  # plot +1 to col arg keeps one box from being black
  par(mar=c(8,5,1,1))
  boxplot(acc,col=(as.numeric(as.factor(rs$methods))+1), las=2,
    ylab=paste("training set cross-validation",metric) )
  legend("bottomright", legend=unique(rs$methods),
    fill=(as.numeric(as.factor(rs$methods))+1) )
}

#get best accuracy
#cat(paste("best model(s): ", maxmodels, "\n"))
#cat(sprintf("%s: %.4f \n", metric, maxacc))
#lapply(maxmodels, function(x) models[[x]])
```

Best model(s): `r maxmodels`

`r sprintf("%s: %.4f", metric, maxacc)`

```{r Output model,echo = FALSE}
lapply(maxmodels, function(x) models[[x]])
```

```{r if LOOCV CV plot ROC,echo = FALSE, fig.width=8, fig.height=10}
if(params$leaveOneOut){
  print(paste("Building ROC Curve for model", maxmodels[[1]]))
  #recursiely subset
  bestmod <- models[[maxmodels[[1]]]] #pick first by default
  modpars <- bestmod$bestTune
  results <- bestmod$pred
  
  #filter predictions on best model
  for(kkk in 1:length(modpars)){
    results <- subset(results, results[,names(modpars)[kkk]]==modpars[,kkk])
  }
  
  print(caret::confusionMatrix(results$pred, reference = results$obs, positive = params$positive_class))
  ROC1 <- roc(response=results$obs,
              predictor = results[,params$positive_class],
              levels = c(setdiff(
                levels(results$obs), params$positive_class),
                         params$positive_class))
  print(ROC1)
  plot(ROC1, print.auc=T, print.auc.y = 0.2, print.auc.x = 0.5)
  # results is a frame with LOOCV data
  
  #best threshold based on sum of sens/spec
  kable(coords(roc=ROC1, x = "best", ret=c('threshold','sens', 'spec')))
  
}
```


## Testing best model on validataion cohort:
```{r experimental cohort, echo = FALSE}

if(!is.null(params$test_csv)){
cat("Reading validation dataset:\n", params$test_csv)
full_data <- read.csv(params$test_csv)

bestModelInputs <- names(bestmod$trainingData[-1])

cat(c("Best model using inputs:", paste(bestModelInputs,collapse="\n"),"\n"))

sub_data <- full_data[,c(bestModelInputs, params$target)]

cat(sprintf("Number of cases ommitted for missing values %d out of %d\n", sum(!complete.cases(sub_data)), nrow(sub_data)))

new_data = sub_data[complete.cases(sub_data),]

preds <- predict(object=bestmod$finalModel,newdata = new_data,type="prob")

kable(preds)
} else {cat("No validation dataset provided")}


```