1. Introduction

This report presents a comprehensive analysis of the Diabetes Risk Factors & Diagnosis Dataset, including data exploration, cleaning, regression modeling, and cluster analysis to identify diabetes risk subgroups.

1.1 Core objectives

Complete the full-process dataset processing: Clarify the basic characteristics (scale, variable types, distribution, etc.) and quality issues (duplicate values, outliers, class imbalance, etc.) of the dataset through exploratory analysis, and optimize data quality through data cleaning to provide a reliable data foundation for subsequent modeling; Build regression models: Based on patients’ demographic characteristics (age, gender, etc.) and clinical indicators (HbA1c, BMI, etc.), compare linear and nonlinear machine learning algorithms to construct a high-precision blood glucose level prediction model and identify key factors affecting blood glucose; Explore diabetes risk subgroups: Identify patient subgroups with different diabetes risk profiles through cluster analysis, and clarify the differences in clinical characteristics and diabetes prevalence among each subgroup; Realize individual risk assessment: Construct an interpretable classification model to accurately determine the diabetes risk level of individual patients, and reveal the core clinical factors driving risk, providing reference for stratified management and early screening of diabetes risk.

2. Data Exploration & Validation

2.1 Overview

The core objective of data exploration is to understand the basic characteristics of the dataset, the distribution patterns of variables, and the quality of the data, and to identify potential problems such as anomalies, imbalances, and rare categories in the data, so as to provide a clear basis for subsequent data cleaning and data analysis.

2.2 Raw Dataset Basic Information

The raw dataset consisted of 100,000 rows × 16 columns, spanning from 2015 to 2022. It was primarily used for two types of supervised learning tasks: classification and regression. In terms of sample size, 100,000 records constitute a large-scale dataset, providing ample sample support for subsequent modeling. The seven-year time span also ensures the temporal representativeness of the data, avoiding the randomness of data from a single year.

2.3 Data Integrity and Uniqueness

Missing Values: There are no NA values in the dataset, indicating that no information was missed during data collection, but there are outliers that need to be addressed during data cleaning. Duplicate Rows: 14 duplicate records were detected, accounting for 0.014%. Although the percentage is extremely low, duplicate records can lead to statistical bias (such as overestimating the proportion of a certain feature) and overfitting in model training. Therefore, they still need to be removed during the cleaning phase.

2.4 Interpretation of the target variable (diabetes status)

The distribution of the target variable diabetes is as follows: 91,500 entries (91.5%) for those without the disease (0) and 8,500 entries (8.5%) for those with the disease (1), exhibiting a highly unbalanced characteristic. The proportion of diabetic patients in the sample is far lower than that of non-patients, which is consistent with the basic characteristics of diabetes prevalence in the general population (the global prevalence of diabetes in adults is approximately 10%), but it also reflects the “natural imbalance” of the data.

# Module 1: Data Exploration Core Code
# Basic dataset validation
cat("Title: Diabetes Risk Factors & Diagnosis Dataset\n")
## Title: Diabetes Risk Factors & Diagnosis Dataset
cat("Purpose: This dataset can be used for supervised learning, e.g.,\n")
## Purpose: This dataset can be used for supervised learning, e.g.,
cat("- Classification: predict diabetes status (0/1)\n")
## - Classification: predict diabetes status (0/1)
cat("- Regression: predict continuous indicators such as HbA1c or blood glucose level\n")
## - Regression: predict continuous indicators such as HbA1c or blood glucose level
# Path compatibility (use current directory if data folder does not exist)
data_path <- "data/diabetes_dataset.csv"
if (!file.exists(data_path)) {
  data_path <- "./diabetes_dataset.csv" # Alternative path
}
diabetes_raw <- read.csv(data_path, stringsAsFactors = FALSE)

cat("Raw Dataset Basic Info\n")
## Raw Dataset Basic Info
cat("Rows:", nrow(diabetes_raw), "\nCols:", ncol(diabetes_raw))
## Rows: 100000 
## Cols: 16
cat("\nFirst 5 rows:\n")
## 
## First 5 rows:
print(head(diabetes_raw, 5))
##   year gender age location race.AfricanAmerican race.Asian race.Caucasian
## 1 2020 Female  32  Alabama                    0          0              0
## 2 2015 Female  29  Alabama                    0          1              0
## 3 2015   Male  18  Alabama                    0          0              0
## 4 2015   Male  41  Alabama                    0          0              1
## 5 2016 Female  52  Alabama                    1          0              0
##   race.Hispanic race.Other hypertension heart_disease smoking_history   bmi
## 1             0          1            0             0           never 27.32
## 2             0          0            0             0           never 19.95
## 3             0          1            0             0           never 23.76
## 4             0          0            0             0           never 27.32
## 5             0          0            0             0           never 23.75
##   hbA1c_level blood_glucose_level diabetes
## 1         5.0                 100        0
## 2         5.0                  90        0
## 3         4.8                 160        0
## 4         4.0                 159        0
## 5         6.5                  90        0
if ("year" %in% names(diabetes_raw)) {
  cat("\nYear range:", min(diabetes_raw$year, na.rm = TRUE), "~", max(diabetes_raw$year, na.rm = TRUE), "\n")
} else {
  cat("\nYear range: (No 'year' column found)\n")
}
## 
## Year range: 2015 ~ 2022
# Target variable (diabetes) distribution
diabetes_dist <- diabetes_raw %>%
  group_by(diabetes) %>%
  summarise(Count = n(), Percentage = round(n()/nrow(diabetes_raw)*100, 1))
cat("\nDiabetes Distribution:\n"); print(diabetes_dist)
## 
## Diabetes Distribution:
## # A tibble: 2 × 3
##   diabetes Count Percentage
##      <int> <int>      <dbl>
## 1        0 91500       91.5
## 2        1  8500        8.5
# Content & Structure: variable types + unique counts
cat("\nContent & Structure: Column Types & Unique Levels\n")
## 
## Content & Structure: Column Types & Unique Levels
col_type <- sapply(diabetes_raw, function(x) class(x)[1])
col_unique <- sapply(diabetes_raw, function(x) length(unique(x)))

structure_table <- data.frame(
  Variable = names(diabetes_raw),
  Type = col_type,
  Unique_Values = as.integer(col_unique),
  stringsAsFactors = FALSE
)
print(structure_table)
##                                  Variable      Type Unique_Values
## year                                 year   integer             7
## gender                             gender character             3
## age                                   age   numeric           102
## location                         location character            55
## race.AfricanAmerican race.AfricanAmerican   integer             2
## race.Asian                     race.Asian   integer             2
## race.Caucasian             race.Caucasian   integer             2
## race.Hispanic               race.Hispanic   integer             2
## race.Other                     race.Other   integer             2
## hypertension                 hypertension   integer             2
## heart_disease               heart_disease   integer             2
## smoking_history           smoking_history character             6
## bmi                                   bmi   numeric          4247
## hbA1c_level                   hbA1c_level   numeric            18
## blood_glucose_level   blood_glucose_level   integer            18
## diabetes                         diabetes   integer             2
# Data quality check (missing + duplicate values)
cat("\nData Quality Check\n")
## 
## Data Quality Check
# Missing values:
missing_summary <- diabetes_raw %>%
  summarise_all(~sum(is.na(.))) %>%
  pivot_longer(everything(), names_to = "Variable", values_to = "Missing_Count") %>%
  filter(Missing_Count > 0)
if (nrow(missing_summary) == 0) {
  cat("1) Missing Values: None detected (NA)\n")
} else {
  cat("1) Missing Values: Found NA in these columns:\n")
  print(missing_summary)
}
## 1) Missing Values: None detected (NA)
# Duplicate values:
duplicate_count <- sum(duplicated(diabetes_raw))
cat("2) Duplicate Rows:", duplicate_count,
    "(", round(duplicate_count/nrow(diabetes_raw)*100, 4), "%)\n")
## 2) Duplicate Rows: 14 ( 0.014 %)

3. Numeric & Categorical Variable Analysis

3.1 Variable Structure and Types

The dataset contains 16 variables, divided into two categories: numerical and categorical/character variables: Numerical variables Covering year, age, 5 ethnic identifiers (race series), hypertension, heart disease, body mass index (BMI), glycated hemoglobin (hbA1c), blood glucose, and the target variable (diabetes); among them, ethnic identifier, hypertension, heart disease, and diabetes are binary variables (0/1 encoding). Categorical/Characteristic variables Gender, location, and smoking history, containing 3, 55, and 6 classes of values respectively.

3.2 Key Variable Statistical Characteristics

