Machine Prediction for Human Resources Analytics

Introduction

Here I work on kaggle data set "Human Resources Analytics", which can be found by url
https://www.kaggle.com/ludobenistant/hr-analytics.
We're supposed to predict which valuable employees will leave next. Fields in the data set include:

• Employee satisfaction level
• Last evaluation
• Number of projects
• Average monthly hours
• Time spent at the company
• Whether they have had a work accident
• Whether they have had a promotion in the last 5 years
• Department
• Salary
• Whether the employee has left

Data Load

I start with loading data set and looking at its properties: dimensions, column names and types of data.

dt=read.csv("HR_comma_sep.csv", stringsAsFactors=F)
str(dt)

## 'data.frame': 14999 obs. of  10 variables:
##  $ satisfaction_level   : num  0.38 0.8 0.11 0.72 0.37 0.41 0.1 0.92 0.89 0.42 ...
##  $ last_evaluation      : num  0.53 0.86 0.88 0.87 0.52 0.5 0.77 0.85 1 0.53 ...
##  $ number_project       : int  2 5 7 5 2 2 6 5 5 2 ...
##  $ average_montly_hours : int  157 262 272 223 159 153 247 259 224 142 ...
##  $ time_spend_company   : int  3 6 4 5 3 3 4 5 5 3 ...
##  $ Work_accident        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ left                 : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ promotion_last_5years: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ sales                : chr  "sales" "sales" "sales" "sales" ...
##  $ salary               : chr  "low" "medium" "medium" "low" ...

sapply(dt, function(x) sum(is.na(x)))

##    satisfaction_level       last_evaluation        number_project 
##                     0                     0                     0 
##  average_montly_hours    time_spend_company         Work_accident 
##                     0                     0                     0 
##                  left promotion_last_5years                 sales 
##                     0                     0                     0 
##                salary 
##                     0

As we see all data except for the last 2 columns are numeric. In addition, here are no missing values.
We can check pairwise plots, correlations and densities.

library(ggplot2)
library(GGally)
ggpairs(data=dt, columns=1:8,
    mapping = aes(color = "navy"),
    axisLabels="show")

Now let us look at what the values in last 2 columns.

unique(dt$salary)

## [1] "low"    "medium" "high"

unique(dt$sales)

##  [1] "sales"       "accounting"  "hr"          "technical"   "support"    
##  [6] "management"  "IT"          "product_mng" "marketing"   "RandD"

Here are 3 classes for salary values and it looks like the column "sales" represents departments. Let us see how uniform is the distribution of the data:

table(dt$sales)

## 
##  accounting          hr          IT  management   marketing product_mng 
##         767         739        1227         630         858         902 
##       RandD       sales     support   technical 
##         787        4140        2229        2720

table(dt$salary)

## 
##   high    low medium 
##   1237   7316   6446

It is not very uniform, but at least each class is not too small. I would like to replace the columns with dummy variable columns. As first salary values:

library(dummies)
dumdt=dummy(dt$salary)
dumdt=as.data.frame(dumdt)
names(dumdt)

## [1] "MachinePredictionForNorthropGrumman.Rhtmlhigh"  
## [2] "MachinePredictionForNorthropGrumman.Rhtmllow"   
## [3] "MachinePredictionForNorthropGrumman.Rhtmlmedium"

Firstly, these names are not good because they are too long. Secondly, I do not need all of them, because they are correlated: having a low salary means not having high or medium. I will remove the low salary column and I will rename the rest of them.

dumdt=dumdt[,-2]
names(dumdt)=c("high_salary", "medium_salary")

Now I will attach my new variables to existing data frame.

dt=cbind(dt,dumdt)
dt$salary=NULL
str(dt)

## 'data.frame': 14999 obs. of  11 variables:
##  $ satisfaction_level   : num  0.38 0.8 0.11 0.72 0.37 0.41 0.1 0.92 0.89 0.42 ...
##  $ last_evaluation      : num  0.53 0.86 0.88 0.87 0.52 0.5 0.77 0.85 1 0.53 ...
##  $ number_project       : int  2 5 7 5 2 2 6 5 5 2 ...
##  $ average_montly_hours : int  157 262 272 223 159 153 247 259 224 142 ...
##  $ time_spend_company   : int  3 6 4 5 3 3 4 5 5 3 ...
##  $ Work_accident        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ left                 : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ promotion_last_5years: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ sales                : chr  "sales" "sales" "sales" "sales" ...
##  $ high_salary          : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ medium_salary        : int  0 1 1 0 0 0 0 0 0 0 ...

