---
title: "Big Data Machine Learning Methods for Metabolomics"
subtitle: Part I
author: 
  - Kangni Alemdjrodo, Department of Statistics, Purdue University
  - Dabao Zhang, Department of Statistics, Purdue University
  - Daniel Raftery, Department of Anesthesiology and Pain Medicine, University of Washington
  - Min Zhang, Department of Statistics, Purdue University
date: "Summer 2022"
output:
  pdf_document:
    number_sections: yes
    toc: yes
    toc_depth: 3
  html_document:
    number_sections: yes
    toc: yes
    toc_depth: 3
---

```{r, setup, include=FALSE}
knitr::opts_chunk$set(message=FALSE, warning=FALSE, results="hide")
```


# Introduction

This notebook aims to give an overview on popular multivariate methods used in metabolomics, such as logistic regression and partial least squares-discriminant analysis (PLS-DA). It also introduces the Extreme Gradient Boosting (XGBoost) method which builds a tree ensemble model, and our newly developed method gPOCRE-based screening (GPS) which can improve metabolomics data analysis by selecting a small subset of features from a large number of candidates.


## Required R packages

Load the necessary R packages for the following data analyses.

```{r, message=FALSE, warning=FALSE}
library(caret)   ## For model tuning using resampling
library(dplyr)   ## For some data preparation
library(Matrix)  ## For sparse matrix
library(missRanger) ## For missing value imputation
library(mixOmics)## For pls-da
library(POCRE)   ## For gps
library(ROCR)    ## For ROC curves
library(SIS)     ## For sis
library(splitTools) ## For stratified samples
library(xgboost) ## For xgboost
```



## The Data 

