Dataset from Velloso, E.; Bulling, A.; Gellersen, H.; Ugulino, W.; Fuks, H. Qualitative Activity Recognition of Weight Lifting Exercises. Proceedings of 4th International Conference in Cooperation with SIGCHI (Augmented Human ’13) . Stuttgart, Germany: ACM SIGCHI, 2013. Read more: http://groupware.les.inf.puc-rio.br/har#ixzz3uJG8yOEC
We start by loading caret and other potentially useful packages
library(caret)## Warning: package 'caret' was built under R version 3.2.2## Loading required package: lattice
## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 3.2.3library(corrplot)And then provide remote urls to get the data from and set the working directory
train_url = "https://d396qusza40orc.cloudfront.net/predmachlearn/pml-training.csv"
test_url = "https://d396qusza40orc.cloudfront.net/predmachlearn/pml-testing.csv"
directory = "."
setwd(directory)
trainfile = "train.csv"
testfile = "test.csv"
train_dest = paste(directory, trainfile,sep ="/")
test_dest = paste(directory, testfile,sep ="/")
# if files do not exist then download them from url
if (!file.exists(train_dest) | !file.exists(test_dest)){
method = "curl"
download.file(train_url,destfile = train_dest,method = method)
download.file(test_url,destfile = test_dest,method = method)
}Read in “raw” datasets
train_raw <- read.csv(train_dest)
test_raw <-read.csv(test_dest)Preprocess datasets
Begin by removing columns with missing values
cols = apply(train_raw,2,function(x) {sum(is.na(x))==0})
train_raw = train_raw[,cols]And do the same on the test set
test_raw = test_raw[,cols]We now proceed on to removing columns containing “#DIV/0!” as it is not a meaningful value and hence will not be useful for prediction
cols = apply(train_raw,2,function(x) {length(grep("#D",x))==0} )
train_raw = train_raw[,cols]
test_raw = test_raw[,cols]At this point, after basic cleaning, the two datasets contain 60 columns. We can further reduce this number by noting that columns containing user names, time stamps and windows will not be useful for prediction and can thus be eliminated
train_raw = train_raw[,-grep("user_name|window|timestamp|^X",names(train_raw))]
test_raw = test_raw[,-grep("user_name|window|timestamp|^X",names(test_raw))]Which brings the number of columns down to 53.
We can now perform a slightly more sofisticated pre-processing by checking if any pairs of variables is significantly correlated in which case we can further reduce the dimensionality of the problem
We thus construct the correlation matrix and set its diagonal to 0
corrmat = abs(cor(train_raw[,-53]))
diag(corrmat) = 0if any pairs of variables is highly correlated then the dimensionality of the set can be further reduced
length(which(corrmat>0.8))## [1] 38Given the presence on non-negligible correlations we perform principal component analysis and retain only those variables contributing to 95% of the variance are retained.
preobj = preProcess(train_raw[,-53],method = "pca", thresh = 0.95)
trainPC = predict(preobj,train_raw[,-53])
testPC = predict(preobj,test_raw[,-53])
train = trainPC
test = testPC
ncol(train)## [1] 25ncol(test)## [1] 25At the end of this pre-processing phase we have obtained a polished training and test set. We now leave the test set aside and proceed to data slicing, setting a seed for reproducibility.
set.seed(20140604)We then split the training dataset in two parts for validation and correspondingly create two “classe” vectors against which validate the model
cvind = createDataPartition(train_raw$classe,p = 0.7,list = FALSE)
train_cv = train[cvind,]
classe_train_cv = train_raw$classe[cvind]
test_cv = train[-cvind,]
classe_test_cv = train_raw$classe[-cvind]Once the data has been sliced we enter the model training part. We begin by choosing the parameters for cross validation
controls1 = trainControl(method = "cv", number = 5)I chose the number of folds for cross validation to be 5 in order for it not be too small (hence limit the bias) but still manageable from a computational point of view. Normally I would have done few trials trying different values for number. Also, having the time, it could be interesting to set number = 10 and compare cross validation against repeated cross validation with 2 repeats and 5 folds both in terms of computational efficiency and performance of the fitted model.
We now proceed to model training. I choose to use a decision-tree based algorithm as they appear to me best suited for multi-variate classification problems. I begin by training a decision tree with maximum tree-depth = 30 (which should be the maximum allowed value)
if (!exists("modFitDecTree")){
modFitDecTree <- train(classe_train_cv ~.,
data = train_cv,
method="rpart2",
trControl = controls1,
tuneGrid = data.frame(.maxdepth = 30))
}## Loading required package: rpartprint(modFitDecTree$finalModel)## n= 13737
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 13737 9831 A (0.28 0.19 0.17 0.16 0.18)
## 2) PC14>=-0.9175744 12375 8509 A (0.31 0.19 0.19 0.17 0.14)
## 4) PC8< 1.347329 10248 6643 A (0.35 0.18 0.22 0.11 0.14)
## 8) PC3< 0.3309807 6444 3520 A (0.45 0.14 0.23 0.086 0.088) *
## 9) PC3>=0.3309807 3804 2914 B (0.18 0.23 0.2 0.16 0.22)
## 18) PC17< 0.6126656 2834 2025 B (0.23 0.29 0.19 0.14 0.15)
## 36) PC10>=-0.8004982 2153 1575 A (0.27 0.22 0.22 0.18 0.12)
## 72) PC14< 0.3835702 1168 721 A (0.38 0.25 0.08 0.13 0.16) *
## 73) PC14>=0.3835702 985 613 C (0.13 0.19 0.38 0.23 0.07) *
## 37) PC10< -0.8004982 681 346 B (0.11 0.49 0.12 0.037 0.24) *
## 19) PC17>=0.6126656 970 545 E (0.029 0.084 0.23 0.22 0.44) *
## 5) PC8>=1.347329 2127 1223 D (0.12 0.23 0.049 0.43 0.17)
## 10) PC25>=0.3589947 444 253 B (0.074 0.43 0.047 0.09 0.36) *
## 11) PC25< 0.3589947 1683 819 D (0.14 0.18 0.05 0.51 0.12) *
## 3) PC14< -0.9175744 1362 608 E (0.029 0.26 0.031 0.13 0.55) *and then a random forest: in order to improve computation speed I decided to set mtry, i.e. the variables tried at each split, to be 18. The reason is that after PCA all variables, being linearly independent, are likely to contribute on a split. On the other hand I wanted to prevent overfitting and gain some computational speed, so trying approx 2/3 of the available variables seemed to me a good compromise.
if (!exists("modFit1")){
modFit1 <- train(classe_train_cv ~.,
data = train_cv,
method="rf",
trControl = controls1,
prox = TRUE,
allowParallel = TRUE, ntrees = 250,
tuneGrid = data.frame(.mtry = 18))
}## Loading required package: randomForest
## randomForest 4.6-12
## Type rfNews() to see new features/changes/bug fixes.
##
## Attaching package: 'randomForest'
##
## The following object is masked from 'package:ggplot2':
##
## marginprint(modFit1$finalModel)##
## Call:
## randomForest(x = x, y = y, mtry = param$mtry, proximity = TRUE, allowParallel = TRUE, ntrees = 250)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 18
##
## OOB estimate of error rate: 3.49%
## Confusion matrix:
## A B C D E class.error
## A 3833 28 23 17 5 0.01868920
## B 54 2537 47 4 16 0.04552295
## C 11 42 2302 33 8 0.03923205
## D 9 9 99 2126 9 0.05595027
## E 3 20 20 22 2460 0.02574257We now check how the fitted models do against the training set
decision tree
predict0dt = predict(modFitDecTree,train_cv)
confusionMatrix(classe_train_cv,predict0dt)## Confusion Matrix and Statistics
##
## Reference
## Prediction A B C D E
## A 3371 108 131 228 68
## B 1208 526 186 306 432
## C 1575 102 372 84 263
## D 706 65 227 864 390
## E 752 324 69 201 1179
##
## Overall Statistics
##
## Accuracy : 0.4595
## 95% CI : (0.4511, 0.4679)
## No Information Rate : 0.5541
## P-Value [Acc > NIR] : 1
##
## Kappa : 0.2914
## Mcnemar's Test P-Value : <2e-16
##
## Statistics by Class:
##
## Class: A Class: B Class: C Class: D Class: E
## Sensitivity 0.4429 0.46756 0.37766 0.5134 0.50557
## Specificity 0.9127 0.83095 0.84128 0.8849 0.88198
## Pos Pred Value 0.8630 0.19789 0.15526 0.3837 0.46693
## Neg Pred Value 0.5686 0.94593 0.94595 0.9287 0.89716
## Prevalence 0.5541 0.08190 0.07170 0.1225 0.16976
## Detection Rate 0.2454 0.03829 0.02708 0.0629 0.08583
## Detection Prevalence 0.2843 0.19349 0.17442 0.1639 0.18381
## Balanced Accuracy 0.6778 0.64926 0.60947 0.6991 0.69378random forest
predict0rf = predict(modFit1,train_cv)
confusionMatrix(classe_train_cv,predict0rf)## Confusion Matrix and Statistics
##
## Reference
## Prediction A B C D E
## A 3906 0 0 0 0
## B 0 2658 0 0 0
## C 0 0 2396 0 0
## D 0 0 0 2252 0
## E 0 0 0 0 2525
##
## Overall Statistics
##
## Accuracy : 1
## 95% CI : (0.9997, 1)
## No Information Rate : 0.2843
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 1
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: A Class: B Class: C Class: D Class: E
## Sensitivity 1.0000 1.0000 1.0000 1.0000 1.0000
## Specificity 1.0000 1.0000 1.0000 1.0000 1.0000
## Pos Pred Value 1.0000 1.0000 1.0000 1.0000 1.0000
## Neg Pred Value 1.0000 1.0000 1.0000 1.0000 1.0000
## Prevalence 0.2843 0.1935 0.1744 0.1639 0.1838
## Detection Rate 0.2843 0.1935 0.1744 0.1639 0.1838
## Detection Prevalence 0.2843 0.1935 0.1744 0.1639 0.1838
## Balanced Accuracy 1.0000 1.0000 1.0000 1.0000 1.0000From examination of the results, we choose random forest as prediction algorithm due to its improved accuracy. However, the fact that the training model is fitted with accuracy 1 can still be interpreted as a warning sign of overfitting despite cross validation?
We now test the performance of the rf model on the part of the training test left out from training
predict1 = predict(modFit1,test_cv)
confusionMatrix(classe_test_cv,predict1)## Confusion Matrix and Statistics
##
## Reference
## Prediction A B C D E
## A 1655 5 6 7 1
## B 14 1114 6 0 5
## C 4 22 980 15 5
## D 4 1 42 909 8
## E 2 6 8 3 1063
##
## Overall Statistics
##
## Accuracy : 0.9721
## 95% CI : (0.9676, 0.9762)
## No Information Rate : 0.2853
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.9647
## Mcnemar's Test P-Value : 0.0004319
##
## Statistics by Class:
##
## Class: A Class: B Class: C Class: D Class: E
## Sensitivity 0.9857 0.9704 0.9405 0.9732 0.9824
## Specificity 0.9955 0.9947 0.9905 0.9889 0.9960
## Pos Pred Value 0.9886 0.9781 0.9552 0.9429 0.9824
## Neg Pred Value 0.9943 0.9928 0.9872 0.9949 0.9960
## Prevalence 0.2853 0.1951 0.1771 0.1587 0.1839
## Detection Rate 0.2812 0.1893 0.1665 0.1545 0.1806
## Detection Prevalence 0.2845 0.1935 0.1743 0.1638 0.1839
## Balanced Accuracy 0.9906 0.9826 0.9655 0.9811 0.9892Having obtained an estimated accuracy of 97% and out-of-sample error of 3%, and taking into account time/computational limitations, I believe the algorithm performs sufficiently good to proceed to prediction. I therefore apply the fitted model to the test set and write the answers in separate files
predict_on_test = predict(modFit1,test)
pml_write_files = function(x){
n = length(x)
for(i in 1:n){
filename = paste0("problem_id_",i,".txt")
write.table(x[i],file=filename,quote=FALSE,row.names=FALSE,col.names=FALSE)
}
}
answers = as.character(predict_on_test)
answers_dest = paste(directory, "prediction",sep ="/")
if (!file.exists(answers_dest)){
dir.create(answers_dest)
}
setwd(answers_dest)
pml_write_files(answers)