binary_analysis.RMD

---
title: "Binary Classification Analysis"
output: pdf_document
params:
  csvPath:            exampledata.csv
  outCsvName:        !r NULL
  target:             NULL
  positive_class:     NULL
  regexp:            !r "(Volume|Mean)" 
  fixed_inputs:      !r NULL
  exclude_inputs:    !r NULL
  nfold:             !r NULL
  foldID:            !r NULL
  leaveOneOut:       !r TRUE
  test_csv:          !r NULL
  rescale:           !r FALSE
  removeCorrelated:  !r TRUE
  cutoff:               0.7
  plot:              !r TRUE
  l1lasso:    !r FALSE
  kClusters:         !r as.numeric(9)
  univariate:        !r TRUE
  unipValThresh:     !r 0.05

---
## Version: BINARY CLASSIFIER

DOC: Usage - Rscript -e " rmarkdown::render( 'binary_analysis.RMD', output_file = './binary_exampleout.pdf', params = list( csvPath='dicom/radiomicsout.csv', outCsvName = 'exampledata_binary_withpredictions.csv',target='fibbin', positive_class='other', regexp='(original_firstorder|original_glcm|original_glrlm|original_glszm)', fixed_inputs=NULL,leaveOneOut=TRUE, nfold=NULL, foldID=NULL, rescale=FALSE, removeCorrelated=TRUE, plot=TRUE, l1lasso=TRUE, univariate=TRUE, unipValThresh=0.05,cutoff=.7 ) )"


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 libs, echo=FALSE, warning=FALSE}
sessionInfo()
# hack for debuggin
myparams = params
myparams 
options("width"=80)
# WIP reduce number of packages
libs <- c("caret", "kernlab", "knitr", "e1071", "magrittr", "rpart", "nnet", "parallel",
"randomForest", "MASS","cluster", "corrplot", "gridExtra","pROC","glmnet")
# "xgboost", "leaps", "ranger", "subselect", 

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



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

# load data
datFull <- read.csv(myparams$csvPath)

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

inputs <- union(myparams$fixed_inputs,
                grep(myparams$regexp, names(datFull),value=T))


print(sprintf("%d rows and %d columns read from csv file: %s", nrow(datFull), ncol(datFull), myparams$csvPath))
print(sprintf("%d inputs found with regexp parameter and fixed inputs", length(inputs)))


# remove excluded inputs if defined
if (!is.null(myparams$exclude_inputs)){
  inputs <- union( setdiff(inputs, grep(myparams$exclude_inputs, names(datFull),value=T)), myparams$fixed_inputs)
  print(sprintf("%d inputs remained after excluded inputs", length(inputs)))
} else {
  print("No excluded_inputs from parameters")
}


# need binary vars
#datFull <- cbind(
#    fibbin = as.factor(ifelse(datFull[,"fibrosis"] <= 1                  , "zeroone", "other")),
#    datFull )
# remove unnecessary columns and clean NAs
datCrop <- datFull[,c(inputs,myparams$target)]


#remove incomplete cases
incompCases <- !complete.cases(datCrop)
print("The following rows are excluded for missing data:")
print(row.names(datFull)[incompCases])

datCompleteCases <- datCrop[!incompCases,]

if(nrow(datCompleteCases) < 5){
  stop("Less than 5 rows had no missing input values")
}

# remove inputs with NA or infinite values (should do nothing if complete cases was applied)
dat <- datCompleteCases[, (colSums(is.na(datCompleteCases)) == 0 ) & (sapply(datCompleteCases, function(x) sum(is.infinite(x)) ) == 0) ]


