Principal Component Analysis

The first step will be to coerce the data into a list of two variables. The first variable will be termed x and it will hold all the features without the ID or diagnosis, and the second variable, y, will hold only the diagnosis.

dat <- list(x = as.matrix(temp_dat[3:32]), y=as.factor(temp_dat$diagnosis))

# Now I will check the class of the two variables and the dimensions of the
# x feature matrix to ensure all the features I wanted to keep have been
# maintained
class(dat$x)

## [1] "matrix" "array"

class(dat$y)

## [1] "factor"

# dimensions
dim(dat$x)[1]

## [1] 569

dim(dat$x)[2]

## [1] 30

The data is now stored in “dat” and the “dat$x” variable has maintained all 569 samples along with the 30 feature observations.

Time to run some calculations.

First up, I would like to know what proportion of the samples have been classified as malignant, and what proportion of samples have been classified as benign?

mean(dat$y == "M")

## [1] 0.3725835

mean(dat$y == "B")

## [1] 0.6274165

I wonder what feature has the highest mean, and what feature has the lowest standard deviation?

which.max(colMeans(dat$x))

## area_worst 
##         24

which.min(colSds(dat$x))

## fractal_dim_se 
##             20

The area worst feature has the highest average measurement, which makes sense. Column 20 corresponds with the standard error of the fractal dimension (“coastline approximation” - 1), and it has the lowest standard deviation.

Next up, I want to center and then scale the features.

x_centered <- sweep(dat$x, 2, colMeans(dat$x))
x_scaled <- sweep(x_centered, 2, colSds(dat$x), FUN = "/")

sd(x_scaled[,1])

## [1] 1

median(x_scaled[,1])

## [1] -0.2148925

Let’s use the x_scaled data to calculate the distance between all samples. Then we can compute the average distance between the first sample, which is malignant, and the other malignant samples; followed by the average distance between the first malignant sample and other benign samples.

# average distance between all samples
d_samples <- dist(x_scaled)

# first malignant to other malignant samples - average distance
dist_MtoM <- as.matrix(d_samples)[1, dat$y == "M"]
mean(dist_MtoM[2:length(dist_MtoM)])

## [1] 9.50695

# first malignant to other benign samples - average distance
dist_MtoB <- as.matrix(d_samples)[1, dat$y == "B"]
mean(dist_MtoB)

## [1] 13.14004

Now, I will make a heat-map of the relationship between features using the scaled matrix.

d_features <- dist(t(x_scaled))
heatmap(as.matrix(d_features),col = RColorBrewer::brewer.pal(11, "Spectral"))

We can clearly see that there are high correlations between some of our features.

Let’s cut the features into 5 groups and examine the hierarchy of the features.

h <- hclust(d_features)
groups <- cutree(h, k = 5)
split(names(groups), groups)

## $`1`
##  [1] "radius_mean"       "perimeter_mean"    "area_mean"        
##  [4] "concavity_mean"    "concave_pts_mean"  "radius_se"        
##  [7] "perimeter_se"      "area_se"           "radius_worst"     
## [10] "perimeter_worst"   "area_worst"        "concave_pts_worst"
## 
## $`2`
## [1] "texture_mean"  "texture_worst"
## 
## $`3`
## [1] "smoothness_mean"   "compactness_mean"  "symmetry_mean"    
## [4] "fractal_dim_mean"  "smoothness_worst"  "compactness_worst"
## [7] "concavity_worst"   "symmetry_worst"    "fractal_dim_worst"
## 
## $`4`
## [1] "texture_se"    "smoothness_se" "symmetry_se"  
## 
## $`5`
## [1] "compactness_se" "concavity_se"   "concave_pts_se" "fractal_dim_se"

We can use these groups to attempt a delineation between benign and malignant tumors. The following code selects one feature from the 1st group and one feature from the 2nd group and makes a scatter plot.

temp_dat %>% ggplot(aes(radius_mean, texture_worst, color = diagnosis)) +
  geom_point(alpha = 0.5)

OK, we are getting closer. However, there is still quite a bit of overlap between benign and malignant tumors.

The PCA below is exploratory — it’s used to understand feature structure and visualize separability. A PCA-based model is fit and compared against the full-feature models later in the Models section.

pca <- prcomp(x_scaled)
pca_sum <- summary(pca)

# show the summary of the principal component analysis
pca_sum$importance 