Age: Minimum 0.08 (abnormal for very young children), median 43 years, maximum 80 years, 1st percentile 1.08 years, indicating the presence of infants and young children (inconsistent with the target population for diabetes studies). Gender: “Female” (58,552 entries, 58.55%), “Male” (41,430 entries, 41.43%), “Other” (18 entries, 0.018%), “Other” is a rare category. BMI: Minimum 10.01, median 27.32, maximum 95.69, 99th percentile 48.79, indicating the presence of extremely high BMI outliers. HbA1c_level: Median 5.8, Mean 5.53, Maximum 9.0, 99th percentile 8.8 (≥6.5 is the clinical diagnostic threshold for diabetes). Blood_glucose_level: Median 140, Mean 138.06, Maximum 300, 99th percentile 280, indicating the presence of abnormally high blood glucose levels.

3.3 Distribution of Categorical Variables

Smoking history: “No Info” (not recorded, 35816 entries, 35.816%) and “never” (never smoked, 35095 entries, 35.095%) accounted for over 70%. The high proportion of the “No Info” category may affect the association analysis. Location: Covering 55 locations, the sample size is evenly distributed across regions (e.g., Hawaii and Iowa each have 2038 entries), indicating good data representativeness.

# Numeric Columns (Range & Quantiles)
# Quantiles: min, 1%, 25%, 50% (median), mean, 75%, 99%, max, standard deviation
cat("\nNumeric Columns (Range & Quantiles)\n")
## 
## Numeric Columns (Range & Quantiles)
num_cols <- names(diabetes_raw)[sapply(diabetes_raw, is.numeric)]
if (length(num_cols) == 0) {
  cat("No numeric columns found.\n")
} else {
  for (col in num_cols) {
    x <- diabetes_raw[[col]]
    cat("\n[", col, "]\n", sep = "")
    print(c(
      Minimum = min(x, na.rm = TRUE),
      P01_1Percent = quantile(x, 0.01, na.rm = TRUE),
      Q1_25Percent = quantile(x, 0.25, na.rm = TRUE),
      Median_50Percent = median(x, na.rm = TRUE),
      Mean_Average = mean(x, na.rm = TRUE),
      Q3_75Percent = quantile(x, 0.75, na.rm = TRUE),
      P99_99Percent = quantile(x, 0.99, na.rm = TRUE),
      Maximum = max(x, na.rm = TRUE),
      SD_StdDev = sd(x, na.rm = TRUE)
    ))
  }
}
## 
## [year]
##           Minimum   P01_1Percent.1%  Q1_25Percent.25%  Median_50Percent 
##       2015.000000       2015.000000       2019.000000       2019.000000 
##      Mean_Average  Q3_75Percent.75% P99_99Percent.99%           Maximum 
##       2018.360820       2019.000000       2019.000000       2022.000000 
##         SD_StdDev 
##          1.345239 
## 
## [age]
##           Minimum   P01_1Percent.1%  Q1_25Percent.25%  Median_50Percent 
##           0.08000           1.08000          24.00000          43.00000 
##      Mean_Average  Q3_75Percent.75% P99_99Percent.99%           Maximum 
##          41.88586          60.00000          80.00000          80.00000 
##         SD_StdDev 
##          22.51684 
## 
## [race.AfricanAmerican]
##           Minimum   P01_1Percent.1%  Q1_25Percent.25%  Median_50Percent 
##         0.0000000         0.0000000         0.0000000         0.0000000 
##      Mean_Average  Q3_75Percent.75% P99_99Percent.99%           Maximum 
##         0.2022300         0.0000000         1.0000000         1.0000000 
##         SD_StdDev 
##         0.4016648 
## 
## [race.Asian]
##           Minimum   P01_1Percent.1%  Q1_25Percent.25%  Median_50Percent 
##         0.0000000         0.0000000         0.0000000         0.0000000 
##      Mean_Average  Q3_75Percent.75% P99_99Percent.99%           Maximum 
##         0.2001500         0.0000000         1.0000000         1.0000000 
##         SD_StdDev 
##         0.4001145 
## 
## [race.Caucasian]
##           Minimum   P01_1Percent.1%  Q1_25Percent.25%  Median_50Percent 
##          0.000000          0.000000          0.000000          0.000000 
##      Mean_Average  Q3_75Percent.75% P99_99Percent.99%           Maximum 
##          0.198760          0.000000          1.000000          1.000000 
##         SD_StdDev 
##          0.399069 
## 
## [race.Hispanic]
##           Minimum   P01_1Percent.1%  Q1_25Percent.25%  Median_50Percent 
##         0.0000000         0.0000000         0.0000000         0.0000000 
##      Mean_Average  Q3_75Percent.75% P99_99Percent.99%           Maximum 
##         0.1988800         0.0000000         1.0000000         1.0000000 
##         SD_StdDev 
##         0.3991595 
## 
## [race.Other]
##           Minimum   P01_1Percent.1%  Q1_25Percent.25%  Median_50Percent 
##          0.000000          0.000000          0.000000          0.000000 
##      Mean_Average  Q3_75Percent.75% P99_99Percent.99%           Maximum 
##          0.199980          0.000000          1.000000          1.000000 
##         SD_StdDev 
##          0.399987 
## 
## [hypertension]
##           Minimum   P01_1Percent.1%  Q1_25Percent.25%  Median_50Percent 
##         0.0000000         0.0000000         0.0000000         0.0000000 
##      Mean_Average  Q3_75Percent.75% P99_99Percent.99%           Maximum 
##         0.0748500         0.0000000         1.0000000         1.0000000 
##         SD_StdDev 
##         0.2631505 
## 
## [heart_disease]
##           Minimum   P01_1Percent.1%  Q1_25Percent.25%  Median_50Percent 
##          0.000000          0.000000          0.000000          0.000000 
##      Mean_Average  Q3_75Percent.75% P99_99Percent.99%           Maximum 
##          0.039420          0.000000          1.000000          1.000000 
##         SD_StdDev 
##          0.194593 
## 
## [bmi]
##           Minimum   P01_1Percent.1%  Q1_25Percent.25%  Median_50Percent 
##         10.010000         14.600000         23.630000         27.320000 
##      Mean_Average  Q3_75Percent.75% P99_99Percent.99%           Maximum 
##         27.320767         29.580000         48.790100         95.690000 
##         SD_StdDev 
##          6.636783 
## 
## [hbA1c_level]
##           Minimum   P01_1Percent.1%  Q1_25Percent.25%  Median_50Percent 
##          3.500000          3.500000          4.800000          5.800000 
##      Mean_Average  Q3_75Percent.75% P99_99Percent.99%           Maximum 
##          5.527507          6.200000          8.800000          9.000000 
##         SD_StdDev 
##          1.070672 
## 
## [blood_glucose_level]
##           Minimum   P01_1Percent.1%  Q1_25Percent.25%  Median_50Percent 
##          80.00000          80.00000         100.00000         140.00000 
##      Mean_Average  Q3_75Percent.75% P99_99Percent.99%           Maximum 
##         138.05806         159.00000         280.00000         300.00000 
##         SD_StdDev 
##          40.70814 
## 
## [diabetes]
##           Minimum   P01_1Percent.1%  Q1_25Percent.25%  Median_50Percent 
##         0.0000000         0.0000000         0.0000000         0.0000000 
##      Mean_Average  Q3_75Percent.75% P99_99Percent.99%           Maximum 
##         0.0850000         0.0000000         1.0000000         1.0000000 
##         SD_StdDev 
##         0.2788831
# Unique levels + top10 levels + rare levels count
cat("\n\nCategorical Columns (Levels & Top Frequencies)\n\n")
## 
## 
## Categorical Columns (Levels & Top Frequencies)
cat_cols <- names(diabetes_raw)[!sapply(diabetes_raw, is.numeric)]
cat_cols <- setdiff(cat_cols, c("diabetes"))
cat_cols_to_show <- head(cat_cols, 10)