# subset based on validation set if provided
if(!is.null(myparams$foldID)){
  set.seed(5)
#TODO: translate validation indices back to full data frame

  if(!("validationFolds" %in% names(datFull)) ){
    datFull$validationFolds <- rep(NA,nrow(datFull))
    validationFolds <- createFolds(dat[,myparams$target],k=myparams$nfold,list=F)

    # add kfold indices to data matrix
    dat$validationFolds <- validationFolds
    datFull[row.names(dat),"validationFolds"] <- validationFolds

  } else {
    print(sprintf("validationFolds column already exists in the csv, reading it instead of generating folds."))
    dat$validationFolds <- datFull[row.names(dat),"validationFolds"]
  }

  validationidx <- which(validationFolds == myparams$foldID) #gives row numbers (not row names), this will not match the original (nonfiltered) data
  print(paste("Leaving out ", myparams$validationID, " fold number ", myparams$foldID,sep=""))

  #bad code: order of these two lines matters
  validationdat <- dat[validationidx,]

  print("removing validation fold from data frame")
  dat <- dat[-validationidx,]
}
print("summary of classes")
print(summary(dat[,myparams$target]))

# Legacy code: old script uses dsraw instead of dat
dsraw <- dat

# 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")
)



```
## 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
# careful with factor levels
print(sprintf("target variable %s ",myparams$target) )  



if(!is.factor(dsraw[,myparams$target])){
  dsraw[myparams$target]<- as.factor(dsraw[,myparams$target])
  print("changed target variable to factor")
  print(sprintf("initial level: %s", levels(dsraw[,myparams$target])))

  if(!is.null(myparams$positive_class)) {
  positive_class <- myparams$positive_class
  } else {
#   stop(myparams$positive_class)
  positive_class <- "case"
  }

  levels(dsraw[,myparams$target]) <- c("control",positive_class)  
  print(sprintf("reassigned target levels to: %s", levels(dsraw[,myparams$target])))
}

if(!is.null(myparams$positive_class)) {
  positive_class <- myparams$positive_class
} else {
  positive_class <- levels(dsraw[,myparams$target])[2]
  print(sprintf("assuming positive class is: %s", positive_class))
}

#check for positive class in factor levels
if(!(positive_class %in% levels(dsraw[,myparams$target])) ){
  stop(sprintf("positive class %s  not found in levels of target variable %s", positive_class, levels(dsraw[,myparams$target])))
}


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

#rename positive class if needed
#levels(dsraw[,myparams$target])<- c("control","case")
#positive_class <- "case"


#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),myparams$target)

#ppMethods <- "zv"
#if(myparams$rescale)          {ppMethods <- c(ppMethods,"center","scale")}
#if(myparams$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[myparams$target])

# reduce input set after correlation reduction
FilteredinputsZeroVar <- names(ds)[names(ds) %in% inputs]
cor_matrix <- cor(ds[FilteredinputsZeroVar]) 
cor_high <- findCorrelation(cor_matrix, myparams$cutoff)
high_cor_remove <- row.names(cor_matrix)[cor_high] 

#get reduced input set/dataset
Filteredinputs <- setdiff(FilteredinputsZeroVar, high_cor_remove )
ds <- cbind(ds[Filteredinputs],ds[myparams$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 Wilcoxon test (P value cutoff `r myparams$unipValThresh`/ `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)  )
   print(" ")
  }

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, "~", myparams$target))
           test <- wilcox.test(formula, ds[, c(var, myparams$target)])
           test$p.value #could use 1-pchisq(test$statistic, df= test$parameter)
       })
names(pvals) <- nums
```

```{r Lasso Section (if TRUE), echo=FALSE, warning=FALSE, fig.width=8, fig.height=8}
# l1 regression
if(myparams$l1lasso){
lassodat = na.omit(ds)
glmnetfit <- cv.glmnet(as.matrix(lassodat[Filteredinputs]) , as.numeric(unlist(lassodat[myparams$target])), alpha = 1.0, nlambda = 100,family = "binomial", intercept = FALSE)
plot(glmnetfit )
lassocoef = as.matrix(coef(glmnetfit , s = "lambda.min"))
lassonames = rownames(lassocoef ) 
# WIP - filter 0 coeff
variableSelections[["lasso"]] <- lassonames[abs(lassocoef) > 0]

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

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

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

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


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

if(variableSelections$univariate %>% length > 1 & myparams$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]], "~", myparams$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[,myparams$target]),myparams$positive_class), 'summary\n')
 print(summary(univariable[which(ds[,myparams$target] != myparams$positive_class ),]))
 cat('\n')
 cat(myparams$positive_class, 'summary\n')
 print(summary(univariable[which(ds[,myparams$target] == myparams$positive_class ),]))
 cat('\n')
}# end boxplot
} #end if >0 variables
} #end if myparams$univariate




