# Question 4.2: Iris clustering with k-means

data(iris)

dim(iris)
names(iris)
head(iris)
table(iris$Species)

# Step 2: I will Keep only the four predictors for k-means
# Species is not included.

iris_predictors <- iris[, 1:4]

head(iris_predictors)
summary(iris_predictors)

# Step 3: Test possible numbers of clusters
# Elbow Method

set.seed(6501)

k_values <- 2:6
within_ss <- numeric(length(k_values))

for (i in seq_along(k_values)) {
  
  model <- kmeans(
    iris_predictors,
    centers = k_values[i],
    nstart = 50
  )
  
  within_ss[i] <- model$tot.withinss
}

data.frame(
  k = k_values,
  total_within_cluster_sum_of_squares = within_ss
)

plot(
  k_values,
  within_ss,
  type = "b",
  xlab = "Number of clusters (k)",
  ylab = "Total within-cluster sum of squares",
  main = "Elbow Method for Iris Data"
)

# Step 4: Run k-means with k = 3 using all four predictors

set.seed(6501)

iris_kmeans_all <- kmeans(
  iris_predictors,
  centers = 3,
  nstart = 50
)

cluster_table_all <- table(
  Cluster = iris_kmeans_all$cluster,
  Species = iris$Species
)

cluster_table_all


# Step 5: Compare all predictor combinations with k = 3

cluster_accuracy <- function(selected_predictors, actual_species) {
  
  set.seed(6501)
  
  model <- kmeans(
    selected_predictors,
    centers = 3,
    nstart = 50
  )
  
  comparison_table <- table(
    Cluster = model$cluster,
    Species = actual_species
  )
  
  sum(apply(comparison_table, 1, max)) / length(actual_species)
}

predictor_names <- names(iris)[1:4]

combination_results <- data.frame(
  predictors = character(),
  accuracy = numeric()
)

for (size in seq_along(predictor_names)) {
  
  combinations <- combn(
    predictor_names,
    size,
    simplify = FALSE
  )
  
  for (variables in combinations) {
    
    current_accuracy <- cluster_accuracy(
      iris[, variables, drop = FALSE],
      iris$Species
    )
    
    combination_results <- rbind(
      combination_results,
      data.frame(
        predictors = paste(variables, collapse = " + "),
        accuracy = current_accuracy
      )
    )
  }
}

combination_results <- combination_results[
  order(-combination_results$accuracy),
]

combination_results


# Step 6: Final best k-means model

best_predictors <- c(
  "Petal.Length",
  "Petal.Width"
)

set.seed(6501)

best_iris_model <- kmeans(
  iris[, best_predictors],
  centers = 3,
  nstart = 50
)

best_cluster_table <- table(
  Cluster = best_iris_model$cluster,
  Species = iris$Species
)

best_accuracy <- sum(
  apply(best_cluster_table, 1, max)
) / nrow(iris)

best_cluster_table
best_accuracy