for (col in cat_cols_to_show) {
  cat("\n[", col, "]\n", sep = "")
  freq <- sort(table(diabetes_raw[[col]]), decreasing = TRUE)
  cat("Unique levels:", length(freq), "\n")
  cat("Top 10 levels:\n")
  print(head(freq, 10))

  # Count rare levels (<0.1% of total observations)
  rare_levels <- sum((freq / sum(freq)) < 0.001)
  cat("Rare levels (<0.1%):", rare_levels, "\n")
}
## 
## [gender]
## Unique levels: 3 
## Top 10 levels:
## 
## Female   Male  Other 
##  58552  41430     18 
## Rare levels (<0.1%): 1 
## 
## [location]
## Unique levels: 55 
## Top 10 levels:
## 
##     Hawaii       Iowa   Kentucky   Nebraska   Arkansas    Florida  Minnesota 
##       2038       2038       2038       2038       2037       2037       2037 
## New Jersey    Alabama   Delaware 
##       2037       2036       2036 
## Rare levels (<0.1%): 0 
## 
## [smoking_history]
## Unique levels: 6 
## Top 10 levels:
## 
##     No Info       never      former     current not current        ever 
##       35816       35095        9352        9286        6447        4004 
## Rare levels (<0.1%): 0
# Check special categories in smoking_history
cat("\nSpecial Category Check (Missing-like categories)\n")
## 
## Special Category Check (Missing-like categories)
N <- nrow(diabetes_raw)
if ("smoking_history" %in% names(diabetes_raw)) {
  noinfo_n <- sum(diabetes_raw$smoking_history == "No Info", na.rm=TRUE)
  cat("smoking_history == 'No Info':", noinfo_n, " (", round(noinfo_n/N*100, 3), "%)\n", sep="")
}
## smoking_history == 'No Info':35816 (35.816%)
if ("gender" %in% names(diabetes_raw)) {
  other_n <- sum(diabetes_raw$gender == "Other", na.rm=TRUE)
  cat("gender == 'Other':", other_n, " (", round(other_n/N*100, 4), "%)\n", sep="")
}
## gender == 'Other':18 (0.018%)

4. Exploratory Data Analysis Visualizations

4.1 Exploratory Visualization Results

Diabetes Label Distribution: Visually presents the imbalance of the target variable (far more non-diabetic samples than diabetic samples).

HbA1c Level Distribution: HbA1c is mainly concentrated in the 4-6 range, consistent with the characteristics of the non-diabetic population.

HbA1c by Diabetes: The median HbA1c in the diabetic group is significantly higher than that in the non-diabetic group, validating the strong correlation between HbA1c and diabetes.

# Module 2: 3 Key Visualizations
base_theme <- theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

cat("\n=== 3 Key Visualizations (EDA) ===\n")
## 
## === 3 Key Visualizations (EDA) ===
# 1) Bar chart: diabetes distribution
# Purpose: Check target variable balance (critical for classification tasks)
if ("diabetes" %in% names(diabetes_raw)) {
  p1 <- ggplot(diabetes_raw, aes(x = factor(diabetes), fill = factor(diabetes))) +
    geom_bar() +
    labs(title = "Diabetes Label Distribution",
         x = "diabetes (0/1)", y = "Count", fill = "diabetes") +
    scale_fill_brewer(palette = "Set2") +
    base_theme
  print(p1)
}

# 2) Histogram: HbA1c/BMI distribution
# Purpose: Inspect distribution shape and potential outliers of key clinical variables
if ("hbA1c_level" %in% names(diabetes_raw)) {
  p2 <- ggplot(diabetes_raw, aes(x = hbA1c_level)) +
    geom_histogram(bins = 40, fill = "darkseagreen3", color = "white") +
    labs(title = "HbA1c Level Distribution",
         x = "HbA1c Level", y = "Count") +
    base_theme
  print(p2)
} else if ("bmi" %in% names(diabetes_raw)) {
  p2 <- ggplot(diabetes_raw, aes(x = bmi)) +
    geom_histogram(bins = 40, fill = "steelblue", color = "white") +
    labs(title = "BMI Distribution",
         x = "BMI", y = "Count") +
    base_theme
  print(p2)
}

# 3) Boxplot: HbA1c/BMI by diabetes status
# Purpose: Visually compare key indicators between diabetes positive/negative groups
if (all(c("diabetes", "hbA1c_level") %in% names(diabetes_raw))) {
  p3 <- ggplot(diabetes_raw, aes(x = factor(diabetes), y = hbA1c_level, fill = factor(diabetes))) +
    geom_boxplot(outlier.alpha = 0.3) +
    labs(title = "HbA1c by Diabetes Status",
         x = "diabetes (0/1)", y = "HbA1c Level", fill = "diabetes") +
    scale_fill_brewer(palette = "Pastel1") +
    base_theme
  print(p3)
} else if (all(c("diabetes", "bmi") %in% names(diabetes_raw))) {
  p3 <- ggplot(diabetes_raw, aes(x = factor(diabetes), y = bmi, fill = factor(diabetes))) +
    geom_boxplot(outlier.alpha = 0.3) +
    labs(title = "BMI by Diabetes Status",
         x = "diabetes (0/1)", y = "BMI", fill = "diabetes") +
    scale_fill_brewer(palette = "Pastel1") +
    base_theme
  print(p3)
}

4.2 Summary

This exploration identified the basic characteristics and potential problems of the dataset: the data was largely complete but contained duplicate rows; the target variable was categorically imbalanced; outliers were present in variables such as age, BMI, and blood glucose; and gender and smoking history included rare/unrecorded categories. These issues need to be addressed specifically in subsequent data cleaning to improve data quality.

5. Data Cleaning Process

5.1 Overview

The core objective of data cleaning is to address quality issues in data, such as duplicate rows, outliers, and rare categories, thereby improving the accuracy, rationality, and usability of the data and providing a high-quality dataset for subsequent supervised learning modeling.

5.2 Background and Objectives

During the data exploration phase, 14 duplicate records were found in the raw data, the age included outliers under 18 years old, the gender included the rare “Other” category, and there were extreme outliers in BMI, blood glucose, and glycated hemoglobin.

The cleanup objectives are to remove completely duplicate records to ensure data uniqueness. Convert key numerical fields such as age, BMI, HbA1c, and blood glucose to numeric types and handle invalid inputs. Remove records where the key modeling variables HbA1c and blood glucose are missing. Filter outliers based on medically reasonable ranges to reduce noise from data entry or measurement errors. Eliminate combinations of HbA1c and blood glucose that clearly contradict each other to avoid data that does not conform to clinical logic and interferes with the model. Standardize the format of categorical variables such as gender and smoking history while retaining categories with practical meaning.

# =========================
# Module 3: Data Cleaning (As Required)
# goals：
# 1) Delete completely duplicate records
# 2) BMI outliers: Set the range to 15-40, and set NA for values exceeding this range, then fill in the gaps using the BMI median.
# 3) Categorical field standardization
# 4) Data validation: No missing values, no BMI out of range
# =========================

library(tidyverse, warn.conflicts = FALSE)
library(stringr)
raw_path <- data_path


# Step 0: Load raw data
diabetes_raw <- read.csv(raw_path, stringsAsFactors = FALSE)

cat("\n=== Data Cleaning Start ===\n")
## 
## === Data Cleaning Start ===
cat("Raw rows:", nrow(diabetes_raw), "| Raw cols:", ncol(diabetes_raw), "\n")
## Raw rows: 100000 | Raw cols: 16
# 1) Delete duplicate records
diabetes_clean <- diabetes_raw %>%
  distinct()

dup_removed <- nrow(diabetes_raw) - nrow(diabetes_clean)
cat("\n1) Remove Duplicates\n")
## 
## 1) Remove Duplicates
cat("Rows:", nrow(diabetes_raw), "->", nrow(diabetes_clean),
    "| Removed:", dup_removed, "duplicate rows\n")
## Rows: 100000 -> 99986 | Removed: 14 duplicate rows
# 2) Column naming normalization
diabetes_clean <- diabetes_clean %>%
  rename_with(~ gsub(":", ".", .x))

cat("\n2) Rename Columns\n")
## 
## 2) Rename Columns
cat("Renamed ':' to '.' in column names.\n")
## Renamed ':' to '.' in column names.
# 3) Outlier handling
cat("\n3) Numeric Clipping\n")
## 
## 3) Numeric Clipping
# Ensure that the key column is numeric.
diabetes_clean <- diabetes_clean %>%
  mutate(
    bmi = as.numeric(bmi),
    hbA1c_level = as.numeric(hbA1c_level),
    blood_glucose_level = as.numeric(blood_glucose_level)
  ) %>%
  mutate(
    bmi = pmin(pmax(bmi, 15), 40),                   
    # BMI limited to [15, 40]
    hbA1c_level = pmin(hbA1c_level, 8.8),            
    # HbA1c upper limit 8.8
    blood_glucose_level = pmin(blood_glucose_level, 280) 
    # Upper limit of blood sugar 280
  )

cat("BMI after clip -> Min:", round(min(diabetes_clean$bmi), 2),
    "| Max:", round(max(diabetes_clean$bmi), 2), "\n")