```


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

#print results and error handle if no significant variables
if(myparams$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 myparams$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(myparams$target,"~."))
dataparams   <- list(form = modelformula,
#                      data = ds[,c(myparams$target,Filteredinputs)],
                      metric="Accuracy", #other option: AUC?
                      trControl=trainControl(allowParallel  = T,
                                             method = ifelse(myparams$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()

```
```{r Modeling part duex, echo=FALSE, warning=FALSE, message=FALSE}
for(jjj in 1:length(variableSelections)){
modeldata  <- ds[,c(myparams$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(myparams$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(myparams$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 = myparams$positive_class))
  ROC1 <- roc(response=results$obs,
              predictor = results[,myparams$positive_class],
              levels = c(setdiff(
                levels(results$obs), myparams$positive_class),
                         myparams$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')))
  
}
```

```{r write csv with validation predictions, echo=FALSE}
#write cross-validation predictions


if(!is.null(myparams$foldID)){
  trainingpredname <- paste("trainindPreds_exceptfold",myparams$foldID,sep="")

  if(trainingpredname %in% names(datFull)){
    print(sprintf("cross-validation predictions for fold %s already found, compare differences below:", myparams$foldID))
    print("existing predictions (may be out of order):")
    kable(data.frame(existing=datFull[rownames(dat),trainingpredname],new=as.character(results$pred)))

  } else {
  datFull[trainingpredname] <- NA
  datFull[rownames(dat)[results$rowIndex],trainingpredname] <- as.character(results$pred)  #trickey line of code since results$rowIndex referes to rows of dat not datFull
    }

} else { #if no validation fold, write all predictions
  trainingpredname <- "trainingPreds"
  datFull[trainingpredname] <- NA
  datFull[rownames(dat)[results$rowIndex],trainingpredname] <- as.character(results$pred)  #trickey line of code since results$rowIndex referes to rows of dat not datFull
}


# write validation predictions
if(!is.null(myparams$foldID)){

predname <- paste("validationPreds_fold",myparams$foldID,sep="")


  if(predname %in% names(datFull)){
    print(sprintf("Validation predictions for fold %s already found, compare differences below:", myparams$foldID))
    print("existing predictions:")
    print(datFull$predname)
  } else {

datFull[predname] <- NA
datFull[which(datFull$validationFolds == myparams$foldID), predname] <- validationpreds
  }
}

print("summary of predictions:")

if(!is.null(myparams$foldID)){
  kable(datFull[,c(myparams$target,predname,trainingpredname)])
} else {
  kable(datFull[,c(myparams$target,trainingpredname)])
}


if(is.null(myparams$outCsvName)){
  print("No output csv file specified")

# Option to auto generate csv name
#  if(!is.null(myparams$foldID)){
#    outCsv <- sub(".csv",sprintf("_Fold%s_validation.csv",myparams$foldID),basename(myparams$csvPath)) 
#  } else {
#    outCsv <- sub(".csv","_no_validation.csv",basename(myparams$csvPath))
#  }

} else {
    outCsv <- myparams$outCsvName
}

print(sprintf("writing csv with validation predictions to %s", outCsv))
write.csv(datFull,outCsv,row.names=F)


```



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

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

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

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

sub_data <- full_data[,c(bestModelInputs, myparams$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")}


```