Background

A large amount of data about personal activity is nowadays available, thanks to devices such as FitBit, Jawbone Up and Nike FuelBand. People are interested in tracking how much they do of an activity, but they also can be interested in knowing how correctly they do the activity. In this project, we use data from accelerometers on the belt, forearm, arm, and dumbell of 6 participants. They were asked to perform barbell lifts correctly and incorrectly in 5 different ways. The goal is to use the provided data to build a model that predicts the quality of barbell lifting.

Executive Summary

In this report, we propose a machine learning model to predict the quality of dumbbell lifting during an exercise. The provided data represent measurements collected from accelerometers of the participants while they were performing Biceps Curl. We first did some pre-processing for the data by removing the mostly-empty columns and keeping the relevant features. We then sliced the data into training and testing sets. Using cross-validation on the training set, several models were implemented and compared. In particular, the ensemble methods (random forest and gradient boosted trees) had the best performance in terms of accuracy. After fine-tuning the more promising models, the final model was chosen to be gradient boosted trees with an out-of-sample accuracy of 99.71%.

Exploration and Pre-Processing of the Dataset

The provided data, which is available here https://d396qusza40orc.cloudfront.net/predmachlearn/pml-training.csv, represent measurements collected from accelerometers on the belt, forearm, arm, and dumbbell of 6 participants and belong to one of the following 5 classes: exactly according to the specification (Class A), throwing the elbows to the front (Class B), lifting the dumbbell only halfway (Class C), lowering the dumbbell only halfway (Class D) and throwing the hips to the front (Class E). More information about the data is available from this source: http://web.archive.org/web/20161224072740/http:/groupware.les.inf.puc-rio.br/har The data used consist of 19622 instances where each instance has 160 features including the category class. However, not all features are relevant. The first few columns consist of participants’ name and timestamps, and some other columns are mostly empty. This is why we first extract the relevant features:

data<-read.csv("pml-training.csv")
pos=grep("magnet|^accel|^(total_accel)|gyros|^pitch|^roll|^yaw",names(data));
pos=append(pos,ncol(data))
new_data=data[,pos]

Each observation now consists of 53 features (the first 52 measurements consist of features on the Euler angles (roll, pitch and yaw), as well as the raw accelerometer, gyroscope and magnetometer readings and the last column is for class type (A,B,C,D or E)). We split the data into training and testing sets (75% for training and 25% for testing). The training set will be used for model training and selection and the testing set will be used to predict the out-of-sample error for the final chosen model.

library(caret)
set.seed(101)
inTrain <- createDataPartition(new_data$classe, p=0.75, list=FALSE)
training<-new_data[inTrain,]
testing<- new_data[-inTrain,]

We now explore the statistics of each feature using the skim_to_wide function, which computes the statistics of each row:

library(skimr)
statTrain <- skim_to_wide(training)
#1:type, 2:variable, 3:missing, 4:complete, 5:n
#9:mean, 10:sd, 11:p0, 13:p50, 15:p100, 16:hist
statTrain[, c(1:5, 9:11, 13, 15:16)]
## # A tibble: 53 x 11
##    type  variable missing complete n     mean  sd    p0    p50   p100 
##    <chr> <chr>    <chr>   <chr>    <chr> <chr> <chr> <chr> <chr> <chr>
##  1 fact~ classe   0       14718    14718 <NA>  <NA>  <NA>  <NA>  <NA> 
##  2 inte~ accel_a~ 0       14718    14718 " -6~ 182.~ -404  -44   435  
##  3 inte~ accel_a~ 0       14718    14718 "  3~ 109.~ -318  14    308  
##  4 inte~ accel_a~ 0       14718    14718 " -7~ 134.~ -636  -48   292  
##  5 inte~ accel_b~ 0       14718    14718 "  -~ " 29~ -83   -15   85   
##  6 inte~ accel_b~ 0       14718    14718 "  3~ " 28~ -65   33    149  
##  7 inte~ accel_b~ 0       14718    14718 " -7~ 100.~ -275  -151  104  
##  8 inte~ accel_d~ 0       14718    14718 " -2~ " 66~ -237  -8    235  
##  9 inte~ accel_d~ 0       14718    14718 "  5~ " 80~ -189  41    315  
## 10 inte~ accel_d~ 0       14718    14718 " -3~ 109.~ -334  -1    318  
## # ... with 43 more rows, and 1 more variable: hist <chr>

No missing values were found and all features are numeric and have normal values (note that we did not show all rows because of space limitation). We then scale and center the features as follows:

preProcess_model <- preProcess(training, method=c('center','scale'))
training<-predict(preProcess_model,newdata=training)
testing<-predict(preProcess_model,newdata=testing)

Models Implemented

The data is now ready for model training. For that, we start with implementing some models with their default tuning parameters. We will then fine-tune the models that look more promising. In order to compare the models and tune their parameters, we choose 5-fold cross-validation as the resampling method:

fitControl <- trainControl(method = "cv",number = 5,savePredictions = "final")