##                             PC1      PC2      PC3      PC4      PC5      PC6
## Standard deviation     3.644394 2.385656 1.678675 1.407352 1.284029 1.098798
## Proportion of Variance 0.442720 0.189710 0.093930 0.066020 0.054960 0.040250
## Cumulative Proportion  0.442720 0.632430 0.726360 0.792390 0.847340 0.887590
##                              PC7       PC8       PC9      PC10      PC11
## Standard deviation     0.8217178 0.6903746 0.6456739 0.5921938 0.5421399
## Proportion of Variance 0.0225100 0.0158900 0.0139000 0.0116900 0.0098000
## Cumulative Proportion  0.9101000 0.9259800 0.9398800 0.9515700 0.9613700
##                             PC12      PC13      PC14      PC15      PC16
## Standard deviation     0.5110395 0.4912815 0.3962445 0.3068142 0.2826001
## Proportion of Variance 0.0087100 0.0080500 0.0052300 0.0031400 0.0026600
## Cumulative Proportion  0.9700700 0.9781200 0.9833500 0.9864900 0.9891500
##                             PC17      PC18      PC19      PC20      PC21
## Standard deviation     0.2437192 0.2293878 0.2224356 0.1765203 0.1731268
## Proportion of Variance 0.0019800 0.0017500 0.0016500 0.0010400 0.0010000
## Cumulative Proportion  0.9911300 0.9928800 0.9945300 0.9955700 0.9965700
##                             PC22      PC23      PC24      PC25      PC26
## Standard deviation     0.1656484 0.1560155 0.1343689 0.1244238 0.0904303
## Proportion of Variance 0.0009100 0.0008100 0.0006000 0.0005200 0.0002700
## Cumulative Proportion  0.9974900 0.9983000 0.9989000 0.9994200 0.9996900
##                              PC27      PC28       PC29       PC30
## Standard deviation     0.08306903 0.0398665 0.02736427 0.01153451
## Proportion of Variance 0.00023000 0.0000500 0.00002000 0.00000000
## Cumulative Proportion  0.99992000 0.9999700 1.00000000 1.00000000

# plot the proportion of variance explained and the principal component index
var_explained <- cumsum(pca$sdev^2/sum(pca$sdev^2))
plot(var_explained, ylab = "Proportion of Variation Explained", xlab = "PC Index")

Here we see that 90% of the variance is explained by the first 7 principal components alone. Let’s make a quick scatter plot of the first two principal components with color representing the tumor type.

data.frame(pca$x[,1:2], type = dat$y) %>%
  ggplot(aes(PC1, PC2, color = type)) +
  geom_point(alpha = 0.5)

Here we see that benign tumors tend to have higher values of PC1; we also observe that benign and malignant tumors have a similar spread of values for PC2. We also achieve much better delineation between classes than the previous analysis achieved. However, there are still some outliers present.

Let’s make box-plots of the first 7 principal components grouped by tumor type, because as we already discovered, the first 7 PC’s explain 90% of the variance.

data.frame(type = dat$y, pca$x[,1:7]) %>%
  gather(key = "PC", value = "value", -type) %>%
  ggplot(aes(PC, value, fill = type)) +
  geom_boxplot()

Now that we have better idea of the predictive challenges ahead, let’s move on to splitting the data before we scale it for modeling, to avoid any leakage from the test/validation sets into our training statistics.

Logistic Regression

train_glm <- train(y~.,data = train, method = "glm")
glm_preds <- predict(train_glm, test)
glm_acc <- mean(glm_preds == test$y)
Accuracy_Results <- tibble(Method = "GLM", Accuracy = glm_acc)
Accuracy_Results %>% knitr::kable(align='l')

Linear Discriminatory Analysis

train_lda <- train(y~., data = train, method = "lda", preProcess = "center")
lda_preds <- predict(train_lda, test)
lda_acc <- mean(lda_preds == test$y)
Accuracy_Results <- bind_rows(Accuracy_Results,
                              tibble(Method = "LDA", Accuracy = lda_acc))
Accuracy_Results %>% knitr::kable(align='l')

Let’s examine the LDA model by making a plot of the predictors and their importance in the model for classifying malignant tumors.

t(train_lda$finalModel$means) %>% 
  data.frame() %>% 
  mutate(predictor_name = rownames(.)) %>% 
  ggplot(aes(predictor_name, M, label=predictor_name)) +
  geom_point()+
  coord_flip()

Interesting, concave_pts_worst was the most important variable for classifying malignant tumors in the LDA model.

Quadratic Discriminant Analysis

train_qda <- train(y~., train, method = "qda", preProcess = "center")
qda_preds <- predict(train_qda, test)
qda_acc <- mean(qda_preds == test$y)
Accuracy_Results <- bind_rows(Accuracy_Results,
                              tibble(Method = "QDA", Accuracy = qda_acc))
