Data_Analysis_Process_Code.Rmd

This report illustrates the codes and explanation of the data analysis process that I've done for Predictive Analysis Competition Project.

# Load Packages
Before I start doing anything, I loaded all the packages needed first, so that I won't accidentally load any packages for multiple times.
```{r}
#data exploration and tidying
library(dplyr)
library(tidyr)
library(stringr)
library(caTools)
library(caret)
library(janitor)
library(skimr)
library(ggcorrplot)
#intersect
library(base)
#create dummy variable
library(vtreat)
#regression tree
library(rpart); library(rpart.plot)
#bag ipred
library(ipred)
#random Forest
library(randomForest)
#ranger
library(ranger)
#gbm
library(gbm)
#xgboost
library(xgboost)
#box cox transformation and svm
library(e1071)
```

# Data Exploration

## Read Data
The first step is to set up the environment and import analysisData and scoringData.
```{r}
setwd("/Users/alicezhou/Documents/Columbia/5200 Applied Analytics Frameworks & Methods I/Notes/PAC")
data = read.csv("analysisData.csv")
sdata = read.csv(file="scoringData.csv")
```

## Variable Overview
Summary and comparison of two data sets are shown below. Here are some highlights and my thoughts for further analysis.

analysisData:
1.  There are 19 variables, with 1 identifier (id) and 1 outcome / dependent variable (rating). The other variables can all be the potential predictors / independent variables to predictive models.
2.  Outcome variable (rating) is continuous. Only models for continuous outcome variable should be considered, while logistics regression and classification trees can be excluded.
3.  Predictors include categorical and continuous variables.
4.  It's easier to explore the numeric and logical data types' characteristics by reading the overall summary. I can get a preliminary understanding, in terms of dispersion (standard deviation, interquartile range and range) for numeric variables and central tendency (mean) for both numeric and logical variables.
5.  The character data type has more uniqueness and is less structured than the others.
    -   song: It has 16,542 unique values, which is a similar amount to number of observations 19,485. Almost each song has its own name. This variable is too unique to be predictive. Therefore, I won't include it as the predictor for future modelling.
    -   genre: It has only 2,937 unique values. However, every value may contain several genres of a song. I am going to remove the special characters, and transform it to new categorical variables, so that each genre type can be represented by a category individually. Besides, since it contains missing value, I will fill the missing value as 0 before data transformation.
    -   performer: It has 6,687 unique values. Around 1/3 of songs have their unique performers. It's reasonable to be considered as a predictive variable now.

Comparison between analysisData and scoringData:
1.  Compared to the variables of analysisData, only outcome variable (rating) is missed for scoringData. All the predictors are in same data type. We can use all the predictors in analysisData to train the predictive model, to predict rating of songs in scoringData.
2.  There are significantly less rows for scoringData than analysisData. Some characters of genre and performer can be different from analysisData. Thus, I will only include the inter-sets of both data sets for these variables to data analysis.
```{r}
str(data)
```
```{r}
skim(data)
```
```{r}
str(sdata)
```
```{r}
skim(sdata)
```

# Data Tidying

## Encode Missing Data
```{r}
data[is.na(data)] = 0
```

## Data Parsing
For character type of data, since some variables contain observations with multiple delimited values, I am going to separate the values and places each one in its own row. Moreover, special character and unnecessary space have to be removed.
```{r}
data1 <- data

#character data - performer
##seperate a collapsed column into multiple rows, so that each row only indicates 1 performer
data1 <- data1 %>% 
  separate_rows(performer, sep = ",") 
##remove special character
data1$performer <- gsub("[^[:alnum:] ]","",data1$performer)

#character data - genre
##seperate a collapsed column into multiple rows, so that each row only indicates 1 genre
data1 <- data1 %>% 
  separate_rows(genre, sep = ",") 
##remove special character
data1$genre <- gsub("[^[:alnum:] ]","",data1$genre)
##remove both leading and trailing whitespace
data1$genre <- trimws(data1$genre)

head(data1)
```
As a result, we can see the observations are separated to multiple rows, and each row only indicates 1 genre and 1 performer. For example, the first observation (id: 94500), has been separated to 4 rows, because it represents 4 genres and 1 performer.

