LAB 6.2

#Q1 Make a single plot that examines the association between radius_mean and radius_sd separately for each diagnosis (hint: aes() should have three arguments).

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.0     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(broom)
biopsies <- 
  read_csv("https://www.dropbox.com/s/0rbzonyrzramdgl/cells.csv?dl=1") |>
  mutate(diagnosis = factor(diagnosis, levels = c("B", "M")))

Rows: 569 Columns: 31
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): diagnosis
dbl (30): radius_mean, texture_mean, perimeter_mean, area_mean, smoothness_m...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

ggplot(biopsies, aes(x = radius_mean, y = radius_sd, color = diagnosis)) +
  geom_point() +
  labs(x = "Radius Mean", y = "Radius SD", color = "Diagnosis") +
  ggtitle("Association between Radius Mean and Radius SD by Diagnosis")

#Q2 Calculate the correlation between these two variables for each diagnosis.

biopsies %>%
  group_by(diagnosis) %>%
  summarize(correlation = cor(radius_mean, radius_sd)) 

# A tibble: 2 × 2
  diagnosis correlation
  <fct>           <dbl>
1 B             -0.0278
2 M              0.639 

#Q3 Give at least a two-sentence interpretation of the results in the last two questions. In particular, comment on: “Is the relationship between radius_mean and radius_sd different for benign biopsies vs. malignant biopsies?” Yes “If so, can you give an explanation for this difference?” There is an upward drift observed in the plot for malignant biopsies, and the bigger the radius standard deviation is, the bigger the radius mean becomes. However, based on the plot, most benign biospies tend to have lower radius standard deviation and smaller variation of radius mean.

#Q4 Split the data set into a roughly 80-20 train-test set split.

set.seed(20)

train_or_test <- sample(x = c("train", "test"), 
                   size = 569, 
                   replace = TRUE, 
                   prob = c(0.8, 0.2))

biopsies <- biopsies |>
    mutate(set_type = train_or_test)

biopsies_train <- biopsies |>
    filter(set_type == "train")

biopsies_test <- biopsies |>
    filter(set_type == "test")

biopsies_train

# A tibble: 439 × 32
   diagnosis radius_mean texture_mean perimeter_mean area_mean smoothness_mean
   <fct>           <dbl>        <dbl>          <dbl>     <dbl>           <dbl>
 1 M                20.6         17.8          133.      1326           0.0847
 2 M                19.7         21.2          130       1203           0.110 
 3 M                11.4         20.4           77.6      386.          0.142 
 4 M                18.2         20.0          120.      1040           0.0946
 5 M                13.7         20.8           90.2      578.          0.119 
 6 M                13           21.8           87.5      520.          0.127 
 7 M                12.5         24.0           84.0      476.          0.119 
 8 M                16.0         23.2          103.       798.          0.0821
 9 M                15.8         17.9          104.       781           0.0971
10 M                19.2         24.8          132.      1123           0.0974
# ℹ 429 more rows
# ℹ 26 more variables: compactness_mean <dbl>, concavity_mean <dbl>,
#   concave.points_mean <dbl>, symmetry_mean <dbl>,
#   fractal_dimension_mean <dbl>, radius_sd <dbl>, texture_sd <dbl>,
#   perimeter_sd <dbl>, area_sd <dbl>, smoothness_sd <dbl>,
#   compactness_sd <dbl>, concavity_sd <dbl>, concave.points_sd <dbl>,
#   symmetry_sd <dbl>, fractal_dimension_sd <dbl>, radius_worst <dbl>, …

biopsies_test

# A tibble: 130 × 32
   diagnosis radius_mean texture_mean perimeter_mean area_mean smoothness_mean
   <fct>           <dbl>        <dbl>          <dbl>     <dbl>           <dbl>
 1 M                18.0         10.4          123.      1001           0.118 
 2 M                20.3         14.3          135.      1297           0.100 
 3 M                12.4         15.7           82.6      477.          0.128 
 4 B                13.5         14.4           87.5      566.          0.0978
 5 M                14.6         21.5           97.4      645.          0.105 
 6 M                18.6         20.2          122.      1094           0.0944
 7 M                17.0         24.0          113.       899.          0.120 
 8 B                13.0         18.4           82.6      524.          0.0898
 9 M                18.2         18.7          120.      1033           0.115 