dumdt=dummy(dt$sales)
dumdt=as.data.frame(dumdt)
names(dumdt)

##  [1] "MachinePredictionForNorthropGrumman.Rhtmlaccounting" 
##  [2] "MachinePredictionForNorthropGrumman.Rhtmlhr"         
##  [3] "MachinePredictionForNorthropGrumman.RhtmlIT"         
##  [4] "MachinePredictionForNorthropGrumman.Rhtmlmanagement" 
##  [5] "MachinePredictionForNorthropGrumman.Rhtmlmarketing"  
##  [6] "MachinePredictionForNorthropGrumman.Rhtmlproduct_mng"
##  [7] "MachinePredictionForNorthropGrumman.RhtmlRandD"      
##  [8] "MachinePredictionForNorthropGrumman.Rhtmlsales"      
##  [9] "MachinePredictionForNorthropGrumman.Rhtmlsupport"    
## [10] "MachinePredictionForNorthropGrumman.Rhtmltechnical"

Clearly I need to go through the similar steps.

dumdt=dumdt[,-10]
names(dumdt)=c("accounting","hr","IT","management","marketing","product_mng","RandD","sales","support")
dt=cbind(dt,dumdt)
dt$sales=NULL
# Look at the new data frame:
str(dt)

## 'data.frame': 14999 obs. of  19 variables:
##  $ satisfaction_level   : num  0.38 0.8 0.11 0.72 0.37 0.41 0.1 0.92 0.89 0.42 ...
##  $ last_evaluation      : num  0.53 0.86 0.88 0.87 0.52 0.5 0.77 0.85 1 0.53 ...
##  $ number_project       : int  2 5 7 5 2 2 6 5 5 2 ...
##  $ average_montly_hours : int  157 262 272 223 159 153 247 259 224 142 ...
##  $ time_spend_company   : int  3 6 4 5 3 3 4 5 5 3 ...
##  $ Work_accident        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ left                 : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ promotion_last_5years: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ high_salary          : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ medium_salary        : int  0 1 1 0 0 0 0 0 0 0 ...
##  $ accounting           : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ hr                   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ IT                   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ management           : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ marketing            : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ product_mng          : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ RandD                : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ sales                : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ support              : int  0 0 0 0 0 0 0 0 0 0 ...

rm(dumdt)

Machine Learning

Now my data is ready for prediction. I split the data set into 2 sets: for training and for testing.

indeces=sample(1:dim(dt)[1], round(.7*dim(dt)[1]))
train=dt[indeces,]
test=dt[-indeces,]

str(train)

## 'data.frame': 10499 obs. of  19 variables:
##  $ satisfaction_level   : num  0.71 0.73 0.74 0.41 0.1 0.88 0.71 0.76 0.85 0.93 ...
##  $ last_evaluation      : num  0.5 0.6 0.55 0.49 0.94 1 0.92 0.8 0.65 0.65 ...
##  $ number_project       : int  4 3 6 2 6 5 3 4 4 4 ...
##  $ average_montly_hours : int  253 137 130 130 255 219 202 226 280 212 ...
##  $ time_spend_company   : int  3 3 2 3 4 6 4 5 3 4 ...
##  $ Work_accident        : int  0 0 0 0 0 1 0 0 1 0 ...
##  $ left                 : int  0 0 0 1 1 1 0 0 0 0 ...
##  $ promotion_last_5years: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ high_salary          : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ medium_salary        : int  1 1 1 0 0 0 0 0 1 1 ...
##  $ accounting           : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ hr                   : int  0 1 0 0 0 0 0 0 0 0 ...
##  $ IT                   : int  0 0 0 0 0 0 0 0 1 1 ...
##  $ management           : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ marketing            : int  0 0 0 1 0 0 0 0 0 0 ...
##  $ product_mng          : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ RandD                : int  1 0 0 0 0 0 0 0 0 0 ...
##  $ sales                : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ support              : int  0 0 1 0 0 0 0 1 0 0 ...