## Data Transformation
It's hard to directly input logical and character data into models, so I am going to change them to either numeric or factor data type. Furthermore, to prepare future data analysis (eg. Principal Components Analysis), I will change integer to numeric data.
```{r}
#data type transformation: logical to numeric
data1$track_explicit <- ifelse(data1$track_explicit==TRUE, 1, 0)
#data type transformation: character to factor
data1$performer <- as.factor(data1$performer)
data1$genre <- as.factor(data1$genre)
#data type transformation: integer to numeric
data1$rating <- as.numeric(data1$rating)
data1$key <- as.numeric(data1$key)
data1$mode <- as.numeric(data1$mode)
data1$time_signature <- as.numeric(data1$time_signature)

skim(data1)
```

Even though character data has been transformed to factor, it cannot be integrated over, summed over, or marginalized for further data analysis. Thus, I am going to create dummy variables so that these variables can be used in many machine learning models.
```{r}
#create dummy variable - genre 
data2 <- data1 %>% 
    mutate(n=1) %>%
    pivot_wider(names_from=genre, values_from=n, names_prefix="genre_", values_fill=list(n=0)) 
#create dummy variable - performer
data2 <- data2 %>% 
    mutate(n=1) %>%    
    pivot_wider(names_from=performer, values_from=n, names_prefix="performer_", values_fill=list(n=0))
```

To extract more specific data sets for further variable selection, I drop the variables (id & song) that has been excluded to data analysis.
```{r}
#delete unused variables
data3 <- data2 %>% select(-id) %>% select(-song)

str(data3)
```

Now, the analysisData contains same rows of observation as original, meanwhile it contains the dummy variables derived by performer and genre. Each dummy variable represents a category of either performer or genre.
I am going to perform the same data tidying process for scoringData, before we check the inter-sets of variables/columns between two data sets.
```{r}
#Encode Missing Data
sdata[is.na(sdata)] = 0

#Data Parsing
sdata1 <- sdata
#character data - performer
##seperate a collapsed column into multiple rows, so that each row only indicates 1 performer
sdata1 <- sdata1 %>% 
  separate_rows(performer, sep = ",") 
##remove special character
sdata1$performer <- gsub("[^[:alnum:] ]","",sdata1$performer)

#character data - genre
##seperate a collapsed column into multiple rows, so that each row only indicates 1 genre
sdata1 <- sdata1 %>% 
  separate_rows(genre, sep = ",") 
##remove special character
sdata1$genre <- gsub("[^[:alnum:] ]","",sdata1$genre)
##remove both leading and trailing whitespace
sdata1$genre <- trimws(sdata1$genre)

#Data Transformation
#data type transformation: logical to numeric
sdata1$track_explicit <- ifelse(sdata1$track_explicit==TRUE, 1, 0)
#data type transformation: character to factor
sdata1$performer <- as.factor(sdata1$performer)
sdata1$genre <- as.factor(sdata1$genre)
#data type transformation: integer to numeric
sdata1$key <- as.numeric(sdata1$key)
sdata1$mode <- as.numeric(sdata1$mode)
sdata1$time_signature <- as.numeric(sdata1$time_signature)

#create dummy variables - genre
sdata2 <- sdata1 %>% 
    mutate(n=1) %>%
    pivot_wider(names_from=genre, values_from=n, names_prefix="genre_", values_fill=list(n=0)) 
#create dummy variables - performer
sdata2 <- sdata2 %>% 
    mutate(n=1) %>%    
    pivot_wider(names_from=performer, values_from=n, names_prefix="performer_", values_fill=list(n=0))

#delete unused variables
sdata3 <- sdata2 %>% select(-id) %>% select(-song)

str(sdata3)
```