## BMI after clip -> Min: 15 | Max: 40
cat("HbA1c after clip -> Max:", round(max(diabetes_clean$hbA1c_level), 2), "\n")
## HbA1c after clip -> Max: 8.8
cat("Glucose after clip -> Max:", round(max(diabetes_clean$blood_glucose_level), 2), "\n")
## Glucose after clip -> Max: 280
# 4) Added glucose_group (grouping by blood glucose level)
cat("\n4) Create glucose_group\n")
## 
## 4) Create glucose_group
diabetes_clean <- diabetes_clean %>%
  mutate(
    glucose_group = case_when(
      blood_glucose_level <= 140 ~ "Low",
      blood_glucose_level <= 200 ~ "Medium",
      blood_glucose_level <= 260 ~ "High",
      TRUE ~ "Very High"
    )
  )

# 5) Data Validation
cat("\n5) Data Validation\n")
## 
## 5) Data Validation
# 5.1 Duplicate row check
dup_after <- sum(duplicated(diabetes_clean))
cat("Duplicated rows after cleaning:", dup_after, "\n")
## Duplicated rows after cleaning: 1
# 5.2 Missing value check
na_total <- sum(is.na(diabetes_clean))
cat("Total NA after cleaning:", na_total, "\n")
## Total NA after cleaning: 0
# 5.3 Range check
bmi_bad <- sum(diabetes_clean$bmi < 15 | diabetes_clean$bmi > 40, na.rm = TRUE)
hba1c_bad <- sum(diabetes_clean$hbA1c_level > 8.8, na.rm = TRUE)
glu_bad <- sum(diabetes_clean$blood_glucose_level > 280, na.rm = TRUE)

cat("BMI out-of-range count:", bmi_bad, "\n")
## BMI out-of-range count: 0
cat("HbA1c > 8.8 count:", hba1c_bad, "\n")
## HbA1c > 8.8 count: 0
cat("Glucose > 280 count:", glu_bad, "\n")
## Glucose > 280 count: 0
# 6) Save the cleaned data
out_path <- "data/diabetes_dataset_cleaned.csv"
write.csv(diabetes_clean, file = out_path, row.names = FALSE, fileEncoding = "UTF-8")

cat("\n=== Save Cleaned Dataset ===\n")
## 
## === Save Cleaned Dataset ===
cat("Saved to:", out_path, "\n")
## Saved to: data/diabetes_dataset_cleaned.csv
cat("Cleaned Dimension:", nrow(diabetes_clean), "rows x", ncol(diabetes_clean), "cols\n")
## Cleaned Dimension: 99986 rows x 17 cols
cat("=== Data Cleaning Completed ===\n")
## === Data Cleaning Completed ===

5.3 Implementation Steps

  1. Loaded the raw diabetes dataset containing 100,000 rows and 16 columns.

  2. Removed completely duplicated records using a full-row comparison, eliminating 14 duplicate rows and ensuring that each observation was unique.

  3. Standardized column names by replacing the character : with . to improve naming consistency and compatibility with data processing functions.

  4. Handled outliers in key numeric variables using a clipping approach:

BMI values were restricted to the range 15–40. HbA1c values were capped at 8.8. Blood glucose levels were capped at 280. No records were removed and no statistical imputation was applied.

  1. Created a new categorical variable, glucose_group, by grouping blood glucose levels into four categories: Low, Medium, High, and Very High.

  2. Performed data validation to confirm that there were no duplicated records, no missing values introduced during processing, and that all numeric variables fell within their predefined ranges.

  3. Saved the cleaned dataset as a CSV file for subsequent analysis.

5.4 Summary

The data cleaning process improved the quality and consistency of the diabetes dataset by removing duplicate records, standardizing variable names, and controlling extreme values in key numeric variables through clipping. A new glucose-level grouping variable was also created to support stratified analysis. After validation, the final dataset contains 99,986 rows with consistent formatting and valid numeric ranges, making it suitable for further statistical analysis and modeling tasks.

6. Regression Modeling (Linear, Tree, Random Forest)

6.1 Overview

In the data analysis stage of this project, we are committed to exploring the complex relationship between key health indicators and patients’ Blood Glucose Level through Regression Analysis. Considering the complexity and potential nonlinear characteristics of medical data, we constructed and compared three different mechanisms of machine learning models: Lasso regression (linear regularization model), random forest (nonlinear integration model) and stacking Ensemble.

After data reprocessing and strict model evaluation, the results show that there are huge differences in the performance of different algorithms. The analysis shows that Random Forest performed well in all indicators, significantly surpassing the linear model. This reveals that there is a significant nonlinear interaction in the data set, so we select the random forest as the final prediction model of this project.

6.2 Research Questions

According to the project requirements, we need to identify the problems to be solved through data processing. For the regression analysis part, we put forward the following two core issues:

RQ1: Based on the patient’s demographic characteristics (such as age, gender) and clinical indicators (such as HbA1c, BMI), can a high-precision model be built to predict blood sugar levels?

RQ2: Compared with traditional linear statistical models, can nonlinear machine learning algorithms, such as random forests, better capture the complex relationship between health indicators and blood sugar?

6.3 Research Objectives

The specific research objectives of this part are as follows:

  1. Use Lasso, random forest and stacking technology to build a prediction model.

  2. By comparing RMSE and R², quantify the generalization ability of different models on the test set.

  3. Identify the most critical factors affecting blood sugar levels through model explanatory tools and verify the necessity of nonlinear modeling.

6.4 Expected Outcomes:

  1. High-precision model: Obtain a robust regression model of R² exceeding 0.7, which can reliably predict the patient’s blood sugar.

  2. Model comparison insight: Prove the superiority of nonlinear models over linear models in processing such physiological data.

  3. Key feature identification: Identify the top three variables,such as HbA1c or age that contribute the most to the prediction results to provide reference for medical screening.