And now we can try a few predicting algorithms.

Logistic Regression

glm_mod=glm(left~., data=train, family=binomial)
options(width=120)
 summary(glm_mod)

## 
## Call:
## glm(formula = left ~ ., family = binomial, data = train)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -2.252360  -0.669234  -0.404160  -0.118092   3.005301  
## 
## Coefficients:
##                           Estimate   Std. Error   z value   Pr(>|z|)    
## (Intercept)            0.513671571  0.152904509   3.35943 0.00078104 ***
## satisfaction_level    -4.102569262  0.116375474 -35.25287 < 2.22e-16 ***
## last_evaluation        0.636893223  0.177888069   3.58030 0.00034320 ***
## number_project        -0.315622781  0.025421761 -12.41546 < 2.22e-16 ***
## average_montly_hours   0.004809019  0.000613495   7.83872 4.5515e-15 ***
## time_spend_company     0.270089407  0.018619140  14.50601 < 2.22e-16 ***
## Work_accident         -1.489204860  0.104868336 -14.20071 < 2.22e-16 ***
## promotion_last_5years -1.268547473  0.277706609  -4.56794 4.9254e-06 ***
## high_salary           -1.903687581  0.153808902 -12.37697 < 2.22e-16 ***
## medium_salary         -0.485601493  0.054435288  -8.92071 < 2.22e-16 ***
## accounting            -0.216824268  0.130317943  -1.66381 0.09615045 .  
## hr                     0.177414039  0.129105528   1.37418 0.16938628    
## IT                    -0.271998822  0.110571123  -2.45994 0.01389585 *  
## management            -0.523601663  0.162647148  -3.21925 0.00128527 ** 
## marketing             -0.098831483  0.127918381  -0.77261 0.43975108    
## product_mng           -0.272539608  0.124254440  -2.19340 0.02827862 *  
## RandD                 -0.599234520  0.144357350  -4.15105 3.3095e-05 ***
## sales                 -0.129873866  0.078066780  -1.66363 0.09618735 .  
## support               -0.048540168  0.089959935  -0.53958 0.58948988    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 11544.393  on 10498  degrees of freedom
## Residual deviance:  9035.048  on 10480  degrees of freedom
## AIC: 9073.048
## 
## Number of Fisher Scoring iterations: 5

In this model being in a specific department is not a reliable predictor except for being in management or R and D. Let us see what is suitable threshold for decision making.

fittedValues=fitted(glm_mod)
# How well is predicted people who stayed
hist(fittedValues[train$left==0])

# And the ones who left.
hist(fittedValues[train$left==1])

Even on train data we do not get good predictions! Nevertheless we can check it on our test data.

glm_predictions=predict(glm_mod,test[,-match("left", names(test))])
# How well is predicted people who stayed, colored red
hist(glm_predictions[test$left==0], breaks =15, col=2 )
# And the ones who stayed, colored blue
hist(glm_predictions[test$left==1],breaks =15,  col=rgb(0,0.7,0.8,1/2),add=T)

Taking into account that to define whose who left is more important it looks like a threshold may be about -1.

One of the best random forest implementations

I was looking for ways to run a parallel version of random forest. Of course we can do it with the widely known "randomForest", "caret" and "doParallel" packages. But my attempt to apply it to Digit Recognizer problem on Kaggle was not satisfactory: when I tried their whole data set my R instance was crushed, and when I took less than a quarter of it then it was running for too long. Here are my code and my notes, assuming that the data set is already here as a data frame "df":

library(caret); library(e1071)
library(kernlab); library(doParallel)
library(randomForest)

partial_1=df[1:10000, ]

# All my 8 cores are detectable.  I will take 7 for calculations.

cl=makeCluster(detectCores()-1)
registerDoParallel(cl)
t=Sys.time()
rf_1=train(label~., data=partial_1, method="rf" )
Sys.time()-t
stopCluster(cl)

# Time difference of 3.000682 hours
# Accuracy 0.9463709
# Accuracy was used to select the optimal model using  the largest value.
# The final value used for the model was mtry = 39.

