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).
Training Data: https://d396qusza40orc.cloudfront.net/predmachlearn/pml-training.csv
Test Data: 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. If you use the document you create for this class for any purpose please cite them as they have been very generous in allowing their data to be used for this kind of assignment.
The goal of your project is to predict the manner in which they did the exercise. This is the “classe” variable in the training set. You may use any of the other variables to predict with. You should create a report describing:
You will also use your prediction model to predict 20 different test cases.
The required libraries are listed and loaded into the environment below. For reproduce-ability, we have also set a seed.
library(tidyverse)
library(data.table)
library(caret)
library(rattle)
library(parallel)
library(doParallel)
set.seed(19930212)Load data into R.
train_url <- "https://d396qusza40orc.cloudfront.net/predmachlearn/pml-training.csv"
test_url <- "https://d396qusza40orc.cloudfront.net/predmachlearn/pml-testing.csv"
raw_train <- fread(train_url)
raw_test <- fread(test_url)Some initial exploratory analysis was performed on the data. This can be found in notebooks repository in several different notebooks. The original data has dimensions of 19622, 160. This many features are likely not required and could result in over-fitting. The next section of this document will outline how the final features were selected for use in the model.
Before selecting any features, we will first clean the data so it is ready for use with the caret package.
# Create a new dataframe for training
training <- as.data.frame(raw_train)
# clean the data
# turn the classe column into a factor
training$classe <- factor(training$classe)
# turn all integers into numeric
for (i in colnames(training)){
class_i <- class(training[[i]])
if(class_i == "integer"){training[[i]] <- as.numeric(training[[i]])}
}
# how many features are there?
dim(training)## [1] 19622 160Features that have little variance will not be useful to include in the model. We can use R to search for features with near zero variance, and remove them from the model.
# check for near zero variance
nsv <- nearZeroVar(training, saveMetrics = TRUE)
nsv[nsv$nzv == TRUE,]# remove the columns with near zero variances
drop_cols <- rownames(nsv[nsv$nzv == TRUE,])
training <- training %>% select(-drop_cols)
# how many features are there?
dim(training)## [1] 19622 124The data has some columns with many missing values. Generally, the columns that look to be summary level statistics are missing values in most cases (e.g. columns starting with max, min, var, etc.)
# total missing values
sum(is.na(training))## [1] 1250007# missing values by column
sapply(training, function(x) sum(is.na(x))/nrow(training))## V1 user_name raw_timestamp_part_1
## 0.0000000 0.0000000 0.0000000
## raw_timestamp_part_2 cvtd_timestamp num_window
## 0.0000000 0.0000000 0.0000000
## roll_belt pitch_belt yaw_belt
## 0.0000000 0.0000000 0.0000000
## total_accel_belt kurtosis_roll_belt kurtosis_picth_belt
## 0.0000000 0.9798186 0.9809398
## skewness_roll_belt skewness_roll_belt.1 max_roll_belt
## 0.9797676 0.9809398 0.9793089
## max_picth_belt max_yaw_belt min_roll_belt
## 0.9793089 0.9798186 0.9793089
## min_pitch_belt min_yaw_belt amplitude_roll_belt
## 0.9793089 0.9798186 0.9793089
## amplitude_pitch_belt var_total_accel_belt avg_roll_belt
## 0.9793089 0.9793089 0.9793089
## stddev_roll_belt var_roll_belt avg_pitch_belt
## 0.9793089 0.9793089 0.9793089
## stddev_pitch_belt var_pitch_belt avg_yaw_belt
## 0.9793089 0.9793089 0.9793089
## stddev_yaw_belt var_yaw_belt gyros_belt_x
## 0.9793089 0.9793089 0.0000000
## gyros_belt_y gyros_belt_z accel_belt_x
## 0.0000000 0.0000000 0.0000000
## accel_belt_y accel_belt_z magnet_belt_x
## 0.0000000 0.0000000 0.0000000
## magnet_belt_y magnet_belt_z roll_arm
## 0.0000000 0.0000000 0.0000000
## pitch_arm yaw_arm total_accel_arm
## 0.0000000 0.0000000 0.0000000
## var_accel_arm gyros_arm_x gyros_arm_y
## 0.9793089 0.0000000 0.0000000
## gyros_arm_z accel_arm_x accel_arm_y
## 0.0000000 0.0000000 0.0000000
## accel_arm_z magnet_arm_x magnet_arm_y
## 0.0000000 0.0000000 0.0000000
## magnet_arm_z kurtosis_roll_arm kurtosis_picth_arm
## 0.0000000 0.9832841 0.9833860
## kurtosis_yaw_arm skewness_roll_arm skewness_pitch_arm
## 0.9798695 0.9832331 0.9833860
## skewness_yaw_arm max_picth_arm max_yaw_arm
## 0.9798695 0.9793089 0.9793089
## min_yaw_arm amplitude_yaw_arm roll_dumbbell
## 0.9793089 0.9793089 0.0000000
## pitch_dumbbell yaw_dumbbell kurtosis_roll_dumbbell
## 0.0000000 0.0000000 0.9795638
## kurtosis_picth_dumbbell skewness_roll_dumbbell skewness_pitch_dumbbell
## 0.9794109 0.9795128 0.9793599
## max_roll_dumbbell max_picth_dumbbell max_yaw_dumbbell
## 0.9793089 0.9793089 0.9795638
## min_roll_dumbbell min_pitch_dumbbell min_yaw_dumbbell
## 0.9793089 0.9793089 0.9795638
## amplitude_roll_dumbbell amplitude_pitch_dumbbell total_accel_dumbbell
## 0.9793089 0.9793089 0.0000000
## var_accel_dumbbell avg_roll_dumbbell stddev_roll_dumbbell
## 0.9793089 0.9793089 0.9793089
## var_roll_dumbbell avg_pitch_dumbbell stddev_pitch_dumbbell
## 0.9793089 0.9793089 0.9793089
## var_pitch_dumbbell avg_yaw_dumbbell stddev_yaw_dumbbell
## 0.9793089 0.9793089 0.9793089
## var_yaw_dumbbell gyros_dumbbell_x gyros_dumbbell_y
## 0.9793089 0.0000000 0.0000000
## gyros_dumbbell_z accel_dumbbell_x accel_dumbbell_y
## 0.0000000 0.0000000 0.0000000
## accel_dumbbell_z magnet_dumbbell_x magnet_dumbbell_y
## 0.0000000 0.0000000 0.0000000
## magnet_dumbbell_z roll_forearm pitch_forearm
## 0.0000000 0.0000000 0.0000000
## yaw_forearm kurtosis_roll_forearm kurtosis_picth_forearm
## 0.0000000 0.9835898 0.9836408
## skewness_roll_forearm skewness_pitch_forearm max_picth_forearm
## 0.9835389 0.9836408 0.9793089
## max_yaw_forearm min_pitch_forearm min_yaw_forearm
## 0.9835898 0.9793089 0.9835898
## amplitude_pitch_forearm total_accel_forearm var_accel_forearm
## 0.9793089 0.0000000 0.9793089
## gyros_forearm_x gyros_forearm_y gyros_forearm_z
## 0.0000000 0.0000000 0.0000000
## accel_forearm_x accel_forearm_y accel_forearm_z
## 0.0000000 0.0000000 0.0000000
## magnet_forearm_x magnet_forearm_y magnet_forearm_z
## 0.0000000 0.0000000 0.0000000
## classe
## 0.0000000# Remove columns that do not have enough data. We will remove any column that does not have data for at least 95% of observations.
greater_than_95 <- sapply(training, function(x) sum(is.na(x))/nrow(training)<0.05)
training <- training[, greater_than_95]
# how many features are there?
dim(training)## [1] 19622 59Lastly, based on our knowledge of the data and problem we are looking to solve, we will remove specific features that are not relevant:
# drop selected columns
drop_cols <- c("V1", "user_name", "raw_timestamp_part_1", "raw_timestamp_part_2", "cvtd_timestamp", "new_window", "num_window")
training <- training %>% select(-one_of(drop_cols))## Warning: Unknown columns: `new_window`# remove columns that start with "total"
training <- training %>% select(-starts_with("total"))
# how many features are there?
dim(training)## [1] 19622 49Our model will use the following features.
str(training)## 'data.frame': 19622 obs. of 49 variables:
## $ roll_belt : num 1.41 1.41 1.42 1.48 1.48 1.45 1.42 1.42 1.43 1.45 ...
## $ pitch_belt : num 8.07 8.07 8.07 8.05 8.07 8.06 8.09 8.13 8.16 8.17 ...
## $ yaw_belt : num -94.4 -94.4 -94.4 -94.4 -94.4 -94.4 -94.4 -94.4 -94.4 -94.4 ...
## $ gyros_belt_x : num 0 0.02 0 0.02 0.02 0.02 0.02 0.02 0.02 0.03 ...
## $ gyros_belt_y : num 0 0 0 0 0.02 0 0 0 0 0 ...
## $ gyros_belt_z : num -0.02 -0.02 -0.02 -0.03 -0.02 -0.02 -0.02 -0.02 -0.02 0 ...
## $ accel_belt_x : num -21 -22 -20 -22 -21 -21 -22 -22 -20 -21 ...
## $ accel_belt_y : num 4 4 5 3 2 4 3 4 2 4 ...
## $ accel_belt_z : num 22 22 23 21 24 21 21 21 24 22 ...
## $ magnet_belt_x : num -3 -7 -2 -6 -6 0 -4 -2 1 -3 ...
## $ magnet_belt_y : num 599 608 600 604 600 603 599 603 602 609 ...
## $ magnet_belt_z : num -313 -311 -305 -310 -302 -312 -311 -313 -312 -308 ...
## $ roll_arm : num -128 -128 -128 -128 -128 -128 -128 -128 -128 -128 ...
## $ pitch_arm : num 22.5 22.5 22.5 22.1 22.1 22 21.9 21.8 21.7 21.6 ...
## $ yaw_arm : num -161 -161 -161 -161 -161 -161 -161 -161 -161 -161 ...
## $ gyros_arm_x : num 0 0.02 0.02 0.02 0 0.02 0 0.02 0.02 0.02 ...
## $ gyros_arm_y : num 0 -0.02 -0.02 -0.03 -0.03 -0.03 -0.03 -0.02 -0.03 -0.03 ...
## $ gyros_arm_z : num -0.02 -0.02 -0.02 0.02 0 0 0 0 -0.02 -0.02 ...
## $ accel_arm_x : num -288 -290 -289 -289 -289 -289 -289 -289 -288 -288 ...
## $ accel_arm_y : num 109 110 110 111 111 111 111 111 109 110 ...
## $ accel_arm_z : num -123 -125 -126 -123 -123 -122 -125 -124 -122 -124 ...
## $ magnet_arm_x : num -368 -369 -368 -372 -374 -369 -373 -372 -369 -376 ...
## $ magnet_arm_y : num 337 337 344 344 337 342 336 338 341 334 ...
## $ magnet_arm_z : num 516 513 513 512 506 513 509 510 518 516 ...
## $ roll_dumbbell : num 13.1 13.1 12.9 13.4 13.4 ...
## $ pitch_dumbbell : num -70.5 -70.6 -70.3 -70.4 -70.4 ...
## $ yaw_dumbbell : num -84.9 -84.7 -85.1 -84.9 -84.9 ...
## $ gyros_dumbbell_x : num 0 0 0 0 0 0 0 0 0 0 ...
## $ gyros_dumbbell_y : num -0.02 -0.02 -0.02 -0.02 -0.02 -0.02 -0.02 -0.02 -0.02 -0.02 ...
## $ gyros_dumbbell_z : num 0 0 0 -0.02 0 0 0 0 0 0 ...
## $ accel_dumbbell_x : num -234 -233 -232 -232 -233 -234 -232 -234 -232 -235 ...
## $ accel_dumbbell_y : num 47 47 46 48 48 48 47 46 47 48 ...
## $ accel_dumbbell_z : num -271 -269 -270 -269 -270 -269 -270 -272 -269 -270 ...
## $ magnet_dumbbell_x: num -559 -555 -561 -552 -554 -558 -551 -555 -549 -558 ...
## $ magnet_dumbbell_y: num 293 296 298 303 292 294 295 300 292 291 ...
## $ magnet_dumbbell_z: num -65 -64 -63 -60 -68 -66 -70 -74 -65 -69 ...
## $ roll_forearm : num 28.4 28.3 28.3 28.1 28 27.9 27.9 27.8 27.7 27.7 ...
## $ pitch_forearm : num -63.9 -63.9 -63.9 -63.9 -63.9 -63.9 -63.9 -63.8 -63.8 -63.8 ...
## $ yaw_forearm : num -153 -153 -152 -152 -152 -152 -152 -152 -152 -152 ...
## $ gyros_forearm_x : num 0.03 0.02 0.03 0.02 0.02 0.02 0.02 0.02 0.03 0.02 ...
## $ gyros_forearm_y : num 0 0 -0.02 -0.02 0 -0.02 0 -0.02 0 0 ...
## $ gyros_forearm_z : num -0.02 -0.02 0 0 -0.02 -0.03 -0.02 0 -0.02 -0.02 ...
## $ accel_forearm_x : num 192 192 196 189 189 193 195 193 193 190 ...
## $ accel_forearm_y : num 203 203 204 206 206 203 205 205 204 205 ...
## $ accel_forearm_z : num -215 -216 -213 -214 -214 -215 -215 -213 -214 -215 ...
## $ magnet_forearm_x : num -17 -18 -18 -16 -17 -9 -18 -9 -16 -22 ...
## $ magnet_forearm_y : num 654 661 658 658 655 660 659 660 653 656 ...
## $ magnet_forearm_z : num 476 473 469 469 473 478 470 474 476 473 ...
## $ classe : Factor w/ 5 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...Partition the training data into a subset so we can better estimate the out of sample error rate.
inTrain <- createDataPartition(y = training$classe, p = 0.6, list = FALSE)
training_train <- training[inTrain,]
training_test <- training[-inTrain,]Random forests was selected as the model of choice due to the high accuracy of the predictions. Originally, a classification tree model was built using “rpart”. While this model was highly interpret-able, the accuracy was not very good. Using rpart I was not able to achieve the 80% accuracy required to pass the week 4 quiz. This lead to the random forest model described below as the final model.
Random forests was selected as model due to its ability to make very accurate predictions. On this large dataset, random forest can be quite slow. To speed up the model tips were used from this blog post on parrallel implementation.
# Step 1: configure parallel processing
cluster <- makeCluster(detectCores() - 1) # leave one core for the OS
registerDoParallel(cluster)
# Step 2: Configure train control object
fitControl <- trainControl(method = "cv", # use cross validation, faster than random forest defaults
number = 5, # using larger number could improve accuracy, but 5 should be suffecient
allowParallel = TRUE)
# Step 3: Develop training model
model_1 <- train(classe ~ .,
data = training_train,
method = "rf",
trControl = fitControl)
# Step 4: De-register the parallel processing cluster
stopCluster(cluster)
registerDoSEQ()# view model
model_1$finalModel##
## Call:
## randomForest(x = x, y = y, mtry = param$mtry)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 2
##
## OOB estimate of error rate: 0.85%
## Confusion matrix:
## A B C D E class.error
## A 3341 4 0 1 2 0.002090800
## B 21 2251 7 0 0 0.012286090
## C 0 19 2032 3 0 0.010710808
## D 0 0 35 1893 2 0.019170984
## E 0 0 1 5 2159 0.002771363First we will assess the accuracy of our model using the same data that was used to build the model.
predict_in_sample <- predict(model_1$finalModel, newdata = training_train, type = "class")
conf_in_sample <- confusionMatrix(predict_in_sample, training_train$classe)
conf_in_sample## Confusion Matrix and Statistics
##
## Reference
## Prediction A B C D E
## A 3348 0 0 0 0
## B 0 2279 0 0 0
## C 0 0 2054 0 0
## D 0 0 0 1930 0
## E 0 0 0 0 2165
##
## 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.0000conf_in_sample$overall[[1]]## [1] 1Our model has predicted the correct classe with an accuracy of 1. This number gives us a good sense of what the out of sample error may be. However we can get a better estimate by using the portion of the training data that we set aside to estimate the out of sample error rate.
To get a better estimate of the true out of sample error rate, we can predict using the training data that was set aside.
predict_out_sample <- predict(model_1$finalModel, newdata = training_test, type = "class")
conf_out_sample <- confusionMatrix(predict_out_sample, training_test$classe)
conf_out_sample## Confusion Matrix and Statistics
##
## Reference
## Prediction A B C D E
## A 2230 13 0 0 0
## B 1 1503 20 0 0
## C 0 2 1346 38 0
## D 0 0 2 1248 6
## E 1 0 0 0 1436
##
## Overall Statistics
##
## Accuracy : 0.9894
## 95% CI : (0.9869, 0.9916)
## No Information Rate : 0.2845
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.9866
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: A Class: B Class: C Class: D Class: E
## Sensitivity 0.9991 0.9901 0.9839 0.9705 0.9958
## Specificity 0.9977 0.9967 0.9938 0.9988 0.9998
## Pos Pred Value 0.9942 0.9862 0.9711 0.9936 0.9993
## Neg Pred Value 0.9996 0.9976 0.9966 0.9942 0.9991
## Prevalence 0.2845 0.1935 0.1744 0.1639 0.1838
## Detection Rate 0.2842 0.1916 0.1716 0.1591 0.1830
## Detection Prevalence 0.2859 0.1942 0.1767 0.1601 0.1832
## Balanced Accuracy 0.9984 0.9934 0.9889 0.9846 0.9978Our model has predicted the correct classe with an accuracy of 0.9894214. This number is a good estimate of the out of sample error rate, as these records were not used to build the model.
The final step is to apply our model to the test data. Before applying the model, we must first make the same transformations to the test data that we made to the training data.
# Create a new data frame
testing <- as.data.frame(raw_test)
# turn all integers into numeric
for (i in colnames(testing)){
class_i <- class(testing[[i]])
if(class_i == "integer"){testing[[i]] <- as.numeric(testing[[i]])}
}
# Keep only the same columns that were used in the training model
keep_cols <- names(training)
keep_cols <- keep_cols[!keep_cols == "classe"] # remove classe, since this does not exist in the testing data
testing <- testing[ ,keep_cols]Our data is now ready to be fed into the model for predictions.
predict_final <- predict(model_1$finalModel, newdata = testing, type = "class")
(predict_final)## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## B A B A A E D B A A B C B A E E A B B B
## Levels: A B C D EThese predictions were submitted into week 4 Coursera project quiz with 100% accuracy.