The primary data (in .zip format) is available [here](https://figshare.com/articles/dataset/Metabolomic_and_Lipidomic_Analysis_of_human_CSF/14816622), which includes the following two .csv files:

* `subject_info.csv`: including the age, disease type (AD, PD or control), and sex of the subjects; 
* `targeted.csv`: including abundances of 108 known metabolites.

The data consist of lipidomic and metabolomic assays from the cerebrospinal fluid of 198 CSF samples, collected from 85 controls, 57 patients of Alzheimer's disease (AD), and 56 patients of Parksinson's disease (PD) (note that our analysis will exclude the PD patients). Each subject was also recorded with the age and sex.



## Preprocessing the Data

For demonstration purpose, we will focus on predicting Alzheimer's disease (AD). Therefore, we only include the healthy controls (including CO, CM, CY) and AD patients in the following analysis.

**Caution:** Before you continue to the next step, you have to download the data, copy the two files `targeted.csv` and `subject_info.csv` to your working directory, and set up your working directory as follows,

```{r}
# setwd("<you_working_directory>")
```

We first load the two data sets, and then exclude AD patients
```{r}
# Load the data sets
targeted  <- read.csv("targeted.csv") 
subject_info <- read.csv("subject_info.csv")

# Cohorts
xtabs(~Type,data=subject_info)

# Reclassify "CM", "CO", "CY" as "Control"
subject_info[subject_info$Type=="CM",]$Type <- "Control"
subject_info[subject_info$Type=="CO",]$Type <- "Control"
subject_info[subject_info$Type=="CY",]$Type <- "Control"

# Exclude PD patients
targeted <- targeted[subject_info$Type!="PD",]
subject_info <- subject_info[subject_info$Type!="PD",]
xtabs(~Type,data=subject_info)
```

**Caution:** The metabolite levels included in `targeted.csv` has already been normalized. Otherwise, you have to normalize the metabolite levels first.



###  Missing Value Imputation

A brief check of the data can show that some missing values represented by `NA` in recorded metabolomic information,

```{r}
#head(targeted)
summary(targeted)
summary(subject_info)
```


While we don't have any missing value in `subject_info`, we can check which metabolite has missing values and how many missing values there are.
```{r}
which(is.na(colSums(targeted)))
sum(is.na(targeted)) 
```

It turns out that there are a total of 6 missing values spreading out in 6 metabolites.

**Caution:** If there are too many missing values (e.g., over 20\%) in a single subject or single metabolite (covariate), we may have to exclude that subject or metabolite (covariate) from the followup analysis. Otherwise, we can impute the missing values via different strategies.

The simplest method of missing value imputation is replacing with mean/median. However, it may be risky, especially in the case that many missing values in observed metabolites may be caused by insufficient abundance of metabolites. A better practice is predicting via statistical models, like KNN and random forest. Here we will impute the missing values via random forest as proposed by [Stekhoven and Buehlmann (2012)](https://academic.oup.com/bioinformatics/article/28/1/112/219101?login=true).


```{r, message=FALSE, warning=FALSE}
set.seed(20220620) ## To replicate the results
targeted <- missRanger(targeted[, -1], num.trees = 100) ## Remove id before imputation
```

We can check if the data is complete now.

```{r}
sum(is.na(targeted)) ## number of missing values
```



### Outliers

An outlier is a value that is extremely large/small or remote compared to the majority of its peers. It can be caused by a measurement error, or may be an incorrectly reported or misplaced value. If the outliers are not appropriately processed before the data analysis, the results may be distorted and misleading.

There are different graphic methods to check possible outliers, such as histogram and boxplot. You may also conduct hypothesis tests on possible outliers, such as Rosner's test for outliers (available as `rosnerTest()` in R package **EnvStats**). However, as there are many metabolites, graphic check on all metabolites may demand too many work. On the other hand, Rosner's test indeed check possible outliers from a Gaussian distribution but it may be unrealistic to assume that a metabolite's levels follow a Gaussian distribution even that we have normalized the data. 

Here we start with calculating the largest possible deviation (with median absolute deviation, i.e., MAD, as unit) from its median value for each metabolite, 

```{r}
mad_t  <- sapply(targeted, mad) ## Calcuate MAD for each metabolite
median_t <- sapply(targeted, median) ## Calcuate median for each metabolite

nMAD_t <- apply(targeted,1,function(x) abs(x-median_t)/mad_t)
max_nMAD_t <- apply(nMAD_t,1,max)
hist(max_nMAD_t)
which(max_nMAD_t>=10) # Metabolites with possible outliers
```

We can check the histograms and boxplots of those metabolites with large deviations,

```{r}
ids <- c(16, 19, 49, 71, 83, 105)
mnames <- colnames(targeted)
for( j in ids ) {
  par(mfrow = c(1, 2))
  hist(as.matrix(targeted[,j]), main=paste(mnames[j],"(",j,"): Histogram",sep=''))
  boxplot(targeted[,j], main = paste(mnames[j],"(",j,"): Boxplot",sep=''))
}
```


While the boxplots show outliers for each of these metabolites, we observe that those extremal values are attached to skewed tails, and may not really be inaccurate values. Therefore, we will keep these values as they are in the following analysis.

We prepare `subject_info` by converting both `Type` and `Sex` into factors with properly defined levels.

```{r}
subject_info <- subject_info[,-1] ## Remove the id variable
# Type: 0-Contro; 1-AD
subject_info$Type <- factor(subject_info$Type,levels=c("Control","AD"),labels=c(0,1))
# Sex: 0-M; 1-F
subject_info$Sex <- factor(subject_info$Sex,levels=c("M","F"),labels=c(0,1))
summary(subject_info)
```



# Predictive Models with Main Effects of Metabolite Levels


Here we intend to build predictive models to predict AD patients using only metabolomic information. We first create `csf_m`, which include both disease type (control vs. AD) and standardized metabolites levels. 

```{r}
subject_info$Type <- as.numeric(as.character(subject_info$Type))
#subject_info$Sex <- as.numeric(as.character(subject_info$Sex))
#subject_info$Age <- as.numeric(subject_info$Age)
csf_m <- data.frame(Type=subject_info$Type, scale(targeted))
```

Here we use the data `csf_m` to build different predictive models with only main effects of 108 metabolites from 85 `Control` and 57 `AD` subjects.


**Machine Learning Methods and Their Performance Evaluation**: Here we will consider three classification methods, i.e., logistic regression, PLS-DA, and XGBoost. In general, we will build and evaluate the model based on cross validation. That is, for a specific method, we will splitting the available data into training and test data sets, build the predictive model based on the training data set, and then evaluate its predictability based on the test data set.

* **Remark 1.** We can calculate the area under the receiver operating characteristic (ROC) curve (AUC) to gauge the performance of the predictive model built by each method. ROC is the plot of the true positive rates against false positive rates for all possible thresholds. The AUC is always between 0 and 1. The closer to 1 the better the built model.
* **Remark 2.** When there is a large sample size, we only have to split the data once into a set of training data and a set of test data. However, when the sample size is small, we may have to take the K-fold cross-validation, which splits the data into K folds (each with the same sample size), and we can take each fold to serve as the test data once while using the other K-1 folds as the training data.
* **Remark 3.** When the sample size is small, K-fold cross-validation may still rely too much on the K folds. For example, we may consider 5-fold cross-validation so each time we include 80\% of the samples in the training set and the other 20\% in the test set. However, the result of AUC may rely on the randomness of the split 5-fold. Therefore, we may better randomly split the data for many times, say 100, each time with training and test sets including 80% and 20% of the samples, respectively.
* **Remark 4.** To keep the data structure, we usually stratify each splitting within `Control` and `AD` subjects to preserve the same ratio of `Control` and `AD` in both training and test sets. When there are other categorical variables (like `Sex`) involved, we should also stratify the splitting based on these categorical variables.




## Predictive Models with All Metabolites

###  Logistic regression

[Logistic regression](https://www.sciencedirect.com/topics/computer-science/logistic-regression) builds a statistical model to describe the probability of a binary event (`Control` vs `AD` here) based on a set of predictors (metabolite levels here). The probability is modeled as a logistic function.

With the built logistic regression model, we can calculate the probability of `AD` for a subject with metabolite profile and may claim the subject as `AD` when the probability is above a threshold (as `Control` otherwise).

```{r, warning=FALSE}
set.seed(20220620) ## To replicate the results
start.time <- Sys.time() # Start the clock to calculate computing time
nsim <- 100 ## Number of times to split the data into training & test sets

csf_m.logit.test.preds <- vector(mode="list", length=nsim)
csf_m.logit.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  train.idx <- partition(csf_m$Type, p=c(train=0.8, test=0.2))
  train <- csf_m[train.idx$train,]
  test <- csf_m[train.idx$test,]
  
  ## Fit LOGISTIC 
  csf_m.logit.train <- glm(Type~., data=train, family=binomial(link="logit"))

  ## Predict
  csf_m.logit.test.preds[[j]] <- predict(csf_m.logit.train, test, type="response")
  csf_m.logit.test.labels[[j]] <- test$Type
}

csf_m.logit.ROCR <- ROCR::prediction(csf_m.logit.test.preds, csf_m.logit.test.labels)
csf_m.logit.roc <- performance(csf_m.logit.ROCR,"tpr","fpr")
csf_m.logit.auc <- performance(csf_m.logit.ROCR, measure = "auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```


Let us plot the ROC curve and get the AUC.

```{r}
#mar.default <- c(5,4,4,2) + 0.1
#par(mar = mar.default + c(0, 2, 0, 0)) 
par(mar = c(5,6,4,2) + 0.1) # (bottom, left, top, right) margins in numbers of lines
plot(csf_m.logit.roc,col="grey50",lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_m.logit.roc,lwd=3,avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="Logistic Regression with All Metabolites")
text(x=0.8, y=0.1,paste("AUC=", round(median(unlist(csf_m.logit.auc@y.values)),2)),cex = 1.5)  # AUC
abline(a=0, b=1) 
par(new=FALSE)
```

While the above plot shows with AUC=0.65, we can further confirm this and calculate its standard error.

```{r}
round(median(unlist(csf_m.logit.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_m.logit.auc@y.values)),2)/10  ## standard error
```



### PLS-DA 

The [PLS-DA](https://analyticalsciencejournals.onlinelibrary.wiley.com/doi/pdfdirect/10.1002/cem.785?casa_token=jRc0OxmFfhoAAAAA:rWNzDs13fIj7fSuhQ4qdt0uqIynAyw3PzYsBbhoPBTx213HELD7OLfGxqWDdsuap6vL3oJW8ZvlnBOE) is the discriminant version of the Partial least squares (PLS). As for the discrimination of principal components, the goal is to seek a model of orthogonal components constructed from linear combinations of centered explanatory variables. PLS-DA replaces the initial space of explanatory variables with a low-dimensional space, implied by a small number of variables called "factors" or "latent variables," which are iteratively built one after the other. These factors constitute the new explanatory variables of the model. They are orthogonal (uncorrelated) and are linear combinations of the initial explanatory variables.

A PLS-DA model is fitted with 10 components to evaluate the performance and the number of components necessary for the final model (see below).

```{r}
csf_m.plsda <- plsda(csf_m[,-1], csf_m$Type, ncomp=10)  # set ncomp to 10
```

The samples are then projected into the subspace spanned by the first two components.

```{r}
plotIndiv(csf_m.plsda, comp=1:2, group=csf_m$Type, ind.names=FALSE,
          ellipse=TRUE, legend=TRUE, title='PLS-DA with All Metabolites')
```

From the sample plot, we observe a better separation of 2 groups compared to an unsupervised PCA sample plot (skipped). Confidence ellipses for each class are plotted to highlight the strength of the discrimination (confidence level set to 95% by default, see argument ellipse.level). 

The prediction area can be visualized by calculating a background surface first, before overlaying the sample plot. See `? background.predict` for more details. The prediction distance needs to be specified:

```{r}
bg.csf_m.plsda <- background.predict(csf_m.plsda, comp.predicted=2, dist="max.dist")
plotIndiv(csf_m.plsda, comp=1:2, group=csf_m$Type, ind.names=FALSE, title="Maximum Distance",
          legend=TRUE,  background=bg.csf_m.plsda)
```



**Performance**

The classification performance of the PLS-DA model is assessed with the `perf` function using 5-fold
cross-validation repeated 50 times. The number of repeats is necessary to ensure a good estimate of the classification error rate (as the CV-folds are determined in a random manner). From the performance results we can decide on the choice of the number of components of the final PLSDA model. The maximal number of components assessed corresponds to the number of components of the `csf_m.plsda` model ran above.

This takes a couple of minutes (2 min) to run.

```{r}
set.seed(20220620) # for reproducibility, only when the `cpus' argument is not used
start.time <- Sys.time() # Start the clock!

csf_m.plsda.perf <- perf(csf_m.plsda, validation="Mfold", folds=5,
                         progressBar=FALSE, auc=TRUE, nrepeat=50)

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken) 
```

**Error rates**
```{r}
plot(csf_m.plsda.perf, col=color.mixo(5:7), sd=TRUE, legend.position="horizontal")
```


From the performance plot, we observe that the centroid error rate decreases slowly from 1 to 6 components. The error rates does not stabilize after 6 components.

The perf function outputs the optimal number of components in `csf_m.plsda.perf$choice.ncomp` based on t-tests that test for a significant difference in the mean error rate between components. The result is output for each performance measure and distance, when applicable.

```{r}
ncomps.csf_m.plsda <- csf_m.plsda.perf$choice.ncomp
ncomps.csf_m.plsda
```

```{r}
csf_m.plsda.perf$error.rate
```

**Conclusion:** The performance displayed suggests that `ncomp = 5`  with centroids distance is sufficient to achieve performance (0.31 error rate).


#### PLS-DA final model performance

```{r}
ncomps.csf_m.plsda = 5 ## Number of component for PLS-DA 
set.seed(20220620) ## To replicate the results
start.time <- Sys.time() # Start the clock!
nsim <- 100 ## Number of repetitions
csf_m.plsda.test.preds <- vector(mode="list", length=nsim)
csf_m.plsda.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  train.idx <- partition(csf_m$Type, p=c(train=0.8, test=0.2))
  train <- csf_m[train.idx$train,]
  test  <- csf_m[train.idx$test,]
  
  ## Fit PLS-DA
  csf_m.plsda.train <- plsda(train[,-1], train$Type, ncomp=ncomps.csf_m.plsda)
  
  ## Predict
  csf_m.plsda.test.preds[[j]] <- predict(csf_m.plsda.train, test[,-1], 
                                         dist="centroids.dist")$predict[,,ncomps.csf_m.plsda][,2]
  csf_m.plsda.test.labels[[j]]<- test$Type
}

csf_m.plsda.ROCR = ROCR::prediction(csf_m.plsda.test.preds, csf_m.plsda.test.labels)
csf_m.plsda.roc = performance(csf_m.plsda.ROCR,"tpr","fpr")
csf_m.plsda.auc = performance(csf_m.plsda.ROCR, measure="auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```


**ROC curve and AUC**

```{r}
#mar.default <- c(5,4,4,2) + 0.1
#par(mar = mar.default + c(0, 2, 0, 0)) 
par(mar = c(5,6,4,2) + 0.1) # (bottom, left, top, right) margins in numbers of lines
plot(csf_m.plsda.roc,col="grey50", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_m.plsda.roc, lwd=3, avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="PLS-DA with All Metabolites")
text(x=0.8, y=0.1, paste("AUC=", round(median(unlist(csf_m.plsda.auc@y.values)),2)), cex=1.5)  # AUC
abline(a=0, b=1)
par(new=FALSE)
```

From the above figure, we see that PLS-DA could build predictive models with all metabolites with AUC at 0.72, which is confirmed in the following with standard error.

```{r}
round(median(unlist(csf_m.plsda.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_m.plsda.auc@y.values)),2)/10  ## standard error
```



### XGBoost

[XGBoost](https://dl.acm.org/doi/pdf/10.1145/2939672.2939785) combines the advantages of bagging (bootstrap aggregating) and boosting and effectively reduces the bias-related error and variance-related error of the model. Recall that bagging and boosting are ensemble learning methods that can integrate decision trees to reduce the model error.

We will take XGBoost to build tree ensemble models (XGBoost-T) and linear-regression ensemble models (XGBoost-L), respectively.



#### XGBoost-T

##### Determine optimal parameters to use with the XGBoost-T


Parameters for Tree Booster

**nrounds(default=100)**
 
* It controls the maximum number of iterations. For classification, it is similar to the number of trees to grow.
* Should be tuned using CV

 	
**eta[default=0.3][range: (0,1)]**

* It controls the learning rate, i.e., the rate at which our model learns patterns in data. After every round, it shrinks the feature weights to reach the best optimum.
* Lower eta leads to slower computation. It must be supported by increase in nrounds.
* Typically, it lies between 0.01 - 0.3
	
**gamma[default=0][range: (0,Inf)]**

* It controls regularization (or prevents overfitting). The optimal value of gamma depends on the data set and other parameter values.
* Higher the value, higher the regularization. Regularization means penalizing large coefficients which don't improve the model's performance. default = 0 means no regularization.
* Tune trick: Start with 0 and check CV error rate. If you see train error >>> test error, bring gamma into action. Higher the gamma, lower the difference in train and test CV. If you have no clue what value to use, use gamma=5 and see the performance. Remember that gamma brings improvement when you want to use shallow (low max_depth) trees.
	
**max_depth[default=6][range: (0,Inf)]**

* It controls the depth of the tree.
*Larger the depth, more complex the model; higher chances of overfitting. There is no standard value for max_depth. Larger data sets require deep trees to learn the rules from data.
*Should be tuned using CV
    
**min_child_weight[default=1][range:(0,Inf)]**

* In regression, it refers to the minimum number of instances required in a child node. In classification, if the leaf node has a minimum sum of instance weight (calculated by second order partial derivative) lower than min_child_weight, the tree splitting stops.
* In simple words, it blocks the potential feature interactions to prevent overfitting. Should be tuned using CV.
    
**subsample[default=1][range: (0,1)]**

* It controls the number of samples (observations) supplied to a tree.
* Typically, its values lie between (0.5-0.8)
    
**colsample_bytree[default=1][range: (0,1)]**

* It control the number of features (variables) supplied to a tree
* Typically, its values lie between (0.5,0.9)
    
**lambda[default=0]**

* It controls L2 regularization (equivalent to Ridge regression) on weights. It is used to avoid overfitting.
    
**alpha[default=1]**

*It controls L1 regularization (equivalent to Lasso regression) on weights. In addition to shrinkage, enabling alpha also results in feature selection. Hence, it's more useful on high dimensional data sets.



**XGBoost-T with default parameters**

```{r}
start.time <- Sys.time() # Start the clock!
set.seed(20220620)

# Split the data
train.idx <- partition(csf_m$Type, p=c(train=0.8, test=0.2))
train <- csf_m[train.idx$train,]
test  <- csf_m[train.idx$test,]
  
### Make sure the dependent variable is numeric 0,1,2... (required by XGboost package)
train.X <- sparse.model.matrix(Type~.-1, data=train)
train.Y <- train$Type

test.X <- sparse.model.matrix(Type~.-1, data=test)
test.Y <- test$Type
  
## Convert the train and test data into xgb matrix type (required by XGboost package)
xgb_train = xgb.DMatrix(data=train.X, label=train.Y)
xgb_test = xgb.DMatrix(data=test.X, label=test.Y)
  
## Defining the parameters
xgb_params <- list("objective"="binary:logistic", ## logistic regression
                   "eval_metric"="logloss") ## multiclass logloss (used in classification)
watchlist <- list(train=xgb_train, test=xgb_test)

xgbc = xgb.train(booster="gbtree",
                 params=xgb_params,
                 data=xgb_train,
                 max_depth = 6,
                 eta = 0.3,
                 gamma=0,
                 nthread = 1,
                 nrounds=100, ## number of XGBoost iterations
                 watchlist=watchlist,
                 min_child_weight=1, ##
                 subsample=1, 
                 colsample_bytree=1, ##
                 verbose=F) ## Print or not

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
# print(xgbc)
```

**Check error rates**

Visually

```{r}
# Training & test error plot
error_xgbc <- data.frame(xgbc$evaluation_log)
plot(error_xgbc$iter, error_xgbc$train_logloss, col='blue', type="l", ylim=c(0,1))
lines(error_xgbc$iter, error_xgbc$test_logloss, col='red')
legend(x=60,y=0.4,c("Train","Test"), col=c('blue','red'), lty=1)
```


Using the Root mean square error(RMSE) 

```{r}
csf_m.xgbtree.test.preds <- as.numeric(predict(xgbc, xgb_test))
csf_m.xgbtree.test.labels <- test.Y
rmse.xgbtree <- sqrt(mean((csf_m.xgbtree.test.preds-csf_m.xgbtree.test.labels)^2)) 
rmse.xgbtree # RMSE
```



**Tune parameters**

There are several ways to tune the hyperparameters. We choose here to use grid search with the **caret** package.

The XGBoost model has many parameters. If all the parameters are optimized, the computer's memory will be solicited, and the optimization time will be significantly increased. We will optimize the following parameters: `nrounds`, `eta`, `min_child_weight`, and `colsample_bytree`.

**Note:* You can break down the tuning process into five steps:

* Fixing learning rate `eta` and number of iterations `nrounds`;
* Maximum depth `max_depth` and child weight `min_child_weight`;
* Setting column `colsample_bytree` and row sampling `subsample`;
* Experimenting with different `gamma` values;
* Reducing the learning rate `eta`. 

Each following step is fed by the best values of the previous step. See [here](https://www.kaggle.com/code/pelkoja/visual-xgboost-tuning-with-caret/report) for detailed information. 


This tunning process can be time consuming (here 2 min)

```{r, message=FALSE, warning=FALSE, results = "hide"}
set.seed(20220620)
start.time <- Sys.time() # Start the clock!

xgbtree_caret <- train(x=as.matrix(train.X) ,
                   y=as.factor(train.Y),
                   method="xgbTree",
                   trControl=trainControl(method="repeatedcv",
                                          number=3,
                                          repeats=2,
                                          verboseIter=FALSE),
                   tuneGrid=expand.grid(nrounds=c(50,100,500),
                                        eta=c(0.01,0.05,0.1),
                                        max_depth=6,
                                        colsample_bytree=c(0.5,1),
                                        subsample=c(0.5,1),
                                        gamma=0,
                                        min_child_weight=c(0.5,1)))

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```

The best parameters

```{r}
xgbtree_caret$bestTune 
```


Run the XGBoost with the tuned parameters

```{r}
set.seed(20220620)
xgbtree_tuned <- xgb.train(booster="gbtree",
                    params=xgb_params,
                    data=xgb_train,
                    max_depth=6,
                    eta=xgbtree_caret$bestTune$eta,
                    gamma=0,
                    nthread=1,
                    nrounds=xgbtree_caret$bestTune$nrounds, ## number of XGBoost iterations
                    watchlist=watchlist,
                    min_child_weight=xgbtree_caret$bestTune$min_child_weight, ##
                    subsample=xgbtree_caret$bestTune$subsample, 
                    colsample_bytree=xgbtree_caret$bestTune$colsample_bytree, ##
                    verbose=F) ## Print or not
print(xgbtree_tuned)
```

**Check error rates**

Visually

```{r}
# Training & test error plot
error_xgbtree_tuned <- data.frame(xgbtree_tuned$evaluation_log)
plot(error_xgbtree_tuned$iter, error_xgbtree_tuned$train_logloss, col='blue', type="l", ylim=c(0,1))
lines(error_xgbtree_tuned$iter, error_xgbtree_tuned$test_logloss, col='red')
legend(x=40,y=0.5,c("Train","Test"), col=c('blue','red'), lty=1)
```

Using the Root mean square error(RMSE) 

```{r}
csf_m.xgbtree.test.preds <- as.numeric(predict(xgbtree_tuned, xgb_test))
csf_m.xgbtree.test.labels <- test.Y
rmse.xgbtree_tuned <- sqrt(mean((csf_m.xgbtree.test.preds-csf_m.xgbtree.test.labels)^2))
rmse.xgbtree_tuned # RMSE
```


**Note**: The RMSE was 0.496 with the default parameters.



##### Final model performance

```{r}
set.seed(20220620) ## To replicate the results
nsim <- 100

csf_m.xgbtree.test.preds <- vector(mode="list", length=nsim)
csf_m.xgbtree.test.labels <- vector(mode="list", length=nsim)

## Extract x - feature and y - label parts.
## To use the XGBoost algorithm we need to create a matrix and use the one-hot encoding.
## This will create a dummy variable for the factor variables: DONE

for(j in 1:nsim) {
  train.idx <- partition(csf_m$Type, p=c(train=0.8, test=0.2))
  
  train <- csf_m[train.idx$train,]
  test  <- csf_m[train.idx$test,]
  
  ### Make sure the dependent variable is numeric 0,1,2... (required by XGboost package)
  train.X <- sparse.model.matrix(Type~.-1, data=train)
  train.Y <- train$Type
  
  test.X <- sparse.model.matrix(Type~.-1, data=test)
  test.Y <- test$Type
  
  ## Convert the train and test data into xgb matrix type (required by XGboost package)
  xgb_train = xgb.DMatrix(data=train.X, label=train.Y)
  xgb_test = xgb.DMatrix(data=test.X, label=test.Y)
  
  ## Defining the parameters
  xgb_params <- list("objective"="binary:logistic", ## logistic regression for binary classification
                     "eval_metric" = "logloss") ## multiclass logloss (used in classification)
  watchlist <- list(train=xgb_train, test=xgb_test)
  
  xgbc <- xgb.train(booster="gbtree",
                    params=xgb_params,
                    data=xgb_train,
                    max_depth=6,
                    eta=xgbtree_caret$bestTune$eta,
                    gamma=0,
                    nthread=1,
                    nrounds=xgbtree_caret$bestTune$nrounds, ## number of XGBoost iterations
                    watchlist=watchlist,
                    min_child_weight=xgbtree_caret$bestTune$min_child_weight, ##
                    subsample=xgbtree_caret$bestTune$subsample, 
                    colsample_bytree=xgbtree_caret$bestTune$colsample_bytree, ##
                    verbose=F) ## Print or not
  ### print(xgbc)

  ## The model is ready and we can predict the test data.
  ## Predicting test data (Hold-out)
  csf_m.xgbtree.test.preds[[j]] <- as.numeric(predict(xgbc, xgb_test))
  csf_m.xgbtree.test.labels[[j]] <- test.Y
}

csf_m.xgbtree.ROCR = ROCR::prediction(csf_m.xgbtree.test.preds, csf_m.xgbtree.test.labels)
csf_m.xgbtree.roc = performance(csf_m.xgbtree.ROCR,"tpr","fpr")
csf_m.xgbtree.auc = performance(csf_m.xgbtree.ROCR, measure="auc")
```


ROC curve and AUC

```{r}
#mar.default <- c(5,4,4,2) + 0.1   # plot margins
#par(mar = mar.default + c(0, 2, 0, 0)) 
par(mar = c(5,6,4,2) + 0.1) # (bottom, left, top, right) margins in numbers of lines
plot(csf_m.xgbtree.roc,col="grey50", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_m.xgbtree.roc, lwd=3, avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="XGBoost-T with All Metabolites")
text(x=0.8, y=0.1, paste("AUC=", round(median(unlist(csf_m.xgbtree.auc@y.values)),2)), cex=1.5) # AUC
abline(a=0, b=1)
par(new=FALSE)
```


```{r}
round(median(unlist(csf_m.xgbtree.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_m.xgbtree.auc@y.values)),2)/10  ## standard error
```



####  XGBoost-L
##### Determine optimal parameters to use with the XGBoost-L

Parameters for Linear Booster

Using linear booster has relatively lesser parameters to tune, hence it computes much faster than gbtree booster.

**nrounds[default=100]**

* It controls the maximum number of iterations (steps) required for gradient descent to converge.
Should be tuned using CV
  
**lambda[default=0]**

* It enables Ridge Regression. Same as above
	
**alpha[default=1]**

* It enables Lasso Regression. Same as above

**XGBoost-L with default parameters**


```{r}
set.seed(20220620)
start.time <- Sys.time() # Start the clock!
train.idx <- partition(csf_m$Type, p=c(train=0.8, test=0.2)) ## library(SpliTools)
  
train <- csf_m[train.idx$train,]
test  <- csf_m[train.idx$test,]

### Make sure the dependent variable is numeric 0,1,2... (required by XGboost package)
train.X <- as.matrix(train[,-1])
train.Y <- train$Type
  
test.X <- as.matrix(test[,-1])
test.Y <- test$Type
  
## Need to convert the train and test data into xgb matrix type (required by XGboost package)
xgb_train <- xgb.DMatrix(data=train.X, label=train.Y)
xgb_test <- xgb.DMatrix(data=test.X, label=test.Y)

## Defining the parameters
#numberOfClasses <- length(unique(csf_m$Type))
xgb_params <- list("objective"="binary:logistic",  ## logistic regression for binary classification
                   "eval_metric"="logloss") ## multiclass logloss (used in classification)
watchlist <- list(train=xgb_train, test=xgb_test)
  
xgbc <- xgb.train(data=xgb_train, booster="gblinear",
                  params=xgb_params,
                  nrounds=100,
                  verbose=0,
                  watchlist=watchlist)

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken) 
```


**Check error rates**

```{r}
csf_m.xgblin.test.preds <- as.numeric(predict(xgbc, xgb_test))
csf_m.xgblin.test.labels <- test.Y
rmse_xgbc <- sqrt(mean((csf_m.xgblin.test.preds-csf_m.xgblin.test.labels)^2))
rmse_xgbc # RMSE .
```


**Tune parameters**

```{r, message=FALSE, warning=FALSE}
set.seed(20220620)
start.time <- Sys.time() # Start the clock!

xgblin_caret <- train(x=as.matrix(train.X) ,
                   y=as.factor(train.Y),
                   method="xgbLinear",
                   trControl=trainControl(method="repeatedcv",
                                          number=5,
                                          repeats=10,
                                          verboseIter=FALSE),
                   tuneGrid=expand.grid(nrounds=c(1,2,5),
                                        lambda =0,
                                        alpha = 1,
                                        eta=0.3
                                          ))

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```
The best parameters

```{r}
xgblin_caret$bestTune 
```


```{r}
xgblin_tuned <- xgb.train(data=xgb_train, booster="gblinear",
                        params=xgb_params,
                        nrounds=xgblin_caret$bestTune$nrounds,
                        verbose=0,
                        watchlist=watchlist)
```

```{r}
csf_m.xgblin.test.preds_tuned <- as.numeric(predict(xgblin_tuned, xgb_test))
csf_m.xgblin.test.labels_tuned <- test.Y
rmse_xgblin_tuned <- sqrt(mean((csf_m.xgblin.test.preds_tuned-csf_m.xgblin.test.labels_tuned)^2)) 
rmse_xgblin_tuned # RMSE
```
The error rate with default parameters was 0.61


#####  Final model performance

```{r}
set.seed(20220620) ## To replicate the results
start.time <- Sys.time() # Start the clock!
nsim <- 100

csf_m.xgblin.test.preds <- vector(mode="list", length=nsim)
csf_m.xgblin.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  ## Extract x - feature and y - label parts.
  ## To use the XGBoost algorithm we need to create a matrix and use the one-hot encoding.
  ## This will create a dummy variable for the factor variables: DONE
  
  train.idx <- partition(csf_m$Type, p=c(train=0.8, test=0.2)) ## library(SpliTools)
  
  train <- csf_m[train.idx$train,]
  test  <- csf_m[train.idx$test,]
  
  train.X <- as.matrix(train[,-1])
  train.Y <- train$Type
  
  test.X <- as.matrix(test[,-1])
  test.Y <- test$Type
  
  ## Convert the train and test data into xgb matrix type (required by XGboost package)
  xgb_train <- xgb.DMatrix(data=train.X, label=train.Y)
  xgb_test <- xgb.DMatrix(data=test.X, label=test.Y)

  ## Defining the parameters
  xgb_params <- list("objective"="binary:logistic",  ## logistic regression
                     "eval_metric"="logloss") ## multiclass logloss (used in classification)
  watchlist <- list(train=xgb_train, test=xgb_test)
  
  xgbc <- xgb.train(data=xgb_train, booster="gblinear",
                    params=xgb_params,
                    nrounds=xgblin_caret$bestTune$nrounds,
                    verbose=0,
                    watchlist=watchlist)

  ## The model is ready and we can predict the test data.
  ## Predicting test data (Hold-out)
  csf_m.xgblin.test.preds[[j]] <- as.numeric(predict(xgbc, xgb_test))
  csf_m.xgblin.test.labels[[j]]<- test.Y
}

csf_m.xgblin.ROCR = ROCR::prediction(csf_m.xgblin.test.preds, csf_m.xgblin.test.labels)

csf_m.xgblin.roc = performance(csf_m.xgblin.ROCR,"tpr","fpr")
csf_m.xgblin.auc = performance(csf_m.xgblin.ROCR, measure="auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken) 
```



**ROC curve and AUC**

```{r}
#mar.default <- c(5,4,4,2) + 0.1
#par(mar = mar.default + c(0, 2, 0, 0)) 
par(mar = c(5,6,4,2) + 0.1) # (bottom, left, top, right) margins in numbers of lines
plot(csf_m.xgblin.roc, col="grey50", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_m.xgblin.roc,lwd=3,avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="XGBoost-L with All Metabolites")
abline(a=0, b=1) 
text(x=0.8, y=0.1,paste("AUC=", round(median(unlist(csf_m.xgblin.auc@y.values)),2)),cex = 1.5)  # Show AUC
par(new=FALSE)
```


```{r}
round(median(unlist(csf_m.xgblin.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_m.xgblin.auc@y.values)),2)/10  ## standard error
```




##  Predictive Models with Selected Metabolites

We can consider two different methods to select a pre-specified number of metabolites to build the predictive models for AD patients. GPS & SIS...

Note that a cross-validation method can be used to optimize the number of metabolites to select.



### Logistic regression with 35 Metabolites Selected by GPS

```{r, warning=FALSE, message=FALSE, results="hide"}
set.seed(20220620) ## To replicate the results
start.time <- Sys.time() # Start the clock!
nsim <- 100 ## Number of repetitions

csf_m35_GPS.logit.test.preds <- vector(mode="list", length=nsim)
csf_m35_GPS.logit.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  train.idx <- partition(csf_m$Type, p = c(train=0.8, test=0.2)) ## library(SpliTools)
  train <- csf_m[train.idx$train,]
  test  <- csf_m[train.idx$test,]
  train.X <- train[,-1]
  test.X <- test[,-1]
  
  m35_GPS <- gps(train$Type, train.X, family='binomial', maxvar=35, maxcmp=7)
  train_GPS <- data.frame(Type=train$Type, train.X[,m35_GPS])
  
  ## Fit LOGISTIC 
  csf_m35_GPS.logit.train <- glm(Type~., data=train_GPS, family=binomial(link="logit"))

  ## Predict
  csf_m35_GPS.logit.test.preds[[j]] <- predict(csf_m35_GPS.logit.train, test.X[,m35_GPS], type="response")
  csf_m35_GPS.logit.test.labels[[j]]<- test$Type
}

csf_m35_GPS.logit.ROCR <- ROCR::prediction(csf_m35_GPS.logit.test.preds, csf_m35_GPS.logit.test.labels)
csf_m35_GPS.logit.roc <- performance(csf_m35_GPS.logit.ROCR,"tpr","fpr")
csf_m35_GPS.logit.auc <- performance(csf_m35_GPS.logit.ROCR, measure = "auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```


Let us plot the ROC curve and get the AUC.

```{r}
#mar.default <- c(5,4,4,2) + 0.1
#par(mar = mar.default + c(0, 2, 0, 0)) 
par(mar = c(5,6,4,2) + 0.1) # (bottom, left, top, right) margins in numbers of lines
plot(csf_m35_GPS.logit.roc,col="grey50",lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_m35_GPS.logit.roc,lwd=3,avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="Logistic Regression with 35 Metabolites Selected by GPS")
text(x=0.8, y=0.1,paste("AUC=", round(median(unlist(csf_m35_GPS.logit.auc@y.values)),2)),cex = 1.5)  # AUC
abline(a=0, b=1) 
par(new=FALSE)
```

While the above plot shows with AUC=0.67, we can further confirm this and calculate its standard error.

```{r}
round(median(unlist(csf_m35_GPS.logit.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_m35_GPS.logit.auc@y.values)),2)/10  ## standard error
```


### Logistic regression with 35 Metabolites Selected by SIS


```{r, warning=FALSE}
set.seed(20220620) ## To replicate the results
start.time <- Sys.time() # Start the clock!
nsim <- 100 ## Number of repetitions

csf_m35_SIS.logit.test.preds <- vector(mode="list", length=nsim)
csf_m35_SIS.logit.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  train.idx <- partition(csf_m$Type, p=c(train=0.8, test=0.2)) ## library(SpliTools)
  train <- csf_m[train.idx$train,]
  test  <- csf_m[train.idx$test,]
  train.X <- train[,-1]
  test.X <- test[,-1]
  
  m35_SIS <- SIS(as.matrix(train.X), train$Type, family='binomial', iter=FALSE, tune='cv', nsis=35)$ix0
  train_SIS <- data.frame(Type=train$Type,train.X[,m35_SIS])

  ## Fit LOGISTIC 
  csf_m35_SIS.logit.train <- glm(Type~., data=train_SIS, family=binomial(link="logit"))

  ## Predict
  csf_m35_SIS.logit.test.preds[[j]] <- predict(csf_m35_SIS.logit.train, test.X[,m35_SIS], type="response")
  csf_m35_SIS.logit.test.labels[[j]]<- test$Type
}

csf_m35_SIS.logit.ROCR <- ROCR::prediction(csf_m35_SIS.logit.test.preds, csf_m35_SIS.logit.test.labels)
csf_m35_SIS.logit.roc <- performance(csf_m35_SIS.logit.ROCR,"tpr","fpr")
csf_m35_SIS.logit.auc <- performance(csf_m35_SIS.logit.ROCR, measure = "auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)

```


Let us plot the ROC curve and get the AUC.

```{r}
#mar.default <- c(5,4,4,2) + 0.1
#par(mar = mar.default + c(0, 2, 0, 0)) 
par(mar = c(5,6,4,2) + 0.1) # (bottom, left, top, right) margins in numbers of lines
plot(csf_m35_SIS.logit.roc,col="grey50",lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_m35_SIS.logit.roc,lwd=3,avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="Logistic Regression with 35 Metabolites Selected by SIS")
text(x=0.8, y=0.1,paste("AUC=", round(median(unlist(csf_m35_SIS.logit.auc@y.values)),2)),cex = 1.5)  # AUC
abline(a=0, b=1) 
par(new=FALSE)
```

While the above plot shows with AUC=0.67, we can further confirm this and calculate its standard error.

```{r}
round(median(unlist(csf_m35_SIS.logit.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_m35_SIS.logit.auc@y.values)),2)/10  ## standard error
```



###  PLS-DA with 35 Metabolites Selected by GPS

**Tune the number of components for the PLS-DA**

```{r}
set.seed(20220620) # for reproducibility, only when the `cpus' argument is not used
start.time <- Sys.time() # Start the clock!

train.idx <- partition(csf_m$Type, p = c(train=0.8, test=0.2)) ## library(SpliTools)
  
train <- csf_m[train.idx$train,]
test  <- csf_m[train.idx$test,]
  
train.X <- as.matrix(train[,-1])
train.Y<- train$Type
  
m35_GPS <- gps(train.Y, train.X, family='binomial', maxvar=35, maxcmp=3)
m35_GPS.plsda <- plsda(train.X[, m35_GPS], train.Y, ncomp=10)  # set ncomp to 10

m35_GPS.plsda.perf <- perf(m35_GPS.plsda, validation="Mfold", folds=5,
                         progressBar=FALSE, auc=TRUE, nrepeat=50)
                         
ncomps.m35_GPS.plsda <- m35_GPS.plsda.perf$choice.ncomp
ncomps.m35_GPS.plsda

m35_GPS.plsda.perf$error.rate

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken) 
```


**Final model**

```{r, message=FALSE, warning=FALSE, results="hide"}
ncomps.csf_m35.plsda <- 1 ## Number of component for PLS-DA 
set.seed(20220620) ## To replicate the results
start.time <- Sys.time() # Start the clock!
nsim <- 100 ## Number of repetitions

csf_m35_GPS.plsda.test.preds <- vector(mode="list", length=nsim)
csf_m35_GPS.plsda.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  train.idx <- partition(csf_m$Type, p = c(train=0.8, test=0.2)) ## library(SpliTools)
  
  train <- csf_m[train.idx$train,]
  test  <- csf_m[train.idx$test,]
  
  train.X <- as.matrix(train[,-1])
  train.Y<- train$Type
  
  test.X <- as.matrix(test[,-1])
  test.Y <- test$Type
  
  ## Select 35 Metabolites by GPS
  #m35_GPS <- gps(train.Y, train.X, family='binomial', maxvar=35, maxcmp=4)
  m35_GPS <- gps(train.Y, train.X, family='binomial', maxvar=35, maxcmp=5)
  
  ## PLS-DA with 35 Metabolites Selected by GPS
  csf_m35_GPS.plsda.train <- plsda(train.X[,m35_GPS], train.Y, ncomp=ncomps.csf_m35.plsda)
  
  ## Predict
  csf_m35_GPS.plsda.test.preds[[j]] <- predict(csf_m35_GPS.plsda.train, test.X[,m35_GPS],
                                               dist="mahalanobis.dist")$predict[,,ncomps.csf_m35.plsda][,2]
  csf_m35_GPS.plsda.test.labels[[j]]<- test.Y
}

csf_m35_GPS.plsda.ROCR = ROCR::prediction(csf_m35_GPS.plsda.test.preds, csf_m35_GPS.plsda.test.labels)
csf_m35_GPS.plsda.roc = performance(csf_m35_GPS.plsda.ROCR,"tpr","fpr")
csf_m35_GPS.plsda.auc = performance(csf_m35_GPS.plsda.ROCR, measure="auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```




```{r}
## ROC curve and AUC
#mar.default <- c(5,4,4,2) + 0.1
#par(mar = mar.default + c(0, 2, 0, 0)) 
par(mar = c(5,6,4,2) + 0.1) # (bottom, left, top, right) margins in numbers of lines
plot(csf_m35_GPS.plsda.roc,col="grey50", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_m35_GPS.plsda.roc, lwd=3, avg="vertical", cex.lab=1.5, cex.main=1.5, xlab="", ylab="",
     main="PLS-DA with 35 Metabolites Selected by GPS")
text(x=0.8, y=0.1,paste("AUC=", round(median(unlist(csf_m35_GPS.plsda.auc@y.values)),2)),cex=1.5) # Show AUC
abline(a=0, b=1) 
par(new=FALSE)
```



```{r}
round(median(unlist(csf_m35_GPS.plsda.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_m35_GPS.plsda.auc@y.values)),2)/10  ## standard error
```



### PLS-DA with 35 Metabolites Selected by SIS

Create a training data with 35 metabolites to tune the number of components for the PLS-DA

```{r}
set.seed(20220620) # for reproducibility, only when the `cpus' argument is not used
start.time <- Sys.time() # Start the clock!

train.idx <- partition(csf_m$Type, p = c(train=0.8, test=0.2)) ## library(SpliTools)
  
train <- csf_m[train.idx$train,]
test  <- csf_m[train.idx$test,]
  
train.X <- as.matrix(train[,-1])
train.Y<- train$Type
  
m35_SIS <- SIS(train.X, train.Y, family='binomial', iter=FALSE, tune='cv', nsis=35)$ix0
train_SIS <- data.frame(Type=train$Type, train.X[, m35_SIS])
  
m35_SIS.plsda <- plsda(train_SIS[,-1], train_SIS$Type, ncomp=10)  # set ncomp to 10

m35_SIS.plsda.perf <- perf(m35_SIS.plsda, validation="Mfold", folds=5,
                         progressBar=FALSE, auc=TRUE, nrepeat=50)
                         
ncomps.m35_SIS.plsda <- m35_SIS.plsda.perf$choice.ncomp
ncomps.m35_SIS.plsda

m35_SIS.plsda.perf$error.rate

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken) 
```


Run the model

```{r, message=FALSE, warning=FALSE}
set.seed(20220620) ## To replicate the results
ncomps.csf_m35.plsda <- 3 ## Number of component for PLS-DA 
start.time <- Sys.time() # Start the clock!
nsim <- 100 ## Number of repetitions

csf_m35_SIS.plsda.test.preds <- vector(mode="list", length=nsim)
csf_m35_SIS.plsda.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  train.idx <- partition(csf_m$Type, p=c(train=0.8, test=0.2)) ## library(SpliTools)
  
  train <- csf_m[train.idx$train,]
  test  <- csf_m[train.idx$test,]
  
  train.X <- as.matrix(train[,-1])
  train.Y<- train$Type
  
  test.X <- as.matrix(test[,-1])
  test.Y <- test$Type
  
  # Use SIS to select variables
  #m35_SIS <- SIS(train.X, train.Y, family='binomial', iter=FALSE, tune='cv', nsis=35)$ix0
  m35_SIS <- SIS(train.X, train.Y, family='binomial', iter=FALSE, tune='cv', nsis=35)$ix0

  ## PLS-DA
  csf_m35_SIS.plsda.train <- plsda(train.X[,m35_SIS], train.Y, ncomp=ncomps.csf_m35.plsda)
  
  ## Predict
  csf_m35_SIS.plsda.test.preds[[j]] <- predict(csf_m35_SIS.plsda.train, test.X[,m35_SIS],
                                               dist="centroids.dist")$predict[,,ncomps.csf_m35.plsda][,2]
  csf_m35_SIS.plsda.test.labels[[j]]<- test.Y
}

csf_m35_SIS.plsda.ROCR <- ROCR::prediction(csf_m35_SIS.plsda.test.preds, csf_m35_SIS.plsda.test.labels)
csf_m35_SIS.plsda.roc <- performance(csf_m35_SIS.plsda.ROCR,"tpr","fpr")
csf_m35_SIS.plsda.auc <- performance(csf_m35_SIS.plsda.ROCR, measure="auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```



```{r}
## ROC curve and AUC
#mar.default <- c(5,4,4,2) + 0.1
#par(mar = mar.default + c(0, 2, 0, 0)) 
par(mar = c(5,6,4,2) + 0.1) # (bottom, left, top, right) margins in numbers of lines
plot(csf_m35_SIS.plsda.roc, col="grey50", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_m35_SIS.plsda.roc, lwd=3, avg="vertical", xlab="", ylab="",
      cex.lab=1.5, cex.main=1.5, main="PLS-DA with 35 Metabolites Selected by SIS")
text(x=0.8, y=0.1,paste("AUC=", round(median(unlist(csf_m35_SIS.plsda.auc@y.values)),2)), cex=1.5)  # Show AUC
abline(a=0, b=1) 
par(new=FALSE)
```

```{r}
round(median(unlist(csf_m35_SIS.plsda.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_m35_SIS.plsda.auc@y.values)),2)/10  ## standard error
```



###  XGBoost-T with 35 Metabolites Selected by GPS

** Tune the parameters for the XGBoost Tree **

```{r, message=FALSE, warning=FALSE, results="hide"}

set.seed(20220620) # for reproducibility, only when the `cpus' argument is not used
start.time <- Sys.time() # Start the clock!

train.idx <- partition(csf_m$Type, p = c(train=0.8, test=0.2)) ## library(SpliTools)
  
train <- csf_m[train.idx$train,]
test <- csf_m[train.idx$test,]
  
train.X <- as.matrix(train[,-1])
train.Y <- train$Type
  
m35_GPS <- gps(train.Y, train.X, family='binomial', maxvar=35, maxcmp=5)

xgbtree_caret <- train(x=as.matrix(train.X[, m35_GPS]) ,
                   y=as.factor(train.Y),
                   method="xgbTree",
                   trControl=trainControl(method="repeatedcv",
                                          number=3,
                                          repeats=2,
                                          verboseIter=FALSE),
                   tuneGrid=expand.grid(nrounds=c(50,100,500),
                                        eta=c(0.01,0.05,0.1),
                                        max_depth=6,
                                        colsample_bytree=c(0.25, 0.5),
                                        subsample=c(0.5,2),
                                        gamma=0,
                                        min_child_weight=c(1,2)))
  
end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken) 

xgbtree_caret$bestTune ## Best parameters 
```



```{r, message=FALSE, warning=FALSE, results="hide"}
set.seed(20220620) # for reproducibility here
start.time <- Sys.time() # Start the clock!
nsim <- 100    # number of simulations, could vary

csf_m35_GPS.xgbtree.test.preds <- vector(mode="list", length=nsim)
csf_m35_GPS.xgbtree.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  # Random Split into training and test data sets
  train.idx <- partition(csf_m$Type, p=c(train=0.8, test=0.2)) ## library(SpliTools)
  
  train <- csf_m[train.idx$train,]
  test <- csf_m[train.idx$test,]
  
  train.X <- as.matrix(train[,-1])
  train.Y <- train$Type
  
  test.X <- as.matrix(test[,-1])
  test.Y <- test$Type
  
  ## Select 35 Metabolites by GPS
  m35_GPS <- gps(train.Y, train.X, family='binomial', maxvar=35, maxcmp=3)

  ## Convert the train and test data into xgb matrix type (required by XGboost package)
  xgb_train <- xgb.DMatrix(data=train.X[,m35_GPS], label=train.Y)
  xgb_test <- xgb.DMatrix(data=test.X[,m35_GPS], label=test.Y)
  
  ## Defining the parameters
  xgb_params <- list("objective"="binary:logistic",   ## logistic regression
                     "eval_metric"="logloss") ## multiclass logloss (used in classification)
  watchlist <- list(train=xgb_train, test=xgb_test)
  
  xgbc <- xgb.train(booster="gbtree",
                    params=xgb_params,
                    data=xgb_train,
                    max_depth=6,
                    eta=xgbtree_caret$bestTune$eta,
                    gamma=0,
                    nthread=1,
                    nrounds=xgbtree_caret$bestTune$nrounds, ## number of XGBoost iterations
                    watchlist=watchlist,
                    min_child_weight=xgbtree_caret$bestTune$min_child_weight, ##
                    subsample=xgbtree_caret$bestTune$subsample, 
                    colsample_bytree=xgbtree_caret$bestTune$colsample_bytree, ##
                    verbose=F) ## Print or not

  ### print(xgbc)
  
  ## The model is ready and we can predict the test data.
  ## Predicting test data (Hold-out)
  csf_m35_GPS.xgbtree.test.preds[[j]] <- as.numeric(predict(xgbc, xgb_test))
  csf_m35_GPS.xgbtree.test.labels[[j]] <- test.Y
}

csf_m35_GPS.xgbtree.ROCR <- ROCR::prediction(csf_m35_GPS.xgbtree.test.preds, csf_m35_GPS.xgbtree.test.labels)
csf_m35_GPS.xgbtree.roc <- performance(csf_m35_GPS.xgbtree.ROCR,"tpr","fpr")
csf_m35_GPS.xgbtree.auc <- performance(csf_m35_GPS.xgbtree.ROCR, measure = "auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```



```{r}
#mar.default <- c(5,4,4,2) + 0.1
#par(mar = mar.default + c(0, 2, 0, 0)) 
par(mar = c(5,6,4,2) + 0.1) # (bottom, left, top, right) margins in numbers of lines
plot(csf_m35_GPS.xgbtree.roc, col="grey", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_m35_GPS.xgbtree.roc, lwd=3, avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="XGBoost-T with 35 Metabolites Selected by GPS")
text(x=0.8, y=0.1,paste("AUC=", round(median(unlist(csf_m35_GPS.xgbtree.auc@y.values)),2)), cex=1.5) # AUC
abline(a=0, b=1) 
par(new=FALSE)
```

```{r}
round(median(unlist(csf_m35_GPS.xgbtree.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_m35_GPS.xgbtree.auc@y.values)),2)/10  ## standard error
```



###  XGBoost-T with 35 Metabolites Selected by SIS

**Tune parameters **

```{r, message=FALSE, warning=FALSE, results="hide"}
set.seed(20220620) # for reproducibility, only when the `cpus' argument is not used
start.time <- Sys.time() # Start the clock!

train.idx <- partition(csf_m$Type, p = c(train=0.8, test=0.2)) ## library(SpliTools)
  
train <- csf_m[train.idx$train,]
test  <- csf_m[train.idx$test,]
  
train.X <- as.matrix(train[,-1])
train.Y <- as.numeric(as.character(train$Type))
  
m35_SIS <- SIS(train.X, train.Y, family='binomial', iter=FALSE, tune='cv', nsis=35)$ix0
  
xgbtree_caret <- train(x=as.matrix(train.X[, m35_SIS]) ,
                   y=as.factor(train.Y),
                   method="xgbTree",
                   trControl=trainControl(method="repeatedcv",
                                          number=3,
                                          repeats=2,
                                          verboseIter=FALSE),
                   tuneGrid=expand.grid(nrounds=c(50,100,500),
                                        eta=c(0.01,0.05,0.1),
                                        max_depth=6,
                                        colsample_bytree=c(0.5,1),
                                        subsample=c(0.5,1),
                                        gamma=0,
                                        min_child_weight=c(1,2)))
  
end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken) 

xgbtree_caret$bestTune ## Best parameters 
```


```{r, message=FALSE, warning=FALSE}
set.seed(20220620) # for reproducibility here
start.time <- Sys.time() # Start the clock!
nsim <- 100    # number of simulations, could vary

csf_m35_SIS.xgbtree.test.preds <- vector(mode="list", length=nsim)
csf_m35_SIS.xgbtree.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  # Randomly split into training & test data sets
  train.idx <- partition(csf_m$Type, p=c(train=0.8, test=0.2)) ## library(SpliTools)

  train <- csf_m[train.idx$train,]
  test <- csf_m[train.idx$test,]
  
  train.X <- as.matrix(train[,-1])
  train.Y <- train$Type
  
  test.X = as.matrix(test[,-1])
  test.Y <- test$Type
  
  # Use SIS to select variables
  m35_SIS <- SIS(train.X, train.Y, family='binomial', iter=FALSE, tune='cv', nsis=35)$ix0

  ## Convert the train and test data into xgb matrix type (required by XGboost package)
  xgb_train <- xgb.DMatrix(data=train.X[,m35_SIS], label=train.Y)
  xgb_test <- xgb.DMatrix(data=test.X[,m35_SIS], label=test.Y)
  
  ## Defining the parameters
  xgb_params <- list("objective"="binary:logistic",  ## logistic regression
                     "eval_metric" = "logloss") ## multiclass logloss (used in classification)
  watchlist <- list(train=xgb_train, test=xgb_test)

  xgbc <- xgb.train(booster="gbtree",
                    params=xgb_params,
                    data=xgb_train,
                    max_depth=6,
                    eta=xgbtree_caret$bestTune$eta,
                    gamma=0,
                    nthread=1,
                    nrounds=xgbtree_caret$bestTune$nrounds, ## number of XGBoost iterations
                    watchlist=watchlist,
                    min_child_weight=xgbtree_caret$bestTune$min_child_weight, ##
                    subsample=xgbtree_caret$bestTune$subsample, 
                    colsample_bytree=xgbtree_caret$bestTune$colsample_bytree, ##
                    verbose=F) ## Print or not
 
  ## print(xgbc)
  
  ## The model is ready and we can predict the test data.
  ## Predicting test data (Hold-out)
  csf_m35_SIS.xgbtree.test.preds[[j]] <- as.numeric(predict(xgbc, xgb_test))
  csf_m35_SIS.xgbtree.test.labels[[j]] <- test.Y
}