# Module 4: Regression Analysis
# ==============================================================================
# Step 1: Check and install any missing packages
# ==============================================================================
required_packages <- c("tidyverse", "caret", "glmnet", "randomForest", "Metrics", "corrplot", "gridExtra")
new_packages <- required_packages[!(required_packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)

library(tidyverse)    
library(caret)        
library(glmnet)       
library(randomForest)
library(Metrics)      
library(corrplot)     
library(gridExtra)    

# ==============================================================================
# Step 2: Data Preparation
# ==============================================================================
cat("\n--- Reading data ---\n")
## 
## --- Reading data ---
df <- read.csv("data/diabetes_dataset_cleaned.csv", stringsAsFactors = FALSE)

# 2.1 Examine the distribution of the target variable
p1 <- ggplot(df, aes(x = blood_glucose_level)) +
  geom_histogram(bins = 30, fill = "steelblue", color = "white") +
  theme_minimal() +
  labs(title = "Target variable distribution: blood glucose level", x = "blood sugar", y = "Frequency")

# 2.2 Correlation analysis of numerical variables
# Filtering numerical columns
numeric_vars <- df %>% select_if(is.numeric)
cor_matrix <- cor(numeric_vars, use = "complete.obs")

# Plotting correlation heatmaps
cat("Generating EDA charts...\n")
## Generating EDA charts...
corrplot(cor_matrix, method = "color", type = "upper",
         tl.col = "black", tl.cex = 0.8,
         title = "Correlation Matrix", mar = c(0,0,1,0))

print(p1) 

# ==============================================================================
# Step 3：Data preprocessing
# ==============================================================================
X <- df %>% select(-blood_glucose_level)
y <- df$blood_glucose_level

# Define feature type
categorical_features <- c("gender", "location", "smoking_history") 
# Ensure that categorical variables are present in the data.
categorical_features <- intersect(categorical_features, colnames(X))
numeric_features <- setdiff(colnames(X), categorical_features)

# Standardized numerical features
preprocess_steps <- preProcess(X[, numeric_features], method = c("center", "scale"))
X_numeric_processed <- predict(preprocess_steps, X[, numeric_features])

# One-hot coding classification features
if(length(categorical_features) > 0){
  dummy <- dummyVars(" ~ .", data = X[, categorical_features])
  X_categorical_encoded <- predict(dummy, newdata = X[, categorical_features])
  X_processed <- cbind(X_numeric_processed, X_categorical_encoded)
} else {
  X_processed <- X_numeric_processed
}

# Split the dataset (80% training, 20% testing)
set.seed(42)
train_index <- createDataPartition(y, p = 0.8, list = FALSE)
X_train <- X_processed[train_index, ]
X_test <- X_processed[-train_index, ]
y_train <- y[train_index]
y_test <- y[-train_index]

# Matrixing (adapted for glmnet)
X_train_matrix <- as.matrix(X_train)
X_test_matrix <- as.matrix(X_test)

# ==============================================================================
# Step 4: Modeling – This corresponds to the theoretical concept of "model fitting".
# ==============================================================================
cat("\n--- Training Lasso to return to normal ---\n")
## 
## --- Training Lasso to return to normal ---
# Model A: Lasso
lasso_cv <- cv.glmnet(x = X_train_matrix, y = y_train, alpha = 1, nfolds = 5)
lasso_model <- glmnet(x = X_train_matrix, y = y_train, alpha = 1, lambda = lasso_cv$lambda.min)

y_pred_lasso_train <- predict(lasso_model, X_train_matrix)
y_pred_lasso_test <- predict(lasso_model, X_test_matrix)

cat("\n--- Random Forest is being trained. ---\n")
## 
## --- Random Forest is being trained. ---
# Model B: Random Forest (Parameter Optimization)
rf_optimized <- randomForest(
  x = X_train, y = y_train,
  ntree = 300, mtry = 10, maxdepth = 15, nodesize = 3,
  importance = TRUE, keep.forest = TRUE
)

y_pred_rf_train <- predict(rf_optimized, X_train)
y_pred_rf_test <- predict(rf_optimized, X_test)

cat("\n--- Training a stacking model. ---\n")
## 
## --- Training a stacking model. ---
# Model C: Stacking
stack_train <- data.frame(lasso = as.vector(y_pred_lasso_train), rf = as.vector(y_pred_rf_train), y = y_train)
stack_test <- data.frame(lasso = as.vector(y_pred_lasso_test), rf = as.vector(y_pred_rf_test), y = y_test)

# Using linear regression as the meta-learner
stack_model <- lm(y ~ ., data = stack_train)
y_pred_stack <- predict(stack_model, stack_test)

# ==============================================================================
# Step 5: Model Diagnostics – “Testing the Hypothesis” of the Corresponding Theory
# ==============================================================================
cat("\n--- Execution model diagnostics (based on stacked models) ---\n")
## 
## --- Execution model diagnostics (based on stacked models) ---
# Calculate residuals
residuals_stack <- y_test - y_pred_stack
diagnostic_df <- data.frame(
  Predicted = y_pred_stack,
  Residuals = residuals_stack
)

# 5.1 Residual vs. Predicted Value 
p_resid <- ggplot(diagnostic_df, aes(x = Predicted, y = Residuals)) +
  geom_point(alpha = 0.5, color = "darkblue") +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  labs(title = "Residual Plot", x = "Predicted", y = "Residual") +
  theme_minimal()

# 5.2 Q-Q plot (check for normality)
p_qq <- ggplot(diagnostic_df, aes(sample = Residuals)) +
  stat_qq(color = "darkblue", alpha = 0.5) +
  stat_qq_line(color = "red") +
  labs(title = "Residual Analysis - Normal Q-Q Plot", x = "Theoretical quantiles", y = "Sample quantiles") +
  theme_minimal()

# 5.3 Actual value vs. predicted value
p_fit <- ggplot(data.frame(Actual = y_test, Predicted = y_pred_stack), aes(x = Actual, y = Predicted)) +
  geom_point(alpha = 0.5, color = "forestgreen") +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed") +
  labs(title = "Predicted vs Actual Plot", x = "Actual", y = "Predicted") +
  theme_minimal()

grid.arrange(p_resid, p_qq, p_fit, ncol = 3)

# ==============================================================================
# Step 6: Model Evaluation and Results Interpretation – Corresponding to the "Evaluation and Interpretation" of the Theory
# ==============================================================================
cat("\n====================== Final Model Evaluation Table ======================\n")
## 
## ====================== Final Model Evaluation Table ======================
# Define the evaluation function
calc_metrics <- function(actual, pred, model_name) {
  r2 <- R2(actual, pred)
  rmse_val <- rmse(actual, pred)
  mae_val <- mae(actual, pred)
  return(data.frame(Model = model_name, R2 = r2, RMSE = rmse_val, MAE = mae_val))
}

# Summary Results
results <- rbind(
  calc_metrics(y_test, as.vector(y_pred_lasso_test), "Lasso"),
  calc_metrics(y_test, y_pred_rf_test, "Random Forest"),
  calc_metrics(y_test, y_pred_stack, "Stacking")
)

print(results)
##           Model        R2     RMSE      MAE
## 1         Lasso 0.1616773 36.62113 30.54724
## 2 Random Forest 0.7441867 20.29623 16.84022
## 3      Stacking 0.7440471 20.32161 17.33117
cat("\n====================== Key features (based on random forest) ======================\n")
## 
## ====================== Key features (based on random forest) ======================
# Feature Importance Visualization
rf_imp <- importance(rf_optimized) %>%
  as.data.frame() %>%
  rownames_to_column("Feature") %>%
  arrange(desc(`%IncMSE`)) %>% # 
  head(10)

p_imp <- ggplot(rf_imp, aes(x = reorder(Feature, `%IncMSE`), y = `%IncMSE`)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  labs(title = "Feature Importance(Top 10)", x = "Feature", y = "%IncMSE") +
  theme_minimal()

print(p_imp)

cat("\nAll steps are complete. Check the pop-up window.\n")
## 
## All steps are complete. Check the pop-up window.

6.5 Results Summary

We use 80% of the data for training and 20% for testing. By comparing the performance of the three models on the test set, the results are shown in the following table:

1)The victory of the nonlinear model: The R² of the random forest reached 0.744, which means that the model successfully explained about 74.4% of the fluctuations in blood glucose data. In contrast, the R²of Lasso regression is only 0.162. This huge performance gap strongly proves that there is a high degree of nonlinear relationship between independent variables and dependent variables.

2)Stacking model vs. Single model: Despite the use of a complex Stacking model, its performance (R²=0.7440) is almost equal to or even slightly lower that of the random forest (R²=0.7442). According to the Oakham razor principle, under the condition of equivalent effect, we give priority to random forests with lower calculation cost and easier to explain as the final model.

3)The error was greatly reduced: the RMSE of the random forest was reduced to 20.30, which was nearly 45% lower than the 36.62 of Lasso. This shows that the prediction accuracy of the model in practical applications has been qualitatively improved.

6.6 Visualization

(1)Correlation Matrix

The linear correlation between variables is shown. It is worth noting that although the linear correlation coefficient between some variables and blood sugar is not high, they are important in random forests, which further confirms the necessity of nonlinear modeling.

(2)Predicted vs. Actual Plot

Draw for random forest models. Unlike the linear model, the chart shows that the data points are closely surrounded by the diagonal, indicating that the predicted value is highly consistent with the real value.

(3)Residual Plot

The residual plot shows that most of the errors are concentrated near zero. Although there is still a small number of extreme “long tail” phenomena, the overall fit has improved significantly compared with before.

(4)Residual Analysis - Normal Q-Q Plot

It shows an obvious S-shaped curve, indicating that the residual distribution deviates from the theoretical normal distribution at both ends and has light tail characteristics. Given that the random forest is a Non-parametric model, it does not strictly rely on the normality hypothesis of residual, and the overall prediction accuracy of the model is high (R² > 0.74), so this distribution pattern is within the acceptable range.

(5)Feature Importance

Based on random forest algorithm extraction. The figure clearly shows HbA1c_level, age and other characteristics, which are the key factors that dominate the change of blood sugar.

6.7 Summary

This regression analysis has achieved remarkable success. By optimizing the data preprocessing process, the performance of the model has been greatly improved. Random forest proved to be the best algorithm to handle the diabetes data set. It not only raises R² to more than 74%, but also effectively reduces the prediction error. The results show that the prediction of blood glucose levels cannot only rely on simple linear weighting, but must take into account the complex interaction between physiological indicators. The current random forest model has a high practical value and can be used to assist the early screening of potential risks.

7. Cluster Analysis for Diabetes Risk Subgroups

7.1 Research Question

Are there distinct subgroups of patients with different diabetes risk profiles, and what clinical factors drive diabetes risk across these subgroups?

7.2 Research Objective

The objective of processing this dataset is to examine whether diabetes risk is heterogeneous across patients and whether meaningful patient subgroups with different risk profiles can be identified. In addition, the analysis aims to determine whether an interpretable model can assess diabetes risk at the individual level and explain why a patient is classified as high risk or low risk. From this dataset, we are interested in understanding how diabetes risk varies across patients and whether this risk can be meaningfully structured into distinct profiles. Specifically, we aim to identify whether patients can be grouped into subgroups with different levels of diabetes risk, to determine how diabetes prevalence differs across these subgroups, and to understand which clinical characteristics contribute most strongly to elevated diabetes risk. In addition, we are interested in assessing whether the identified risk patterns can be reliably used to evaluate diabetes risk at the individual level.