10 M                14.2         22.2           96.4      646.          0.105 
# ℹ 120 more rows
# ℹ 26 more variables: compactness_mean <dbl>, concavity_mean <dbl>,
#   concave.points_mean <dbl>, symmetry_mean <dbl>,
#   fractal_dimension_mean <dbl>, radius_sd <dbl>, texture_sd <dbl>,
#   perimeter_sd <dbl>, area_sd <dbl>, smoothness_sd <dbl>,
#   compactness_sd <dbl>, concavity_sd <dbl>, concave.points_sd <dbl>,
#   symmetry_sd <dbl>, fractal_dimension_sd <dbl>, radius_worst <dbl>, …

#Q5 Using the training data, fit a simple logistic regression model that predicts the diagnosis using the mean of the texture index.

m1 <- glm(diagnosis ~ texture_mean, data = biopsies_train, family = "binomial")
m1 |>
  tidy()|>
  select(term,estimate, std.error)

# A tibble: 2 × 3
  term         estimate std.error
  <chr>           <dbl>     <dbl>
1 (Intercept)    -5.26     0.599 
2 texture_mean    0.240    0.0297

#Q6 Using a threshold of .5, What would your model predict for a biopsy with a mean texture of 15? What probability does it assign to that outcome?

new_data <- data.frame(texture_mean = 15)
prediction <- predict(m1, newdata = new_data, type = "response")
print(prediction)

       1 
0.160454 

predicted_diagnosis <- ifelse(prediction > 0.5, "M", "B")
print(predicted_diagnosis)

  1 
"B" 

#Q7 Calculate and report two misclassification rates for your simple model: first on the training data and then on the testing data.

biopsies <- biopsies_train |>
  mutate(p_hat_m1 = predict(object = m1,
                            newdata = biopsies_train, 
                            type = "response"),
         y_hat_m1 = ifelse(test = p_hat_m1> .5,
                           yes = "M",
                           no = "B"))
biopsies |>
  mutate(misclassified = (diagnosis !=y_hat_m1)) |>
  summarise(MCR =mean(misclassified))

# A tibble: 1 × 1
    MCR
  <dbl>
1 0.298

biopsies <- biopsies_test |>
  mutate(p_hat_m1 = predict(object = m1,
                            newdata = biopsies_test, 
                            type = "response"),
         y_hat_m1 = ifelse(test = p_hat_m1> .5,
                           yes = "M",
                           no = "B"))

biopsies |>
  mutate(misclassified = (diagnosis !=y_hat_m1)) |>
  summarise(MCR = mean(misclassified))

# A tibble: 1 × 1
    MCR
  <dbl>
1   0.3

#Q8 Build a more complex model to predict the diagnosis using five predictors of your choosing.

complex_model <- glm(diagnosis ~ radius_mean + texture_mean + compactness_mean + symmetry_mean + smoothness_mean, 
                     family = binomial(),
                     data = biopsies_train)
summary(complex_model)

Call:
glm(formula = diagnosis ~ radius_mean + texture_mean + compactness_mean + 
    symmetry_mean + smoothness_mean, family = binomial(), data = biopsies_train)

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -42.90841    6.00907  -7.141 9.29e-13 ***
radius_mean        1.39055    0.19534   7.119 1.09e-12 ***
texture_mean       0.39979    0.07115   5.619 1.92e-08 ***
compactness_mean  18.66494    8.66875   2.153  0.03131 *  
symmetry_mean     22.08002   11.76892   1.876  0.06064 .  
smoothness_mean   88.85644   27.51082   3.230  0.00124 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 578.1  on 438  degrees of freedom
Residual deviance: 125.6  on 433  degrees of freedom
AIC: 137.6

Number of Fisher Scoring iterations: 8

