RandomForestModel.py

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import joblib
import os
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score, classification_report, confusion_matrix, roc_auc_score
from sklearn.model_selection import train_test_split
from imblearn.over_sampling import SMOTE, RandomOverSampler

# Define directories
data_path = r'C:\Users\jsntg\OneDrive\Desktop\3YP\code\data\good.csv'
graph_path = r'C:\Users\jsntg\OneDrive\Desktop\3YP\code\graph'
result_path = r'C:\Users\jsntg\OneDrive\Desktop\3YP\code\result'

class SepsisRandomForest:
    def __init__(self):
        self.data_path = data_path
        self.graph_path = graph_path
        self.result_path = result_path
        self.model = None
        self.feature_names = None
        self.X_train = None
        self.X_test = None
        self.y_train = None
        self.y_test = None

    def load_data(self):
        print(f"Loading data from: {self.data_path}")
        if not os.path.exists(self.data_path):
            raise FileNotFoundError(f"Error! Data file not found: {self.data_path}")
        df = pd.read_csv(self.data_path)
        print(f"Data shape: {df.shape}")
        print(f"Columns: {df.columns.tolist()}")
        return df

    def preprocess_data(self, df, window_size=2, test_size=0.2, random_state=42):
        print(f"Preprocessing data (window_size={window_size})...")

        # Stratified split by patient (using max sepsis label per patient)
        patients = df['patient_id'].unique()
        patient_labels = df.groupby('patient_id')['sepsis'].max()  # 0 or 1 per patient
        train_patients, test_patients = train_test_split(
            patients, test_size=test_size, random_state=random_state, stratify=patient_labels
        )
        print(f"Training patients: {len(train_patients)}, Testing patients: {len(test_patients)}")

        train_df = df[df['patient_id'].isin(train_patients)]
        test_df = df[df['patient_id'].isin(test_patients)]

        # Create features for each set
        X_train, y_train, _, feature_names = self._create_sliding_window_features(train_df, window_size)
        X_test, y_test, _, _ = self._create_sliding_window_features(test_df, window_size)

        self.X_train = X_train
        self.X_test = X_test
        self.y_train = y_train
        self.y_test = y_test
        self.feature_names = feature_names

        print(f"Training set shape: {X_train.shape}, Testing set shape: {X_test.shape}")
        print(f"Number of features: {len(feature_names)}")
        return X_train, X_test, y_train, y_test, feature_names

    def _create_sliding_window_features(self, df, window_size=2):
        features_list = []
        labels_list = []
        patient_ids_list = []

        feature_columns = [
            'respiratory_rate', 'oxygen_saturation',
            'temperature', 'systolic_bp', 'diastolic_bp', 'heart_rate',
            'consciousness'
        ]
        # Remove any missing columns
        missing = [c for c in feature_columns if c not in df.columns]
        if missing:
            print(f"Warning: missing columns {missing}, removing them")
            feature_columns = [c for c in feature_columns if c not in missing]

        time_col = 'timestamp'
        df = df.copy()
        df[time_col] = pd.to_datetime(df[time_col])

        for pid in df['patient_id'].unique():
            patient_data = df[df['patient_id'] == pid].sort_values(time_col)
            values = patient_data[feature_columns].values
            labels = patient_data['sepsis'].values
            times = patient_data[time_col].values

            n = len(values)
            for i in range(1, n):          # start from second record
                current_time = times[i]
                current_vals = values[i]
                current_label = labels[i]

                # Find record with timestamp closest to 0.5h before current_time
                best_j = None
                min_diff = float('inf')
                for j in range(i):
                    time_diff = (current_time - times[j]).astype('timedelta64[s]').astype(float) / 3600.0
                    diff_abs = abs(time_diff - 0.5)
                    if diff_abs < min_diff:
                        min_diff = diff_abs
                        best_j = j

                if best_j is not None:
                    prev_vals = values[best_j]
                    changes = current_vals - prev_vals
                    window_feat = np.concatenate([current_vals, changes])
                    features_list.append(window_feat)
                    labels_list.append(current_label)
                    patient_ids_list.append(pid)

        # Build feature names: current values + changes
        feature_names = []
        for feat in feature_columns:
            feature_names.append(f"{feat}_t")
        for feat in feature_columns:
            feature_names.append(f"{feat}_change")

        return np.array(features_list), np.array(labels_list), np.array(patient_ids_list), feature_names

    def train_model(self, n_estimators=200, max_depth=15, min_samples_split=2,
                    min_samples_leaf=1, class_weight='balanced', use_smote=True):
        """
        Train a Random Forest classifier, optionally applying SMOTE to balance classes.
        """
        print("\nTraining Random Forest model...")
        X_train, y_train = self.X_train, self.y_train

        # Handle class imbalance with SMOTE (or fallback to random oversampling)
        if use_smote:
            unique, counts = np.unique(y_train, return_counts=True)
            if len(unique) == 2:
                min_class_count = counts[unique == 1][0]
                if min_class_count < 6:
                    print(f"Warning: minority class has only {min_class_count} samples. Using RandomOverSampler instead of SMOTE.")
                    oversampler = RandomOverSampler(random_state=42)
                else:
                    oversampler = SMOTE(random_state=42, k_neighbors=min(5, min_class_count-1))
            else:
                oversampler = RandomOverSampler(random_state=42)

            X_train, y_train = oversampler.fit_resample(X_train, y_train)
            print(f"After oversampling – positive: {np.sum(y_train)}, negative: {len(y_train)-np.sum(y_train)}")

        self.model = RandomForestClassifier(
            n_estimators=n_estimators,
            max_depth=max_depth,
            min_samples_split=min_samples_split,
            min_samples_leaf=min_samples_leaf,
            class_weight=class_weight,
            random_state=42,
            n_jobs=-1
        )
        self.model.fit(X_train, y_train)
        print(f"Number of trees: {len(self.model.estimators_)}")
        return self.model

    def evaluate_model(self):
        if self.model is None:
            raise ValueError("Model not trained yet.")

        y_train_pred = self.model.predict(self.X_train)
        y_test_pred = self.model.predict(self.X_test)

        train_acc = accuracy_score(self.y_train, y_train_pred)
        test_acc = accuracy_score(self.y_test, y_test_pred)
        print(f"Training accuracy: {train_acc:.4f} ({train_acc*100:.2f}%)")
        print(f"Testing accuracy: {test_acc:.4f} ({test_acc*100:.2f}%)")

        print("\nTraining classification report:")
        print(classification_report(self.y_train, y_train_pred))
        print("\nTesting classification report:")
        print(classification_report(self.y_test, y_test_pred))

        # ROC‑AUC on test set
        y_prob = self.model.predict_proba(self.X_test)[:, 1]
        auc = roc_auc_score(self.y_test, y_prob)
        print(f"Test AUC: {auc:.4f}")

        results = {
            'train_accuracy': train_acc,
            'test_accuracy': test_acc,
            'test_auc': auc,
            'train_report': classification_report(self.y_train, y_train_pred, output_dict=True),
            'test_report': classification_report(self.y_test, y_test_pred, output_dict=True),
            'confusion_matrix_train': confusion_matrix(self.y_train, y_train_pred),
            'confusion_matrix_test': confusion_matrix(self.y_test, y_test_pred)
        }
        return results

    def plot_feature_importance(self, top_n=20):
        if self.model is None:
            raise ValueError("Model not trained.")

        importances = self.model.feature_importances_
        indices = np.argsort(importances)[::-1]

        if top_n < len(importances):
            display_idx = indices[:top_n]
        else:
            display_idx = indices

        plt.figure(figsize=(12, 8))
        plt.bar(range(len(display_idx)), importances[display_idx])
        plt.xticks(range(len(display_idx)),
                   np.array(self.feature_names)[display_idx],
                   rotation=45, ha='right')
        plt.xlabel('Feature')
        plt.ylabel('Importance')
        plt.title(f'Random Forest Feature Importance (Top {len(display_idx)})')
        plt.tight_layout()

        save_path = os.path.join(self.graph_path, 'feature_importance_RandomForestModeln.png')
        plt.savefig(save_path, dpi=300, bbox_inches='tight')
        plt.show()
        print(f"Feature importance plot saved to {save_path}")

        print(f"\nTop {len(display_idx)} important features:")
        for i, idx in enumerate(display_idx):
            print(f"{i+1:2d}. {self.feature_names[idx]:30s} : {importances[idx]:.4f}")

    def plot_confusion_matrices(self):
        if self.model is None:
            raise ValueError("Model not trained.")

        y_train_pred = self.model.predict(self.X_train)
        y_test_pred = self.model.predict(self.X_test)
        cm_train = confusion_matrix(self.y_train, y_train_pred)
        cm_test = confusion_matrix(self.y_test, y_test_pred)

        fig, axes = plt.subplots(1, 2, figsize=(14, 6))
        sns.heatmap(cm_train, annot=True, fmt='d', cmap='Blues', ax=axes[0])
        axes[0].set_title('Training Confusion Matrix')
        axes[0].set_xlabel('Predicted')
        axes[0].set_ylabel('True')

        sns.heatmap(cm_test, annot=True, fmt='d', cmap='Blues', ax=axes[1])
        axes[1].set_title('Testing Confusion Matrix')
        axes[1].set_xlabel('Predicted')
        axes[1].set_ylabel('True')

        plt.tight_layout()
        save_path = os.path.join(self.graph_path, 'confusion_matrices_RandomForestModeln.png')
        plt.savefig(save_path, dpi=300, bbox_inches='tight')
        plt.show()
        print(f"Confusion matrices saved to {save_path}")

    def save_model(self, model_path=None):
        """Save the trained model and feature names."""
        if self.model is None:
            raise ValueError("No model to save.")

        if model_path is None:
            model_path = os.path.join(self.result_path, 'RandomForestModeln.pkl')
        joblib.dump(self.model, model_path)
        print(f"Model saved to {model_path}")

        # Save feature names for later reference
        feat_path = os.path.join(self.result_path, 'feature_names_RandomForestModeln.txt')
        with open(feat_path, 'w') as f:
            for name in self.feature_names:
                f.write(name + '\n')
        print(f"Feature names saved to {feat_path}")

def main():
    model = SepsisRandomForest()
    try:
        df = model.load_data()
        model.preprocess_data(df, window_size=2)
        model.train_model(n_estimators=100, max_depth=10, use_smote=False)
        model.evaluate_model()
        model.plot_feature_importance(top_n=15)
        model.plot_confusion_matrices()
        model.save_model()
    except Exception as e:
        print(f"Error: {e}")
        import traceback
        traceback.print_exc()

if __name__ == "__main__":
    main()