Adding new notebook

python/Testing real data with logit_2_Zahra.ipynb

@@ -0,0 +1,153 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8f50efc2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import seaborn as sns\n",
+    "from scipy import io\n",
+    "from sklearn.linear_model import LogisticRegression\n",
+    "from sklearn.preprocessing import StandardScaler  # , PolynomialFeatures\n",
+    "from sklearn.model_selection import StratifiedKFold, cross_val_score\n",
+    "from sklearn.metrics import roc_curve, auc\n",
+    "\n",
+    "# ---------------------------\n",
+    "# Load the data from .mat files\n",
+    "# ---------------------------\n",
+    "data_path = \"../data/\"\n",
+    "suv_file = data_path + \"suv_percentilesSLOthenUWM.mat\"\n",
+    "flags_file = data_path + \"flags_combined.mat\"\n",
+    "normal_range_file = data_path + \"normal_range.mat\"\n",
+    "\n",
+    "# Load .mat files\n",
+    "suv_dict = io.loadmat(suv_file)\n",
+    "flags_dict = io.loadmat(flags_file)\n",
+    "normal_range_dict = io.loadmat(normal_range_file)\n",
+    "\n",
+    "# Extract relevant data\n",
+    "# Adjust the slicing as needed based on your data dimensions.\n",
+    "suv = suv_dict['lung_SUVperc_COMBINED'][0:58, :, :]\n",
+    "flags = flags_dict['flags'][0:58, 3]\n",
+    "\n",
+    "# ---------------------------\n",
+    "# Feature Engineering\n",
+    "# ---------------------------\n",
+    "# Choose the percentile of interest. For lung and bowel, p = 94 (0-indexed)\n",
+    "p = 94\n",
+    "\n",
+    "# Compute a set of summary statistics for the chosen percentile.\n",
+    "feature_max = np.nanmax(suv[:, :, p], axis=1)\n",
+    "feature_mean = np.nanmean(suv[:, :, p], axis=1)\n",
+    "feature_std = np.nanstd(suv[:, :, p], axis=1)\n",
+    "\n",
+    "# Stack the features so that each row corresponds to one sample and\n",
+    "# each column is a different feature.\n",
+    "X = np.column_stack((feature_max, feature_mean, feature_std))\n",
+    "\n",
+    "# Handle missing values by selecting only rows with complete data.\n",
+    "valid_indices = ~np.isnan(X).any(axis=1)\n",
+    "X = X[valid_indices]\n",
+    "flags = flags[valid_indices]\n",
+    "\n",
+    "# ---------------------------\n",
+    "# Preprocessing: Feature Scaling\n",
+    "# ---------------------------\n",
+    "scaler = StandardScaler()\n",
+    "X_scaled = scaler.fit_transform(X)\n",
+    "\n",
+    "# --------------\n",
+    "# (Optional) Experiment with Polynomial Features for non-linearity\n",
+    "# --------------\n",
+    "# Uncomment below to include non-linear interactions\n",
+    "# from sklearn.preprocessing import PolynomialFeatures\n",
+    "# poly = PolynomialFeatures(degree=2, include_bias=False)\n",
+    "# X_scaled = poly.fit_transform(X_scaled)\n",
+    "\n",
+    "# ---------------------------\n",
+    "# Model: Regularized Logistic Regression\n",
+    "# ---------------------------\n",
+    "# We use L2 regularization; you can change penalty to 'l1' to test Lasso.\n",
+    "clf = LogisticRegression(penalty='l2', solver='liblinear', random_state=42)\n",
+    "clf.fit(X_scaled, flags)\n",
+    "\n",
+    "# Display model coefficients (note: these are in the scaled feature space)\n",
+    "print(\"Intercept:\", clf.intercept_)\n",
+    "print(\"Coefficients:\", clf.coef_)\n",
+    "\n",
+    "# ---------------------------\n",
+    "# Model Evaluation: Cross-Validation Accuracy\n",
+    "# ---------------------------\n",
+    "cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)\n",
+    "accuracy_scores = cross_val_score(clf, X_scaled, flags, cv=cv, scoring='accuracy')\n",
+    "print(\"Cross-validated accuracy scores:\", accuracy_scores)\n",
+    "print(\"Mean cross-validated accuracy:\", np.mean(accuracy_scores))\n",
+    "\n",
+    "# ---------------------------\n",
+    "# Additional Evaluation: ROC Curve and AUC\n",
+    "# ---------------------------\n",
+    "# Get predicted probabilities for the positive class\n",
+    "y_scores = clf.predict_proba(X_scaled)[:, 1]\n",
+    "\n",
+    "# Compute ROC curve\n",
+    "fpr, tpr, thresholds = roc_curve(flags, y_scores)\n",
+    "roc_auc = auc(fpr, tpr)\n",
+    "print(\"AUC:\", roc_auc)\n",
+    "\n",
+    "# Plot ROC curve\n",
+    "plt.figure(figsize=(8, 6))\n",
+    "plt.plot(fpr, tpr, color='darkorange', lw=2, label='ROC curve (AUC = %0.2f)' % roc_auc)\n",
+    "plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--')\n",
+    "plt.xlim([0.0, 1.0])\n",
+    "plt.ylim([0.0, 1.05])\n",
+    "plt.xlabel('False Positive Rate')\n",
+    "plt.ylabel('True Positive Rate')\n",
+    "plt.title('Receiver Operating Characteristic (ROC)')\n",
+    "plt.legend(loc=\"lower right\")\n",
+    "plt.show()\n",
+    "\n",
+    "# ---------------------------\n",
+    "# Plotting Logistic Curve for Primary Feature (maxSUV)\n",
+    "# ---------------------------\n",
+    "# For visualization, we pick the feature corresponding to maxSUV (first column).\n",
+    "feature_index = 0\n",
+    "x_range = np.linspace(X[:, feature_index].min(), X[:, feature_index].max(), 100)\n",
+    "\n",
+    "# For the additional features (mean and std), we set them to their average values.\n",
+    "avg_mean = np.mean(X[:, 1])\n",
+    "avg_std = np.mean(X[:, 2])\n",
+    "X_plot = np.column_stack((x_range, np.repeat(avg_mean, 100), np.repeat(avg_std, 100)))\n",
+    "\n",
+    "# Scale the plot data using the same scaler.\n",
+    "X_plot_scaled = scaler.transform(X_plot)\n",
+    "\n",
+    "# Compute predicted probabilities along the x_range.\n",
+    "pred_probs = clf.predict_proba(X_plot_scaled)[:, 1]\n",
+    "\n",
+    "plt.figure(figsize=(10, 6))\n",
+    "# Scatter plot of observed data along maxSUV (using flags to indicate classes).\n",
+    "plt.scatter(X[flags == 0, feature_index], flags[flags == 0], label=\"NC\", color='blue', alpha=0.6)\n",
+    "plt.scatter(X[flags == 1, feature_index], flags[flags == 1], label=\"AE\", color='red', alpha=0.6)\n",
+    "# Plot the logistic regression prediction curve.\n",
+    "plt.plot(x_range, pred_probs, color='green', lw=2, label=\"Regularized Logistic Regression\")\n",
+    "plt.title(\"Prediction Probability of AE Depending on maxSUV\")\n",
+    "plt.xlabel(\"maxSUV\")\n",
+    "plt.ylabel(\"P(AE)\")\n",
+    "plt.legend(loc='best')\n",
+    "plt.grid(True)\n",
+    "plt.show()\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "language_info": {
+   "name": "python"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}