csf_m35_SIS.xgbtree.ROCR <- ROCR::prediction(csf_m35_SIS.xgbtree.test.preds, csf_m35_SIS.xgbtree.test.labels)
csf_m35_SIS.xgbtree.roc <- performance(csf_m35_SIS.xgbtree.ROCR,"tpr","fpr")
csf_m35_SIS.xgbtree.auc <- performance(csf_m35_SIS.xgbtree.ROCR, measure="auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```




```{r}
#mar.default <- c(5,4,4,2) + 0.1
#par(mar = mar.default + c(0, 2, 0, 0)) 
par(mar = c(5,6,4,2) + 0.1) # (bottom, left, top, right) margins in numbers of lines
plot(csf_m35_SIS.xgbtree.roc, col="grey", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_m35_SIS.xgbtree.roc, lwd=3, avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="XGBoost-T with 35 Metabolites Selected by SIS")
text(x=0.8, y=0.1,paste("AUC=", round(median(unlist(csf_m35_SIS.xgbtree.auc@y.values)),2)), cex=1.5) # AUC
abline(a=0, b=1) 
par(new=FALSE)
```


```{r}
round(median(unlist(csf_m35_SIS.xgbtree.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_m35_SIS.xgbtree.auc@y.values)),2)/10  ## standard error

```




### XGBoost-L with 35 Metabolites Selected by GPS


**Tune parameters**

```{r, message=FALSE, warning=FALSE}
set.seed(20220620)
start.time <- Sys.time() # Start the clock!

train.idx <- partition(csf_m$Type, p=c(train=0.8, test=0.2)) ## library(SpliTools)
  
train <- csf_m[train.idx$train,]
test <- csf_m[train.idx$test,]
  
train.X <- as.matrix(train[,-1])
train.Y <- train$Type
  
m35_GPS <- gps(train.Y, train.X, family='binomial', maxvar=35, maxcmp=5)
xgblin_caret <- train(x=as.matrix(train.X[,m35_GPS]) ,
                   y=as.factor(train.Y),
                   method="xgbLinear",
                   trControl=trainControl(method="repeatedcv",
                                          number=5,
                                          repeats=10,
                                          verboseIter=FALSE),
                   tuneGrid=expand.grid(nrounds=c(1,2,5),
                                        lambda =0,
                                        alpha = 1,
                                        eta=0.3))

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```

The best parameters

```{r}
xgblin_caret$bestTune 

```

<!---
```{r}
## Convert the train and test data into xgb matrix type (required by XGboost package)
xgb_train <- xgb.DMatrix(data=train.X[,m35_GPS], label=train.Y)
xgblin_tuned <- xgb.train(data=xgb_train, booster="gblinear",
                        params=xgb_params,
                        nrounds=xgblin_caret$bestTune$nrounds,
                        verbose=0,
                        watchlist=watchlist)
```
--->



```{r, message=FALSE, warning=FALSE, results="hide"}
set.seed(20220620) # for reproducibility here
start.time <- Sys.time() # Start the clock!
nsim <- 100    # number of simulations, could vary

csf_m35_GPS.xgblin.test.preds <- vector(mode="list", length=nsim)
csf_m35_GPS.xgblin.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  # Randomly split into training & test data sets
  train.idx <- partition(csf_m$Type, p=c(train=0.8, test=0.2)) ## library(SpliTools)
  
  train <- csf_m[train.idx$train,]
  test <- csf_m[train.idx$test,]
  
  train.X <- as.matrix(train[,-1])
  train.Y <- train$Type
  
  test.X <- as.matrix(test[,-1])
  test.Y <- test$Type
  
  ## Select 35 Metabolites by GPS
  m35_GPS <- gps(train.Y, train.X, family='binomial', maxvar=35, maxcmp=3)

  ## Convert the train and test data into xgb matrix type (required by XGboost package)
  xgb_train <- xgb.DMatrix(data=train.X[,m35_GPS], label=train.Y)
  xgb_test <- xgb.DMatrix(data=test.X[,m35_GPS], label=test.Y)
  
  ## Defining the parameters
  xgb_params <- list("objective"="binary:logistic",  ## logistic regression
                     "eval_metric" = "logloss") ## multiclass logloss (used in classification)
  watchlist <- list(train=xgb_train, test=xgb_test)

  xgbc <- xgb.train(data=xgb_train, booster="gblinear",
                    params=xgb_params,
                    nrounds=xgblin_caret$bestTune$nrounds ,
                    verbose=0,
                    watchlist=watchlist)
  ### print(xgbc)
  
  ## The model is ready and we can predict the test data.
  ## Predicting test data (Hold-out)
  csf_m35_GPS.xgblin.test.preds[[j]] <- as.numeric(predict(xgbc, xgb_test))
  csf_m35_GPS.xgblin.test.labels[[j]] <- test.Y
  
}

csf_m35_GPS.xgblin.ROCR <- ROCR::prediction(csf_m35_GPS.xgblin.test.preds, csf_m35_GPS.xgblin.test.labels)
csf_m35_GPS.xgblin.roc <- performance(csf_m35_GPS.xgblin.ROCR,"tpr","fpr")
csf_m35_GPS.xgblin.auc <- performance(csf_m35_GPS.xgblin.ROCR, measure="auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```


```{r}
#mar.default <- c(5,4,4,2) + 0.1
#par(mar = mar.default + c(0, 2, 0, 0)) 
par(mar = c(5,6,4,2) + 0.1) # (bottom, left, top, right) margins in numbers of lines
plot(csf_m35_GPS.xgblin.roc, col="grey", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_m35_GPS.xgblin.roc, lwd=3, avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="XGBoost-L with 35 Metabolites Selected by GPS")
text(x=0.8, y=0.1, paste("AUC=", round(median(unlist(csf_m35_GPS.xgblin.auc@y.values)),2)), cex=1.5) # AUC
abline(a=0, b=1) 
par(new=FALSE)
```


```{r}
round(median(unlist(csf_m35_GPS.xgblin.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_m35_GPS.xgblin.auc@y.values)),2)/10  ## standard error
```




### XGBoost-L with 35 Metabolites Selected by SIS

**Tune parameters**

```{r, message=FALSE, warning=FALSE}
set.seed(20220620)
start.time <- Sys.time() # Start the clock!

train.idx <- partition(csf_m$Type, p = c(train=0.8, test=0.2)) ## library(SpliTools)
  
train <- csf_m[train.idx$train,]
test <- csf_m[train.idx$test,]
  
train.X <- as.matrix(train[,-1])
train.Y <- train$Type
  
m35_SIS <- SIS(train.X, train.Y, family='binomial', iter=FALSE, tune='cv', nsis=35)$ix0

xgblin_caret <- train(x=as.matrix(train.X) ,
                   y=as.factor(train.Y),
                   method="xgbLinear",
                   trControl=trainControl(method="repeatedcv",
                                          number=5,
                                          repeats=10,
                                          verboseIter=FALSE),
                   tuneGrid=expand.grid(nrounds=c(1,2,5),
                                        lambda =0,
                                        alpha = 1,
                                        eta=0.3))

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```

The best parameters

```{r}
xgblin_caret$bestTune 
```


**Final model**

```{r, message=FALSE, warning=FALSE}
set.seed(20220620) # for reproducibility here
start.time <- Sys.time() # Start the clock!
nsim <- 100    # number of simulations, could vary

csf_m35_SIS.xgblin.test.preds <- vector(mode="list", length=nsim)
csf_m35_SIS.xgblin.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  # Randomly split into training & test data sets
  train.idx <- partition(csf_m$Type, p=c(train=0.8, test=0.2)) ## library(SpliTools)
  
  train <- csf_m[train.idx$train,]
  test <- csf_m[train.idx$test,]
  
  train.X <- as.matrix(train[,-1])
  train.Y <- train$Type
  
  test.X <- as.matrix(test[,-1])
  test.Y <- test$Type
  
  # Use SIS to select variables
  m35_SIS <- SIS(train.X, train.Y, family='binomial', iter=FALSE, tune='cv', nsis=35)$ix0

  ## Convert the train and test data into xgb matrix type (required by XGboost package)
  xgb_train <- xgb.DMatrix(data=train.X[,m35_SIS], label=train.Y)
  xgb_test <- xgb.DMatrix(data=test.X[,m35_SIS], label=test.Y)
  
  ## Defining the parameters
  xgb_params <- list("objective"="binary:logistic",  ## logistic regression
                     "eval_metric" = "logloss") ## multiclass logloss (used in classification)
  watchlist <- list(train=xgb_train, test=xgb_test)
  
  xgbc <- xgb.train(data=xgb_train, booster="gblinear",
                    params=xgb_params,
                    nrounds=xgblin_caret$bestTune$nrounds ,
                    verbose = 0,
                    watchlist=watchlist)
  ### print(xgbc)
  
  ## The model is ready and we can predict the test data.
  ## Predicting test data (Hold-out)
  csf_m35_SIS.xgblin.test.preds[[j]] <- as.numeric(predict(xgbc, xgb_test))
  csf_m35_SIS.xgblin.test.labels[[j]] <- test.Y
}

csf_m35_SIS.xgblin.ROCR = ROCR::prediction(csf_m35_SIS.xgblin.test.preds, csf_m35_SIS.xgblin.test.labels)
csf_m35_SIS.xgblin.roc = performance(csf_m35_SIS.xgblin.ROCR,"tpr","fpr")
csf_m35_SIS.xgblin.auc = performance(csf_m35_SIS.xgblin.ROCR, measure="auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```



```{r}
#mar.default <- c(5,4,4,2) + 0.1
#par(mar = mar.default + c(0, 2, 0, 0)) 
par(mar = c(5,6,4,2) + 0.1) # (bottom, left, top, right) margins in numbers of lines
plot(csf_m35_SIS.xgblin.roc, col="grey", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_m35_SIS.xgblin.roc, lwd=3, avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="XGBoost-L with 35 Metabolites Selected by SIS")
text(x = 0.8, y = 0.1,paste("AUC=", round(median(unlist(csf_m35_SIS.xgblin.auc@y.values)),2)),cex = 1.5)  # Add text element
#mtext("B", side=3, line=1, cex=1.5, adj=-0.05)
abline(a=0, b=1) 
par(new=FALSE)
```



```{r}
round(median(unlist(csf_m35_SIS.xgblin.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_m35_SIS.xgblin.auc@y.values)),2)/10  ## standard error
```



* AUC summary table for the **csf_m** study.

```{r}
Summary_AUC <- data.frame(Variable =c('All','35 with GPS','35 with SIS'),
               Logistic=c(0.64, 0.65, 0.65), PLSDA=c(0.70, 0.73, 0.71),
               XGBoostT=c(0.71, 0.76, 0.74), XGBoostL=c(0.75, 0.75, 0.74))
Summary_AUC
```