We also use Accuracy as the performance metric.

Linear Models

We first start with the following linear models: multinomial logistic regression and linear SVM (SVM: support vector machine, which similarly to logistic regression, classifies the observations by finding a linear boundary between classes. For SVM, the linear boundary is found by maximizing the margins between classes):

##Logistic regression
glmFit<-train(classe~.,data=training,method="multinom",trControl=fitControl,verbose=FALSE,trace=FALSE)
##SVM 
svmFit<-train(classe~.,data=training,method="svmLinear",trControl=fitControl,verbose=FALSE)

The resulting accuracy of each model is:

max(glmFit$results$Accuracy)
## [1] 0.7333876
max(svmFit$results$Accuracy)
## [1] 0.783665

Non-Linear Models

We now implement non-linear models to see if the obtained performance with linear models can be improved. We use the models of naive-Bayes, tree and SVM with radial kernel:

##tree
treeFit<-train(classe~.,data=training,method="rpart2",trControl=fitControl,tuneGrid=expand.grid(maxdepth=c(10,20,30)))
##Naive Bayes 
library(klaR)
nbFit<-train(classe~.,data=training, method="nb",trControl=fitControl,verbose=FALSE)
## Radial SVM
svmRadFit<-train(classe~.,data=training,method="svmRadial",trControl=fitControl,verbose=FALSE)

The resulting accuracies are as follows:

max(treeFit$results$Accuracy)
## [1] 0.7303279
max(nbFit$results$Accuracy)
## [1] 0.752613
max(svmRadFit$results$Accuracy)
## [1] 0.9230869

We see that with SVM with radial kernel, there is a great enhancement in performance.

Ensemble Methods

We now try ensemble methods such as random forests and gradient boosted trees:

rfFit<-train(classe~.,data=training, method="rf",trControl=fitControl)
boostFit<-train(classe~.,data=training,method="xgbTree",trControl=fitControl)

and the obtained accuracies are as follows.

max(rfFit$results$Accuracy)
## [1] 0.9917112
max(boostFit$results$Accuracy)
## [1] 0.9944282

The ensemble methods and the SVM with radial kernel seem more promising. We next are going to fine-tune the paramaters of each of these promising models.

Parameters Tuning

We now tune the parameter mtry (number of variables randomly sampled as candidates at each split) for the model of random forest.

set.seed(102)
rfGrid<- expand.grid(mtry=c(seq(6,10,2),15))
rfFit2<-train(classe~.,data=training, method="rf",trControl=fitControl,tuneGrid=rfGrid)

We obtain the following results:

rfFit2
## Random Forest 
## 
## 14718 samples
##    52 predictor
##     5 classes: 'A', 'B', 'C', 'D', 'E' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 11774, 11775, 11774, 11775, 11774 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##    6    0.9937491  0.9920922
##    8    0.9944965  0.9930379
##   10    0.9937491  0.9920923
##   15    0.9934775  0.9917488
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 8.

We now tune the parameters for the model of boosted trees:

set.seed(103)
boostGrid<-expand.grid(max_depth=6,nrounds=150,eta=c(0.3,0.4),gamma=0,colsample_bytree=c(0.6,0.8,1),subsample=c(0.5,0.75,1),min_child_weight=1)
boostFit2<-train(classe~.,data=training,method="xgbTree",trControl=fitControl,tuneGrid=boostGrid)

We obtain the following results:

boostFit2
## eXtreme Gradient Boosting 
## 
## 14718 samples
##    52 predictor
##     5 classes: 'A', 'B', 'C', 'D', 'E' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 11775, 11776, 11774, 11774, 11773 
## Resampling results across tuning parameters:
## 
##   eta  colsample_bytree  subsample  Accuracy   Kappa    
##   0.3  0.6               0.50       0.9966708  0.9957887
##   0.3  0.6               0.75       0.9966028  0.9957027
##   0.3  0.6               1.00       0.9963990  0.9954451
##   0.3  0.8               0.50       0.9959233  0.9948432
##   0.3  0.8               0.75       0.9963990  0.9954450
##   0.3  0.8               1.00       0.9966708  0.9957888
##   0.3  1.0               0.50       0.9963311  0.9953590
##   0.3  1.0               0.75       0.9959912  0.9949293
##   0.3  1.0               1.00       0.9961272  0.9951014
##   0.4  0.6               0.50       0.9962631  0.9952729
##   0.4  0.6               0.75       0.9963990  0.9954450
##   0.4  0.6               1.00       0.9963310  0.9953592
##   0.4  0.8               0.50       0.9955838  0.9944138
##   0.4  0.8               0.75       0.9961953  0.9951874
##   0.4  0.8               1.00       0.9965347  0.9956167
##   0.4  1.0               0.50       0.9961951  0.9951872
##   0.4  1.0               0.75       0.9961951  0.9951873
##   0.4  1.0               1.00       0.9961271  0.9951012
## 
## Tuning parameter 'nrounds' was held constant at a value of 150
## 
## Tuning parameter 'gamma' was held constant at a value of 0
## 
## Tuning parameter 'min_child_weight' was held constant at a value of 1
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were nrounds = 150, max_depth = 6,
##  eta = 0.3, gamma = 0, colsample_bytree = 0.6, min_child_weight = 1
##  and subsample = 0.5.

