---
title: "logistic-regression"
output: html_document
date: "2026-02-22"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## R Markdown

```{r}
library(dplyr)
library(tidyr)
library(ggplot2)
library(caret)
library(pROC)

ecommerce_data <- read_csv("df_clustered.csv")

# Prepare numeric features
numeric_features <- c("Tenure", "CityTier", "WarehouseToHome", "HourSpendOnApp",
                       "NumberOfDeviceRegistered", "SatisfactionScore", "NumberOfAddress",
                       "Complain", "OrderAmountHikeFromlastYear", "CouponUsed",
                       "OrderCount", "DaySinceLastOrder", "CashbackAmount")

# Scale features
ecommerce_scaled <- ecommerce_data %>%
  mutate(across(all_of(numeric_features), scale))

# Fit global logistic regression
global_lr <- glm(Churn ~ Tenure + CityTier + WarehouseToHome + HourSpendOnApp +
                   NumberOfDeviceRegistered + SatisfactionScore + NumberOfAddress +
                   Complain + OrderAmountHikeFromlastYear + CouponUsed +
                   OrderCount + DaySinceLastOrder + CashbackAmount,
                 data = ecommerce_scaled, family = binomial)

summary(global_lr)

# Extract and plot coefficients
global_coefs <- broom::tidy(global_lr) %>%
  filter(term != "(Intercept)") %>%
  arrange(desc(abs(estimate)))

ggplot(global_coefs, aes(x = reorder(term, estimate), y = estimate, fill = estimate)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  scale_fill_gradient2(low = "darkseagreen2", mid = "white", high = "firebrick", midpoint = 0) +
  labs(title = "Appendix D: Global Regression Coefficients (Standardized)",
       x = "Feature", y = "Coefficient (standardized)", fill="Estimate") +
  theme_minimal()
```

```{r}
#Chi-squared test 
churn_cluster_table <- table(df_clustered$cluster, df_clustered$Churn)
chi_test <- chisq.test(churn_cluster_table)
print(chi_test)
```

```{r}
library(purrr)

# Scale within each cluster and fit LR
cluster_models <- df_clustered %>%
  mutate(across(all_of(numeric_features), ~ as.numeric(scale(.))), .by = cluster) %>%
  group_by(cluster) %>%
  group_map(~ {
    model <- glm(Churn ~ Tenure + CityTier + WarehouseToHome + HourSpendOnApp +
                   NumberOfDeviceRegistered + SatisfactionScore + NumberOfAddress +
                   Complain + OrderAmountHikeFromlastYear + CouponUsed +
                   OrderCount + DaySinceLastOrder + CashbackAmount,
                 data = .x, family = binomial)
    broom::tidy(model) %>% mutate(cluster = .y$cluster)
  }) %>%
  bind_rows()

# Combine with global coefficients for heatmap
all_coefs <- cluster_models %>%
  filter(term != "(Intercept)") %>%
  select(term, estimate, cluster) %>%
  mutate(cluster = paste0("Cluster ", cluster)) %>%
  bind_rows(
    global_coefs %>% select(term, estimate) %>% mutate(cluster = "Global")
  )
```

```{r}
# Heatmap of coefficients
ggplot(all_coefs, aes(x = cluster, y = reorder(term, estimate), fill = estimate)) +
  geom_tile(color = "white") +
  geom_text(aes(label = round(estimate, 2)), size = 3) +
  scale_fill_gradient2(low = "darkseagreen2", mid = "white", high = "firebrick", midpoint = 0) +
  labs(title = "Figure Churn Drivers Vary by Customer Segment",
       subtitle = "Logistic Regression Coefficients (Standardized)",
       x = "Model", y = "Feature", fill = "Coefficient") +
  theme_minimal(base_size = 12)

# Figure 16: Segment-specific logistic regression coefficients (standardised)
library(ggplot2)
library(dplyr)
library(forcats)

plot_df <- all_coefs %>%
  mutate(
    cluster = as.factor(cluster),
    term = fct_reorder(term, abs(estimate), .fun = mean, .desc = FALSE)
  )

ggplot(plot_df, aes(x = cluster, y = term, fill = estimate)) +
  geom_tile(color = "white", size = 0.6) +
  geom_text(aes(label = sprintf("%.2f", estimate)),
            size = 3.2, colour = "#111111") +
  scale_fill_gradient2(
    low = "#2F6C8C",   
    mid = "white",
    high = "#B23A48",  
    midpoint = 0,
    name = "Coefficient"
  ) +
  labs(
    title = "Figure 16: Segment-Specific Churn Drivers",
    subtitle = "Standardised logistic regression coefficients comparing segment-level and global effects..",
    x = "Segment",
    y = "Feature"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0),
    plot.subtitle = element_text(hjust = 0),
    plot.title.position = "plot",
    panel.grid = element_blank(),
    axis.title.x = element_text(margin = ggplot2::margin(t = 8)),
    axis.title.y = element_text(margin = ggplot2::margin(r = 8)),
    axis.text.x = element_text(face = "bold"),
    axis.text.y = element_text(size = 10),
    legend.position = "right",
    legend.title = element_text(face = "bold"),
    plot.margin = ggplot2::margin(10, 10, 10, 10, "pt")
  )
```

