bio2019_script.Rmd

---
title: "Invent the Future"
output: html_notebook
---

Authors: Raquel Aoki and Alice kang Yue


# Data

We will work with RNAseq data to try predict the gestational age in weeks. 

## Input Files:

- features: (samples x 32830/925032 gene/probeset) RNAseq counts (+extracted features); data has been batch and count normalized
- Value: (367 train sample) gestational age 8-42 weeks

## Output Files:
- Value: (368 test sample) gestational age 8-42 weeks rounded to 1 decimal place

# Data Analysis 

##1. Preprocessing  

In the first part we will work on the environment we will use. We will clean the environment, declare the project directory and install/load libraries.  


```{r, warning=FALSE, message=FALSE, error = FALSE  , results = 'hide' }
rm(list=ls(all=T)) # clean the environment
set.seed(10)


## Defining the work directory of the project
root = '~/GitHub/2019_bio'
#root = '~/Documents/GitHub/2019_bio'
setwd(root)


## Creating sub directories to save data, features, model and results 
input_dir = paste0(root,"/00_input") # raw data directory
feat_dir = paste0(root,"/01_features") # feature directory
model_dir = paste0(root, "/02_models") # model directory
result_dir = paste0(root, "/03_results") # stats/plots directory
sapply(c(input_dir,feat_dir, model_dir, result_dir), 
       function(x) dir.create(x, showWarnings=F))

## load packages; need to fix according to what model we'll be using
pkgs = c("Rfast", "stringr", "plyr", "dplyr", "Matrix", # var, str_, llply, etc
         "lattice",  # barplot, plots
         'ggplot','gridExtra','grid',
         "foreach", "doMC", # parallel back-end
         "caret", "e1071", "ranger", "ANN2", "randomForest",
         "elasticnet", "fastICA", "foba", "glmnet","kernlab", 
         "KRLS", "lars", "leaps", "nnls", "nodeHarvest", 
         "partDSA", "pls", "plsRglm", "rpart", "rqPen",
         "RSNNS", "spikeslab", "xgboost", "Metrics",'caretEnsemble') # ml

pkgs_ui = setdiff(pkgs, rownames(installed.packages()))
if (length(pkgs_ui) > 0) install.packages(pkgs_ui, verbose=F)
sapply(pkgs, require, character.only=T)

## script options
#no_cores = detectCores()-1 # number of cores to use in parallel
#registerDoMC(no_cores)
overwrite = F 


#load workspace for time consuming tasks 
load(paste(input_dir,'script.RData',sep='/'))

```

Next, we will load the 3 input files:

1. Sample Annotation file: 

* SampleID: unique identifier of the sample (matching the name of the .CEL file in HTA20 folder, except for extension .CEL);
* GA: gestational age as determined by the last menstrual period and or ultrasound;
* Batch: the batch identifier;
* Set: name of the source dataset;
* Train: 1 for samples to be used for training, 0 for samples to be used for test;

2. RNASEQ Data: each row is a gene and each column is a sample 

*  probeset: gene ID

3. Submission template for the competition

```{r,warning=FALSE, message=FALSE, error = FALSE}
## load input files

#Sample Annotation file
meta=  read.csv(paste0(input_dir,"/anoSC1_v11_nokey.csv"))

# RNASEQ data
data0 = t(get(load(paste0(input_dir,"/HTA20_RMA.RData"))))


# submission template
submission = read.csv(paste0(input_dir,"/TeamX_SC1_prediction.csv")) 

```


Data Exploration: 

```{r}
cat('Range\n');range(data0) ; cat('Genes IDs\n');gid = colnames(data0); head(gid) ;cat('Patients IDs\n');pid = rownames(data0); head(pid)
cat('Meta Data head and shape\n'); head(meta); dim(meta)
cat('RNASEQ head and shape\n'); head(data0[,c(1:5)]); dim(data0)
cat('Submission head and shape\n'); head(submission); dim(submission)

```

Now we will explore the data with some plots. First we need to prepare the data and summarize the informatoin we will use in the plots. 