Accuracy_Results %>% knitr::kable(align='l')

Let’s again make a plot of the predictors and their importance to the QDA model in classifying malignant tumors.

t(train_qda$finalModel$means) %>% data.frame() %>%
  mutate(predictor_name = rownames(.)) %>%
  ggplot(aes(predictor_name, M, label=predictor_name)) +
  geom_point()+
  coord_flip()

So, we see the level of feature importance for the QDA model is the same as the LDA model.

GamLOESS

set.seed(5, sample.kind = "Rounding")
train_loess <- train(y~., data = train, method = "gamLoess")
loess_preds <- predict(train_loess, test)
loess_acc <- mean(loess_preds == test$y)
Accuracy_Results <- bind_rows(Accuracy_Results,
                              tibble(Method = "gamLoess", Accuracy = loess_acc))
Accuracy_Results %>% knitr::kable(align='l')

K nearest neighbors

set.seed(7, sample.kind="Rounding")
train_knn <- train(y~.,
                   data = train,
                   method = "knn", 
                   tuneGrid = knn_grid)

# show the best tune
plot(train_knn)

knn_preds <- predict(train_knn, test)
knn_acc <- mean(knn_preds == test$y)
Accuracy_Results <- bind_rows(Accuracy_Results,
                              tibble(Method = "KNN", Accuracy = knn_acc))
Accuracy_Results %>% knitr::kable(align='l')

Random Forest

set.seed(9, sample.kind="Rounding")
train_rf <- train(y~., 
                  data = train,
                  method = "rf",
                  tuneGrid = rf_grid,
                  trControl = trainControl(method = "cv",
                                           number = 10),
                  importance = TRUE)

We can examine the tuning with a plot.

plot(train_rf)

Now that we know we have a reasonably tuned random forest model, let’s examine the importance of the features to the classification process.

plot(varImp(train_rf))

OK, we find that many of the most important features are the “worst” features. We also find that the importance of features for the Random Forest model is different than what we found with the LDA and QDA models.

Time to run the predictions and report the accuracy.

rf_preds <- predict(train_rf, test)
rf_acc <- mean(rf_preds == test$y)
Accuracy_Results <- bind_rows(Accuracy_Results,
                              tibble(Method = "RF", Accuracy = rf_acc))
Accuracy_Results %>% knitr::kable(align='l')

Neural Network with repeated Cross Validation

set.seed(2976, sample.kind = "Rounding")
# I will use cross validation in the train set in order to find the optimal 
# hidden layers and decay.
tc <- trainControl(method = "repeatedcv", number = 10, repeats=3)
train_nn <- train(y~., 
                  data=train, 
                  method='nnet', 
                  linout=FALSE, 
                  trace = FALSE, 
                  trControl = tc,
                  tuneGrid = nn_grid) 

Plot the different models to see the effect on accuracy

plot(train_nn)

Show the best Tune

train_nn$bestTune

##    size decay
## 46    9  0.11

nn_preds <- predict(train_nn, test)
nn_acc <- mean(nn_preds == test$y)
Accuracy_Results <- bind_rows(Accuracy_Results,
                              tibble(Method = "NNet", Accuracy = nn_acc))
Accuracy_Results %>% knitr::kable(align='l')

Recursive Partitioning - Rpart

set.seed(10, sample.kind = "Rounding")
train_rpart <- train(y ~ ., 
                     data = train,
                     method = "rpart",
                     trControl = trainControl(method = "cv", number = 10),
                     tuneGrid = rpart_grid)

# Show the best tune
train_rpart$bestTune

##       cp
## 11 0.023

# plot the resulting decision tree from the rpart model.
fancyRpartPlot(train_rpart$finalModel, yesno = 2)

rpart_preds <- predict(train_rpart, test)
rpart_acc <- mean(rpart_preds == test$y)
Accuracy_Results <- bind_rows(Accuracy_Results,
                              tibble(Method = "RPart", Accuracy = rpart_acc))
Accuracy_Results %>% knitr::kable(align='l')

All of the models we have run so far have an accuracy over 90%. We know from examining the feature importance in the models that they are weighting the features differently for their use in classifying the tumor samples. Therefore, we should expect an increase in sensitivity by creating an ensemble.

First, let’s examine the most accurate model so far with a confusion matrix.

