2.1 Identifier variable

The sample_code_number variable is an identifier and not relevant in a machine learning approach.
However, we have 699 patients, but only 645 unique identifiers.

# check number of unique identifiers => 645
nlevels(cancer$sample_code_number)

## [1] 645

Which observations have the same identifier ?

# extract multiple observations with the same identifier
multi_obs <- cancer %>% 
  count(sample_code_number) %>% 
  arrange(desc(n)) %>% 
  filter(n > 1)

multi_obs %>% 
  head(10) %>% 
  kable() %>% 
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")

Some observations occur up to 6 times using the same identifier !!
Since we don’t have any information related to the exact time of evaluation, we can assume that some patients were followed over several months and were evaluated several times during the campaign.

However, a more concerning issue is that some observations (8 in total) are redundant in the database, which could bias our analysis.

# redundant observations
duplicated_obs <- cancer %>% 
  filter(duplicated(.))

# check duplicated observations
cancer %>%
  filter(sample_code_number %in% duplicated_obs$sample_code_number) %>%
  arrange(sample_code_number) %>% 
  kable() %>% 
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left") %>% 
  scroll_box(width = "800px")

As we don’t know why these data are redundant but could bias the analysis, we will delete them.

# remove redundant observations and identifier column
cancer <- cancer %>% 
  filter(!duplicated(.)) %>% 
  select(-sample_code_number)

glimpse(cancer)

## Observations: 691
## Variables: 10
## $ clump_thickness             <dbl> 5, 5, 3, 6, 4, 8, 1, 2, 2, 4, 1, 2, …
## $ uniformity_of_cell_size     <dbl> 1, 4, 1, 8, 1, 10, 1, 1, 1, 2, 1, 1,…
## $ uniformity_of_cell_shape    <dbl> 1, 4, 1, 8, 1, 10, 1, 2, 1, 1, 1, 1,…
## $ marginal_adhesion           <dbl> 1, 5, 1, 1, 3, 8, 1, 1, 1, 1, 1, 1, …
## $ single_epithelial_cell_size <dbl> 2, 7, 2, 3, 2, 7, 2, 2, 2, 2, 1, 2, …
## $ bare_nuclei                 <dbl> 1, 10, 2, 4, 1, 10, 10, 1, 1, 1, 1, …
## $ bland_chromatin             <dbl> 3, 3, 3, 3, 3, 9, 3, 3, 1, 2, 3, 2, …
## $ normal_nucleoli             <dbl> 1, 2, 1, 7, 1, 7, 1, 1, 1, 1, 1, 1, …
## $ mitoses                     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 5, 1, 1, 1, …
## $ class                       <fct> 2, 2, 2, 2, 2, 4, 2, 2, 2, 2, 2, 2, …

2.2 Response variable

The class variable is our variable of interest, and can be recoded as a factor using appropriate labels.

# recode 'class' variable
cancer <- cancer %>% 
  mutate(class = as.factor(class),
         class = fct_recode(class, "benign" = "2", "malignant" = "4"))

# Plot distribution
ggplot(cancer, aes(x = class, fill = class)) +
  geom_bar(show.legend = FALSE) +
  geom_label(aes(label = paste(..count.., " (", round(prop.table(..count..) * 100, 2), "%)")), 
             stat = "count", 
             position = position_stack(vjust = 0.75),
             show.legend = FALSE) +
  labs(title = "Tumor type distribution", x = "Class", y = "Frequency")

Even though the data is clearly imbalanced, we wont’ treat this issue in this analysis.

2.3 Handling missing values

summary(cancer)