7.3 Data Analysis Approach

The analysis integrates both unsupervised and supervised learning techniques to explore diabetes risk patterns in the dataset. First, analyze data and visualization, including histograms and summary plots, are used to examine the distribution and variability of key clinical variables. Next, K-means clustering is applied to standardized clinical features to identify patient subgroups that represent different combinations of metabolic and cardiovascular risk factors. The clustering results are further examined using principal component analysis (PCA) for visualization, and diabetes prevalence is compared across the identified subgroups to assess differences in risk profiles. Finally, an interpretable logistic regression classification model is used to predict diabetes risk at the individual level. Model performance is evaluated using accuracy and ROC-AUC, and model outputs are interpreted using Odds Ratios to explain how individual clinical features and subgroup membership contribute to diabetes risk assessment. The code implementation is provided in Module 5.

# Module 5: Cluster Analysis & Diabetes Risk Classification
# Research Question: Are there distinct subgroups of patients with different diabetes risk profiles that may require differentiated clinical attention?
# Step 0 — Initialize Environment
rm(list = ls()[!ls() %in% c("out_path")]); gc() # Preserve key variables
##           used  (Mb) gc trigger   (Mb)  max used   (Mb)
## Ncells 2757878 147.3    7547902  403.2   7547902  403.2
## Vcells 4854386  37.1  374436057 2856.8 465129859 3548.7
# Step 1 — Load Cleaned Data, Standardize Column Names & Set Target Label
dat <- read_csv(out_path) |>
  clean_names()

stopifnot(is.data.frame(dat), "diabetes" %in% names(dat))
dat$diabetes <- factor(dat$diabetes, levels = c(0,1), labels = c("No","Yes"))

# Fix HbA1c column name (if needed)
if (!("hba1c_level" %in% names(dat)) && ("hb_a1c_level" %in% names(dat))) {
  names(dat)[names(dat)=="hb_a1c_level"] <- "hba1c_level"
}

cat("Dataset Dimensions - Rows:", nrow(dat), "Cols:", ncol(dat), "\n")
## Dataset Dimensions - Rows: 99986 Cols: 17
print(table(dat$diabetes))
## 
##    No   Yes 
## 91486  8500
# Step 2 — Visualization: Diabetes Outcome Distribution
barplot(table(dat$diabetes),
        main="Diabetes Outcome Distribution",
        xlab="Diabetes Status", ylab="Patient Count")

png("01_outcome_distribution.png", width=900, height=600, res=150)
barplot(table(dat$diabetes),
        main="Diabetes Outcome Distribution",
        xlab="Diabetes Status", ylab="Patient Count")
dev.off()
## png 
##   2
# Step 3 — Visualization: Key Clinical Variable Distributions
key_vars <- c("age","bmi","hba1c_level","blood_glucose_level")

# stopifnot(all(key_vars %in% names(dat)))
#
# par(mfrow=c(2,2))
# for(v in key_vars){
#   hist(dat[[v]], main=paste("Histogram:", v), xlab=v, col="grey")
# }
# par(mfrow=c(1,1))
#
# png("02_key_variables_hist.png", width=1200, height=800, res=150)
# par(mfrow=c(2,2))
# for(v in key_vars){
#   hist(dat[[v]], main=paste("Histogram:", v), xlab=v, col="grey")
# }
# par(mfrow=c(1,1))
# dev.off()


# Step 4 — Prepare Clustering Input: Select Variables & Standardize
# Select clinically meaningful variables for clustering
cluster_vars <- c("age","bmi","hba1c_level","blood_glucose_level","hypertension","heart_disease")
stopifnot(all(cluster_vars %in% names(dat)))

x_cluster <- dat[, cluster_vars]
x_scaled  <- scale(x_cluster) # Standardize variables to equalize scale impact

# Step 5 — Determine Optimal Cluster Number: Elbow Method (WSS)
# WSS = Within-Cluster Sum of Squares
set.seed(123)
ks <- 2:6 # Test cluster numbers from 2 to 6

wss <- sapply(ks, function(k){
  kmeans(x_scaled, centers=k, nstart=10, iter.max=20)$tot.withinss
})

plot(ks, wss, type="b", pch=19,
     main="Elbow Method for Optimal k (WSS)", xlab="Number of Clusters (k)", ylab="Total Within-Cluster Sum of Squares")

png("03_elbow.png", width=900, height=600, res=150)
plot(ks, wss, type="b", pch=19,
     main="Elbow Method for Optimal k (WSS)", xlab="Number of Clusters (k)", ylab="Total Within-Cluster Sum of Squares")
dev.off()
## png 
##   2
# Step 6 — Determine Optimal Cluster Number: Silhouette Score
# Silhouette Score ranges from -1 to 1; higher values indicate better cluster separation
set.seed(123)
n_all <- nrow(x_scaled)
m_sub <- min(3000, n_all) # Use subset to reduce computational cost
idx_sub <- sample.int(n_all, m_sub)
x_sub <- x_scaled[idx_sub, , drop=FALSE]

sil <- sapply(ks, function(k){
  km_tmp <- kmeans(x_sub, centers=k, nstart=10, iter.max=20)
  ss <- silhouette(km_tmp$cluster, dist(x_sub))
  mean(ss[,3])
})

plot(ks, sil, type="b", pch=19,
     main=paste("Silhouette Score (subset n =", m_sub, ")"),
     xlab="Number of Clusters (k)", ylab="Average Silhouette Score")

png("04_silhouette.png", width=900, height=600, res=150)
plot(ks, sil, type="b", pch=19,
     main=paste("Silhouette Score (subset n =", m_sub, ")"),
     xlab="Number of Clusters (k)", ylab="Average Silhouette Score")
dev.off()
## png 
##   2
# Step 7 — Run K-means Clustering with Optimal k
k <- 2 # Selected based on Elbow and Silhouette methods
set.seed(123)
km <- kmeans(x_scaled, centers=k, nstart=20, iter.max=30)

# Assign cluster labels to original dataset
dat$cluster <- factor(km$cluster)
print(table(dat$cluster))
## 
##     1     2 
## 89475 10511
# Step 8 — Visualization: PCA Projection of Clusters (Subset)
# Reduce dimensionality to 2D using PCA for visualization
set.seed(123)
m_pca <- min(5000, nrow(x_scaled)) # Use subset for faster plotting
idx_pca <- sample.int(nrow(x_scaled), m_pca)

pca <- prcomp(x_scaled[idx_pca, , drop=FALSE], center=FALSE, scale.=FALSE)

plot(pca$x[,1], pca$x[,2],
     col=as.numeric(dat$cluster[idx_pca]), pch=19,
     xlab="Principal Component 1", ylab="Principal Component 2",
     main=paste("PCA Projection of Patient Clusters (subset n =", m_pca, ")"))
legend("topright", legend=levels(dat$cluster),
       col=1:length(levels(dat$cluster)), pch=19, title="Cluster")

png("05_cluster_pca.png", width=900, height=650, res=150)
plot(pca$x[,1], pca$x[,2],
     col=as.numeric(dat$cluster[idx_pca]), pch=19,
     xlab="Principal Component 1", ylab="Principal Component 2",
     main=paste("PCA Projection of Patient Clusters (subset n =", m_pca, ")"))
legend("topright", legend=levels(dat$cluster),
       col=1:length(levels(dat$cluster)), pch=19, title="Cluster")
dev.off()
## png 
##   2
# Step 9 — Cluster Profile Analysis: Clinical Characteristics & Diabetes Prevalence
# Calculate mean clinical characteristics per cluster
cluster_profile <- aggregate(dat[, cluster_vars],
                             by=list(cluster=dat$cluster),
                             FUN=function(x) mean(x, na.rm=TRUE))

# Calculate diabetes prevalence per cluster
cluster_prev <- aggregate(as.numeric(dat$diabetes=="Yes"),
                          by=list(cluster=dat$cluster),
                          FUN=mean)
names(cluster_prev)[2] <- "diabetes_rate"

# Add cluster size information
cluster_n <- as.data.frame(table(dat$cluster))
names(cluster_n) <- c("cluster","n")
cluster_prev <- merge(cluster_prev, cluster_n, by="cluster")