confusionMatrix(as.factor(qda_preds), as.factor(test$y))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  B  M
##          B 56  2
##          M  1 32
##                                           
##                Accuracy : 0.967           
##                  95% CI : (0.9067, 0.9931)
##     No Information Rate : 0.6264          
##     P-Value [Acc > NIR] : 8.87e-15        
##                                           
##                   Kappa : 0.9291          
##                                           
##  Mcnemar's Test P-Value : 1               
##                                           
##             Sensitivity : 0.9825          
##             Specificity : 0.9412          
##          Pos Pred Value : 0.9655          
##          Neg Pred Value : 0.9697          
##              Prevalence : 0.6264          
##          Detection Rate : 0.6154          
##    Detection Prevalence : 0.6374          
##       Balanced Accuracy : 0.9618          
##                                           
##        'Positive' Class : B               
## 

We will compare these results to a confusion matrix of an ensemble approach. But first, let’s also check whether reducing to a handful of principal components costs us much accuracy compared to using all 30 scaled features.

PCA-based comparison model

# Fit PCA on TRAIN only, then project test onto the same rotation.
# This keeps the comparison honest — no leakage from test into the PCA basis.
train_x <- train %>% dplyr::select(-y)
test_x  <- test  %>% dplyr::select(-y)

pca_train <- prcomp(train_x)
n_pcs <- which(cumsum(pca_train$sdev^2 / sum(pca_train$sdev^2)) >= 0.90)[1]

train_pcs <- as.data.frame(pca_train$x[, 1:n_pcs]) %>% cbind(y = train$y)
test_pcs  <- as.data.frame(predict(pca_train, test_x)[, 1:n_pcs]) %>% cbind(y = test$y)

train_glm_pca <- train(y ~ ., data = train_pcs, method = "glm")
glm_pca_preds <- predict(train_glm_pca, test_pcs)
glm_pca_acc <- mean(glm_pca_preds == test$y)

Accuracy_Results <- bind_rows(Accuracy_Results,
                              tibble(Method = paste0("GLM (", n_pcs, " PCs)"),
                                     Accuracy = glm_pca_acc))
Accuracy_Results %>% knitr::kable(align='l')

Ensemble Method

I will create a new data frame that holds all of the benign predictions from the previous generative models.

ensemble <- cbind(glm = glm_preds == "B", lda = lda_preds == "B", 
                  qda = qda_preds == "B", loess = loess_preds == "B", 
                  rf = rf_preds == "B", nn = nn_preds == "B", 
                  rp = rpart_preds == "B", knn = knn_preds == "B")

I will say that if more than half of the algorithms predicted benign then the ensemble will predict benign. This also means that if less than half predict benign then the ensemble will predict malignant.

ensemble_preds <- ifelse(rowMeans(ensemble) > 0.5, "B", "M")

As an alternative to a simple majority vote, I will also fit a stacked ensemble: a logistic regression meta-learner trained on the base models’ predictions, so it can learn how much to trust each one rather than treating them all equally.

stack_train_df <- as.data.frame(ensemble) %>%
  mutate(across(everything(), as.numeric)) %>%
  cbind(y = test$y)

train_stack <- train(y ~ ., data = stack_train_df, method = "glm")

# Inspect learned weights per base model
summary(train_stack$finalModel)$coefficients %>% knitr::kable(align = 'l')

Let’s check the accuracy of the ensemble method against the test set and present the findings by attaching it to our running table.

ensemble_acc <- mean(ensemble_preds == test$y)
Accuracy_Results <- bind_rows(Accuracy_Results,
                              tibble(Method = "Ensemble", 
                                     Accuracy = ensemble_acc))
Accuracy_Results %>% knitr::kable(align='l')

Now, let’s examine a confusion matrix to see if we have better results than the QDA model had on it’s own.

confusionMatrix(as.factor(ensemble_preds), as.factor(test$y))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  B  M
##          B 57  3
##          M  0 31
##                                           
##                Accuracy : 0.967           
##                  95% CI : (0.9067, 0.9931)
##     No Information Rate : 0.6264          
##     P-Value [Acc > NIR] : 8.87e-15        
##                                           
##                   Kappa : 0.9283          
##                                           
##  Mcnemar's Test P-Value : 0.2482          
##                                           
##             Sensitivity : 1.0000          
##             Specificity : 0.9118          
##          Pos Pred Value : 0.9500          
##          Neg Pred Value : 1.0000          
##              Prevalence : 0.6264          
##          Detection Rate : 0.6264          
##    Detection Prevalence : 0.6593          
##       Balanced Accuracy : 0.9559          
##                                           
##        'Positive' Class : B               
## 

OK! The ensemble method has improved the sensitivity. We also see that the negative prediction value is 1. This means that every time (100% of the time) the ensemble method predicted malignancy, the tumor was indeed malignant.

Let’s repeat the above steps with the data from the first split, and run the final validation on the hold out set.