I watched how much RAM has been used and it was a bit under 9 Gb. I've googled if there are any improvements to this, and turned out they are, in new packages: better memory usage and optimized algorithms. In particular, a package "ranger" is in a native for R "cran" repository. It has a built-in parallel option.  I tried it on the whole data set with the same options as before, with a number of trees equals to 500 (default option for "randomForest" and "ranger"), "mtry=39" (from the caret model above) and 7 cores ("num.threads=7"), and it was a breeze. Computation time was less than 1 min 15 sec,  accuracy was higher than 94.6% as well. The used RAM was a little bit above 3 Gb, which was pleasant, too. A default option for "verbose" is TRUE, so if you want to enjoy watching  computations, you can do it. A formula was written beforehand as "frml".  You can see my commands below.

library(ranger)
t=Sys.time()
frml=as.formula(paste0("label~", paste(names(dt)[-1], sep="", collapse='+')))
ranger(frml, data=df, mtry=39, num.threads=7, verbose=FALSE, classification=TRUE)
Sys.time()-t
# Time difference of 1.240249 mins
# prediction error:         3.20 % 

After this I've explored different values of "mtry", since it was so easy, and changed a number of trees to 5000. My final choice for "mtry" was 32. Peaks for amount of RAM were between 9 and a bit above 10 Gb. Times used were a around 10-11 minutes for a model.
There are 2 options for prediction with the package. First, you can do it in one move with the "ranger" function. Or second, you should use their option "write.forest" for the "ranger" function, and then you can use "predict" function as usual. My last steps follow, assuming that the required test data frame is here under name "testdf":

t=Sys.time()
rf_mod=ranger(frml, data=df, mtry = 32, num.threads=7, 
               num.trees =5000, verbose= FALSE,
               write.forest=T, classification=TRUE)
Sys.time()-t
# Time difference of 10.48414 mins
# prediction error:         3.06 %

Computing final predictions and writing a file for kaggle.com submission as a required table in csv format:

Label=predict(rf_mod, testdf)
Label=Label$predictions
submission=data.frame(id=testdf$id, outcome=Label)
write.csv(submission,"submission.csv", row.names=F)

Improving performance of random forests for a particular value of outcome by selecting better features

Summary

Choosing features to improve a performance of a particular algorithm is a difficult question. Currently here is PCA, which is hard to understand (although it can be used out-of-the-box), is not easy to interpret and requires centralizing and scaling of features. In addition, it does not allow to improve prediction performance for a particular outcome (if its accuracy is lower than for others or it has a particular importance). My method enables to use features without preprocessing. Therefore a resulting prediction is easy to explain. Plus it can be used to improve a accuracy prediction of a specified outcome value. It based on comparison of feature densities and has a good visual interpretation, which does not require thorough knowledge of linear algebra or calculus.

Application Example

Here is a worked out example of Choosing features for random forests algorithm with R code. It is supplemented with choosing additional features to improve prediction for a particular value of outcome. The method for comparing densities is described in detail there: Computing Ratio of Areas.

I will use for my computations data for Human Activity Recognition. The short description of the problem follows:

Human Activity Recognition - HAR - is a supervised classification problem, whose training data is obtained via an experiment having human subjects perform different activities. It focuses on the data collection from multiple accelerometers and possibly other sensors placed on a human body. There are many potential applications for HAR, like: elderly monitoring, life log systems for monitoring energy expenditure and for supporting weight-loss programs, and digital assistants for weight lifting exercises.

Data

The training data for this project are available here:

https://d396qusza40orc.cloudfront.net/predmachlearn/pml-training.csv

The test data are available here: https://d396qusza40orc.cloudfront.net/predmachlearn/pml-testing.csv

The data for this project come from this source: http://groupware.les.inf.puc-rio.br/har.

In this project we use data from accelerometers on the belt, forearm, arm, and dumbbell of 6 participants. They were asked to perform barbell lifts correctly (marked as “A”) and incorrectly in 4 different ways.

A - exactly according to the specification

B - throwing the elbows to the front

C - lifting the dumbbell only halfway

D - lowering the dumbbell only halfway

E - throwing the hips to the front

An outcome column with the letters is called “classe”. Read more, with pictures: http://groupware.les.inf.puc-rio.br/har#ixzz3jfosRia6