```{r}
# plot stats: mean count, pearson/spearman corr

#Using only the train dataset to calculate the correlation, mean and variance 
data1 = data.frame('SampleID' = rownames(data0), data0)
data1 = merge(meta[,c(1,2,5)],data1,by.x = 'SampleID',by.y = 'SampleID', all = T)

#splitting training set and testing set 
data1 = subset(data1, Train == 1) 
data1 = subset(data1, select = -c(Train,SampleID))

#train set 
data_cor = data.frame(col = names(data1),
                      corr_p = apply(data1, 2, cor, x = data1$GA, method = 'pearson'))
data_cor = data.frame(data_cor, corr_s = apply(data1, 2, cor, x = data1$GA, method = 'spearman'))
rownames(data_cor) = NULL 


#Dataset with all information 
data_cor = data.frame(data_cor, variance = c(var(data1$GA),colVars(data0)), mean = c(mean(data1$GA),colMeans(data0))) 
data_cor = data_cor[-1,] # removing GA from the dataset



```

To make the plots, we will use a library called ggplot2. Here are shown some plots we can make using this library. 

```{r}
# plot stats
p1 <- ggplot(data_cor, aes(x=mean)) + 
  geom_histogram(fill = 'lightgreen') + 
  xlab('Average Expression') + labs(title='(a)')

p2 <- ggplot(data_cor, aes(x=mean, y = variance)) + 
  geom_point(color = 'lightgreen')+
  xlab('Average Expression') + 
  ylab('Variance')+ labs(title='(b)')

p3 <- ggplot(data_cor, aes(x=corr_p)) + 
  geom_histogram(fill = 'lightgreen')+
  xlab('Pearson Correlation between Genes and Gestacional Age')+ labs(title='(c)')

p4 <- ggplot(data_cor, aes(x=corr_p, y = corr_s)) + 
  geom_point(color = 'lightgreen')+ labs(title='(d)')+
  xlab('Pearson Correlation') + ylab('Spearman Correlation')

grid.arrange(p1, p2, p3, p4, nrow = 2)
```

### Exercise 1: 
1. What can we conclude from the plots above? 
2. Why do we not see correlation between gene expression and GA for the test set? 
3. Can you think of any other intesresting plots or analysis? 

```{r}
#insert here a new plot or analysis you think is instering

```


##2. Features Extraction 

The original dataset has 32830 genes. However, from the correlation plot, we know that some of these genes aren't associate with GA. In addition, genes with low variance might not contribute towards prediction of gestational age. Also, large datasets contain noise that may affect the efficacy of our machine leanring models. 

Therefore, we will use 5 methods to extract features.
1. Elimination of 30% of genes with lowest variance; 
2. Elimination of 30% of genes with lowest correlation with GA; 
3. PCA (principal component analysis)
4. Random Forest 
5. Autoenconder

```{r, results='hide'}

## 1 and 2) Removing genes with low variance and low absolute correlation
data_cor$corr_p = abs(data_cor$corr_p)
keep = subset(data_cor, variance>quantile(data_cor$variance, 0.5) & corr_p>quantile(data_cor$corr_p,0.5))

data2 = subset(data0, select = keep$col)
data2 = scale(data2)

## Saving for future use
write.csv(data2, paste0(feat_dir,'/features_raw.csv'), row.names=T)

```

```{r}
## 3) PCA
PrePCA = preProcess(data2,method="pca")
feat.pca = predict(PrePCA,data2)
write.csv(feat.pca, paste0(feat_dir,'/features_pca.csv'), row.names=T)

```

```{r}
## 4) Random Forest 
metric = "Accuracy"
GA = data1$GA
data1a = subset(data1, select =keep$col)
data1 = data.frame(GA, data1)
if(!exists('PreRF')){
  PreRF = caret::train(y=data1$GA, x=data1[,-1],  method='ranger',importance='impurity')
  PreRF.i = varImp(PreRF)$importance
}else{
  PreRF.i = varImp(PreRF)$importance
}

PreRF.i= data.frame(col = rownames(PreRF.i),PreRF.i)
rownames(PreRF.i) = NULL
PreRF.i = PreRF.i[-1,]
PreRF.i = PreRF.i[order(PreRF.i$Overall,decreasing = T),]

feat.ra = data0[,PreRF.i$col[0:500]]
write.csv(feat.ra, paste0(feat_dir,'/features_ra.csv'), row.names = T)

```


