Practical Learning Project
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. In this project, your goal will be to 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. More information is available from the website here: http://groupware.les.inf.puc-rio.br/har (see the section on the Weight Lifting Exercise Dataset).
We are given 2 sets, a training set of 19622 obs of 160 variables, and a testing set with only 20 obs of 160 variables.
most rows have a lot of NAs for a lot of variables. what happens, is that at some intervals, the devices were making a “full checkup” and collecting more info (maybe average of values during an interval ?) theres 406 such complete observations in the set. Outside those 406 obs, a lot of variables are NA. For this study, i removed the variables with NA, and/or invalid values such as “”. some variables had incorrectly their format as “factor” and i converted it to numeric. I’m ending up with “only”" 84 variables.
str(aa[,1:7])## 'data.frame': 19622 obs. of 7 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ user_name : Factor w/ 6 levels "adelmo","carlitos",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ raw_timestamp_part_1: int 1323084231 1323084231 1323084231 1323084232 1323084232 1323084232 1323084232 1323084232 1323084232 1323084232 ...
## $ raw_timestamp_part_2: int 788290 808298 820366 120339 196328 304277 368296 440390 484323 484434 ...
## $ cvtd_timestamp : Factor w/ 20 levels "02/12/2011 13:32",..: 9 9 9 9 9 9 9 9 9 9 ...
## $ new_window : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ num_window : int 11 11 11 12 12 12 12 12 12 12 ...among the first 7 variables are user_name: who did the exercise timestamps new_window: a boolean variable. “yes” if we are among the 406 obs described above num_window
the last is
str(aa[,84])## Factor w/ 5 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...class: a factor variable taking values A,B,C,D,E which is the value we want to predict.
Now i will show you how easy it was to make a very simple model with 100% accuracy, and i insist its 100%, not 95%, not 99% not 99.5%, not 99.9% !
Before using the brute force of the CPU, let’s do some exploratory data analysis to find the variables with the most explanating power. It turns out, you need only ONE variable, and that variable is called… num_window !
answers = rep("A", 20) #initialize vector
testing_hack <- testing[7]
index <- c(74, 431, 439, 194, 235, 504, 485, 440, 323, 664, 859, 461, 257, 408, 779, 302, 48, 361, 72, 255)
hack_table <- aa[,c(7,84)]
hack_table1 <- table(hack_table)
head(hack_table1)## classe
## num_window A B C D E
## 1 0 0 0 0 20
## 2 0 0 0 0 21
## 3 0 0 0 0 3
## 4 0 0 0 0 28
## 5 0 0 0 0 21
## 6 0 0 0 0 21yes to num_window correspond only one classe ! no even overlapping problem ! the testing data does provide a num_window and simply making a lookup for the num_window in the table will provide the correct classe of course !!! because the authors (or the teaching staff ?) didn’t take the elementary precaution of using testing data OUTSIDE the training data.
hack_table2 <- cbind(rownames(hack_table1),colnames(hack_table1[,apply(hack_table1,1,which.max)]))
answers <- hack_table2[match(index,hack_table2[,1]),2]
# print(answers) try it on your PCI’ve got 20/20 predictions in a few minutes, before even running a single model…
Now, call me a cheater… well nobody said that you should submit the 20 values output by a ML model ! and definitely i can call my procedure a model ! I’m ROFL when i see people in the forum complaining they can’t find the right values whereas they have a so-called model with 99.5% accuracy.
I feel I could end the assignment here, job done, but I’ll play ball, and try to think what can we do if we DO NOT use this num_window. Lets continue some exploratory data analysis.
head(aa[,c(1:5),7])## X user_name raw_timestamp_part_1 raw_timestamp_part_2 cvtd_timestamp
## 1 1 carlitos 1323084231 788290 05/12/2011 11:23
## 2 2 carlitos 1323084231 808298 05/12/2011 11:23
## 3 3 carlitos 1323084231 820366 05/12/2011 11:23
## 4 4 carlitos 1323084232 120339 05/12/2011 11:23
## 5 5 carlitos 1323084232 196328 05/12/2011 11:23
## 6 6 carlitos 1323084232 304277 05/12/2011 11:23Look at the timestamp variable together with the username and classe. If we browse through the data, we can see that the people who participated in the experience did the whole thing in a short time frame. its not like they went one day and did exercise A, and went back 1 day later to re-do it again. No, they did a set of exercise A for a short period of time, then exercise B, … that makes the whole procedure questionable, really… even if any model was successful, it wouldn’t prove it would work for even the SAME persons on another day (testing and validation test taken from the same data…) !!!
Of course, using the timestamp or any of the header variables would be cheating as much as using just num_window. As Tomek Pyda said on the forum https://class.coursera.org/predmachlearn-014/forum/thread?thread_id=129: “Lets look at yaw_belt: [he plots yaw_belt vs num_window]
It gives you much information about the time. Looks like some curl lifts had different starting values than others. Other variables carry a lot of information about num_window too (see yaw_belt,magnet_dumbbell_y etc)
I think what we are so effectively predicting is: “to which particular exercise the data belongs to”. And the classe as a result. Not how well somebody did it. That is why we can predict classes C and D.“[end of Tomek’s quote]
one can also suspect some of the instruments were not reset when participants changed and thats altering the data. Unfortunately, by lack of time and information, I will ignore that. I will just have a didactic approach for the rest.
Ok, for this project, as I’ve started a bit late, I will just focus on the RCART model, a Random Forest may be more appropriate. I was afraid of performance problems for RF, and using CARET for the RCART was already quite slow. Actually, reading the forum afterwards, it looks like randomForest using the package randomForest is quite fast, but I’m running out of time and will leave it for another day. Anyway, the rest of the discussion is more about methodology and RCART is enough for discussion about Tree modelling.
Let’s create as always some training sets and validation sets we will use cross validation with K-folds (k=10) and 10 repeats for our CART models.
library(caret)
inTrain <- createDataPartition(y=aa$classe,p=0.7,list=FALSE)
set1 <- aa[inTrain,]
training1 <- set1[,8:84]
set1v <- aa[-inTrain,]
valid1 <- set1v[,8:84]
FitControl <- trainControl(method="repeatedcv",
number=10,
repeats=10
)the first reference model I will use is on the training1 set, and tuneLength=10 to run it in a reasonable time.
library(caret)
modFit1 <- train(classe~.,method="rpart",data=training1,trControl=FitControl,tuneLength=10)
print(confusionMatrix(predict(modFit1,newdata=valid1),valid1$classe))## Confusion Matrix and Statistics
##
## Reference
## Prediction A B C D E
## A 1367 211 16 101 34
## B 35 560 67 40 66
## C 57 140 777 152 132
## D 190 185 139 636 128
## E 25 43 27 35 722
##
## Overall Statistics
##
## Accuracy : 0.6902
## 95% CI : (0.6782, 0.702)
## No Information Rate : 0.2845
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.6087
## Mcnemar's Test P-Value : < 2.2e-16
##
## Statistics by Class:
##
## Class: A Class: B Class: C Class: D Class: E
## Sensitivity 0.8166 0.49166 0.7573 0.6598 0.6673
## Specificity 0.9140 0.95617 0.9010 0.8695 0.9729
## Pos Pred Value 0.7906 0.72917 0.6176 0.4977 0.8474
## Neg Pred Value 0.9261 0.88685 0.9462 0.9288 0.9285
## Prevalence 0.2845 0.19354 0.1743 0.1638 0.1839
## Detection Rate 0.2323 0.09516 0.1320 0.1081 0.1227
## Detection Prevalence 0.2938 0.13050 0.2138 0.2172 0.1448
## Balanced Accuracy 0.8653 0.72392 0.8292 0.7646 0.8201see ? without much effort and without any pre-processing, without even trying randomForest, just being dumb, we’ve already got an honourable 70% accuracy. We can see also some patterns: we have already a very high specificity for all classes A to E (95%), but the sensitivity to B, C, and D is quite lower than for A (sitting) and E (walking). As a result, trying to predict the testing set, and comparing to the answers, we get only 8/20. far from the “predicted”" 14/20 (70%). It seems natural that sitting being the most stable and E the most instable they get detected more easily than B (sitting down) C(standing) and D(standing up)
the model tree uses 12 variables.
I wanted to check if those 12 variables are very random. I mean, if we take off those 12 variables and re-run a CART model, theres still more than 60 variables to play with, and can we build a model with similar accuracy ? It turns out i’ve tried it, and for the same command line, the accuracy of such a model drops to 0.5643 (not enough space to publish all the results, try it on your PC)
# library(caret)
# testing without 12 variables found in the model 1
# training2 <- training1[,-c(1,59,56,58,57,3,53,2,19,71,39,48)]
# modFit2 <- train(classe~.,method="rpart",data=training2,trControl=FitControl,tuneLength=10)
# accuracy drops to 0.5643 with modFit2so here are those 12 variables (the number is the index they have in the training1 object) roll_belt 1
pitch_forearm 59
magnet_dumbbell_y 56 <–[we’ll talk a bit about that one later]
roll_forearm 58
magnet_dumbbell_z 57
yaw_belt 3 <– [remember what Tomek said about it]
accel_dumbbell_y 53
pitch_belt 2
magnet_belt_z 19
accel_forearm_x 71
roll_dumbbell 39
total_accel_dumbbell 48
Let’s put these variables aside for later and try to explain why they seem important.
Let’s keep on and move on to tuneLength=30
library(caret)
modFit1pro <- train(classe~.,method="rpart",data=training1,trControl=FitControl,tuneLength=30)
print(confusionMatrix(predict(modFit1pro,newdata=valid1),valid1$classe))## Confusion Matrix and Statistics
##
## Reference
## Prediction A B C D E
## A 1495 130 30 52 24
## B 61 796 49 56 60
## C 34 122 878 57 58
## D 56 59 46 764 65
## E 28 32 23 35 875
##
## Overall Statistics
##
## Accuracy : 0.817
## 95% CI : (0.8069, 0.8268)
## No Information Rate : 0.2845
## P-Value [Acc > NIR] : < 2e-16
##
## Kappa : 0.7684
## Mcnemar's Test P-Value : 3.9e-15
##
## Statistics by Class:
##
## Class: A Class: B Class: C Class: D Class: E
## Sensitivity 0.8931 0.6989 0.8558 0.7925 0.8087
## Specificity 0.9440 0.9524 0.9442 0.9541 0.9754
## Pos Pred Value 0.8637 0.7789 0.7641 0.7717 0.8812
## Neg Pred Value 0.9569 0.9295 0.9688 0.9591 0.9577
## Prevalence 0.2845 0.1935 0.1743 0.1638 0.1839
## Detection Rate 0.2540 0.1353 0.1492 0.1298 0.1487
## Detection Prevalence 0.2941 0.1737 0.1952 0.1682 0.1687
## Balanced Accuracy 0.9185 0.8256 0.9000 0.8733 0.8921final_res <- predict(modFit1pro$finalModel,newdata=testing)
print(colnames(final_res[,apply(final_res,1,which.max)]))## [1] "B" "A" "C" "C" "A" "E" "D" "D" "A" "A" "C" "C" "B" "A" "E" "E" "A"
## [18] "D" "D" "B"accuracy is improved to over 80% (0.817). now only class B suffers from a somewhat lower sensitivity. we got 16/20 out of the answers of the test, more in line with the (predicted) accuracy. B was wrongly predicted twice (once to A and once to C), E was predicted wrongly once to C, and B was wrongly predicted to A once. again B causes problem.
Let’s look what the final Tree looks like
# print(modFit1pro$finalModel) try it on your PC. not enough space hereif you get a headache when trying to read this, its normal !
So, we can’t do a writeup without at least a nice plot so here it is,
We can see that the RED leaves corresponding to E on the right, and DARK GREEN leaves corresponding to A on the left are easily found. the Roll Belt detects very well if someone is walking. while the pitch arm and yaw belt are very useful to detect someone sitting. for B,C,D its more difficult to interpret, there are more variables than in the system with accuracy 70. but we see again that the variables magnet_dumbbell (x,y,z), total accel_dumbbell, accel_dumbell_y and the rest of the earlier 12 variables used in the less complex system appear quite a lot in the tree
So more tuning and more CPU time, surely would get us an even better model ?
wait a minute. lets try something first. I did the following experiment. create a training set without info on pedro, and test it on the subset exclusively composed of pedro’s data. (i tried with tunelength = 10 because of lack of cpu time.)
## Warning: package 'caret' was built under R version 3.1.3## Loading required package: lattice
## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 3.1.3## Loading required package: rparthere’s an extract of the result.
Confusion Matrix and Statistics
Reference <br>Prediction A B C D E
A 0 0 0 0 0
B 192 159 157 134 38
C 0 0 0 0 0
D 9 0 0 1 0
E 0 0 0 0 130
Overall Statistics
Accuracy : 0.3537 <br>
The accuracy dropped by half from 0.704 to 0.354 !!! the model kept predicting always B !! Such a result puts again into perspective these algos. Surely in the study we took 6 people, its because it was assumed it may be enough to predict do HAR on ANYBODY ? Looking at our procedure, we can highly suspect that if we took a testing set from a 7th person, the prediction of the model would look ugly. Remember i was already suspicious that even if the SAME person took the test again, the results would look much more ugly. Really theres 19622 obs, but theres only 6 persons doing 5 exercices, so in some sense, only 30 experiments.
Also, we should be suspicious because to have a model which gets “only” to 80% accuracy like I did its already impossible to make much sense of the tree result and values.
I tried fitting the same kind of CART trees to the subset of 406 observations. the accuracy was significantly below the models i showed here for the same complexity
Now I’ll show something quite interesting that I found through data exploration.
## [,1] [,2]
## [1,] 207 396
## [2,] 283 435
## [3,] 269 408
## [4,] -66 437
## [5,] 281 487
## [6,] -586 -477
## [7,] 211 589
## [8,] -3600 632
## [9,] 330 633
## [10,] 147 425
## [11,] 234 617
## [12,] -730 -462This is an extract of the range (min and max) of this series for class A (lines 1 to 6) and B (lines 7 to 12), each with 6 observations for each participant in the experiment. notably one participant, PEDRO (lines 6 and 12)! yes ! the one who made a failure of the test above, now you can see why perhaps ? that important variable that we see on the tree, has NEGATIVE values for PEDRO. has he inverted the dumbbells ? did he do his exercise upside down ? theres some big outliers in this series for other participants as well. maybe thats also why my little experiment above with pedro went SOOOO bad ! ;-)
In a last effort before the deadline, i wanted to try randomForest as i took courage from the fact in the forum that it could be run in reasonable time after all.
Here are the astonishing results:
## Warning: package 'randomForest' was built under R version 3.1.3## randomForest 4.6-10
## Type rfNews() to see new features/changes/bug fixes.## Confusion Matrix and Statistics
##
## Reference
## Prediction A B C D E
## A 1672 10 0 0 0
## B 2 1127 10 0 0
## C 0 2 1015 6 0
## D 0 0 1 957 0
## E 0 0 0 1 1082
##
## Overall Statistics
##
## Accuracy : 0.9946
## 95% CI : (0.9923, 0.9963)
## No Information Rate : 0.2845
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.9931
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: A Class: B Class: C Class: D Class: E
## Sensitivity 0.9988 0.9895 0.9893 0.9927 1.0000
## Specificity 0.9976 0.9975 0.9984 0.9998 0.9998
## Pos Pred Value 0.9941 0.9895 0.9922 0.9990 0.9991
## Neg Pred Value 0.9995 0.9975 0.9977 0.9986 1.0000
## Prevalence 0.2845 0.1935 0.1743 0.1638 0.1839
## Detection Rate 0.2841 0.1915 0.1725 0.1626 0.1839
## Detection Prevalence 0.2858 0.1935 0.1738 0.1628 0.1840
## Balanced Accuracy 0.9982 0.9935 0.9938 0.9963 0.999999.5% accuracy with default settings !
it gets better, 20/20 on the testing set !
So what are the important variables ?
## roll_belt yaw_belt pitch_forearm magnet_dumbbell_z
## 872.5302 629.9239 552.7861 539.3314
## pitch_belt magnet_dumbbell_y roll_forearm magnet_dumbbell_x
## 480.0703 475.9672 429.1212 332.2248
## roll_dumbbell magnet_belt_y accel_dumbbell_y magnet_belt_z
## 294.1145 281.8044 269.5255 267.6885very interestingly, 10 out of the 12 most important variables of the RF model are the same than in the simpler CART model we documented above, who had only 69% accuracy.
So, is computer science so great these days that a total noob can just press a button and get a 99.5% accurate model ? (or rather, a computer can learn everything by itself very fast ? )
the answer to both questions, i fear, is YES.
However, my hunch feeling is that we are overfitting data. I highly doubt the model works as it is with 99% accuracy ITRL (in the real life). The reason is , i think the protocol of the experience is bad, so that the data is not independent (always a key assumption in statistical modelling). Actually it is almost continuous as we have seen. There are also artefacts (pedro) in the data, which means the more complex models had a turnaround to get along with the fact that Pedro was upside down (?), electrods inverted or whatever…Its great to have an intelligent system able to recognize errors made by users, but this time it may also be an error by the searchers? have they noticed such things? was it voluntary from their side ?
As the data is continuous, theres in reality much less observations than 19,000. So in the validation set, theres already some of the data used in the training set, if not ALL of it ! and I’ve demonstrated that its the same for the testing set.
to test that: (sorry no time for that): take only 1 every 3 or 4 row for example and see the loss of accuracy…
I’d be more impressed if the RF worked with a sample of 10,000 people or 1,000 people doing 10 times the exercices at different times with data taken only every second…and I’m sure the methods used here would work fine.
Also, a truly intelligent system would not look just at a photo of instant values in the instruments, but would look at how they evolve through time. thats not the approach taken by the searchers there with a testing set being instant values.