Background

Using devices such as Jawbone Up, Nike FuelBand, and Fitbit, it is now possible to collect a large amount of data about personal activity relatively inexpensively. These type of devices are part of the quantified self movement – a group of enthusiasts who take measurements about themselves regularly to improve their health, to find patterns in their behavior, or because they are tech geeks. One thing that people regularly do is quantify how much of a particular activity they do, but they rarely quantify how well they do it.

The goal of this project was to predict to the quality of exercise using the Weight Lifting Exercise Dataset from http://groupware.les.inf.puc-rio.br/har. This dataset consists of measurements taken from six participants who performed barbell lifts correctly (class ‘A’) and incorrectly (classes ‘B’, ‘C’, ‘D’, and ‘E’). Measurements included accelerometer, gyroscope, and magnetometer readings from the participants’ armband, glove, belt, and barbell. We combined these measurements to predict the type of lift the participants performed (A-F).

Feature selection

In this section, we show how to reduce the number of features from 159 down to 34.

We start by loading the training set given to us.

training <- read.csv("pml-training.csv")
dim(training)
## [1] 19622   160
colnames(training)[160]
## [1] "classe"

We see that the training set has 19622 measurements and 159 features (the last column is the class label). Let’s check if any of these features are empty.

ef <- lapply(lapply(training,is.na),sum)
ef <- 100*as.numeric(ef)/nrow(training)
table(ef>0)
## 
## FALSE  TRUE 
##    93    67

We see that 67 features are almost entirely empty, i.e., full of NAs. We remove these from further consideration.

empty_features <- ef>0
reduced_set <- training[,!empty_features]

Next, let’s find features that are nearly zero. These features have very low variance, and are unlikely to be useful.

library(caret)
## Loading required package: lattice
nzv_features <- nearZeroVar(reduced_set)

These near-zero features contain only a small number of unique values (< 450) relative to the size of feature (19622). Therefore, we can remove them from further consideration.

Finally, we look at features that are highly correlated with each other.

reduced_set <- reduced_set[,-nzv_features] # Remove near zero features
reduced_set_M <- data.matrix(reduced_set[,7:58]) # Remove factor variables user_name, timestamps, new_window, and classe
corrMatrix <- cor(reduced_set_M)
highCorr <- findCorrelation(corrMatrix, cutoff=0.75)

As an example, we plot three of these features (accel_belt_x,accel_belt_y, and accel_belt_z) against each other. It is easy to see that they are highly correlated. Therefore, it makes sense to remove such features from further consideration.

# Remove corelated features, and retain some factor variables we had previously ignored
keep_features <- c("user_name","num_window",colnames(reduced_set_M)[-highCorr])
feature_names <- colnames(training)
selectedFeatures <- pmatch(keep_features,feature_names)
length(keep_features)
## [1] 34
keep_features
##  [1] "user_name"            "num_window"           "yaw_belt"            
##  [4] "gyros_belt_x"         "gyros_belt_y"         "gyros_belt_z"        
##  [7] "magnet_belt_x"        "magnet_belt_z"        "roll_arm"            
## [10] "pitch_arm"            "yaw_arm"              "total_accel_arm"     
## [13] "gyros_arm_x"          "gyros_arm_z"          "magnet_arm_x"        
## [16] "magnet_arm_z"         "roll_dumbbell"        "pitch_dumbbell"      
## [19] "yaw_dumbbell"         "total_accel_dumbbell" "gyros_dumbbell_y"    
## [22] "gyros_dumbbell_z"     "magnet_dumbbell_z"    "roll_forearm"        
## [25] "pitch_forearm"        "yaw_forearm"          "total_accel_forearm" 
## [28] "gyros_forearm_x"      "gyros_forearm_y"      "accel_forearm_x"     
## [31] "accel_forearm_z"      "magnet_forearm_x"     "magnet_forearm_y"    
## [34] "magnet_forearm_z"

Split training data

Now that we have chosen the features, let us partition the training data using 2-fold cross validation. We consider two equal sets (trainSet and validSet). We will train our model on one set and evaluate on the other to estimate the out-of-sample error.

NOTE: To speed up the execution, we will only consider half the training set here for purposes of demonstration. On a faster machine, we would consider the full training set.

set.seed(1234)
# To speed things up, only consider half the dataset
inConsideration <- createDataPartition(y=training$classe,p=0.5,list=FALSE)
trainSmall <- training[inConsideration,]
# Split further into training and validation sets
inTrain <- createDataPartition(y=trainSmall$classe,p=0.5,list=FALSE)
trainSet <- trainSmall[inTrain,]
valSet <- trainSmall[-inTrain,]