# Feature Selection 

## Variable Inter-set
From the summary of data3 and sdata3, we can see the numbers of variable are different (7,880 vs. 3,500). The reason is that dummy variables derived from the values of performer and genre can be different between two data sets. To make sure the variables to be analysed are aligned in both data sets, I am going to take their inter-sets of the column names to be the predictors.
```{r}
variable1 <- intersect(names(data3), names(sdata3)) 
data4 <- data3[c(variable1,"rating")] %>% relocate(rating, .before=track_duration)
sdata4 <- sdata3[variable1]

#tidy up the column names, especially the dummy variables derived from the value of performer and genre 
data4 <- clean_names(data4)
sdata4 <- clean_names(sdata4)

str(data4)
str(sdata4)
```

## Remove Near Zero Variance
The variables with zero variance or near zero variance have little predictable value, which should be removed from predictors.
```{r}
removeZeroVar <- data4[, sapply(data4, var) != 0]
removenearZeroVar <- nearZeroVar(removeZeroVar, names=TRUE, freqCut = 9999/1, uniqueCut = 1)
data5 <- data4[, setdiff(names(data4), removenearZeroVar)]

s_removeZeroVar <- sdata4[, sapply(sdata4, var) != 0]
s_removenearZeroVar <- nearZeroVar(s_removeZeroVar, names=TRUE, freqCut = 9999/1, uniqueCut = 1)
sdata5 <- sdata4[, setdiff(names(sdata4), s_removenearZeroVar)]

#repeat variable inter-set after removing near zero variance
variable2 <- intersect(names(data5), names(sdata5))
data6 <- data5[c(variable2,"rating")] %>% relocate(rating, .before=track_duration)
sdata6 <- sdata5[variable2]

str(data6)
str(sdata6)
```

# Split Data
I split data5 into train and test sample, so that I can build the model with train sample and check the model performance with test sample later.
```{r}
set.seed(617)
split=sample.split(data5$rating, SplitRatio = 0.7)
train=data6[split,]
test=data6[!split,]
```

To be safe, I remove zero variance and tidy up the column names again after splitting the data.
```{r}
train <- train[ , which(apply(train, 2, var) != 0)]
test <- test[ , which(apply(test, 2, var) != 0)]
train <- clean_names(train)
test <- clean_names(test)
str(train)
str(test)
```

Some dummy variables can be lost from removing zero variance. To make sure the variables to be analysed are aligned in both train and test samples, I am going to take their inter-sets of the column names to be the predictors.
```{r}
variable3 <- intersect(names(train), names(test))
train <- train[c(variable3)] %>% relocate(rating, .before=track_duration)
test <- test[c(variable3)] %>% relocate(rating, .before=track_duration)

#same treatment to scoringData
sdata6 <- clean_names(sdata6)
sdata6 <- sdata6[variable3[-1]]
#prepare a set of scoringData with id, to prepare the final step of result extraction 
sdata6_id <- clean_names(sdata2)[c(variable3[-1],"id")]

str(train)
str(test)
str(sdata6)
str(sdata6_id)
```

## Principal Components Analysis
Considering the large amount of predictors (1,664), it's heavily time-consuming for subset selection or skrinkage. It's also too complicated if we are going to use filter methods, and check the relencancy and non-redundancy among cariables one by one.
However, dimension reduction (Principal Components Analysis) would be an efficient feature selection approach here. With Principal Components Analysis ("PCA"), 1,664 predictors will be reduced to a smaller number (e.g., 70%) of components based on a measure of similarity (e.g., correlation). I am going to use the reduced number of components to predict the the outcome instead of the original set of predictors.
```{r}
trainPredictors = train[,-1]
pca = prcomp(trainPredictors,scale. = T)
train_components = data.frame(rating = train$rating, cbind(pca$x[,1:(1664*0.7)])) 

testPredictors = test[,-1]
test_pca = predict(pca,newdata=testPredictors)
test_components = data.frame(rating = test$rating, cbind(test_pca[,1:(1664*0.7)]))

#scoringdata
sdata6_pca = predict(pca,newdata=sdata5)
sdata6_components = data.frame(sdata5_pca[,1:(1664*0.7)])

str(train_components)
str(test_components)
str(sdata6_components)
```

