using CSV, DataFrames, GLM, StatsPlots, MLDataUtils  # Install these packages if you haven't already

# 1. Load the data
df = CSV.read("pima-indians-diabetes.csv", DataFrame, header=false) # Replace with your file path

# Add column names (the dataset often doesn't have them)
names!(df, [:Pregnancies, :Glucose, :BloodPressure, :SkinThickness, :Insulin, :BMI, :DiabetesPedigreeFunction, :Age, :Outcome])

# 2. Data Preprocessing (Important!)

# Handle missing values (if any).  In this dataset, zeros in some columns are likely placeholders for missing data.
# Several strategies exist.  Here's a simple imputation with the median for Glucose, BloodPressure, SkinThickness, Insulin, and BMI.  You might explore other methods.
for col in [:Glucose, :BloodPressure, :SkinThickness, :Insulin, :BMI]
    median_val = median(filter(x -> x > 0, df[!, col])) # Median of non-zero values
    df[!, col] = map(x -> x == 0 ? median_val : x, df[!, col]) # Replace 0s with median
end

# Split into features (X) and target (y)
X = Matrix(df[:, Not(:Outcome)])
y = df.Outcome

# Split into training and testing sets (essential for evaluating model performance)
X_train, X_test, y_train, y_test = splitobs(X, y, 0.8, shuffle=true) # 80% train, 20% test

# 3. Train the Logistic Regression Model

# Create the model.  `Binomial()` specifies logistic regression. `LogitLink()` is the standard link function.
model = glm(@formula(Outcome ~ 1 + Pregnancies + Glucose + BloodPressure + SkinThickness + Insulin + BMI + DiabetesPedigreeFunction + Age), df, Binomial(), LogitLink())

# Or, if you've already split the data (recommended):

model_train = glm(@formula(Outcome ~ 1 + Pregnancies + Glucose + BloodPressure + SkinThickness + Insulin + BMI + DiabetesPedigreeFunction + Age), DataFrame(hcat(X_train, y_train), names(df)), Binomial(), LogitLink())

# 4. Evaluate the Model

# Make predictions on the test set
y_pred_prob = predict(model_train, DataFrame(X_test, names(df)[1:end-1]))  # Probabilities
y_pred = ifelse.(y_pred_prob .>= 0.5, 1, 0) # Classifications (threshold at 0.5)

# Calculate performance metrics
accuracy = mean(y_pred .== y_test)
println("Accuracy: ", accuracy)

# Confusion Matrix (for more detailed analysis)
confusion_matrix = @. y_test' * (y_pred == 1) + (1 - y_test)' * (y_pred == 0)
println("Confusion Matrix:\n", confusion_matrix)

# ROC Curve and AUC (Area Under the Curve)
roc_curve = roc(y_test, y_pred_prob)
plot(roc_curve, xlabel="False Positive Rate", ylabel="True Positive Rate", title="ROC Curve")
auc_val = auc(roc_curve)
println("AUC: ", auc_val)

# 5. Interpret Model Coefficients (Optional)

println("Coefficients:\n", coef(model_train)) # See the impact of each feature

# 6. Further Improvements (Optional)

# * Feature Scaling/Normalization: Can improve model performance.
# * Regularization (L1 or L2):  Can prevent overfitting.
# * Hyperparameter Tuning: Optimize model parameters.
# * Cross-validation: A more robust way to evaluate model performance.
# * Different Imputation Strategies: Explore other ways to handle missing values.