Apply random forest classifier

1st fold

Next, we fit a random forest classifier to the selected features in trainSet.

trainfeatures <- trainSet[,selectedFeatures] 
trainresponse <- trainSet$classe
modfit <- train(trainfeatures,trainresponse,method="rf",prox=TRUE,ntree=10)
modfit
## Random Forest 
## 
## 4907 samples
##   34 predictor
##    5 classes: 'A', 'B', 'C', 'D', 'E' 
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## 
## Summary of sample sizes: 4907, 4907, 4907, 4907, 4907, 4907, ... 
## 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa      Accuracy SD  Kappa SD   
##    2    0.9328023  0.9149584  0.007500722  0.009466817
##   18    0.9696993  0.9616631  0.006214723  0.007853551
##   34    0.9680542  0.9595954  0.007537980  0.009494613
## 
## Accuracy was used to select the optimal model using  the largest value.
## The final value used for the model was mtry = 18.

Having trained the model on trainSet, let us evaluate the in-sample errors. The overall in-sample accuracy is 99.98%.

pred <- predict(modfit,trainfeatures)
## Loading required package: randomForest
## randomForest 4.6-10
## Type rfNews() to see new features/changes/bug fixes.
confusionMatrix(pred,trainresponse)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    A    B    C    D    E
##          A 1395    0    0    0    0
##          B    0  950    0    0    0
##          C    0    0  855    0    0
##          D    0    0    1  804    0
##          E    0    0    0    0  902
## 
## Overall Statistics
##                                      
##                Accuracy : 0.9998     
##                  95% CI : (0.9989, 1)
##     No Information Rate : 0.2843     
##     P-Value [Acc > NIR] : < 2.2e-16  
##                                      
##                   Kappa : 0.9997     
##  Mcnemar's Test P-Value : NA         
## 
## Statistics by Class:
## 
##                      Class: A Class: B Class: C Class: D Class: E
## Sensitivity            1.0000   1.0000   0.9988   1.0000   1.0000
## Specificity            1.0000   1.0000   1.0000   0.9998   1.0000
## Pos Pred Value         1.0000   1.0000   1.0000   0.9988   1.0000
## Neg Pred Value         1.0000   1.0000   0.9998   1.0000   1.0000
## Prevalence             0.2843   0.1936   0.1744   0.1638   0.1838
## Detection Rate         0.2843   0.1936   0.1742   0.1638   0.1838
## Detection Prevalence   0.2843   0.1936   0.1742   0.1641   0.1838
## Balanced Accuracy      1.0000   1.0000   0.9994   0.9999   1.0000

But what about the out-of-sample errors? We run the model trained on trainSet on the held out validationSet to estimate the out-of-sample error. We see that the overall out-of-sample accuracy is 98.12%!

valfeatures <- valSet[,selectedFeatures]  
valresponse <- valSet$classe
pred <- predict(modfit,valfeatures)
confusionMatrix(pred,valresponse)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    A    B    C    D    E
##          A 1391   20    0    4    2
##          B    2  909   11    1    1
##          C    1    8  835   10    0
##          D    1    3    7  787    8
##          E    0    9    2    2  891
## 
## Overall Statistics
##                                          
##                Accuracy : 0.9812         
##                  95% CI : (0.977, 0.9849)
##     No Information Rate : 0.2844         
##     P-Value [Acc > NIR] : < 2.2e-16      
##                                          
##                   Kappa : 0.9763         
##  Mcnemar's Test P-Value : 0.0002218      
## 
## Statistics by Class:
## 
##                      Class: A Class: B Class: C Class: D Class: E
## Sensitivity            0.9971   0.9579   0.9766   0.9789   0.9878
## Specificity            0.9926   0.9962   0.9953   0.9954   0.9968
## Pos Pred Value         0.9817   0.9838   0.9778   0.9764   0.9856
## Neg Pred Value         0.9989   0.9900   0.9951   0.9959   0.9973
## Prevalence             0.2844   0.1935   0.1743   0.1639   0.1839
## Detection Rate         0.2836   0.1853   0.1702   0.1604   0.1817
## Detection Prevalence   0.2889   0.1884   0.1741   0.1643   0.1843
## Balanced Accuracy      0.9949   0.9770   0.9860   0.9871   0.9923