```{r}
## 5) Autoencoder 
if(!exists('preA')){
  preA = autoencoder(data2, hidden.layers = c(1000, 500, 1000))
  feat.A = encode(preA, data2)  
}else{
  feat.A = encode(preA, data2) 
}
rownames(feat.A) = rownames(feat.ra)

feat.A = read.csv(paste0(feat_dir,'/features_a.csv'), header=T, sep = ',')
#write.csv(feat.A, paste0(feat_dir,'/features_a.csv'), row.names = T)
```

### Exercise 2: 
1. What parameters do you think you could change? 
2. Choose between PCA and Random Forest. Try to save a new features set with a different number of features. 

```{r}

a = read.csv(paste0(feat_dir,'/features_pca.csv'), header=T, sep = ',')
b = read.csv(paste0(feat_dir,'/features_ra.csv'), header=T, sep = ',')
d = read.csv(paste0(feat_dir,'/features_a.csv'), header=T, sep = ',')

#CHANGE THESE NUMBERS BELLOW, OPTIONS ---> a:0-556,b:0-501, d:0-501 
feat.unique = data.frame(a[,c(0:100)],b[,c(1:10)], d[,c(100:200)] )

#COPY HOW TO DO IT FROM CODE ABOVE AND BE CAREFUL TO NOT REPLACE OTHER FEATURES
write.csv(feat.unique, paste0(feat_dir,'/features_unique.csv'), row.names = T)

```


##3. Machine Learning models

In this section we will work with machine learning models. In the last section, we created 5 different sets of features (RAW, PCA, RA, A, UNIQUE). We will apply these features onto 3 models: Linear Regression, RVM (relevance vector machines), and Random Forest. The best feature+model combination will be selected based on RMSE (root mean square errer) scores, or the difference between the values we observed on the training set for GA (Gestacional Age) and the predicted GA values i.e. the lower the better.

Before we begin to test the models, we need to set up an experimental framework so that we can evaluate our results once they come out. One such technique is a resampling procedure called ```cvn```-fold cross validation (here we set ```cvn``` to 10).

Typically, data sets are split up into a train and test set; we extract these indices into the ```train_index_``` and ```train_index_val``` variables respectively. The former set of samples is used to train our model, which is later tested on the test set, a set of samples our model has not seen before. However, since we assume we do not have the GA for our test set, we would not know how our model performs on the test set.

Therefore, we further split the train set into ```cvn```=10 equal sized chunks. Since we know the Ga to all our train samples, we can evaluate the model a total of 10 times, each time training the model 9 of the chunks and testing the model on 1 of the chunks. Thereafter, these metrics can be combined (e.g. via a mean) to produce an appropriate evaluation of the model.

### Exercise 3:
1. Why do we want to perform cross validation? Why would we cross validating the model ```cvn``` times over just 1 time?



```{r}
## 1) prep cvn-fold cross validation & rmse function
cvinds_path = paste0(root,"/cvinds.Rdata")
if(file.exists(cvinds_path)){
  load(cvinds_path)
} else {
  cvn = 10
  train_index = which(meta$Train==1) #selecting only the training examples
  test_index = which(meta$Train==0)
  train_index_val = sample(train_index, ceiling(length(train_index)/11)) #cross validation set
  train_index_ = sample(train_index[!train_index%in%train_index_val]) #removing cross validation set from training set 
  ga_val = as.numeric(meta$GA[train_index_val])
  ga_ = as.numeric(meta$GA[train_index_]) 
  save(cvn,train_index,train_index_val,train_index_,ga_val,ga_, file=cvinds_path)
}

require(Metrics)

```

The RMSE (root mean squared error) is a metric that indicates how different our models' predicted GA are compared to the true GA. In other words, the smaller the RMSE, the better.This will be the metric used to compare the models we will be working on. 

###3.1 Linear Regression 