#Q9 Calculate and report two misclassification rates for your complex model: first on the training data and then on the testing data.

biopsies <- biopsies_train |>
  mutate(p_hat_m2 = predict(object = complex_model,
                            newdata = biopsies_train,
                            type = "response"),
         y_hat_m2 = ifelse(test = p_hat_m2 > .5,
                           yes = "M",
                           no = "B"))

biopsies |>
  mutate(misclassified = (diagnosis !=y_hat_m2)) |>
  summarise(MCR = mean(misclassified))

# A tibble: 1 × 1
     MCR
   <dbl>
1 0.0683

biopsies <- biopsies_test |>
  mutate(p_hat_m2 = predict(object = complex_model,
                            newdata = biopsies_test,
                            type = "response"),
         y_hat_m2 = ifelse(test = p_hat_m2 > .5,
                           yes = "M",
                           no = "B"))

biopsies |>
  mutate(misclassified = (diagnosis !=y_hat_m2)) |>
  summarise(MCR = mean(misclassified))

# A tibble: 1 × 1
     MCR
   <dbl>
1 0.0846

#Q10 Is there any evidence that your model is overfitting? Explain in at least two sentences. Yes, because the misclassification rate on the testing data (9.23%) is higher than that on the training data (5.69%), indicating that the model may have learned the noise in the training data too well and is struggling to generalize to unseen data. This discrepancy suggests that the model’s performance might degrade when applied to new data.

#Q11 Move back to your simple model for the next few questions.Report the total number of false negatives in the test data set.

test_predictions_simple <- predict(m1, newdata = biopsies_test, type = "response")
test_predicted_classes <- ifelse(test_predictions_simple > 0.5, "M", "B")
results <- data.frame(biopsies_test$diagnosis, test_predicted_classes)
results |>
  group_by(biopsies_test$diagnosis, test_predicted_classes) |>
  summarise(n = n())

`summarise()` has grouped output by 'biopsies_test$diagnosis'. You can override
using the `.groups` argument.

# A tibble: 4 × 3
# Groups:   biopsies_test$diagnosis [2]
  `biopsies_test$diagnosis` test_predicted_classes     n
  <fct>                     <chr>                  <int>
1 B                         B                         70
2 B                         M                         10
3 M                         B                         29
4 M                         M                         21

#Q12 What can you change about your classification rule to lower the number of false negatives? To lower the number of false negatives, one approach is to adjust the classification threshold used to determine predictions. Decreasing the threshold increases the sensitivity of the model, making it more likely to classify malignant cases correctly.

#Q13 Make the change you identified in the previous question and calculate the new number of false negatives.

test_predictions_simple <- predict(m1, newdata = biopsies_test, type = "response")
test_predicted_classes <- ifelse(test_predictions_simple > 0.3, "M", "B")
results <- data.frame(biopsies_test$diagnosis, test_predicted_classes)
results |>
  group_by(biopsies_test$diagnosis, test_predicted_classes) |>
  summarise(n = n())

`summarise()` has grouped output by 'biopsies_test$diagnosis'. You can override
using the `.groups` argument.

# A tibble: 4 × 3
# Groups:   biopsies_test$diagnosis [2]
  `biopsies_test$diagnosis` test_predicted_classes     n
  <fct>                     <chr>                  <int>
1 B                         B                         44
2 B                         M                         36
3 M                         B                         12
4 M                         M                         38

There are 12 false negatives.

#Q14 Calculate the testing misclassification rate using your new classification rule.

results |>
  mutate(misclassified = (biopsies_test$diagnosis != test_predicted_classes)) |>
  summarise(MCR = mean(misclassified))

        MCR
1 0.3692308

#Q15 Did your misclassification rule go up or down? Answer this question and explain why it went up or down in at least two sentences.

It went up from 0.3 to 0.369. By lowering the threshold, we skew the error, which leads to a decrease in specificity, resulting in more false positives and, consequently, a higher misclassification rate. This occurs because the model becomes more liberal in classifying cases as positive, potentially increasing both true positives and false positives.