The goal of the project is to predict the manner in which participants did the exercise. In particular, we are to figure out which features to use to reach our goal.

I will load only the “training” file, because in the “testing” file there are no marks for exercise correctness.

training=read.csv("pml-training.csv", stringsAsFactors=F)

Now we can look at the “training” file more closely. It is easy to establish that first columns contain names of subjects, days and times of recording and other classifiers which do not represent physical movements. For example the very first “X” column is used to enumerate rows and the last column “classe” contains the letters which mark the performance quality. We can look at the data table dimensions, first 10 column names, first 20 values of first column and values of last “classe” column.

dim(training)

## [1] 19622   160

names(training)[1:10]

##  [1] "X"                    "user_name"            "raw_timestamp_part_1"
##  [4] "raw_timestamp_part_2" "cvtd_timestamp"       "new_window"          
##  [7] "num_window"           "roll_belt"            "pitch_belt"          
## [10] "yaw_belt"

training$X[1:20]

##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

unique(training$classe)

## [1] "A" "B" "C" "D" "E"

We see users’ names and so on. They should be removed for a predicting model.

Cleaning data

The next step is less obvious: by checking out the “training” data I’ve discovered a lot of features which have mostly undefined values. Some of them are read as logical with nonexistent values and some as characters. We can check the amount of columns with undefined values and the amount of numeric features in the test set and see that the number of useful features is fewer than 40%.

sum(colSums(is.na(training)) !=0)

## [1] 67

sum(sapply(training, is.numeric))

## [1] 123

I will remove non numeric features from the training set for my work together with the first 7 columns. In addition some of the columns with sparse data have mostly undefined values (NA), and we should get rid of them as well. The new data set with numeric features will be called “trai”. Note that the classifier column “classe” is removed after first command line as one of character columns, and I need to put it back.

prVec=sapply(training, is.numeric)
prVec[1:7]=F
trai=training[ , prVec]
trai <- trai[, colSums(is.na(trai)) == 0] 
trai$classe=as.factor(training$classe)
dim(trai)

## [1] 19622    53

We are left with 53 columns. The last is our outcome. Now let us split the data for cross-validation. I will load required packages, choose a number for random sampling and then put 70% of the data frame into a training set and the rest of it into a validation set.

library(caret); library(e1071)

## Loading required package: lattice

## Loading required package: ggplot2

set.seed(11)
inTrain = createDataPartition(y=trai$classe, p=0.7, list=F)
dataFrame=trai[inTrain,]
valdTrai=trai[-inTrain, ]

Choosing features

We need to find ways to distinguish types of performance labeled by letters from each other.

Here are 52 features in our data frame. It is still a lot for my laptop, say unparalleled “random forests” algorithm took 2 hours. Granted my laptop is not a powerful one, but what if our prediction model is supposed to work across different gadgets? It makes sense to make it less dependable on resources and to reduce the number of features. Let us visualize the results and figure out features which show more difference from others.

library(ggplot2)
library(Rmisc)

## Loading required package: plyr

p1=qplot(dataFrame[, 1], color=classe, data=dataFrame, geom="density", xlab="First feature")
p2=qplot(dataFrame[, 2], color=classe, data=dataFrame, geom="density", xlab="Second feature")
p3=qplot(dataFrame[, 3], color=classe, data=dataFrame, geom="density", xlab="Third feature")
p4=qplot(dataFrame[, 4], color=classe, data=dataFrame, geom="density", xlab="Fourth feature")
p5=qplot(dataFrame[, 5], color=classe, data=dataFrame, geom="density", xlab="Fifth feature")
p6=qplot(dataFrame[, 6], color=classe, data=dataFrame, geom="density", xlab="Sixth feature")
multiplot(p1,p2,p3,p4,p5,p6, cols=2)

## removing objects I do not need anymore
rm("p1","p2","p3","p4","p5","p6")

We can see that data behaviors are complicated. Some of features are bimodal and even multimodal. These characteristics could be caused by participants’ different sizes or training levels, but we do not have enough information to check it out. We can pick up features which diverge the most on density graphs by considering difference of densities. The “density” function is the same which is used in R for graphing. My rule for picking up a feature for final prediction data frame is following: \[ \frac{\textrm{an area between two density curves}}{\textrm{an area under one of the curves}} > 0.75 \] When it is true for at least one of pairs, then the feature is picked for prediction. The method itself is described in detail in my post: Computing Ratio of Areas. The threshold .75 is picked up by trial and error.