Remember $y=ax+b$? That's linear regression in a nutshell. Linear regression assumes that there is a linear relationship between the given multidimensional RNAseq train data $x$ and the corresponding variable we want to predict GA $y$. $a$ represents the slope and $b$ represents the y-intercept. In the real world though, there often isn't a perfect linear relationship, so below, we try to estimate the best line (i.e. $a$ and $b$ parameters) given $x$ and $y$.

Linear regression is fast and intuitive, but there are many data sets that do not conform with this linear assumption.

```{r, warning=FALSE}
feature_type = 'pca' #options are 'pca', 'raw', 'ra' for random forest and 'a' for autoencoder
## 0) load feature 
features = read.csv(paste(feat_dir,"/features_",feature_type,'.csv', sep = ''))
rownames(features) = features[,1]
features = as.matrix(features[,-1])

#SPliting features of training and validation set
features_ = features[train_index_,]
features_val = features[train_index_val,]
train_ = data.frame(ga_,features_)

#Linear REgression model
model1 = lm(ga_ ~ .,data = train_)
#Predictions
predictions_m1 = predict(model1, newdata = train_[,-1])
predictions_m1[predictions_m1<8] = 8
predictions_m1[predictions_m1>42] = 42
#Error from the observed values and fitted values
rmse(ga_, predictions_m1)

```

```{r}
plotdata = data.frame(ga_,predictions_m1)
ggplot(plotdata, aes(x = ga_,y = predictions_m1))+
  geom_point(color = 'lightgreen')+
  xlab('Linear Regression Model - Training Set')+
  ylab('Predicted Values')

```

### Exercise 4: 
Calculate the predictions for ```features_val``` and ```ga_val``` using linear regression. With the predicted values, make a plot ```predictions_m1_val``` and ```ga_val```
```{r, warning=FALSE}
predictions_m1_val = predict(model1, data.frame(features_val))
rmse(ga_val,predictions_m1_val)
```
```{r}
plotdata = data.frame(ga_val,predictions_m1_val)
ggplot(plotdata, aes(x = ga_val,y = predictions_m1_val))+geom_point(color = 'blue')+xlab('')+ylab('')

```



###3.2 RVM

RVM (relevance vector machines) is an application of the Bayesian treatment of general linear models to SVMs (support vector machines). 

SVMs is a model that draws support vectors (i.e. a line) to separate data in an effort to predict their corresponding discrete classes (e.g. early GA, late GA) as opposed to predicting continuous values (e.g. 10 week GA, 32 week GA). SVMs extract this line from the train data by first calculating a distance between all the data points. These distances are used to define a best line, one that not only separates the data of different classes, but are also the largest distance from data points of both classes (i.e. the line sites in the middle of the space between the two classes, rather than being closer to one or the other).

RVM is functionally identicle to SVM. The Bayesian treatment is a fancy way of saying that instead of categorizing all data points on one side of the line as a single class, let's give them continuous class values that indicates how far the data points are from the line.

In both cases, the user defined formua that calculates distances between data points are called "kernels".

RVMs are fast to train and easy to use, but it can become less effective if the data set used has too many features because it depends on possibly dimension sensative kernel functions.

```{r}

model2 = rvm(x = features_, ga_, type="regression")
predictions_m2 = predict(model2,data = as.data.frame(features_))
rmse(ga_,predictions_m2)


plotdata = data.frame(ga_,predictions_m2)
ggplot(plotdata, aes(x = ga_,y = predictions_m2))+
  geom_point(color = 'lightgreen')+
  ylab('Predicted Values')+
  xlab('RVM - Training Set')

```
### Exercise 5: 
Calculate the predictions for ```features_val``` and ```ga_val``` using rvm. With the predicted values, make a plot ```predictions_m2_val``` and ```ga_val```
```{r, warning=FALSE}
predictions_m2_val = predict(model2,features_val)
rmse(ga_val,predictions_m2_val)
```


```{r}
plotdata = data.frame(ga_val,predictions_m2_val)
ggplot(plotdata, aes(x = ga_val,y = predictions_m2_val))+geom_point(color = 'red')+xlab('')+ylab('')

```


###3.3 Random Forest

