Getting and cleaning data

Data is coming from Human Activity Recognition project and can be downloaded from the web

download.file('https://d396qusza40orc.cloudfront.net/predmachlearn/pml-training.csv', 'pml-training.csv')
download.file('https://d396qusza40orc.cloudfront.net/predmachlearn/pml-testing.csv', 'pml-testing.csv')

plm.training <- read.csv('pml-training.csv', na.strings = c('', 'NA', '#DIV/0!'))
plm.test_cases <- read.csv('pml-testing.csv', na.strings = c('', 'NA', '#DIV/0!'))

As the training dataset contains a total of 19622, which is big enough, we can split it into training and validation datasets in order to have an estimation of the errors.

set.seed(52534)
inTrain = createDataPartition(plm.training$classe, p = 3/4)[[1]]
training = plm.training [ inTrain,]
testing = plm.training [-inTrain,]

Now we can work with the training sample and see an extract of the training data below:

kable(training[training$X %in% 238:243,1:20], row.names = F)
X user_name raw_timestamp_part_1 raw_timestamp_part_2 cvtd_timestamp new_window num_window roll_belt pitch_belt yaw_belt total_accel_belt kurtosis_roll_belt kurtosis_picth_belt kurtosis_yaw_belt skewness_roll_belt skewness_roll_belt.1 skewness_yaw_belt max_roll_belt max_picth_belt max_yaw_belt
238 pedro 1323094971 796282 05/12/2011 14:22 no 48 129 27.8 1.95 21 NA NA NA NA NA NA NA NA NA
240 pedro 1323094971 948362 05/12/2011 14:22 no 48 128 27.9 1.61 20 NA NA NA NA NA NA NA NA NA
241 pedro 1323094971 980388 05/12/2011 14:22 yes 48 128 27.8 1.52 21 -0.668237 0.988872 NA -0.638799 -1.717824 NA 2.7 21 -0.7
242 pedro 1323094972 32296 05/12/2011 14:22 no 49 128 27.8 1.49 21 NA NA NA NA NA NA NA NA NA
243 pedro 1323094972 48322 05/12/2011 14:22 no 49 128 27.8 1.47 20 NA NA NA NA NA NA NA NA NA

First thing we see is that some columns are empty for most of the observations:

  • Complete cases: 0 out of 14718
  • Columns with low observations: 100 out of 160

The columns with low data seem to be related to an “window aggregation”, as they are only filled in when new_window is yes. We can discard them as we will only work with the raw data, not the aggregations.

Additionally, the user_name is also included in the dataset. As the model is intended to be executer over any user (and not only the training ones), we must discard this column as well.

Finally, the columns related to timestamp, order (X) and window should be discarded as well, since they should not be used in the model.

With this considerations, a new training dataset removing the neccesary columns can be created.

training2 <- training %>%
                select(-one_of("X", "user_name")) %>%
                select(-contains("timestamp")) %>%
                select(-contains("window")) %>%
                select(-starts_with("kurtosis_")) %>%
                select(-starts_with("skewness_")) %>%
                select(-starts_with("max_")) %>%
                select(-starts_with("min_")) %>%
                select(-starts_with("amplitude_")) %>%
                select(-starts_with("var_")) %>%
                select(-starts_with("avg_")) %>%
                select(-starts_with("stddev_"))

It contains now a total of 52 predictors and 14718 observations.

Exploratory data analysis

Once the data set is cleaned, it’s the turn to make some exploratory data analysis.

Since there are many variables and it’s not clear their relationship, one option is to run a tree partition and see what are the firts variables in the tree.

set.seed(12432)
# Apply the model
model_tree <- train(classe ~ ., training2, method='rpart')
# Calculate the most important variables
imp <- varImp(model_tree)
most_imp <- imp$importance %>% mutate(Variable = row.names(imp$importance)) %>% arrange(desc(Overall)) %>% select(Variable, Importance = Overall) %>% slice(1:3)
# Show in a table
kable(most_imp)
Variable Importance
pitch_forearm 100.00000
roll_forearm 72.44185
roll_belt 70.75412

Now we can plot how the classe depends on those variables

ggplot(data=training2, aes(x=classe)) +
  geom_violin(aes(y = pitch_forearm, fill = 'pitch_forearm'), alpha = 0.3) +
  geom_violin(aes(y = roll_forearm, fill = 'roll_forearm'), alpha = 0.3) +
  geom_violin(aes(y = roll_belt, fill = 'roll_belt'), alpha = 0.3) +
  ylab('Value')

And we can see how effectively the figures are different and, depending on the classe, They have diferent intervals where they are bigger than the others.

Model building

For the model, it can be used a combined prediction model using stacking. The models used for building the combined models will be random forest, linear discriminator analysis and generalized boosted regression modeling