My code for it follows.

nmsPred = dim(dataFrame)[2]
classe = unique(dataFrame[, "classe"])
densities=data.frame(rep(0,1000),rep(0,1000),rep(0,1000),rep(0,1000),rep(0,1000))
names(densities)=classe
prefPred = rep(FALSE, nmsPred )
prefPred[53] = TRUE # to keep the "classe" variate
set.seed(11)
for (k in 1:(nmsPred-1) ) {
  lower.limit=min(dataFrame[ , k])
  upper.limit=max(dataFrame[ , k])
     for (lett in classe) {
        den <-density(dataFrame[dataFrame$classe==lett, k],
        kernel="gaussian", n=1000, 
        from = lower.limit, to = upper.limit)
        densities[, lett] <- den$y
     }
     ind=FALSE
     for (l1 in 1:length(classe)) {
        for (l2 in classe[-l1]) {
        ind0 <- ((sum(abs(densities[,l1]-densities[,l2]))/sum(densities[,l1]))>.75)
        ind = (ind | ind0)
              }
       }
           prefPred[k] <- prefPred[k] | (ind)
}
## Gathering the features together:
workdf=dataFrame[ , prefPred]

Random Forests Algorithm Results and Cross-Validation

I will use “random forests” algorithm option from “caret” package. Since it uses random numbers, I need to set up my “random number seed”. Then we can look at the model and its accuracy:

library(caret); library(e1071)
library(randomForest)
library(kernlab); library(doParallel)
set.seed(11)
cl <- makeCluster(detectCores())
registerDoParallel(cl)
modrf=train(classe~., data=workdf, method="rf", allowParallel=T)
stopCluster(cl)
registerDoSEQ()
# Looking at the resulting model
modrf

## Random Forest 
## 
## 13737 samples
##    19 predictors
##     5 classes: 'A', 'B', 'C', 'D', 'E' 
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 13737, 13737, 13737, 13737, 13737, 13737, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##    2    0.9855892  0.9817652
##   10    0.9853784  0.9814997
##   19    0.9728235  0.9656130
## 
## Accuracy was used to select the optimal model using  the largest value.
## The final value used for the model was mtry = 2.

Then we compare predictions on the validation set:

# Forming the testing subset:
valddf=valdTrai[, prefPred]
last=length(prefPred)
# Computing predictions and checking accuracies:
predictRFonVald=predict(modrf, valddf[, -last])
confusionMatrix(valdTrai$classe, predictRFonVald)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    A    B    C    D    E
##          A 1671    2    0    0    1
##          B    6 1130    3    0    0
##          C    0   14 1006    5    1
##          D    0    0   13  950    1
##          E    0    0    4    1 1077
## 
## Overall Statistics
##                                           
##                Accuracy : 0.9913          
##                  95% CI : (0.9886, 0.9935)
##     No Information Rate : 0.285           
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.989           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: A Class: B Class: C Class: D Class: E
## Sensitivity            0.9964   0.9860   0.9805   0.9937   0.9972
## Specificity            0.9993   0.9981   0.9959   0.9972   0.9990
## Pos Pred Value         0.9982   0.9921   0.9805   0.9855   0.9954
## Neg Pred Value         0.9986   0.9966   0.9959   0.9988   0.9994
## Prevalence             0.2850   0.1947   0.1743   0.1624   0.1835
## Detection Rate         0.2839   0.1920   0.1709   0.1614   0.1830
## Detection Prevalence   0.2845   0.1935   0.1743   0.1638   0.1839
## Balanced Accuracy      0.9979   0.9921   0.9882   0.9954   0.9981

As we see on the model data and in the table the resulting accuracy is around 99%, which is quite good. In the same time we have the lowest accuracy for a outcome value C: it is 98.8%, and we may want to increase it to make accuracy for different outcomes more uniform. Let us try to add some features to improve it by using the above method, but only for one outcome value. A threshold in this case is picked up to be 0.65.