A random forest is an aggregate of decision trees formed via model training. A decision tree is a decision making template in the form of "if this then that"; in our case, this decision would be "if this gene is expressed often, then the woman is in her 32nd week of gestation". 

Decision trees and the rules they contain are interpretable and easy to understand. As well combining many decision trees together yield more robust results. However, it can be slow if there are too many features and it is more suited for predicting categorical results.

```{r}
model3 = ranger(ga_~.,data = train_)
predictions_m3 = predict(model3, data = as.data.frame(features_)) 
rmse(ga_,predictions_m3$predictions)

plotdata = data.frame(ga_,predictions_m3$predictions)
ggplot(plotdata, aes(x = ga_,y = predictions_m3.predictions))+
  geom_point(color = 'lightgreen')+
  ylab('Predicted Values')+
  xlab('Random Forest - Training Set')

```
### Exercise 6: 
Calculate the predictions for ```features_val``` and ```ga_val``` using Random Forest. With the predicted values, make a plot ```predictions_m1_val``` and ```ga_val```
```{r, warning=FALSE}
predictions_m3_val = predict(model3, features_val)
rmse(ga_val,predictions_m3_val$predictions)
```


```{r}
plotdata = data.frame(ga_val,predictions_m3_val$predictions)
ggplot(plotdata, aes(x = ga_val,y = predictions_m3_val.predictions))+geom_point(color = 'blue')+xlab('')+ylab('')

```


### Exercise 7:
Repeat Linear Regression, RVM, and Random Forest with 2 other sets of features: 'raw', 'ra' for random forest or 'a' for autoencoder and the 'unique' set of features you made in Exercise 2. From the set of features you chose, which one had the best result? Why?

SET OF FEATURES 1
```{r, warning=FALSE}
#----------------------   CHANGE HERE
feature_type = 'ra' #options are 'pca', 'raw', 'ra' for random forest and 'a' for autoencoder

## 0) load feature 
features = read.csv(paste(feat_dir,"/features_",feature_type,'.csv', sep = ''))
rownames(features) = features[,1]
features = as.matrix(features[,-1])
features_ = features[train_index_,]
features_val = features[train_index_val,]
train_ = data.frame(ga_,features_)

#----------------------   Linear Regression
model1 = lm(ga_ ~ .,data = train_)
predictions_m1 = predict(model1, newdata = train_[,-1])
predictions_m1[predictions_m1<8] = 8
predictions_m1[predictions_m1>42] = 42

rmse(ga_, predictions_m1)

#PLOT 
plotdata = data.frame(ga_,predictions_m1)
ggplot(plotdata, aes(x = ga_,y = predictions_m1))+
  geom_point(color = 'lightgreen')+
  xlab('Linear Regression Model - Training Set')+
  ylab('Predicted Values')

#----------------------   CHANGE HERE: ADD LINEAR REGRESSION FOR THE VALIDATION SET 
predictions_m1_val = predict(model1, newdata = data.frame(features_val))
predictions_m1_val[predictions_m1_val<8] = 8
predictions_m1_val[predictions_m1_val>42] = 42

rmse(ga_val,predictions_m1_val)

#----------------------   ADD SVM FOR TRAINING AND VALIDATION SET 


model2 = rvm(x = features_, ga_, type="regression")
predictions_m2 = predict(model2,data = as.data.frame(features_))
predictions_m2[predictions_m2<8] = 8
predictions_m2[predictions_m2>42] = 42

rmse(ga_,predictions_m2)

predictions_m2_val = predict(model2, data.frame(features_val))
predictions_m2_val[predictions_m2_val<8] = 8
predictions_m2_val[predictions_m2_val>42] = 42
rmse(ga_val,predictions_m2_val)

#----------------------   ADD RAM FOR TRAINING AND VALIDATION SET 


model3 = ranger(ga_~.,data = train_)
predictions_m3 = predict(model3, data = as.data.frame(features_)) 
predictions_m3$predictions[predictions_m3$predictions<8] =8
predictions_m3$predictions[predictions_m3$predictions>42] =42

rmse(ga_,predictions_m3$predictions)

predictions_m3_val = predict(model3, features_val)
predictions_m3_val$predictions[predictions_m3_val$predictions<8] =8
predictions_m3_val$predictions[predictions_m3_val$predictions>42] =42
rmse(ga_val,predictions_m3_val$predictions)

```