We now tune the parameters for the model of radial SVM:

grid_radial <- expand.grid(sigma = c(0.1, 0.2),C = c(1.5,2,5,10))
svmRadFit2<-train(classe~.,data=training,method="svmRadial",trControl=fitControl,tuneGrid=grid_radial,verbose=FALSE)

We obtain the following results:

svmRadFit2
## Support Vector Machines with Radial Basis Function Kernel 
## 
## 14718 samples
##    52 predictor
##     5 classes: 'A', 'B', 'C', 'D', 'E' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 11775, 11774, 11775, 11775, 11773 
## Resampling results across tuning parameters:
## 
##   sigma  C     Accuracy   Kappa    
##   0.1     1.5  0.9836933  0.9793693
##   0.1     2.0  0.9859355  0.9822069
##   0.1     5.0  0.9900122  0.9873665
##   0.1    10.0  0.9919825  0.9898590
##   0.2     1.5  0.9827421  0.9781694
##   0.2     2.0  0.9836254  0.9792878
##   0.2     5.0  0.9854599  0.9816094
##   0.2    10.0  0.9862751  0.9826409
## 
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.1 and C = 10.

Comparison of Models

To compare between the models, we plot the dotplot of the results obtained on the 5 folds.

# Compare model performances using resample()
models_compare <- resamples(list(RF=rfFit2, SVM=svmRadFit2, boost=boostFit2))
# Summary of the models performances
summary(models_compare)
## 
## Call:
## summary.resamples(object = models_compare)
## 
## Models: RF, SVM, boost 
## Number of resamples: 5 
## 
## Accuracy 
##            Min.   1st Qu.    Median      Mean   3rd Qu.      Max. NA's
## RF    0.9925246 0.9942255 0.9942255 0.9944965 0.9949049 0.9966021    0
## SVM   0.9908257 0.9915053 0.9918451 0.9919825 0.9918506 0.9938859    0
## boost 0.9952462 0.9955812 0.9972826 0.9966708 0.9972826 0.9979613    0
## 
## Kappa 
##            Min.   1st Qu.    Median      Mean  3rd Qu.      Max. NA's
## RF    0.9905433 0.9926949 0.9926953 0.9930379 0.993554 0.9957018    0
## SVM   0.9883968 0.9892554 0.9896838 0.9898590 0.989692 0.9922670    0
## boost 0.9939866 0.9944101 0.9965626 0.9957887 0.996563 0.9974214    0
# Draw dor plots to compare models
scales <- list(x=list(relation="free"), y=list(relation="free"))
dotplot(models_compare, scales=scales)

We see from the average performance on the 5 folds and the performance on each fold of each method, that the gradient boosted tree has the best performance. Therefore, we choose the model of gradient boosted trees.

Testing the final model

We now test the model of gradient boosted trees on the testing set in order to predict the out-of-sample accuracy.

pred=predict(boostFit2,testing)
confusionMatrix(data=pred,reference=testing$classe)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    A    B    C    D    E
##          A 1395    2    0    0    0
##          B    0  947    2    0    0
##          C    0    0  851    5    0
##          D    0    0    2  799    3
##          E    0    0    0    0  898
## 
## Overall Statistics
##                                           
##                Accuracy : 0.9971          
##                  95% CI : (0.9952, 0.9984)
##     No Information Rate : 0.2845          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.9964          
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: A Class: B Class: C Class: D Class: E
## Sensitivity            1.0000   0.9979   0.9953   0.9938   0.9967
## Specificity            0.9994   0.9995   0.9988   0.9988   1.0000
## Pos Pred Value         0.9986   0.9979   0.9942   0.9938   1.0000
## Neg Pred Value         1.0000   0.9995   0.9990   0.9988   0.9993
## Prevalence             0.2845   0.1935   0.1743   0.1639   0.1837
## Detection Rate         0.2845   0.1931   0.1735   0.1629   0.1831
## Detection Prevalence   0.2849   0.1935   0.1746   0.1639   0.1831
## Balanced Accuracy      0.9997   0.9987   0.9970   0.9963   0.9983

Therefore the expected out-of-sample accuracy is 99.71%.

Importance of Featurres

We now plot the importnace of each feacture, to see which feature contribute most to predicting the quality of dumbbell lifting.

plot(varImp(boostFit2))

We see that the two variables: roll_belt and yaw_belt are the top two predictors.

Conclusion

We have implemented various models on the provided data to recognize the quality of dumbbell lifting. We have seen that the ensemble method based on gradient boosted trees performed best in terms of accuracy and has an out-of-sample accuracy of 99.71%.