# Data Analysis - Modeling
Now, everything is ready for predictive modeling. Some models usually have higher flexibility and accuracy (eg. Bagging, Boosting, Random Forest, Support Vector Machine), while the other models have higher interpretability (eg. Linear Regression). Since our goal is to improve the predictive accuracy, I prefer the former.
I will still try various models with default model parameters first. Each predictive model's RMSE (root-mean-square error) will be calculated to measure the model's accuracy for prediction. I will compare these models' RMSE, and then pick the models with lowest RMSE (best accuracy) for further parameters turning.

## Multiple Regression
Linear Multiple Regression gets a 15.25 test RMSE.
```{r}
lm = lm(rating~.,train_components)

pred_train_lm=predict(lm)
rmse_train_lm=sqrt(mean((pred_train_lm-train_components$rating)^2)); rmse_train_lm
pred_test_lm=predict(lm, newdata=test_components)
rmse_test_lm=sqrt(mean((pred_test_lm-test_components$rating)^2)); rmse_test_lm
```

## Regression Tree
Regression Tree gets a 15.33 test RMSE.
```{r}
model_tree = rpart(rating~., data=train_components, method = 'anova')

pred_train_tree=predict(model_tree)
rmse_train_tree=sqrt(mean((pred_train_tree-train_components$rating)^2)); rmse_train_tree
pred_test_tree=predict(model_tree, newdata=test_components)
rmse_test_tree=sqrt(mean((pred_test_tree-test_components$rating)^2)); rmse_test_tree
```

## Random Forest
Bagging (Random Forest) gets a 14.79 test RMSE.
```{r}
set.seed(1031)
rf = randomForest(rating~., 
                   data=train_components, 
                   mtry = 12,
                   ntree = 1000)
pred_train_rf = predict(rf)
rmse_train_rf = sqrt(mean((pred_train_rf - train_components$rating)^2)); rmse_train_rf
pred_test_rf = predict(rf, newdata=test_components)
rmse_test_rf = sqrt(mean((pred_test_rf - test_components$rating)^2)); rmse_test_rf
```

## Ranger
Random Forest (Ranger) gets a 14.75 test RMSE.
```{r}
set.seed(1031)
cv_forest_ranger = ranger(rating ~ .,
                          data=train_components,
                          num.trees = 1000)

#test ranger rmse
pred_train = predict(cv_forest_ranger, data = train_components, num.trees = 1000)
rmse_train_cv_forest_ranger = sqrt(mean((pred_train$predictions - train_components$rating)^2)); rmse_train_cv_forest_ranger
pred_test = predict(cv_forest_ranger, data = test_components, num.trees = 1000)
rmse_test_cv_forest_ranger = sqrt(mean((pred_test$predictions - test_components$rating)^2)); rmse_test_cv_forest_ranger
```

## XGBoost
Boosting (XGBoost) gets a 15.85 test RMSE.
```{r}
xgboost = xgboost(data=as.matrix(train_components[,-1]), 
                  label = train_components$rating,
                  nrounds=10000,
                  verbose = 0,
                  early_stopping_rounds = 100)
xgboost$best_iteration

#test xgboost rmse
pred_train = predict(xgboost, 
               newdata=as.matrix(train_components[,-1]))
rmse_train_xgboost = sqrt(mean((pred_train - train_components$rating)^2)); rmse_train_xgboost
pred_test = predict(xgboost, 
               newdata=as.matrix(test_components[,-1]))
rmse_test_xgboost = sqrt(mean((pred_test - test_components$rating)^2)); rmse_test_xgboost
```