SET OF FEATURES unique
```{r, warning=FALSE}
#----------------------   CHANGE HERE
#feature_type =  #options are 'pca', 'raw', 'ra' for random forest and 'a' for autoencoder

## 0) load feature 
features = read.csv(paste(feat_dir,"/features_",feature_type,'.csv', sep = ''))
rownames(features) = features[,1]
features = as.matrix(features[,-1])
features_ = features[train_index_,]
features_val = features[train_index_val,]
train_ = data.frame(ga_,features_)

#----------------------   Linear Regression
model1 = lm(ga_ ~ .,data = train_)
predictions_m1 = predict(model1, newdata = train_[,-1])
predictions_m1[predictions_m1<8] = 8
predictions_m1[predictions_m1>42] = 42
rmse(ga_, predictions_m1)

#PLOT 
plotdata = data.frame(ga_,predictions_m1)
ggplot(plotdata, aes(x = ga_,y = predictions_m1))+
  geom_point(color = 'lightgreen')+
  xlab('Linear Regression Model - Training Set')+
  ylab('Predicted Values')

#----------------------   CHANGE HERE: ADD LINEAR REGRESSION FOR THE VALIDATION SET 

#----------------------   ADD SVM FOR TRAINING AND VALIDATION SET 

#----------------------   ADD RAM FOR TRAINING AND VALIDATION SET 






```

###Exercise 8: 
Create a plot to identify which model or set of features is the best so far. 

```{r, warning=FALSE}
#Saving the results of the predictions when the .unique features were used 
output = data.frame(
  model = c('LR','RVM','RF'),
  feature = c(rep('ra',3) ),
  rmse_t = c(rmse(ga_, predictions_m1),rmse(ga_,predictions_m2),rmse(ga_,predictions_m3$predictions) ),
  rmse_val = c(rmse(ga_val, predictions_m1_val),rmse(ga_val,predictions_m2_val),rmse(ga_val,predictions_m3_val$predictions) )
)
output
#output = read.table(paste0(input_dir,'/output_basic_models.csv'),sep=',', header=T)
#output = rbind(output,output_unique)
```


```{r}
p7<- ggplot(data = output, aes(x = , y = )) +geom_bar(stat='identity')+xlab('Training set')
p8<- ggplot(data = output, aes(x = , y = )) +geom_bar(stat='identity')+xlab('Validation set')
p7 
p8

```



## 3.4 Extra Models
Basides the models shown in the prevous section, there are many other models or variations of that models we could use to explore this problem. However, some of them are more time consuming, thus we trained some extra models and we saved their results, so now all we need is make the predictions and check their RMSE. 


First we define the features and the model we want to explore:
```{r}
# overwrite model?
overwrite = F

# features
feature_name = "features_a" # ra, a, pca, raw

##load the features
features = read.csv(paste0(feat_dir,"/", feature_name,'.csv'))
rownames(features) = features[,1]
features = as.matrix(features[,-1])

#spliting the validation and testing set
features_ = features[train_index_,]
features_val = features[train_index_val,]
```