set.seed(2341)
# Calculate the three models and their predictions
model_rf <- train(classe ~ ., training2, method = 'rf')
model_lda <- train(classe ~ ., training2, method = 'lda')
model_gbm <- train(classe ~ ., training2, method = 'gbm')
pt_rf <- predict(model_rf)
pt_lda <- predict(model_lda)
pt_gbm <- predict(model_gbm)

# Combine them
pt_comb <- data.frame(classe = training2$classe, pt_rf, pt_lda, pt_gbm)
model_comb <- train(classe ~ ., pt_comb, method = 'rf')

Finally, we can apply this prediction model over the testing dataset to evaluate the accuracy for every method and the combination (using confusion matrix)

# Predictions
pt_rf <- predict(model_rf, newdata = testing)
pt_lda <- predict(model_lda, newdata = testing)
pt_gbm <- predict(model_lda, newdata = testing)

pt_comb <- data.frame(pt_rf, pt_lda, pt_gbm)
prediction_testing <- predict(model_comb, newdata = pt_comb)

# Confusion matrixes
conf_rf <- confusionMatrix(pt_rf, testing$classe)
conf_lda <- confusionMatrix(pt_lda, testing$classe)
conf_gbm <- confusionMatrix(pt_gbm, testing$classe)
conf_comb <- confusionMatrix(prediction_testing, testing$classe)

# Plots
plot_rf <- ggplot(as.data.frame(conf_rf$table)) + 
      geom_tile(aes(x=Reference, y=Prediction, fill=Freq)) +
      geom_text(aes(x=Reference, y=Prediction, label=Freq)) +
      scale_fill_gradient(low = 'grey', high='green') +
      ggtitle('Random Forest')

plot_lda <- ggplot(as.data.frame(conf_lda$table)) + 
      geom_tile(aes(x=Reference, y=Prediction, fill=Freq)) +
      geom_text(aes(x=Reference, y=Prediction, label=Freq)) +
      scale_fill_gradient(low = 'grey', high='green') +
      ggtitle('Linear Discriminant Analysis')

plot_gbm <- ggplot(as.data.frame(conf_gbm$table)) + 
      geom_tile(aes(x=Reference, y=Prediction, fill=Freq)) +
      geom_text(aes(x=Reference, y=Prediction, label=Freq)) +
      scale_fill_gradient(low = 'grey', high='green') +
      ggtitle('Generalized Boosted Regression')

plot_comb <- ggplot(as.data.frame(conf_comb$table)) + 
      geom_tile(aes(x=Reference, y=Prediction, fill=Freq)) +
      geom_text(aes(x=Reference, y=Prediction, label=Freq)) +
      scale_fill_gradient(low = 'grey', high='green') +
      ggtitle('Combined Model')

grid.arrange(plot_rf, plot_lda, plot_gbm, plot_comb)

# Finally add accuracy
accuracy_df <- data.frame(Method = c('Random Forest', 'Linear Discriminant Analysis', 'Generalized Boosted Regression', 'Combined Model'),
                          Accuracy = 100*c(conf_rf[['overall']]['Accuracy'], conf_lda[['overall']]['Accuracy'], conf_gbm[['overall']]['Accuracy'], conf_comb[['overall']]['Accuracy']))

kable(accuracy_df, col.names=c('Method', 'Accuracy (%)'))
Method Accuracy (%)
Random Forest 99.34747
Linear Discriminant Analysis 70.96248
Generalized Boosted Regression 70.96248
Combined Model 99.34747

We can see how the accuracy for random forest and the combined model is more than 99%, while lda and gbm is about 70%. Therefore, the combined model will be used for this exercise.

Finally, we can also check the confidence interval of the accuracy coming from the confusion matrix (using the combined model in this case). The 95% confidence interval for the accuracy is [0.9908006, 0.9955326], so around 99% accuracy is expected for the combined model.

kable(as.data.frame(conf_comb[['overall']]), col.names=c('Value'))
Value
Accuracy 0.9934747
Kappa 0.9917462
AccuracyLower 0.9908006
AccuracyUpper 0.9955326
AccuracyNull 0.2844617
AccuracyPValue 0.0000000
McnemarPValue NaN

Apply over test dataset

Finally, we can apply the model over the testing dataset, for which we do not have the observed classe. The predicted values for the 20 tests are:

pt_rf <- predict(model_rf, newdata=plm.test_cases)
pt_lda <- predict(model_lda, newdata=plm.test_cases)
pt_gbm <- predict(model_gbm, newdata=plm.test_cases)

pt_comb <- data.frame(pt_rf, pt_lda, pt_gbm)

exercise_pred <- predict(model_comb, pt_comb)

kable(data.frame(problem_id = plm.test_cases$problem_id, predicted_classe = exercise_pred))
problem_id predicted_classe
1 B
2 A
3 B
4 A
5 A
6 E
7 D
8 B
9 A
10 A
11 B
12 C
13 B
14 A
15 E
16 E
17 A
18 B
19 B
20 B