cat("\nCluster Clinical Profile (Mean Values):\n")
## 
## Cluster Clinical Profile (Mean Values):
print(cluster_profile)
##   cluster      age      bmi hba1c_level blood_glucose_level hypertension
## 1       1 39.35981 26.75319    5.490445            136.5055    0.0000000
## 2       2 63.38959 29.80256    5.830644            150.0059    0.7121111
##   heart_disease
## 1     0.0000000
## 2     0.3750357
cat("\nDiabetes Prevalence by Cluster:\n")
## 
## Diabetes Prevalence by Cluster:
print(cluster_prev)
##   cluster diabetes_rate     n
## 1       1    0.06150321 89475
## 2       2    0.28512986 10511
# Save cluster analysis results
write_csv(cluster_profile, "cluster_profile.csv")
write_csv(cluster_prev, "cluster_prevalence.csv")

# Step 10 — Visualization: Diabetes Prevalence by Cluster
barplot(cluster_prev$diabetes_rate,
        names.arg = cluster_prev$cluster,
        main="Diabetes Prevalence by Patient Cluster",
        xlab="Cluster", ylab="Diabetes Prevalence Rate")

png("06_cluster_prevalence.png", width=900, height=600, res=150)
barplot(cluster_prev$diabetes_rate,
        names.arg = cluster_prev$cluster,
        main="Diabetes Prevalence by Patient Cluster",
        xlab="Cluster", ylab="Diabetes Prevalence Rate")
dev.off()
## png 
##   2
# Step 11 — Train-Test Split for Classification Model (80/20)
set.seed(123)
n <- nrow(dat)
train_idx <- sample.int(n, size=floor(0.8*n))
train <- dat[train_idx, ]
test  <- dat[-train_idx, ]

# Step 12 — Train Interpretable Classification Model: Logistic Regression
# Include cluster membership as a predictor to capture subgroup effects
race_cols <- grep("^race_", names(dat), value=TRUE)

# Convert categorical variables to factor type
train$gender <- as.factor(train$gender); test$gender <- as.factor(test$gender)
train$smoking_history <- as.factor(train$smoking_history); test$smoking_history <- as.factor(test$smoking_history)
train$cluster <- as.factor(train$cluster); test$cluster <- as.factor(test$cluster)

# Define model variables
model_vars <- c("age","bmi","hba1c_level","blood_glucose_level",
                "hypertension","heart_disease",
                "gender","smoking_history",
                race_cols,
                "cluster")
model_vars <- intersect(model_vars, names(dat)) # Ensure variables exist in dataset

# Build model formula
formula_str <- paste("diabetes ~", paste(model_vars, collapse=" + "))
m <- glm(as.formula(formula_str), data=train, family=binomial())

# Generate predictions on test set
prob <- predict(m, newdata=test, type="response")
pred <- factor(ifelse(prob>=0.5, "Yes","No"), levels=c("No","Yes"))

# Step 13 — Numeric Model Evaluation: Confusion Matrix, Accuracy, AUC
# Confusion Matrix
cm <- table(Predicted=pred, Actual=test$diabetes)
cat("\nConfusion Matrix (Test Set):\n")
## 
## Confusion Matrix (Test Set):
print(cm)
##          Actual
## Predicted    No   Yes
##       No  18130   648
##       Yes   147  1073
# Calculate accuracy and AUC
accuracy <- sum(diag(cm))/sum(cm)
roc_obj <- roc(test$diabetes, prob)
auc_val <- auc(roc_obj)

cat("\nModel Performance Metrics (Test Set):\n")
## 
## Model Performance Metrics (Test Set):
cat("Accuracy =", round(as.numeric(accuracy),4), "\n")
## Accuracy = 0.9602
cat("AUC      =", round(as.numeric(auc_val),4), "\n")
## AUC      = 0.9604
# Save evaluation results
write_csv(data.frame(metric=c("Accuracy","AUC"),
                     value=c(as.numeric(accuracy), as.numeric(auc_val))),
          "evaluation_summary.csv")

# Step 14 — Evaluation Plot: ROC Curve
# ROC Curve shows trade-off between Sensitivity and Specificity
plot(roc_obj,
     main=paste0("ROC Curve (AUC = ", round(as.numeric(auc_val),3), ")"),
     lwd=2)
abline(a=0, b=1, lty=2, col="gray") # Reference line (random guess)

png("09_roc_curve.png", width=800, height=600, res=150)
plot(roc_obj,
     main=paste0("ROC Curve (AUC = ", round(as.numeric(auc_val),3), ")"),
     lwd=2)
abline(a=0, b=1, lty=2, col="gray")
dev.off()
## png 
##   2
# Step 15 — Evaluation Plot: Predicted Probability Distributions
# Compare probability distributions between true positive and true negative cases
hist(prob[test$diabetes=="No"], breaks=30, xlim=c(0,1),
     main="Predicted Probability Distribution by True Diabetes Status",
     xlab="Predicted Probability of Diabetes", col="lightblue")
hist(prob[test$diabetes=="Yes"], breaks=30, add=TRUE, col="lightcoral")
legend("topright", legend=c("No Diabetes","Diabetes"), fill=c("lightblue","lightcoral"), bty="n")

png("10_predicted_probability_distribution.png", width=900, height=600, res=150)
hist(prob[test$diabetes=="No"], breaks=30, xlim=c(0,1),
     main="Predicted Probability Distribution by True Diabetes Status",
     xlab="Predicted Probability of Diabetes", col="lightblue")
hist(prob[test$diabetes=="Yes"], breaks=30, add=TRUE, col="lightcoral")
legend("topright", legend=c("No Diabetes","Diabetes"), fill=c("lightblue","lightcoral"), bty="n")
dev.off()
## png 
##   2
# Step 16 — Evaluation Plot: Confusion Matrix Heatmap
cm_mat <- as.matrix(cm)

image(1:ncol(cm_mat), 1:nrow(cm_mat), t(cm_mat),
      axes=FALSE, xlab="Actual Diabetes Status", ylab="Predicted Diabetes Status",
      main="Confusion Matrix Heatmap")
axis(1, at=1:ncol(cm_mat), labels=colnames(cm_mat))
axis(2, at=1:nrow(cm_mat), labels=rownames(cm_mat))
for (i in 1:nrow(cm_mat)) for (j in 1:ncol(cm_mat)) text(j, i, cm_mat[i,j])

png("11_confusion_matrix_heatmap.png", width=800, height=650, res=150)
image(1:ncol(cm_mat), 1:nrow(cm_mat), t(cm_mat),
      axes=FALSE, xlab="Actual Diabetes Status", ylab="Predicted Diabetes Status",
      main="Confusion Matrix Heatmap")
axis(1, at=1:ncol(cm_mat), labels=colnames(cm_mat))
axis(2, at=1:nrow(cm_mat), labels=rownames(cm_mat))
for (i in 1:nrow(cm_mat)) for (j in 1:ncol(cm_mat)) text(j, i, cm_mat[i,j])
dev.off()
## png 
##   2
# Step 17 — Identify Key Diabetes Risk Drivers: Odds Ratios (OR)
# Convert logistic regression coefficients to Odds Ratios for interpretability
coef_est <- coef(m)
or <- exp(coef_est)
ci <- suppressMessages(confint(m))
or_ci <- exp(ci)

# Create Odds Ratio table with 95% Confidence Intervals
or_table <- data.frame(
  term=names(or),
  OR=as.numeric(or),
  CI_low=as.numeric(or_ci[,1]),
  CI_high=as.numeric(or_ci[,2])
)
or_table <- or_table[or_table$term!="(Intercept)", ]
or_table <- or_table[order(or_table$OR, decreasing=TRUE), ]

cat("\nTop 15 Diabetes Risk Drivers (Odds Ratios with 95% CI):\n")
## 
## Top 15 Diabetes Risk Drivers (Odds Ratios with 95% CI):
print(head(or_table, 15))
##                     term         OR    CI_low    CI_high
## 4            hba1c_level 10.3269128 9.5649455 11.1706585
## 20              cluster2  1.6255366 1.2106706  2.1848133
## 7          heart_disease  1.4717010 1.1392650  1.8978990
## 6           hypertension  1.3957743 1.0599960  1.8356866
## 8             genderMale  1.2502490 1.1550129  1.3532736
## 15 race_african_american  1.1338408 1.0031535  1.2817232
## 3                    bmi  1.1233347 1.1151538  1.1316092
## 16            race_asian  1.0713444 0.9467367  1.2124417
## 18         race_hispanic  1.0579636 0.9347348  1.1975102
## 2                    age  1.0469266 1.0443390  1.0495477
## 5    blood_glucose_level  1.0341969 1.0331125  1.0352995
## 17        race_caucasian  1.0018931 0.8846612  1.1346848
## 10   smoking_historyever  0.9798216 0.7986695  1.1997395
## 11 smoking_historyformer  0.9517365 0.8165210  1.1100683
## 12  smoking_historynever  0.8705310 0.7621313  0.9956972
write_csv(or_table, "logistic_or_drivers.csv")