```{r}
# Global vs Cluster-Specific Regression Comparison
set.seed(42)
n_runs <- 20
results <- data.frame()

for (i in 1:n_runs) {
  train_idx <- createDataPartition(df_clustered$Churn, p = 0.8, list = FALSE)
  train_data <- df_clustered[train_idx, ]
  test_data <- df_clustered[-train_idx, ]

  # Global model
  train_scaled <- train_data %>% mutate(across(all_of(numeric_features), scale))
  test_scaled <- test_data %>% mutate(across(all_of(numeric_features), scale))

  global_fit <- glm(Churn ~ ., data = train_scaled %>%
                      select(Churn, all_of(numeric_features)), family = binomial)
  global_probs <- predict(global_fit, test_scaled, type = "response")
  global_preds <- ifelse(global_probs > 0.5, 1, 0)

  # Cluster-specific models
  cluster_probs <- numeric(nrow(test_data))
  cluster_preds <- numeric(nrow(test_data))

  for (cl in unique(df_clustered$cluster)) {
    train_cl <- train_data %>% filter(cluster == cl) %>%
      mutate(across(all_of(numeric_features), scale))
    test_cl_idx <- which(test_data$cluster == cl)
    test_cl <- test_data[test_cl_idx, ] %>%
      mutate(across(all_of(numeric_features), scale))

    cl_fit <- glm(Churn ~ ., data = train_cl %>%
                    select(Churn, all_of(numeric_features)), family = binomial)
    cl_prob <- predict(cl_fit, test_cl, type = "response")
    cluster_probs[test_cl_idx] <- cl_prob
    cluster_preds[test_cl_idx] <- ifelse(cl_prob > 0.5, 1, 0)
  }

  # Compute metrics
  results <- bind_rows(results, data.frame(
  run  = i,
  model = c("Global", "Cluster-Specific"),
  AUC  = c(auc(roc(test_data$Churn, global_probs)),
           auc(roc(test_data$Churn, cluster_probs))),
  F1   = c(confusionMatrix(factor(global_preds), factor(test_data$Churn),
                           positive = "1")$byClass["F1"],
           confusionMatrix(factor(cluster_preds), factor(test_data$Churn),
                           positive = "1")$byClass["F1"]),
  Recall = c(confusionMatrix(factor(global_preds), factor(test_data$Churn),
                             positive = "1")$byClass["Recall"],
             confusionMatrix(factor(cluster_preds), factor(test_data$Churn),
                             positive = "1")$byClass["Recall"]),
  Precision = c(confusionMatrix(factor(global_preds), factor(test_data$Churn),
                                positive = "1")$byClass["Precision"],
               confusionMatrix(factor(cluster_preds), factor(test_data$Churn),
                                positive = "1")$byClass["Precision"])
))
}

summary_tbl <- results %>%
  group_by(model) %>%
  summarise(
    AUC_mean      = mean(AUC),
    F1_mean       = mean(F1),
    Recall_mean   = mean(Recall),
    Precision_mean = mean(Precision)
  ) %>%
  select(model, AUC_mean, F1_mean, Recall_mean, Precision_mean)

summary_tbl

summary_wide <- summary_tbl %>%
  pivot_longer(cols = -model, names_to = "Metric", values_to = "Mean") %>%
  mutate(Metric = gsub("_mean", "", Metric)) %>%
  pivot_wider(names_from = model, values_from = Mean)

summary_wide
```

```{r}
plot_data <- summary_wide %>%
  pivot_longer(cols = c(`Cluster-Specific`, Global),
               names_to = "Model",
               values_to = "Score")

ggplot(plot_data,
       aes(x = Metric, y = Score, fill = Model)) +
  geom_col(position = position_dodge(width = 0.7), width = 0.6) +
  geom_text(aes(label = round(Score, 3)),
            position = position_dodge(width = 0.7),
            vjust = -0.4, size = 3) +
  scale_y_continuous(limits = c(0, 1),
                     breaks = seq(0, 1, by = 0.2)) +
  labs(title = "Cluster-Specific Models Improve Churn Detection",
       subtitle = "Average performance over 20 train/test splits",
       x = "Metric",
       y = "Score",
       fill = "Model") +
  theme_minimal() +
  theme(legend.position = "bottom")

plot_data <- summary_wide %>%
  pivot_longer(cols = c(`Cluster-Specific`, Global),
               names_to = "Model",
               values_to = "Score")

ggplot(plot_data, aes(x = Metric, y = Score, fill = Model)) +
  geom_col(position = position_dodge(width = 0.7), width = 0.6) +
  geom_text(aes(label = round(Score, 3)),
            position = position_dodge(width = 0.7),
            vjust = -0.4, size = 3.2, colour = "#111111", fontface = "bold") +
  scale_y_continuous(limits = c(0, 1),
                     breaks = seq(0, 1, by = 0.2)) +
  scale_fill_manual(values = c("Cluster-Specific" = "#2F6C8C", 
                                "Global" = "#B23A48")) +
  labs(title = "Figure 17: Global vs Segment-Specific Model Performance",
       subtitle = "Average performance over 20 train/test splits",
       x = "Metric",
       y = "Score",
       fill = "Model Type") +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0),
    plot.subtitle = element_text(hjust = 0),
    #plot.title.position = "plot",
    panel.grid.major.x = element_blank(),
    panel.grid.minor = element_blank(),
    axis.title.x = element_text(margin = ggplot2::margin(t = 8)),
    axis.title.y = element_text(margin = ggplot2::margin(r = 8)),
    axis.text.x = element_text(face = "bold"),
    legend.position = "bottom",
    legend.title = element_text(face = "bold"),
    plot.margin = ggplot2::margin(10, 10, 10, 10, "pt")
  )
```
