miniML.Rmd

---
title: "ML_10min"
author: "Eva Kleingeld"
date: "August 4, 2016"
output: pdf_document
---
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)


rm(list=ls())
  
#install required packages here

#install.packages("Metrics")
#install.packages("rattle")
#install.packages("caret", dependencies = c("Imports", "Depends", "Suggests"))
#install.packages("rpart.plot")
#install.packages("DMwR")
#install.packages("RSNNS")
#install.packages("devtools")
#install.packages("kernlab")
```


The following code reads in the 10 min subset and environmental data. 
All "suspect" data is removed.
The two datasets are then merged.
The merged dataset is split in a test and training set. 

```{r}
## Read in subset & environmental data

#For 5 min: load("/usr/people/kleingel/R/Data_subsets/Data_5min.Rda")
#Data_5min<-Data_5min[Data_5min$QUALITY=="valid", ]

load("/usr/people/kleingel/Projects/MLProject/Data_10min.Rda")
Data_10min <- Data_10min[Data_10min$QUALITY == "valid", ]

## Drop the quality column
Data_10min <- Data_10min[ ,-7]

load("/usr/people/kleingel/Projects/MLProject/Env_Data.Rda")

## Merge subset and environmental data
data_GMS<-merge(Data_10min,Env_Data_4,by.x=c("LOCATION","SENSOR"),by.y=c("MISD","SENSOR"))

## Split data_GMS into a training and test set
Time_1 <- as.POSIXct("2016-06-28 00:05:00", format = "%Y-%m-%d %H:%M:%S", tz = "GMT")
Time_2 <- as.POSIXct("2016-06-28 00:10:00", format = "%Y-%m-%d %H:%M:%S", tz = "GMT")  

Train_data <- subset(data_GMS, data_GMS$TIMESTAMP == Time_1) 
Test_data <- subset(data_GMS, data_GMS$TIMESTAMP == Time_2)

## Drop the LOC/Sensor/Timestamp data because this will not be used for building the models
Train_data <- Train_data[ ,-(1:3)]
Test_data <- Test_data[ ,-(1:3)]


```


Center and scale

```{r}
library(caret)

# Test centering and scaling on a dataset without dummy vars.
new_data <- cbind(Train_data[ ,1:3], Train_data[5:7])

# Plot histograms of the new_data columns
par(mfrow = c(1, length(new_data)))
for (i in seq_along(colnames(new_data))){
  
  hist(new_data[ , i], main =  paste("Histogram of" , colnames(new_data)[i]), xlab = colnames(new_data)[i])
}

# Apply centering and scaling
# You can also remove NA values here
CenterScale <- preProcess(new_data, method = c("center", "scale"), na.remove = TRUE)
New_data_processed <- predict(CenterScale, new_data)

# Plot histograms of the centered and scaled data
par(mfrow = c(1, length(New_data_processed)))
for (i in seq_along(colnames(New_data_processed))){
  
  hist(New_data_processed[ , i], main =  paste("Histogram of" , colnames(New_data_processed)[i]), xlab = colnames(New_data_processed)[i])
}
dev.off()

# Now apply centering and scaling to the entire data_GMS dataset
CenterScale_GMS <- preProcess(x = data_GMS, method = c("center", "scale"), na.remove = TRUE)
GMS_processed <- predict(CenterScale_GMS, data_GMS)

hist(GMS_processed$Unix_Time)
hist(GMS_processed$Brug_of_Viaduct_0)

# Apply scaling and centering to the train and test data
# You don't want to center or scale the dummy variables, so you transform only the first part of the data
# Because we use only 5 min datasets it is useful to first drop the Unix_Time column: This value doesn't change because you only have 5 min data (:one measurement). For the bigger dataset you should not drop Unix_time!!!!

Train_data <- Train_data[ ,-4]
Test_data <- Test_data[ ,-4]

# center and scale train and test data
CenterScale_Train <- preProcess(x = Train_data[ ,(1:6)], method = c("center", "scale"), na.remove = TRUE)
Train_data[ ,(1:6)] <- predict(CenterScale_Train, Train_data[ ,(1:6)])
Test_data[ ,(1:6)] <- predict(CenterScale_Train, Test_data[ , (1:6)])