# Step 18 — Visualization: Top Odds Ratio Features
topN <- 10
top_or <- head(or_table, topN)

# Adjust plot margins for long variable names
op <- par(mar = c(5, 14, 4, 2))
barplot(
  rev(top_or$OR),
  names.arg = rev(top_or$term),
  horiz = TRUE,
  las = 1,
  cex.names = 0.9,
  main = "Top 10 Diabetes Risk Drivers (Odds Ratios)",
  xlab = "Odds Ratio (OR)",
  col = "steelblue"
)
abline(v=1, lty=2, col="red") # OR=1 means no effect

par(op)

# Save plot
dev.copy(png, filename="07_top_or_drivers.png", width=1200, height=700, res=150)
## png 
##   3
dev.off()
## png 
##   2
# Step 19 — Driver Heterogeneity: Key Clinical Variables by Cluster (Boxplots)
# Compare distribution of key variables across clusters
par(mfrow=c(2,2))
for(v in key_vars){
  boxplot(dat[[v]] ~ dat$cluster,
          main=paste(v, "by Cluster"),
          xlab="Cluster", ylab=v,
          col=c("lightblue","lightcoral"))
}

par(mfrow=c(1,1))

png("08_drivers_by_cluster.png", width=1200, height=800, res=150)
par(mfrow=c(2,2))
for(v in key_vars){
  boxplot(dat[[v]] ~ dat$cluster,
          main=paste(v, "by Cluster"),
          xlab="Cluster", ylab=v,
          col=c("lightblue","lightcoral"))
}
par(mfrow=c(1,1))
dev.off()
## png 
##   2
# Save complete workspace for future reference
saveRDS(list(dat=dat, x_scaled=x_scaled, km=km, k=k, ks=ks, wss=wss, sil=sil,
             cluster_profile=cluster_profile, cluster_prev=cluster_prev,
             train=train, test=test, m=m, prob=prob, pred=pred, cm=cm,
             roc_obj=roc_obj, auc_val=auc_val, or_table=or_table),
        "workspace_final.rds")

cat("\nCluster Analysis & Classification Modeling Completed Successfully!\n")
## 
## Cluster Analysis & Classification Modeling Completed Successfully!
cat("All output files and plots have been saved to the working directory.\n")
## All output files and plots have been saved to the working directory.

7.4 Results

7.4.1 Evidence of Heterogeneous Diabetes Risk

Exploratory visualizations show substantial variability in key clinical variables such as BMI, HbA1c level, and blood glucose level. These variables display skewed distributions, indicating that a subset of patients exhibits markedly higher metabolic risk. This suggests that diabetes risk is not uniformly distributed across the population.

7.4.2 Identification of Distinct Patient Subgroups

Clustering analysis identifies two distinct patient subgroups based on standardized clinical features. Visualization using principal component analysis (PCA) shows that these subgroups occupy different regions in the feature space, indicating that the clustering captures structured differences rather than random variation.

7.4.3 Differences in Diabetes Risk Between Subgroups

The diabetes prevalence differs substantially between the identified subgroups. One subgroup exhibits a markedly higher diabetes prevalence than the other, directly demonstrating that diabetes risk varies across patient subgroups

7.4.4 Individual-Level Risk Prediction and Model Performance

An interpretable classification model achieves strong predictive performance on the test set, with an accuracy of 0.9602 and a ROC-AUC of 0.9604, indicating excellent overall discrimination between diabetic and non-diabetic patients. The ROC curve further demonstrates robust discriminative ability across different decision thresholds. In addition, the predicted probability distribution shows a clear separation between the two classes: most non-diabetic patients receive predicted probabilities close to 0, while diabetic patients tend to receive much higher predicted probabilities, often approaching 1. This indicates that the model can effectively distinguish low-risk individuals from high-risk individuals at the individual level. The confusion matrix further supports this finding. Out of all test samples, 18,130 non-diabetic patients are correctly classified, while 1,073 diabetic patients are correctly identified. Misclassifications are relatively limited, with 648 false negatives and 147 false positives. The concentration of observations along the diagonal indicates that the model’s predictions are stable and reliable. Overall, these evaluation results demonstrate that the model can accurately assess diabetes risk for individual patients in the test set. Despite class imbalance, the model maintains strong discriminative performance, as reflected by the high ROC-AUC value.

7.4.5 Explanation of Diabetes Risk Drivers

Odds Ratio analysis identifies HbA1c level, blood glucose level, BMI, and age as key drivers of diabetes risk. Importantly, subgroup (cluster) membership is also a significant predictor, indicating that diabetes risk is influenced not only by individual variables but also by combinations of clinical features.

7.4.6 Interpretation of Risk at the Subgroup Level

Boxplots comparing key clinical variables across subgroups show that the high-risk subgroup consistently exhibits higher levels of HbA1c, blood glucose, BMI, and age. These systematic differences explain why this subgroup has a higher diabetes prevalence and provide a clinically meaningful interpretation of subgroup risk profiles.

7.5 Role of the Model in Answering the Question

The model answers the research question at both the subgroup and individual levels. Clustering identifies patient subgroups that represent different combinations of risk factors, while the classification model assesses diabetes risk for individual patients in the test set. For each patient, the model determines whether the patient belongs to a higher-risk or lower-risk group and explains this classification using both individual clinical features and subgroup membership. This demonstrates that diabetes risk is driven by aggregated clinical characteristics rather than single variables alone.

7.6 Conclusion

In conclusion, the analysis provides clear evidence that distinct subgroups of patients with different diabetes risk profiles exist. These subgroups differ in clinical characteristics and diabetes prevalence. The model is able to identify whether an individual patient in the test set belongs to a high-risk or low-risk group and to explain this classification using interpretable clinical factors. Therefore, the findings support the research question and suggest that differentiated clinical attention may be more appropriate than a uniform approach to diabetes risk managem

8 Conclusion

This project conducts a systematic and in-depth analysis of the Diabetes Risk Factors & Diagnosis Dataset, covering the full lifecycle of data processing from exploration and cleaning to modeling and interpretation, and successfully achieves the preset research objectives. In terms of data processing, the project first completed a comprehensive exploration of the 100,000-row × 16-column dataset (spanning 2015–2022), clarifying its variable structure, distribution characteristics, and quality issues such as duplicate records, outliers, and target variable imbalance. Through targeted data cleaning—including removing 14 duplicate rows, standardizing categorical variables, and handling BMI outliers based on medical rationality—a high-quality dataset of 99,986 rows was obtained, laying a solid foundation for subsequent modeling. In the analytical modeling phase, the project effectively addressed core research questions through multiple technical approaches. For the regression task, three models (Lasso regression, random forest, and stacking ensemble) were constructed to predict patients’ blood glucose levels. The results confirmed that the random forest model outperformed linear models significantly, achieving an R² of 0.744 and an RMSE of 20.30, and identified HbA1c level and age as key influencing factors. This verified the existence of complex nonlinear relationships between health indicators and blood glucose, providing a reliable tool for blood glucose prediction. For the exploration of diabetes risk subgroups, K-means clustering combined with PCA visualization successfully identified two distinct patient subgroups with significant differences in diabetes prevalence. The subsequent interpretable logistic regression model achieved an accuracy of 0.9602 and a ROC-AUC of 0.9604 on the test set, effectively realizing individual diabetes risk assessment and confirming that HbA1c level, blood glucose level, BMI, age, and subgroup membership are core drivers of diabetes risk. The findings of this project have important practical significance. The high-precision blood glucose prediction model and diabetes risk assessment tool can assist in clinical early screening and risk stratification, while the identification of risk subgroups provides a basis for personalized diabetes management strategies—suggesting that differentiated clinical interventions may be more effective than uniform management approaches. However, the project also has certain limitations. For example, the “No Info” category in the smoking history variable may have affected the comprehensiveness of association analysis, and the model’s generalization ability may be restricted by the dataset’s regional and demographic coverage. Future research could supplement more detailed lifestyle-related variables (such as diet and exercise habits) and expand the dataset’s sample scope to further improve model robustness and application breadth. Overall, this project not only meets the technical requirements of data science programming (including data processing, multi-model construction, and result visualization) but also provides valuable insights for clinical diabetes risk management.