Preproc_Round2.Rmd

---
title: "Data Preprocessing including dummy variables"
author: "Eva Kleingeld"
date: "October 2, 2016"
output: pdf_document
---
 
 
In this script the test and train set are:
* subsetted so that they only include valid data
* split into test predictors/target variables and train predictors/target variables
* Altitude is heightened by 10m
* Temperatures are put in Kelvin instead of Celsius
* Data is preprocessed using the preProcess function, the transformations that are applied are:
+ Zero variance removal
+ BoxCox transform
+ centering
+ scaling
+ principal component analysis
* Next, some optional code is included for analyzing the remaining PC's 
* Finally, the resulting predictor sets and target variables are saved to .Rdata files

In between optional code is provided. This code can be used to check the distribution of data/missing station numbers, etc. 

```{r}
# Empty environment
rm(list=ls())

# Install packages
# install.packages("corrplot")
# install.packages("GGally")
# install.packages("lubridate")

```

Read in the environmental dataset, train set and test set.
All sets are given a different name here and the old data frames are removed.
This way you can easily apply the script to different test/train sets.

The load_obj function was taken from this stackoverflow post:
[link](http://stackoverflow.com/questions/5577221/how-can-i-load-an-object-into-a-variable-name-that-i-specify-from-an-r-data-file)

```{r}
# Function to read in, rename an .Rdata file and remove the old .Rdata dataframe in one step
load_obj <- function(f)
{
    env <- new.env()
    nm <- load(f, env)[1]
    env[[nm]]
}

# Read in the environmental dataset
#Env_Data <- load_obj("F:/KNMI/MLProject/Env_Data.Rda")
#Env_Data <- load_obj("/run/media/kleingel/Lexar/KNMI/Env_Data.Rda")
Env_Data <- load_obj("/usr/people/kleingel/Projects/MLProject/Env_Data.Rda")

# Read in the train set
#Train_set <- load_obj("F:/KNMI/Train_data_R2_BIG.Rda")
#Train_set <- load_obj("/run/media/kleingel/Lexar/KNMI/Train_data_R2.Rda")
#Train_set <- load_obj("/usr/people/kleingel/Projects/MLProject/Train_data_R2_BIG.Rda")
Train_set <- load_obj("/usr/people/kleingel/Projects/MLProject/Train_data_3D.Rda")

# Read in the test set
#Test_set <- load_obj("F:/KNMI/Test_data_R2_BIG.Rda")
#Test_set <- load_obj("/run/media/kleingel/Lexar/KNMI/Test_data_R2.Rda")
Test_set <- load_obj("/usr/people/kleingel/Projects/MLProject/Test_data_3D.Rda")

```

### Optional code: Test the distribution of station numbers in the train and test sets

The histograms below show that the distribution of LOCATION (: stations) in the train and test sets is very similar. 

We import the Six_Days data set below to test this. 

```{r}
# Test the distribution of station numbers in test and train
# Load the Six_Days dataset and plot LOCATION distribution next to train/test
# NOTE: This WILL make your computer run a lot more slowly
# To test another predictor simply replace LOCATION with the name of the predictor who's distribution you wish to test

#Six_Days <- load_obj("F:/KNMI/Six_Days.Rda")
Six_Days <- load_obj("/usr/people/kleingel/Projects/MLProject/Six_Days.Rda")

par(mfrow= c(1,3))

hist(Six_Days$LOCATION, breaks = seq(1:1505),main = "Six_Days", freq = FALSE)
hist(Train_set$LOCATION, breaks = seq(1:1505), main = "Train",  freq = FALSE)
hist(Test_set$LOCATION, breaks = seq(1:1505), main = "Test",  freq = FALSE)

par(mfrow= c(1,1))
```

### Optional code: Test the number of missing stations

For preprocessing it is important to know how many stations are recording on each day that was selected for the dataset. Thefore the follwing code was included. The code below tests how many stations were recording for each of the six days in the dataset. The stations that were not recording are referred to as 'missing stations'. 

```{r}
# Test if all stations are recording on all six days

# 1. Get unique station numbers for each of the six days
Stat_DOY <- tapply(Six_Days$LOCATION, Six_Days$DOY_Days, unique)

#2. Test if lengths are equal
for (i in 1:6) {
  print(i)
  print(length(Stat_DOY[[i]]))
}

# Not every day has the same amount of stations
# It is likely that on some days a few stations were down. It is also possible that some stations have been added over the years, so that later years have a higher unique number of sensors. 

# Now get a vector with all the station numbers that occur in the data
All_Stat <- unique(Six_Days$LOCATION)

# Which stations are NOT working on a DOY (: Which stations are in the list of ALL station numbers, but not in the list of station numbers on a certain day)
All_Stat[!(All_Stat %in% Stat_DOY[["9"]])] # DOY 9
All_Stat[!(All_Stat %in% Stat_DOY[[2]])]   # DOY 68
All_Stat[!(All_Stat %in% Stat_DOY[[3]])]   # DOY 83
All_Stat[!(All_Stat %in% Stat_DOY[[4]])]   # DOY 328
All_Stat[!(All_Stat %in% Stat_DOY[[5]])]   # DOY 337
All_Stat[!(All_Stat %in% Stat_DOY[[6]])]   # DOY 362

# So, on the first day (9), we should not have stations (1117, 1153, 1236, 1235, 1002 and 354). 
# This is true
Day_One <- subset(Six_Days, Six_Days$DOY_Days == 9)
any(Day_One$LOCATION == 1117)
any(Day_One$LOCATION == 1153)
any(Day_One$LOCATION == 1236)
any(Day_One$LOCATION == 1235)
any(Day_One$LOCATION == 1002)
any(Day_One$LOCATION == 354)


# Which stations are in ALL DOY subsets 
In_All_DOY <- Reduce(intersect, Stat_DOY)

# Now get the stations that are not in all DOY subsets
Missing_Stat <- setdiff(All_Stat, In_All_DOY)

# Select the missing stations per DOY
Missing_DOY <- list()

for (i in 1:6) {
  print(unique(Six_Days$DOY_Days)[i])
  Missing_DOY[[i]] <- (setdiff(All_Stat, Stat_DOY[[i]]))
}

# Plot the number of missing stations per DOY as a barplot
Nr_missing <- as.data.frame(lapply(Missing_DOY, length))
colnames(Nr_missing) <- unique(Six_Days$DOY_Days)
library(reshape2)
Nr_missing <- melt(Nr_missing)
colnames(Nr_missing) <- c("DOY", "Missing")

library(ggplot2)
ggplot(data= Nr_missing, aes(x= DOY, y= Missing)) + geom_bar(stat="identity") + ylab("Number of missing stations")
``` 


## Build the training set

Here the 'suspect' data is removed: The data that did not make the quality check that was developed by Marieke Dirksen. The GMS train set and environmental data set are merged. Target variable is separated from the train set.

```{r}
# Remove al data which is not labeled as 'valid' 
Train_set <- Train_set[Train_set$QUALITY == "valid", ]

# Drop the QUALITY column by name 
To_drop <- c("QUALITY")
Train_set <- Train_set[ ,!(names(Train_set) %in% To_drop)]

# Merge the train set and the environmental dataset
Train_set <-merge(Train_set,Env_Data,by.x=c("LOCATION","SENSOR"),by.y=c("MISD","SENSOR"))

# Split into input (predictors) and output (target variable)
# Target variable
Target_Train <- Train_set$TEMP

# Drop LOCATION/SENSOR/TIMESTAMP and TEMP
# For now, do not drop "LOCATION"
# If you are removing missing stations, do not remove TEMP here
To_drop <- c("SENSOR", "TIMESTAMP", "TEMP")
Train_set <- Train_set[ ,!(names(Train_set) %in% To_drop)]

# Place all predictors in one data frame
#Predictors_Train <- Train_set


```

## Build the test set

Here the 'suspect' data is removed: The data that did not make the quality check that was developed by Marieke Dirksen. The GMS test set and environmental data set are merged. Target variable is separated from the test set.

```{r}
# Build the test set  -----------------------------------------------------

# Read in the test dataset and remove all data which is not labeled as 'valid' 
Test_set <- Test_set[Test_set$QUALITY == "valid", ]

# Drop the QUALITY column by name 
To_drop <- c("QUALITY")
Test_set <- Test_set[ ,!(names(Test_set) %in% To_drop)]

# Merge the test set and the environmental dataset
Test_set <-merge(Test_set,Env_Data,by.x=c("LOCATION","SENSOR"),by.y=c("MISD","SENSOR"))

# Split into input (predictors) and output (target variable)
# Output
Target_Test <- Test_set$TEMP

# Drop LOCATION/SENSOR/TIMESTAMP and TEMP
# For now do not drop "LOCATION"
# If you are removing missing stations, do not remove TEMP here
To_drop <- c("SENSOR", "TIMESTAMP", "TEMP")
Test_set <- Test_set[ ,!(names(Test_set) %in% To_drop)]

# Place all predictors in one data frame
#Predictors_Test <- Test_set

```


### Optional code: Remove small categories

```{r}
# You can use table(Dataframe$Column) to test how many occurences of each unique value there are
table(Test_set$Rijbaan_7)

# The following categories contain very little 1 (:yes) values
# Rijbaan 5/6/7
# Rijstrook 0/4/5
# Spitsstrook
# These columns are here removed from the dataset

To_drop <- c("Rijbaan_5", "Rijbaan_6", "Rijbaan_7",
             "Rijstrook_0", "Rijstrook_4", "Rijstrook_5", 
             "Spitsstrook")
Train_set <- Train_set[ ,!(names(Train_set) %in% To_drop)]
Test_set <- Test_set[ ,!(names(Test_set) %in% To_drop)]


```


### Optional code: Missing stations in test/train sets

We already know how many stations are missing in the original Six_Days dataset. However, because the test/train sets have had NA values & suspect data removed the test/train sets can be expected to have more missing stations. The optional code below tests which stations are missing for the DOY's in the test/train sets. 

```{r}
##### Now test how many unique stations there are in test/train per DOY

# Get all station numbers in test/train sets
All_Stat_TRAIN <- unique(Train_set$LOCATION)
All_Stat_TEST <- unique(Test_set$LOCATION)

# Station numbers per DOY for train/test
Stat_DOY_TRAIN <- tapply(Train_set$LOCATION, Train_set$DOY_Days, unique)
Stat_DOY_TEST <- tapply(Test_set$LOCATION, Test_set$DOY_Days, unique)


# Calculate missing stations per DOY for train/test

#First, build a function that selects the missing days 
Get_Missing <- function(All_Stations, Stations_per_Day){
Missing_Days <- list()

for (i in seq_along(Stations_per_Day)) {
  Missing_Days[[i]] <- (setdiff(All_Stations, Stations_per_Day[[i]]))}

return(Missing_Days)
}

# Missing days train and test
Missing_TRAIN <- Get_Missing(All_Stat_TRAIN, Stat_DOY_TRAIN)
Missing_TEST <- Get_Missing(All_Stat_TEST, Stat_DOY_TEST)

# Test if the stations that miss per DOY in the train set also miss for the same DOY in the test set
for (i in seq_along(Missing_TRAIN)) {
  print(all(Missing_TRAIN[[i]] %in% Missing_TEST[[i]]))
}

# Okay, so how many missing stations do we have in total?
length(unique(unlist(Missing_TRAIN)))
# And which stations are missing?
Missing_TestTrain <- unique(unlist(Missing_TRAIN))

# If you were to remove all missing stations, compared to the original Six_Days set, how many of the original 335 stations would remain?
Remaining_Train <- 335 - 
  (335 - (length(unique(Train_set$LOCATION))) + length(unique(unlist(Missing_TRAIN))))
```

### Optional code: Remove missing stations

The missing stations miss completely on at least one of the days in the test/train sets. Removing all misssing stations may on the one hand improve the model (maybe it will be easier to recognise patterns). On the other hand removing stations reduces the available spatial data points, thereby making it more difficult to build a spatial model. 

```{r}
# # Test what happens to the PCA when you remove all the missing stations from the test/train sets
# Train_set <- subset(Train_set, !(Train_set$LOCATION %in% Missing_TestTrain))
# Test_set <- subset(Test_set, !(Test_set$LOCATION %in% Missing_TestTrain))
# 
# # Get the corresponding target variable (: To do this, do not remove TEMP at the build test/train sections)
# Target_Train <- Train_set$TEMP
# Target_Test <- Test_set$TEMP
# 
# # Drop the TEMP variable
# To_drop <- c("TEMP")
# Train_set <- Train_set[ ,!(names(Train_set) %in% To_drop)]
# Test_set <- Test_set[ ,!(names(Test_set) %in% To_drop)]

```

## Drop the LOCATION variable
```{r}
# Of course, you now still need to drop the LOCATION parameter
To_drop <- c("LOCATION")
Train_set <- Train_set[ ,!(names(Train_set) %in% To_drop)]
Test_set <- Test_set[ ,!(names(Test_set) %in% To_drop)]

```


## Change data format 

Add 10m to all ALT values to ensure there are no negative heights. 
Change all temperatures from Celsius to Kelvin.

**Note: here we also change the target variables to kelvin, even though this is probably not necessary**   

```{r}
# Add 10m to ALT for train and test data
Train_set$ALT <- Train_set$ALT + 10
Test_set$ALT <- Test_set$ALT + 10

# Convert TD, TL in Predictors to Kelvin

# For train
Train_set$TL <- Train_set$TL + 273.15
Train_set$TD <- Train_set$TD + 273.15

# For test
Test_set$TL <- Test_set$TL + 273.15
Test_set$TD <- Test_set$TD + 273.15

# Convert target variables to Kelvin
Target_Train <- Target_Train + 273.15
Target_Test <- Target_Test + 273.15

# Save these sets of predictors 
# save(Train_set, file = "Train_set_unProcessed.Rda")
# save(Test_set, file = "Test_set_unProcessed.Rda")
```

### Optional code: Test the distribution of predictors in test/train

The difference in summary statistics (min, quartile 1, median, mean, quartile 3, max) between the train/test sets is calculated below. 

```{r}
DiffSum <- function(Predictor){
  Summary1 <- summary(Train_set[, Predictor])
  Summary2 <- summary(Test_set[, Predictor])
  
  The_Diff <- Summary1 - Summary2
  
  return(The_Diff)
}

for (i in colnames(Train_set)){
  print(i)
  print(DiffSum(i))
}
```

# Optional code: Remove A LOT of the predictors

Also convert TD to RH 
[link](http://journals.ametsoc.org/doi/pdf/10.1175/BAMS-86-2-225)
```{r}
# Calculate the relative humidity using the formula in the link above. 
#Train_set$RH <- 100 - (5*(Train_set$TL - Train_set$TD))

# What are the lowest/highest TL/TD values?
#lengthTL <- Train_set$TL
#test_1 <- sort(Train_set$TL, decreasing = TRUE)

# Drop a MASSIVE amount of columns
To_drop <- c("TD",  "Rijbaan_3", "Rijbaan_4", "Rijbaan_5", "Rijbaan_6", "Rijbaan_7", "Rijstrook_0", 
             "Rijstrook_5", "Spitsstrook", "Oprit", "Afrit", "Brug_of_Viaduct_0", "Brug_of_Viaduct_1", "Brug_of_Viaduct_2")
Train_set <- Train_set[ ,!(names(Train_set) %in% To_drop)]
Test_set <- Test_set[ ,!(names(Test_set) %in% To_drop)]

```




# PreProcessing

Here we remove variables with zero variance ("zv"), perform a BoxCox transform to resolve skewness ("BoxCox"), center and scale the data ("center", "scale") and perform a principal component analysis ("pca"). 

We remove zero variance predictors because these predictors cannot help the ML algorithms predict the target variable. Centering and scaling the data will speed up the model building process. A pca will reduce the number of predictors and will resolve multicollinearity, thereby speeding up the model building process. In order to perform a pca data needs to be centered and scaled (which we have already done) and the data much not be skewed. With a BoxCox transform we resolve skewness. 

Dummy variables are included in the pca.

```{r}
library(caret)

# Build the object containing the transform information
xTrans <- preProcess(Train_set, method = c("zv", "center", "scale"),
                     na.remove = TRUE)

# Summary of transform information
print(xTrans)

# Number of PC's 
xTrans$numComp

# BoxCox transformed predictors
xTrans$method$BoxCox

# Transform train set
Train_set <- predict(xTrans, Train_set)

# Transform test set
Test_set <- predict(xTrans, Test_set)

```

The total number of principal components is 'r xTrans$numComp'. As the summary of xTrans shows, all predictors were centered, scaled and pca transformed. Only 'r length(xTrans$method$BoxCox)' predictors were BoxCox transformed, namely: 'r xTrans$method$BoxCox'. 


# Analysis of PC's

This section of the code contains an analysis of the PC's in the test and train sets. 
If you wish to cut down the run time of the script this section of the code can be commented/not run without changing the test/train sets. 

First, we look at the variable loadings. 
The rotation column of the xTrans object stores the variable loadings. 
Each principal component after the PCA is a linear combination of the original predictors.The coefficient for each predictor is called loading. A variable loading close to 0 indicates that a predictor did not contribute much to the principal component

```{r}
xTrans$rotation

```

Next, a scree plot is produced which shows how much of the total variance is explained by the principal components. To build the scree plot I adapted a method from: [link](https://www.analyticsvidhya.com/blog/2016/03/practical-guide-principal-component-analysis-python/) 
Important to note!!! By applying the method to the PC's that you get after transformation you can only see the percentage variance explained for the PC's that came out of the preProcessing prediction. In other words: Only the PC's that together explain 95% of the variance are shown. Unfortunately, you cannot extract the rest of the PC's when you use caret to perform a PCA. 

```{r}
# Scree plot of the explained variance

# First, we build a function that makes scree plots of the PC's
Plot_Scree <-function(The_set, TitlePlot){
  
  # Get the standard deviation per PC 
  PC_stdev <- apply(The_set, 2, sd, na.rm = TRUE)
  
  # Get the variance per PC by taking the square of the standard deviation per PC
  PC_var <- PC_stdev^2
  
  # Calculate proportion of total variance explained per PC by dividing PC_var by the total variance explained
  PC_prop_var <- PC_var/(sum(PC_var))
  
  # Make the scree plot
  plot(y = PC_prop_var, x = 1:length(The_set), type = "b", xlab = "Principal Components", ylab = "Proportion of variance explained", main = TitlePlot)
}

# Make a scree plot of the train set PC's 
Plot_Scree(Train_set, "_CatRemoved, train")

# Make a scree plot of the test set PC's 
Plot_Scree(Test_set, "_CatRemoved, test")


```

To get an idea of what these PC's look like you can plot them in a pairwise plot (Also called a scatterplot matrix).

The PC's show that outliers may throw off your pca. 
For more information on this topic, see: 
[link](http://www.math.umn.edu/~lerman/Meetings/SIAM2012_Sujay.pdf) **Note: Do we need to solve this? (Probably yes, but...)** 

Here we make scatterplot matrices for the train set, but the same can be done for the test set.

**Note: Code takes a while to run, so is left commented. Uncomment to plot.**

```{r}
# Scatterplot matrices for the train set, 5x5 plots
#pairs(Train_set[, 1:5])
#pairs(Train_set[ , 6:10]) 
#pairs(Train_set[ , 11:15])
#pairs(Train_set[ , 16:20])
#pairs(Train_set[ , 20:25])

# The large plots stored in the plotting device can slow down your PC.
# In order to empty the plotting device type: dev.off()
```

Next, to check if the pca indeed removed multicollinearity a correlation plot is made. 
If the pca went well the PC's should not be correlated. 

```{r}
library(corrplot)

# Correlation plot for the train set
corTrain_1 <- cor(Train_set, use = "complete.obs")
png(filename = "/usr/people/kleingel/Projects/MLProject/Corrplot_Train_3D_noPCA")
corrplot(corTrain_1, method = "circle")
dev.off()

# Correlation plot for the test set
corTest_1 <- cor(Test_set, use = "complete.obs")
png(filename = "/usr/people/kleingel/Projects/MLProject/Corrplot_Test_3D_noPCA")
corrplot(corTest_1, method = "circle")
dev.off()

``` 

Next, we test for (near) zero variance in the test/train sets
(This takes a while to run)
```{r}
# Test for (near) zero variance in the train set
NearZero_Train <- nearZeroVar(Train_set, names = TRUE, saveMetrics = TRUE)
head(NearZero_Train)

# As expected, the are no zero variance variables, because those are removed during preProcessing
any(NearZero_Train$zeroVar)
# After pca there are also near zero variables
any(NearZero_Train$nzv)


# Test for (near) zero variance in the test set
NearZero_Test <- nearZeroVar(Test_set, names = TRUE, saveMetrics = TRUE)
head(NearZero_Test)

# As expected, the are no zero variance variables, because those are removed during preProcessing
any(NearZero_Test$zeroVar)
# After pca there are also near zero variables
any(NearZero_Test$nzv)

```

### Add the target variable (road temperature) to the pre-processed test/train sets

```{r}
# Test if the length of the target variable is the same as that of the columns in test/train
length(Train_set$PC1) == length(Target_Train)
length(Test_set$PC1) == length(Target_Test)

# If so, we can now add the target variable back to the train/test sets
Train_set$TRoad <- Target_Train
Test_set$TRoad <- Target_Test
```


# Final step: Save sets

To use these sets in the scripts which build TW models we need to save the sets. 

```{r}
# The sets are saved to your working directory as Rdata files
# Your current working directory is:
getwd()
#setwd("C:/Users/Eva/Documents/Stage/ML_project")
#setwd("/usr/people/kleingel/Projects/MLProject")

# Save the test/train predictor sets 
save(Train_set, file = "Train_3D_noPCA.Rda")
save(Test_set, file = "Test_3D_noPCA.Rda")

# Save the corresponding target variables
save(Target_Train, file = "Train_Target_BIG_noM.Rda")
save(Target_Test, file = "Test_Target_BIG_noM.Rda")

```