Then we run this chunk to load or run the model with the features and model specified.
```{r}
# model options:
#enet, foba, gaussprPolyg, gaussprRadial,glmnet, icr, kernelpls, krlsRadial, lars2, lasso, leapBackward, leapForward, leapSeq, nnls, partDSA, pls , plsRglm, rbf, rpart, rqlasso, rvmPoly, rvmRadial, simpls, spikeslab, spls, svmPoly, svmRadial, svmRadialCost, svmRadialSigma, widekernelpls, rqnc, nodeHarvest, mlpML, xgbDART
#see https://topepo.github.io/caret/available-models.html for more information about each model

model = "enet"
parsi = expand.grid(lambda=10^runif(5, min=-5, 1), fraction=runif(5, min=0, max=1)) # parsi = NULL; if no parameters need to be tested
fitcv = trainControl(method="cv", number=cvn)


fname = paste0(model_dir,"/",feature_name,"/",model,".Rdata")
if (!file.exists(fname) | overwrite) { try ({
  t2i = NULL
  if (!is.null(parsi)) {
    t2i = caret::train(y=ga_, x=features_, model, trControl=fitcv, tuneGrid=parsi)
  } else {
    t2i = caret::train(y=ga_, x=features_, model, trControl=fitcv)
  }
  if (!is.null(t2i)) save(t2i, file=fname)
}) }
load(fname)

# results to data frame
df = data.frame(
  rmse=t2i$results$RMSE[which.min(t2i$results$RMSE)],
  time=as.numeric(t2i$times$everything[3]),
  model_=t2i$modelInfo$label, 
  feature=feature_name, model=model, 
  par=paste0( paste0(names(t2i$bestTune), collapse="_"), ": ", paste0(t2i$bestTune, collapse="_") )
  , stringsAsFactors=F)
df

```
Now we will explore the predictions of these models+features:

```{r, results='hide', message=FALSE}
##get test prediction results from models ----------------------
pred = predict(t2i,newdata = as.data.frame(features_))
pred_val = predict(t2i,newdata = as.data.frame(features_val))
title = paste(feature_name,model, sep = '-')

# plot graph to compare models
plot1 = data.frame(ga_,pred)
plot2 = data.frame(ga_val, pred_val)
pred1 <- ggplot(plot1, aes(x = ga_,y=pred)) + geom_point(color='orange')+
  ylab('Predicted')+xlab('Training set')+labs(title=title)

pred2 <- ggplot(plot2, aes(x = ga_val,y=pred_val)) + geom_point(color='orange')+
  ylab('Predicted')+xlab('Validation set')

grid1 <- grid.arrange(pred1, pred2, nrow = 1)
ggsave(filename = paste(result_dir,'/','predictions-',title,'.png',sep=''),grid1)

```


### Exercise 09: 
Explore more combinations of models and features. Keep the results of the top 3 models you found. Discuss with your colleagues their results. Answer: 
* Are they similar? How you defined which models were the best? 
* Can you think in a more efficient way to compare all the models available? 


####Model 1 
```{r}
# features
feature_name = "features_a" # ra, a, pca, raw

##load the features
features = read.csv(paste0(feat_dir,"/", feature_name,'.csv'))
rownames(features) = features[,1]
features = as.matrix(features[,-1])

#spliting the validation and testing set
features_ = features[train_index_,]
features_val = features[train_index_val,]

# model options:
#enet, foba, gaussprPolyg, gaussprRadial,glmnet, icr, kernelpls, krlsRadial, lars2, lasso, leapBackward, leapForward, leapSeq, nnls, partDSA, pls , plsRglm, rbf, rpart, rqlasso, rvmPoly, rvmRadial, simpls, spikeslab, spls, svmPoly, svmRadial, svmRadialCost, svmRadialSigma, widekernelpls, rqnc, nodeHarvest, mlpML, xgbDART
#see https://topepo.github.io/caret/available-models.html for more information about each model

model = "enet"
parsi = expand.grid(lambda=10^runif(5, min=-5, 1), fraction=runif(5, min=0, max=1)) # parsi = NULL; if no parameters need to be tested
fitcv = trainControl(method="cv", number=cvn)

fname = paste0(model_dir,"/",feature_name,"/",model,".Rdata")
if (!file.exists(fname) | overwrite) { try ({
  t2i = NULL
  if (!is.null(parsi)) {
    t2i = caret::train(y=ga_, x=features_, model, trControl=fitcv, tuneGrid=parsi)
  } else {
    t2i = caret::train(y=ga_, x=features_, model, trControl=fitcv)
  }
  if (!is.null(t2i)) save(t2i, file=fname)
}) }
load(fname)

# results to data frame
df = data.frame(
  rmse=t2i$results$RMSE[which.min(t2i$results$RMSE)],
  time=as.numeric(t2i$times$everything[3]),
  model_=t2i$modelInfo$label, 
  feature=feature_name, model=model, 
  par=paste0( paste0(names(t2i$bestTune), collapse="_"), ": ", paste0(t2i$bestTune, collapse="_") )
  , stringsAsFactors=F)
df

##get test prediction results from models ----------------------
pred = predict(t2i,newdata = as.data.frame(features_))
pred_val = predict(t2i,newdata = as.data.frame(features_val))
title = paste(feature_name,model, sep = '-')

# plot graph to compare models
plot1 = data.frame(ga_,pred)
plot2 = data.frame(ga_val, pred_val)
pred1 <- ggplot(plot1, aes(x = ga_,y=pred)) + geom_point(color='orange')+
  ylab('Predicted')+xlab('Training set')+labs(title=title)

pred2 <- ggplot(plot2, aes(x = ga_val,y=pred_val)) + geom_point(color='orange')+
  ylab('Predicted')+xlab('Validation set')

grid1 <- grid.arrange(pred1, pred2, nrow = 1)
grid1
ggsave(filename = paste(result_dir,'/','predictions-',title,'.png',sep=''),grid1)

```