# remove rows with NA values

Train_data <- Train_data[complete.cases(Train_data), ]
Test_data <- Test_data[complete.cases(Test_data), ]

```

Build different inputs and outputs for the ML models

```{r}

#Build different inputs and outputs for the ML models

## Input (: predictors), without TL and TD
Train_X_1 <- Train_data[4:(length(colnames(Train_data)))]
Test_X_1 <- Test_data[4:(length(colnames(Test_data)))]

## Input (: predictors), with TL and TD
Train_X_2 <- Train_data[2:(length(colnames(Train_data)))]
Test_X_2 <- Test_data[2:(length(colnames(Test_data)))]

## Output (:Target variable)
Target_Train <- Train_data[, 1]
Target_Test <- Test_data[, 1]

```







Now, we start building models based on different ML algorithms. 

First, we build a linear model 

```{r}
library(caret)

n<-colnames(Train_data)[2:(length(colnames(Train_data)))]
f <- as.formula(paste("TEMP ~", paste(n[!n %in% "TEMP"], collapse = " + ")))

#LinearModel <- train(form = f, data = Train_data, method = "lm")
LinearModel2 <- train(x = Train_X_2, y = Target_Train, method = "lm")

# Get the lm summary
# The summary indicates that five variables are perfectly collinear (Coefficients: 5 not defined because of singularities)
summary(LinearModel2)
 
# Calculate RMSE 
Linear_Residuals <- residuals(LinearModel2)
sqrt(mean((Linear_Residuals)^2))

# Caret RMSE
LinearModel2$results

# Predict with the lm model
Linear_Predict <- extractPrediction(list(LinearModel2), testX = Test_X_2, testY = Target_Test)

# Plot observed vs predicted
plotObsVsPred(Linear_Predict)
ggplot(data = Linear_Predict) + geom_point(aes(x = obs, y = pred))+ geom_abline(aes(slope =1, intercept = 0))
 


```




Here, a neuralnet neural network is built. 
```{r}

# Build a neural network with neuralnet (from neuralnet package)
library(neuralnet)

# First build a neural network without TL and TD: set 4 in colname selection n
# If you include TL and TD you may have to remove some stations because these variables can be NA -> then NN won't run

n<-colnames(Train_data)[2:(length(colnames(Train_data)))]
f <- as.formula(paste("TEMP ~", paste(n[!n %in% "TEMP"], collapse = " + ")))
nn<-neuralnet(f,data=Train_data,hidden=c(5),linear.output=T)
plot(nn)

pr.nn <- compute(nn,Test_X_2) # select only relevant variables

# Calculate the MSE
library(Metrics)
library(ggplot2)

mse(as.vector(Target_Test), as.vector(pr.nn$net.result))
rmse(as.vector(Target_Test), as.vector(pr.nn$net.result))

## plot the actual values v.s. the predicted values
ggplot() + geom_point(aes(x = as.vector(Target_Test), y = as.vector(pr.nn$net.result))) +
           labs(x = "Measured temperature", y = "Predicted temperature")+ geom_abline(aes(slope =1, intercept = 0))

```

Here a multilayer perceptron is built (: a neural network - like algorithm)

```{r}
library(RSNNS)
library(caret)

# Build the mlp model
MLPModel <-  train(x  = Train_X_2, y = Target_Train, method = "mlp")
compiler::setCompilerOptions(optimize = 1)

# Extract a summary
summary(MLPModel$finalModel)
print(MLPModel$finalModel)

# Plot MLP model
# Import the plot.nnet function from Github. (There is a newer version with additional options, you can also download this function from GitHub, see research log)
# This function can be used to plot multilayer perceptrons
library(devtools)
source_url('https://gist.githubusercontent.com/fawda123/7471137/raw/466c1474d0a505ff044412703516c34f1a4684a5/nnet_plot_update.r')