Given the impressive out-of-sample accuracy, it appears that our random forest classifier has learnt something useful. Let’s see which features the model found most useful.

plot(varImp(modfit))

2nd fold

Let us now swap the trainSet and valSet roles, i.e., train on valSet and test on trainSet.

modfit2 <- train(valfeatures,valresponse,method="rf",prox=TRUE,ntree=10)
modfit2
## Random Forest 
## 
## 4905 samples
##   34 predictor
##    5 classes: 'A', 'B', 'C', 'D', 'E' 
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## 
## Summary of sample sizes: 4905, 4905, 4905, 4905, 4905, 4905, ... 
## 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa      Accuracy SD  Kappa SD   
##    2    0.9273233  0.9079862  0.007065571  0.008967323
##   18    0.9694598  0.9613443  0.005055558  0.006374225
##   34    0.9688497  0.9605668  0.004304966  0.005465781
## 
## Accuracy was used to select the optimal model using  the largest value.
## The final value used for the model was mtry = 18.

Let us evaluate in-sample and out-of-sample errors for this model. We see that the in-sample accuracy is 99.98% and out-of-sample accuracy is 98.43%. These numbers are similar to what we observed in the first fold.

# In-sample error
pred <- predict(modfit2,valfeatures) 
confusionMatrix(pred,valresponse)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    A    B    C    D    E
##          A 1395    0    0    0    0
##          B    0  949    1    0    0
##          C    0    0  854    0    0
##          D    0    0    0  804    0
##          E    0    0    0    0  902
## 
## Overall Statistics
##                                      
##                Accuracy : 0.9998     
##                  95% CI : (0.9989, 1)
##     No Information Rate : 0.2844     
##     P-Value [Acc > NIR] : < 2.2e-16  
##                                      
##                   Kappa : 0.9997     
##  Mcnemar's Test P-Value : NA         
## 
## Statistics by Class:
## 
##                      Class: A Class: B Class: C Class: D Class: E
## Sensitivity            1.0000   1.0000   0.9988   1.0000   1.0000
## Specificity            1.0000   0.9997   1.0000   1.0000   1.0000
## Pos Pred Value         1.0000   0.9989   1.0000   1.0000   1.0000
## Neg Pred Value         1.0000   1.0000   0.9998   1.0000   1.0000
## Prevalence             0.2844   0.1935   0.1743   0.1639   0.1839
## Detection Rate         0.2844   0.1935   0.1741   0.1639   0.1839
## Detection Prevalence   0.2844   0.1937   0.1741   0.1639   0.1839
## Balanced Accuracy      1.0000   0.9999   0.9994   1.0000   1.0000
# Out-of-sample error
pred <- predict(modfit2,trainfeatures) 
confusionMatrix(pred,trainresponse)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    A    B    C    D    E
##          A 1389   10    0    3    0
##          B    4  931   10    5    2
##          C    0    7  844   17    0
##          D    2    2    1  775    9
##          E    0    0    1    4  891
## 
## Overall Statistics
##                                           
##                Accuracy : 0.9843          
##                  95% CI : (0.9804, 0.9876)
##     No Information Rate : 0.2843          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.9801          
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: A Class: B Class: C Class: D Class: E
## Sensitivity            0.9957   0.9800   0.9860   0.9639   0.9878
## Specificity            0.9963   0.9947   0.9941   0.9966   0.9988
## Pos Pred Value         0.9907   0.9779   0.9724   0.9823   0.9944
## Neg Pred Value         0.9983   0.9952   0.9970   0.9930   0.9973
## Prevalence             0.2843   0.1936   0.1744   0.1638   0.1838
## Detection Rate         0.2831   0.1897   0.1720   0.1579   0.1816
## Detection Prevalence   0.2857   0.1940   0.1769   0.1608   0.1826
## Balanced Accuracy      0.9960   0.9873   0.9900   0.9803   0.9933

Let’s see which features this model found most useful. We see that this set of features agrees with that of the other models (previous plot).

plot(varImp(modfit2))

Conclusions

We built a random forest classifier to predict the quality of exercise in the Weight Lifting Exercise Dataset. We evaluated the classifier using two-fold cross validation, and estimated the out-of-sample error rate to be 100 - (98.43 + 98.12)/2 = 1.725%. Note that this value is much less than chance error rate 100*5/6 = 83.33%. Therefore, we expect the classifier performance to generalize well to new samples from the Weight Lifting Exercise Dataset.