####model 2 
```{r}

```

####model 3
```{r}

```


#Prediction for competition submission

As with all Data Science/Machine Learning project, we have spent most of our time on procesing the data, exploring the dataset, creating/selecting features, and exploring different models to figure out which pipeline is the best at answering our main question: predict the GA of our testing set. 

Now, we will create the final predictions for our test set.

```{r}
#RUN HERE ONCE MORE THE MODEL YOU THINK IS THE BEST
#REMEMBER TO LOAD THE FEATURES AND TRAIN THE MODEL USING THE TRAINING SET
```



```{r}
features_pred = features[test_index,]

#INSERT HERE MODEL
#predictions = predict(-MODEL-, data = as.data.frame(features_pred)) 

# submission template
class_final = read.csv(paste0(input_dir,"/TeamX_SC1_prediction.csv")) 
class_final$GA = round(predictions,1)
class_final$GA[class_final$GA<8] = 8 
class_final$GA[class_final$GA>42] = 42 
write.csv(class_final, file=paste0(result_dir,"/submission_file.csv"),row.names=F)
```




#Bonus Round: Ensembles

You can also try to ensemble, or combine results from multiple models to see if you can get a better result, try it out below!

```{r}
load(cvinds_path)
fitcv = trainControl(method="cv", number=10)
rmse = function(x,y) sqrt(mean((x-y)^2))

# load features
xi = "features_pca" # EDIT FEATURE HERE
feature = read.csv(paste0(feat_dir,"/", xi,".csv"))
rownames(feature) = feature[,1]; feature = as.matrix(feature[,-1])
mtr = feature[train_index_,]
mte = feature[train_index_val,]
mte0 = feature[test_index,]

# train models; note some are classification models too
# exercise: pick and choose a few models to ensemble
models_reg = c( # EDIT MODEL(S) HERE
  "enet", # 8.5
  "glmnet", # 8.5
  "kernelpls", # 8.5
  "krlsRadial", # 8.5
  "lars2", # 8.5
  "nnls", # 8.5
  "rbf", # 8.4
  "rqlasso", # 8.4
  "rvmPoly", # 8.5
  "rvmRadial", # 8.5
  "simpls", "spikeslab", # "spls", # 8.5; spls takes a bit longer 950
  "rqnc" # 8.4
)
model_list = caretList(
  y=ga_tr, x=mtr, trControl=fitcv, metric="RMSE",
  methodList=models_reg, continue_on_fail=T)

# weighted linear combination of model predictions
ensemble = caretEnsemble(
  model_list, 
  metric="RMSE",
  trControl=fitcv)
summary(ensemble)

# predict test set GA; get rmse
# exercise: do ensembling methods perform better? why do you think that is?
ens_preds = predict(ensemble, newdata=mte)
rmse(ens_preds, ga_val)

# save as submission
ens_predfinal = predict(ensemble, newdata=mte0)
submission$GA = round(ens_predfinal,1)
submission$GA[submission$GA<8] = 8
submission$GA[submission$GA>42] = 42
write.csv(submission, file=paste0(result_dir,"/submission_ensemble.csv"), row.names=F)
```