##  clump_thickness  uniformity_of_cell_size uniformity_of_cell_shape
##  Min.   : 1.000   Min.   : 1.00           Min.   : 1.000          
##  1st Qu.: 2.000   1st Qu.: 1.00           1st Qu.: 1.000          
##  Median : 4.000   Median : 1.00           Median : 1.000          
##  Mean   : 4.427   Mean   : 3.13           Mean   : 3.201          
##  3rd Qu.: 6.000   3rd Qu.: 5.00           3rd Qu.: 5.000          
##  Max.   :10.000   Max.   :10.00           Max.   :10.000          
##                                                                   
##  marginal_adhesion single_epithelial_cell_size  bare_nuclei    
##  Min.   : 1.000    Min.   : 1.000              Min.   : 1.000  
##  1st Qu.: 1.000    1st Qu.: 2.000              1st Qu.: 1.000  
##  Median : 1.000    Median : 2.000              Median : 1.000  
##  Mean   : 2.825    Mean   : 3.211              Mean   : 3.538  
##  3rd Qu.: 4.000    3rd Qu.: 4.000              3rd Qu.: 6.000  
##  Max.   :10.000    Max.   :10.000              Max.   :10.000  
##                                                NA's   :16      
##  bland_chromatin  normal_nucleoli     mitoses             class    
##  Min.   : 1.000   Min.   : 1.000   Min.   : 1.000   benign   :453  
##  1st Qu.: 2.000   1st Qu.: 1.000   1st Qu.: 1.000   malignant:238  
##  Median : 3.000   Median : 1.000   Median : 1.000                  
##  Mean   : 3.436   Mean   : 2.883   Mean   : 1.593                  
##  3rd Qu.: 5.000   3rd Qu.: 4.000   3rd Qu.: 1.000                  
##  Max.   :10.000   Max.   :10.000   Max.   :10.000                  
## 

We have 16 missing values from the bare_nuclei variable. How do they relate to the response variable ?

p1 <- cancer %>% 
  filter(is.na(bare_nuclei)) %>% 
  ggplot(aes(x = class, fill = class)) +
    geom_bar(show.legend = FALSE) +
    labs(title = "Tumor type distribution", x = "Class",
         y = "Frequency", subtitle = "(missing values)")

p2 <- ggplot(cancer, aes(x = bare_nuclei, color = class)) +
  geom_density(show.legend = FALSE) +
  labs(title = "'Bare Nuclei' distribution \n per tumor class")

gridExtra::grid.arrange(p1, p2, ncol = 2)

## Warning: Removed 16 rows containing non-finite values (stat_density).

We can see that the bare_nuclei distribution is very different for each class value (benign or malignant). So using a basic mean imputation may not be a good solution.

library(mice)
# run multiple imputation ; 
# do not into account the 'id' variable (not relevant for imputation) ; 
# set the seed for reproducibility
cancer_impute <- mice(data = cancer, m = 8, seed = 42, printFlag = FALSE)

# have a a look at all imputed values for the 5 imputations
cancer_impute$imp$bare_nuclei

##      1  2  3  4 5  6 7  8
## 24  10  3 10 10 7 10 8 10
## 41   5  7  2  5 5  5 7  5
## 140  1  1  1  1 1  1 1  1
## 146  1  1  1  1 1  5 1  2
## 159  1  1  1  1 2  2 1  1
## 165  1  1  1  1 1  4 1  1
## 235  1  5  5  1 1  1 1  1
## 249  1  1  1  1 8  1 1  1
## 271  1  1  8  1 1  2 2  1
## 288  6 10  3  6 3 10 9  1
## 290  1  1  5  1 1  1 1  1
## 293  1  1  3  1 1  1 3  1
## 311  8  2  2  1 8  1 7  1
## 317  1  1  1  1 8  1 1  2
## 406  1  1  1  1 1  1 1  1
## 611  1  1  1  1 1  1 1  3

# extract one of the imputed datasets, and add a tracking column (evaluated on rownames)
set.seed(42)
cancer_complete <- mice::complete(data = cancer_impute, action = sample(x = 1:cancer_impute$m, size = 1))
cancer_complete$bare_nuclei_miss <- ifelse(rownames(cancer_complete) %in% rownames(cancer_impute$imp$bare_nuclei), 1, 0)

sum(is.na(cancer_complete))

## [1] 0

# delete observations with missing values
cancer_complete <- cancer[complete.cases(cancer), ]

2.4 Other variables

Quick look at other variables’ distribution.

cancer_complete %>% 
  gather(key = "variable", value = "value", clump_thickness:mitoses) %>% 
  ggplot(aes(x = value, fill = class, color = class)) +
    geom_density(alpha = 0.3) +
    facet_wrap(~variable, scales = "free_y")

Some variables, such as bare_nuclei or clump_thickness, seem to create a clear separation between benign and malignant cases.

3.1 Training/Testing datasets

In order to fine tune parameters’ models, a validation set is needed, in addition to training and testing sets.
Here, we won’t need a validation dataset. Indeed, random forest has a convenient tuning wrapper function tuneRF() that already split the data into training and validation set and perform, automatically, the cross-validation (10 fold) of the data.

# create indices for a 80/20 split
set.seed(31415)
idx <- sample(x = 1:2, size = nrow(cancer_complete), replace = TRUE, prob = c(0.8, 0.2))