# The code below shows that you cannot plot the $finalModel, but you can plot if you use mlp directly instead of via train
plot.nnet(MLPModel$finalModel)
MLP_test <- mlp(x  = Train_X_2, y = Target_Train)
plot.nnet(MLP_test)

# Even though 
class(MLP_test)
class(MLPModel$finalModel)

# Further inspection reveals that MLP$finalModel contains more members (: list items) than the model made directly with the mlp package
# Remove these members
MLP_Adapted <- MLPModel$finalModel
MLP_Adapted[16:19] <- NULL
plot.nnet(MLP_Adapted)

plot.nnet(MLP_test)

# Plot activation map (what is this?)
plotActMap(Train_X_1)

# Plot the iterative error
plotIterativeError(MLPModel$finalModel)

# Extract predictions for MLP
MLP_Predict <- extractPrediction(list(MLPModel), testX = Test_X_2, testY = Target_Test)

# Plot observerd versus predicted
# Deze plots laten zien dat het MLP op het moment maar 1 waarde voorspelt: ongeveer de gemiddelde waarde van Train_data$Temp
# Meer data/normaliseren kan de oplossing zijn
plotObsVsPred(MLP_Predict) 
ggplot(data = MLP_Predict) + geom_point(aes(x = obs, y = pred))+ geom_abline(aes(slope =1, intercept = 0))


```



Build a Radial Basis Function Neural Network (rbf in short)
```{r}
# Takes very long time to run: 10.5 minutes! (commented for now)
# This probably has to do with the default settings. Or the model tries to reduce the bias and this bias is very high right now.
# Or input needs to be scaled first
# RBFModel <- train(x = Train_data[7:(length(colnames(data_GMS)))], y = Train_data$TEMP, method = "rbfDDA")
compiler::setCompilerOptions(optimize = 1)

# If you run from RSNNS package, the RBF is built quite quickly:
RBFModel2 <- rbf(x = Train_X_2, y = Target_Train)
compiler::setCompilerOptions(optimize = 1)

# High RMSE
plot(RBFModel2)
head(summary(RBFModel2))

# Model 2 has high SSE
plotIterativeError(RBFModel2)

# Model from train doesn't predict (pred = 0)
RBF_Predict <- extractPrediction(list(RBFModel), testX = Test_X_2, testY = Target_Test)

# predicting for the other model:
RBF_Predict2 <- predict(RBFModel2, Test_X_2)

# For the other model, observed vs. predicted:
ggplot() + geom_point(aes(x = Test_data$TEMP, y = RBF_Predict2))+ geom_abline(aes(slope =1, intercept = 0))


```



Build a decision tree
```{r}
# Building a tree with the caret package
library(caret)
library(rpart)
library(rattle)

# Build a tree model with the train data
# There are multiple ways to do this: see below
# Warning message rpart, gone but still bugs when using rpart2
#TreeModel <- train(form = f, data = Train_data, method = "rpart")
TreeModel2 <- train(x  = Train_X_2, y = Target_Train, method = "rpart")
compiler::setCompilerOptions(optimize = 1)


# Plot the decision tree you have built (:TreeModel$finalModel)
fancyRpartPlot(TreeModel$finalModel)
fancyRpartPlot(TreeModel2$finalModel)

# In order to plot a basic tree
plot(TreeModel$finalModel, uniform = TRUE)
text(TreeModel$finalModel, use.n = TRUE, all = TRUE)

# To test variable importance you can use
TreeModel2_Imp <- varImp(TreeModel2)
plot(varImp(TreeModel2))

# Predict with test set based on TreeModel
# The predict() function predicts TW values (:numeric)
# You can write both TreeModel$finalModel and TreeModel
Tree_Predict <- predict(TreeModel2, newdata = Test_X_2)
plot(Tree_Predict)

# The extractPrediction() function gives you more information
# Put your model in a list, because otherwise you can get an error
Tree_Predict2 <- extractPrediction(list(TreeModel2), testX = Test_X_2, testY = Target_Test)

# Plot observed vs. predicted 
plotObsVsPred(Tree_Predict2)