## gbm
Boosting (gbm) gets a 14.66 test RMSE.
```{r}
set.seed(1031)
cvboost = gbm(rating ~ .,
              data=train_components,
              distribution="gaussian",
              n.trees=500)

#test gbm rmse
pred_train = predict(cvboost, n.trees = 500)
rmse_train_cv_boost = sqrt(mean((pred_train - train_components$rating)^2)); rmse_train_cv_boost
pred_test = predict(cvboost, newdata = test_components, n.trees = 500)
rmse_test_cv_boost = sqrt(mean((pred_test - test_components$rating)^2)); rmse_test_cv_boost
```

## Radial SVM
Support Vector Machine (Radial) gets a 15.13 test RMSE.
```{r}
svmRadial = svm(rating~.,data = train_components,kernel='radial')
summary(svmRadial)

#test svm rmse
pred_train = predict(svmRadial)
rmse_train_svm = sqrt(mean((pred_train - train_components$rating)^2)); rmse_train_svm
pred_test = predict(svmRadial,newdata=test_components)
rmse_test_svm = sqrt(mean((pred_test - test_components$rating)^2)); rmse_test_svm
```

## Comparison
Among all the models built above, I am going to compare their RMSE with visualization approach.
```{r}
data.frame(
  id = 1:7,
  model = c('multiple regression','regression tree','bag - randomForest','forest - ranger', 'xgboost','boost - gbm','radial svm' ),
  rmse_train = c(rmse_train_lm, rmse_train_tree, rmse_train_rf, rmse_train_cv_forest_ranger, rmse_train_xgboost, rmse_train_cv_boost, rmse_train_svm),
  rmse = c(rmse_test_lm, rmse_test_tree, rmse_test_rf, rmse_test_cv_forest_ranger, rmse_test_xgboost, rmse_test_cv_boost, rmse_test_svm))%>%
  rename('train RMSE' = rmse_train, 'test RMSE' = rmse)%>%
  pivot_longer(cols = 3:4,names_to = 'Sample', values_to = 'RMSE')%>%
  ggplot(aes(x=reorder(model,desc(id)), y = RMSE))+
  geom_col(fill = 'cadetblue')+
  xlab('')+
  coord_flip()+
  theme_bw()+
  facet_wrap(~Sample)
```

Interestingly, eventhough XGBoost get a very low train RMSE (0.37), the test RMSE is the highest (15.85). XGBoost derives the most extreme over-fitting issue here.
We can clearly see that Random Forest, Ranger, and gbm models work the best (get the lowest test RMSE) for test samples. Therefore, we pick these three models for further parameters tuning.

# Data Analysis - Model Tuning
I turned several times for Random Forest, Ranger, and gbm models. The best results and the critical steps for the iterative (decreasing test RMSE and improving predictive accuracy of models) process are shown below.

## Random Forest - Tuned Model Parameters
Tuned Bagging (Random Forest) gets a 14.75 test RMSE.
```{r}
set.seed(1031)
bag_rf = randomForest(rating~., 
                   data=train_components, 
                   mtry = 36,
                   ntree = 1000)
pred_train_rf = predict(bag_rf)
rmse_train_rf = sqrt(mean((pred_train_rf - train_components$rating)^2)); rmse_train_rf
pred_test_rf = predict(bag_rf, newdata=test_components)
rmse_test_rf = sqrt(mean((pred_test_rf - test_components$rating)^2)); rmse_test_rf
```

