Overview

The goal of this course project is to explore the Weight Lifting Exercise Dataset which is available from the website here: http://groupware.les.inf.puc-rio.br/har, and build a model for predicting classe variable which corresponds to the quality of exercising a weight lifting activity. Such a model could be practically used for recognition of user’s activity and providing feedback about quality of exercising.

In this course project the model will be used for predicting classe variable in 20 different test cases for grading purpose. This imposes requirements on the model accuracy which must be hight enough for predicting correctly at least 80% of test cases (passing grade). If we want to be 80% sure about predicting correctly all test cases, the minimum required model accuracy would be:

0.81/20 = 0.9889048

Obtaining and cleaning data

The data that can be used for building a model is available here: https://d396qusza40orc.cloudfront.net/predmachlearn/pml-training.csv

We first load data from this file and remove all rows with new_window = “yes”, because these rows contain malformed values (e.g. “#DIV/0!”). There are only 2.07% of such rows and we hope that removing them won’t affect much the accuracy.

  # Read original data
  data <- read.csv("pml-training.csv", na.strings = c("NA", ""))
  
  # Filter out rows with new_window = "yes" (they contain garbage)
  data <- data[data$new_window == "no",]

We save filtered data to temp file and re-read it again, this time parsing all 152 predictor columns into numeric format:

  tmpFile <- tempfile("tmp_data_", tmpdir = getwd(), fileext = ".csv")
  write.csv(data, file = tmpFile, row.names = FALSE)
  
  # Re-read filtered data from temp file
  data <- read.csv(tmpFile,
                   colClasses = c(
                       "integer",   # X
                       "factor",    # user_name
                       "integer",   # raw_timestamp_part_1
                       "integer",   # raw_timestamp_part_2
                       "character", # cvtd_timestamp
                       "factor",    # new_window
                       "integer",   # num_window
                       rep("numeric", 152),  # all 152 predictors are numeric
                       "factor"     # classe
                   ),
                   na.strings = c("NA"))
  
  # Delete temp file
  unlink(tmpFile)

Then we remove some columns which contain technical information about how and when the data were recorded and are not relevant to the nature of problem. These columns are: X, raw_timestamp_part_1, raw_timestamp_part_2, cvtd_timestamp, new_window and num_window:

    data$X <- NULL
    data$raw_timestamp_part_1 <- NULL
    data$raw_timestamp_part_2 <- NULL
    data$cvtd_timestamp <- NULL
    data$new_window <- NULL
    data$num_window <- NULL

Finally we notice that many columns contain no data at all (all values are NA). We remove these columns too:

   for (col in names(data)) {
      # Actually all column with at least one NA have NAs in all rows
      if (sum(is.na(data[,col])) > 0) {
          data[,col] <- NULL
      }
   }

Remaining data frame has 54 columns (including the class column which is the outcome variable):

##  [1] "user_name"            "roll_belt"            "pitch_belt"          
##  [4] "yaw_belt"             "total_accel_belt"     "gyros_belt_x"        
##  [7] "gyros_belt_y"         "gyros_belt_z"         "accel_belt_x"        
## [10] "accel_belt_y"         "accel_belt_z"         "magnet_belt_x"       
## [13] "magnet_belt_y"        "magnet_belt_z"        "roll_arm"            
## [16] "pitch_arm"            "yaw_arm"              "total_accel_arm"     
## [19] "gyros_arm_x"          "gyros_arm_y"          "gyros_arm_z"         
## [22] "accel_arm_x"          "accel_arm_y"          "accel_arm_z"         
## [25] "magnet_arm_x"         "magnet_arm_y"         "magnet_arm_z"        
## [28] "roll_dumbbell"        "pitch_dumbbell"       "yaw_dumbbell"        
## [31] "total_accel_dumbbell" "gyros_dumbbell_x"     "gyros_dumbbell_y"    
## [34] "gyros_dumbbell_z"     "accel_dumbbell_x"     "accel_dumbbell_y"    
## [37] "accel_dumbbell_z"     "magnet_dumbbell_x"    "magnet_dumbbell_y"   
## [40] "magnet_dumbbell_z"    "roll_forearm"         "pitch_forearm"       
## [43] "yaw_forearm"          "total_accel_forearm"  "gyros_forearm_x"     
## [46] "gyros_forearm_y"      "gyros_forearm_z"      "accel_forearm_x"     
## [49] "accel_forearm_y"      "accel_forearm_z"      "magnet_forearm_x"    
## [52] "magnet_forearm_y"     "magnet_forearm_z"     "classe"

Data partitioning

Before we can begin with data analysis we split the cleaned data into training (75%) and testing (25%) sets for validation purposes. The training set will be used for building the model, and the testing set will be only used for validation and estimation of out-of-sample error of the final model.

set.seed(239237)

inTrain <- createDataPartition(data$class, p = 3/4, list = FALSE)

training <- data[inTrain,]
testing <- data[-inTrain,]

Selecting predictor variables

We choose random forest algorithm for the model for its accuracy. However, the number of predictor variables being large we should try to eliminate some of them and use only the most important predictors for the final model. For this, we first train a simple random forest model with only 20 trees using all 53 predictors:

small_forest_model <- randomForest(classe ~ ., data = training, ntree = 20, importance = TRUE)

and get the information about importance of predictors from that model:

varImpPlot(small_forest_model)

This plot shows two sets of predictor variables orderd by their importance according to two measures. We extract the top 12 predictors variables using the following code:

    # get importance information from the model
    importanceMatrix <- importance(small_forest_model)
    # convert it to data frame with VarName column containing the name of variable
    important_predictors <- as.data.frame(importanceMatrix)
    important_predictors$VarName <- rownames(important_predictors)
    # assign two scores to each variable according to two importance measures
    important_predictors$Score1 <- order(order(important_predictors$MeanDecreaseAccuracy))
    important_predictors$Score2 <- order(order(important_predictors$MeanDecreaseGini))
    # order variables by they total score
    ordered_predictors <- important_predictors$VarName[order(-important_predictors$Score1 - important_predictors$Score2)]
    # and finally select top 12 predictor variables
    selected_predictors <- ordered_predictors[seq(1, min(12, length(ordered_predictors)))]
    
    print(selected_predictors)
##  [1] "roll_belt"            "yaw_belt"             "magnet_dumbbell_y"   
##  [4] "roll_dumbbell"        "pitch_belt"           "pitch_forearm"       
##  [7] "magnet_dumbbell_z"    "roll_forearm"         "roll_arm"            
## [10] "total_accel_dumbbell" "magnet_forearm_z"     "accel_belt_z"

Let’s take look at correlation matrix of selected variables on this heat map plot:

corMatrix <- cor(training[,selected_predictors])
diag(corMatrix) <- 0
corrplot(corMatrix)

We see that two variables are highly correlated: yaw_belt and roll_belt. However, both of them are the most important predictor variables according to the small random forest model, so we keep them in the final model. Another highly correlated pair of variables are accel_belt_z and roll_belt. We will remove accel_belt_z from the final model as the least important of these two.

selected_predictors <- selected_predictors[!selected_predictors %in% c("accel_belt_z")]

Modeling

Finally we have 11 selected predictor variables and we can build the final large random forest model with 500 trees. We use train function from caret package which allows to specify training control parameter to use 3-fold cross validation. Choosing small number of folds allows to reduce computation time.

# train full random forest model (500 trees) using selected predictors
model <- train(x = training[,selected_predictors], 
               y = training$classe,
               method = "rf",
               trControl = trainControl(method = "cv", number = 3))

Validation

We validate the final model using the testing set (remaining 25% part of pml-training.csv file not used for training). The confusion matrix give an idea about the accuracy:

confusionMatrix(predict(model, testing), testing$class)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    A    B    C    D    E
##          A 1362    2    0    0    0
##          B    5  922    3    0    2
##          C    0    5  827    4    3
##          D    0    0    8  781    1
##          E    0    0    0    1  876
## 
## Overall Statistics
##                                           
##                Accuracy : 0.9929          
##                  95% CI : (0.9901, 0.9951)
##     No Information Rate : 0.2847          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.991           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: A Class: B Class: C Class: D Class: E
## Sensitivity            0.9963   0.9925   0.9869   0.9936   0.9932
## Specificity            0.9994   0.9974   0.9970   0.9978   0.9997
## Pos Pred Value         0.9985   0.9893   0.9857   0.9886   0.9989
## Neg Pred Value         0.9985   0.9982   0.9972   0.9988   0.9985
## Prevalence             0.2847   0.1935   0.1745   0.1637   0.1837
## Detection Rate         0.2836   0.1920   0.1722   0.1626   0.1824
## Detection Prevalence   0.2840   0.1941   0.1747   0.1645   0.1826
## Balanced Accuracy      0.9979   0.9949   0.9919   0.9957   0.9965
accuracy <- mean(predict(model, testing) == testing$class)

Accuracy of the final model is 0.9929196 and out of sample error is 0.0070804.

We can calculate the probability of getting 80% or 100% correct predictions out of 20 test cases in the final quiz using this model:

p = accuracy
p_80 <- p^20 + p^19*(1-p)*choose(20,19) + p^18*(1-p)^2*choose(20,18) + p^17*(1-p)^3*choose(20,17) + p^16*(1-p)^4*choose(20,16)
p_100 <- p^20

c(chance_80=p_80, chance_100=p_100)
##  chance_80 chance_100 
##  0.9999997  0.8675246