for (k in 1:(nmsPred-1) ) {
  lower.limit=min(dataFrame[ , k])
  upper.limit=max(dataFrame[ , k])
     for (lett in classe) {
        den <-density(dataFrame[dataFrame$classe==lett, k],
        kernel="gaussian", n=1000, 
        from = lower.limit, to = upper.limit)
        densities[, lett] <- den$y
     }
     ind=FALSE
     for (l1 in classe[-3]) { # working with 3d classifier, "C""
         ind0 <- ((sum(abs(densities[,3]-densities[,l1]))/sum(densities[,3]))>.65)
         ind = (ind | ind0)
     }
           prefPred[k] <- prefPred[k] | (ind)
}
## Gathering the features together:
workdf=dataFrame[ , prefPred]

Running the random forests algorithm again:

set.seed(11)
cl <- makeCluster(detectCores())
registerDoParallel(cl)
modrf=train(classe~., data=workdf, method="rf", allowParallel=T)
stopCluster(cl)
registerDoSEQ()
# Looking at the resulting model
modrf

## Random Forest 
## 
## 13737 samples
##    23 predictors
##     5 classes: 'A', 'B', 'C', 'D', 'E' 
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 13737, 13737, 13737, 13737, 13737, 13737, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa      Accuracy SD  Kappa SD   
##    2    0.9845870  0.9804951  0.001692352  0.002145187
##   12    0.9855843  0.9817603  0.001927464  0.002441153
##   23    0.9732220  0.9661184  0.004619732  0.005841083
## 
## Accuracy was used to select the optimal model using  the largest value.
## The final value used for the model was mtry = 12.

Again comparing results on our validation set we can see that accuracy for the value C is improved.

# Forming the testing subset:
valddf=valdTrai[, prefPred]
last=length(prefPred)
# Computing predictions and checking accuracies:
predictRFonVald=predict(modrf, valddf[, -last])
confusionMatrix(valdTrai$classe, predictRFonVald)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    A    B    C    D    E
##          A 1671    2    0    0    1
##          B    7 1126    4    2    0
##          C    0   12 1010    4    0
##          D    0    0    5  957    2
##          E    0    0    6    0 1076
## 
## Overall Statistics
##                                           
##                Accuracy : 0.9924          
##                  95% CI : (0.9898, 0.9944)
##     No Information Rate : 0.2851          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.9903          
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: A Class: B Class: C Class: D Class: E
## Sensitivity            0.9958   0.9877   0.9854   0.9938   0.9972
## Specificity            0.9993   0.9973   0.9967   0.9986   0.9988
## Pos Pred Value         0.9982   0.9886   0.9844   0.9927   0.9945
## Neg Pred Value         0.9983   0.9971   0.9969   0.9988   0.9994
## Prevalence             0.2851   0.1937   0.1742   0.1636   0.1833
## Detection Rate         0.2839   0.1913   0.1716   0.1626   0.1828
## Detection Prevalence   0.2845   0.1935   0.1743   0.1638   0.1839
## Balanced Accuracy      0.9976   0.9925   0.9910   0.9962   0.9980

Accuracy for other outcome values is slightly changed, too, but not much.

Principal Component Analysis

The principal component analysis is a standard method for reducing number of features. We can compare results with the method used above and ascertain how many components we would be required to use by the method for accuracy of at least 99%.

prep=prcomp(dataFrame[-53], scale=T)
v99<-summary(prep)$importance[3,]>.988
summary(prep)$importance[3, v99]

##    PC35    PC36    PC37    PC38    PC39    PC40    PC41    PC42    PC43 
## 0.98948 0.99102 0.99223 0.99329 0.99434 0.99512 0.99582 0.99649 0.99709 
##    PC44    PC45    PC46    PC47    PC48    PC49    PC50    PC51    PC52 
## 0.99762 0.99814 0.99865 0.99905 0.99943 0.99967 0.99984 0.99996 1.00000

As we see to get the same accuracy as for the random forest model the principal component method advises to use 36 or 37 principal components (PC36 or PC37), which are prepossessed features calculated from original 52 ones. In addition it does not tell us how to improve prediction performance for one of outcome values if we would want it.

Other Applications

The method can be used to evaluate consistency of feature differences during boosting or cross-validation as well. In addition, it could be used to optimize and speed up existing tree based algorithms.