## Ranger - Tuned Model Parameters
Tuned Random Forest (Ranger) gets a 14.65 test RMSE.
```{r}
set.seed(627143)
cv_forest_ranger = ranger(rating ~ .,
                          data=train_components,
                          num.trees = 1000, 
                          mtry=36,
                          min.node.size = 30,
                          splitrule = "extratrees")

#test ranger rmse
pred_train = predict(cv_forest_ranger, data = train_components, num.trees = 1000)
rmse_train_cv_forest_ranger = sqrt(mean((pred_train$predictions - train_components$rating)^2)); rmse_train_cv_forest_ranger
pred_test = predict(cv_forest_ranger, data = test_components, num.trees = 1000)
rmse_test_cv_forest_ranger = sqrt(mean((pred_test$predictions - test_components$rating)^2)); rmse_test_cv_forest_ranger
```

## gbm - Tuned Model Parameters
Tuned Boosting (gbm) gets a 14.53 test RMSE.
```{r}
set.seed(1031)
cvboost = gbm(rating ~ .,
              data=train_components,
              distribution="gaussian",
              n.trees=1000,
              interaction.depth=10,
              shrinkage=0.02,
              n.minobsinnode = 1)

#test gbm rmse
pred_train = predict(cvboost, n.trees = 1000)
rmse_train_cv_boost = sqrt(mean((pred_train - train_components$rating)^2)); rmse_train_cv_boost
pred_test = predict(cvboost, newdata = test_components, n.trees = 1000)
rmse_test_cv_boost = sqrt(mean((pred_test - test_components$rating)^2)); rmse_test_cv_boost
```

## gbm - Changed PCA Parameters
I changed PCA to higher variation (80%) of predictor selection, and the tuned Boosting (gbm) gets a 14.51 test RMSE.
```{r}
#feature selection - pca changed to 80%
#analysisdata
trainPredictors = train[,-1]
pca = prcomp(trainPredictors,scale. = T)
train_components = data.frame(rating = train$rating, cbind(pca$x[,1:1664*0.8])) 

testPredictors = test[,-1]
test_pca = predict(pca,newdata=testPredictors)
test_components = data.frame(rating = test$rating, cbind(test_pca[,1:1664*0.8]))

#scoringdata
sdata6_pca = predict(pca,newdata=sdata6)
sdata6_components = data.frame(sdata6_pca[,1:1664*0.8])


#gbm
set.seed(1031)
cvboost = gbm(rating ~ .,
              data=train_components,
              distribution="gaussian",
              n.trees=1000,
              interaction.depth=10,
              shrinkage=0.02,
              n.minobsinnode = 1)

#test gbm rmse
pred_train = predict(cvboost, n.trees = 1000)
rmse_train_cv_boost = sqrt(mean((pred_train - train_components$rating)^2)); rmse_train_cv_boost
pred_test = predict(cvboost, newdata = test_components, n.trees = 1000)
rmse_test_cv_boost = sqrt(mean((pred_test - test_components$rating)^2)); rmse_test_cv_boost
```

## gbm - Changed Seed
I changed PCA to higher variation (80%) of predictor selection, and then changed seed before training the model. Furthermore, I tuned the model parameters again.
Tuned Boosting (gbm) gets a 14.44 test RMSE, which is the best accuracy I get from the models.
```{r}
set.seed(627143)
cvboost = gbm(rating ~ .,
              data=train_components,
              distribution="gaussian",
              n.trees=500,
              interaction.depth=12,
              shrinkage=0.013,
              n.minobsinnode = 1)

#test gbm rmse
pred_train = predict(cvboost, n.trees = 500)
rmse_train_cv_boost = sqrt(mean((pred_train - train_components$rating)^2)); rmse_train_cv_boost
pred_test = predict(cvboost, newdata = test_components, n.trees = 500)
rmse_test_cv_boost = sqrt(mean((pred_test - test_components$rating)^2)); rmse_test_cv_boost
```

# Prediction
Finally, I fit this model into scoringData, to generate the prediction data set for final submission.
```{r}
pred = predict(cvboost, newdata = sdata6_components, n.trees = 500)
submissionFile=data.frame(id=sdata6_id$id,rating=pred)
write.csv(submissionFile, file="submission_17_20.csv", row.names=F)
```