# Plot manually (ggplot2) 
ggplot() + geom_point(aes(x = Tree_Predict2$obs, y = Tree_Predict2$pred))

```

Build a Random Forest (RF)
``` {r}
library(randomForest)

# Buil a Random Forest (RF) with the training data set
# The RF is much slower to build than the decision tree
# Can be useful to set importance to TRUE, does increase computation time
# You can train the RF model in two different ways.
# The first requires you to specify a formula and a dataframe (comp time: 220.6 seconds)
# The second requires you to specify the predictor values and response variable (comp time: 219.5 seconds)

#RFModel <- train(form = f, data = Train_data, method = "parRF", importance = TRUE)
#compiler::setCompilerOptions(optimize = 1)

RFModel2 <- train(x = Train_X_2, y = Target_Train, method = "parRF", importance = TRUE)
compiler::setCompilerOptions(optimize = 1)

# Plot the RF model you have built: How many trees until the error decreases? 
plot(RFModel$finalModel,plotType="level")

# Find a way of beautifully plotting a RF

# Importance of input variables/features for RF:
# Importance measure = 1 = mean decrease in accuracy = IncMSE
# This can only be calculated if importance = TRUE in train() (part of RandomForest package)
Importance_RF_1 <- importance(RFModel$finalModel, type = 1)
dotplot(Importance_RF_1)

# Importance measure = 2 = decrease in node impurity
Importance_RF_2 <- importance(RFModel$finalModel, type = 2)
dotplot(Importance_RF_2)

# With varImp <- preferable, uses %IncMSE
Importance_RF_3 <- varImp(RFModel, scale = FALSE)
plot(Importance_RF_3)

# Make a partial plot
# A partial plot does not make sense for dummy variables. Relevant columns are 7 until 10. Partial plots are only made for relevant columns. 
# Train_data argument in partialPlot is set to a subset of the data because partialPlot cannot handle NA. But because f excludes columns up until Train_data[ ,7] this should not pose a problem
# The for loop uses seq_along() and vector[i] structure because partialPlot throws an error otherwise

Rel_vars <- colnames(Train_data[ ,2:4])

par(mfrow = c(1, length(Rel_vars)))

for (i in seq_along(Rel_vars)){
  print(i)
  partialPlot(RFModel2$finalModel, pred.data = Train_X_2, x.var = Rel_vars[i], xlab = Rel_vars[i], ylab = "Model predicted temperature", main = paste("Partial plot of", Rel_vars[i]))
  
}

dev.off()

# Predict with your RF
# I suspect that the first RF_predict won't run because the input for the extractPrediction function differs from the input that the train function was given. Therefore it may be a good idea to run the train function with the RFModel2 settings. 
# However, TreeModel does not seem to have this issue?
RF_Predict <-  extractPrediction(list(RFModel), testX = Test_X_1, testY = Target_Test)
RF_Predict2 <-  extractPrediction(list(RFModel2), testX = Test_X_2, testY = Target_Test)

# Plot observations versus predictions
# Why does plotObsVsPred have two the same graphs?
plotObsVsPred(RF_Predict2)
ggplot(data = RF_Predict2) +geom_point(aes(x = pred, y = obs))+ geom_abline(aes(slope =1, intercept = 0))



```



Here, we test a Support Vector Machine (in short: svm) with a Radial Basis Function kernel 
```{r}
# Build the SVM model
SVMradModel <- train(x  = Train_X_2, y = Target_Train, method = "svmRadialCost", scale = FALSE)
compiler::setCompilerOptions(optimize = 1)

# summary SVM
summary(SVMradModel$finalModel)
print(SVMradModel$finalModel)

# Get the RMSE
SVMradModel$results

# Predict with the SVM
SVMrad_Predict <- extractPrediction(list(SVMradModel), testX = Test_X_2, testY = Target_Test)

# Plot predicted vs. observed
plotObsVsPred(SVMrad_Predict)
ggplot(data = SVMrad_Predict, aes(x = obs, y = pred)) + geom_point() + geom_abline(aes(slope =1, intercept = 0))

```