# create train/valid/test sets
cancer_train <- cancer_complete[idx == 1, ]
cancer_test <- cancer_complete[idx == 2, ]

Check distribution in all new sets.

prop.table(table(cancer_train$class))

## 
##    benign malignant 
## 0.6534091 0.3465909

prop.table(table(cancer_test$class))

## 
##    benign malignant 
## 0.6394558 0.3605442

3.2 Scaling

Since all variables are on the same 1-10 scale, we don’t need to scale the data the training variables.

3.3 Logistic regression

model_glm <- glm(class ~ ., data = cancer_train, family = "binomial")
summary(model_glm)

## 
## Call:
## glm(formula = class ~ ., family = "binomial", data = cancer_train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.4016  -0.1252  -0.0643   0.0280   2.5307  
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 -10.07249    1.30683  -7.708 1.28e-14 ***
## clump_thickness               0.56222    0.15897   3.537 0.000405 ***
## uniformity_of_cell_size       0.01942    0.21580   0.090 0.928304    
## uniformity_of_cell_shape      0.24103    0.24071   1.001 0.316657    
## marginal_adhesion             0.30872    0.12484   2.473 0.013401 *  
## single_epithelial_cell_size   0.16143    0.16399   0.984 0.324908    
## bare_nuclei                   0.33035    0.09378   3.523 0.000427 ***
## bland_chromatin               0.44280    0.17722   2.499 0.012467 *  
## normal_nucleoli               0.19318    0.11598   1.666 0.095800 .  
## mitoses                       0.58425    0.38337   1.524 0.127509    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 681.448  on 527  degrees of freedom
## Residual deviance:  91.238  on 518  degrees of freedom
## AIC: 111.24
## 
## Number of Fisher Scoring iterations: 8

As seen in graphs, clump_thickness and bare_nuclei seem to be significant variables, with higher values of these parameters leading to a higher probability of malignant tumor.

3.4 Random Forest

Random forest has a convenient wrapper tuning function to determine the best parameters.
We will tune the number of variables in each tree mtry using tuneRF() function ; then the number of trees in each model will be estimated ‘by hand’.

set.seed(42)
tune_rf <- tuneRF(x = cancer_train[, 1:9],
                  y = cancer_train$class,
                  ntreeTry = 500)

## mtry = 3  OOB error = 2.84% 
## Searching left ...
## mtry = 2     OOB error = 3.22% 
## -0.1333333 0.05 
## Searching right ...
## mtry = 6     OOB error = 3.98% 
## -0.4 0.05

We can see that the classification error is minimized when mtry is equal to 3.
Let’s build a model using this optimized value of mtry, and visualize the Out-Of-Bag error versus the number of trees.

set.seed(42)
model_rf <- randomForest(class ~ ., data = cancer_train, mtry = 3, ntree = 2000)
plot(model_rf)

We can see that over around 700 trees, we don’t have any significant improvement in the classification error value. We will pick this value as the optimized one.

set.seed(42)
rf_opt <- randomForest(class ~ ., data = cancer_train, mtry = 3, ntree = 700)

Plot the variable importance.

varImpPlot(rf_opt)

Contrary to the logistic regression model, uniformity_of_cell_size and uniformity_of_cell_shape seem to be the more relevant variables in the random forest model.

3.5 Models’ evaluation

Let’s predict values on unseen data (testing set) for each of the 2 created models.

# logistic regression
pred_logis <- predict(model_glm, newdata = cancer_test[, 1:9], type = "response")
pred_logis <- ifelse(pred_logis < 0.5, "benign", "malignant")

# optimized RF
pred_rf <- predict(rf_opt, newdata = cancer_test[, 1:9])

# build a dataframe with accuracy of each model
accs_df <- tibble(model = c("logis_reg", "rf_opt"),
                  metric = c(accuracy(actual = cancer_test$class, 
                                      predicted = pred_logis),
                             accuracy(actual = cancer_test$class, 
                                      predicted = pred_rf)))

ggplot(accs_df, aes(x = fct_reorder(model, metric), y = metric)) +
  geom_col(width = 0.6) +
  geom_text(aes(label = scales::percent(metric)), 
            position = position_stack(vjust = 0.94)) +
  coord_flip() +
  labs(title = "Model's accuracy", x = "Model", y = "Accuracy")

It seems that in our case, both models perform similarly.