# -*- coding: utf-8 -*-
"""Appels_et_al_2025_Manual_behavioural_coded_script

Automatically generated by Colab.

Original file is located at
    https://colab.research.google.com/drive/1EuvjfY3aITF_fNJdHW3YnMzQs4jYs4nv

# **Load libraries & import data**
"""

#Import the libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.tree import plot_tree
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score, classification_report
import xgboost as xgb
from xgboost import plot_importance
from sklearn.tree import DecisionTreeClassifier
from xgboost import XGBClassifier

#Import the dataset
from google.colab import drive
drive.mount('/content/drive')
!ls "/content/drive/My Drive"
firework_df = pd.read_excel("/content/drive/My Drive/Fireworks_Video_Data_cleaned.xlsx")
firework_df.head(5)

"""# **Model performance: RF, XGB & DT**"""

firework_df['Condition'] = firework_df['Condition'].replace({'fireworks': 0, 'control': 1})

# Split the data into features (X) and target (y)
X = firework_df[['Ear 1 (low)','Ear 2','Ear 3','Ear 4','Ear 5','Tail 1','Tail 2','Tail 3','Tail 4','Tail 5', 'Posture 1', 'Posture 2', 'Posture 3', 'Posture 4', 'Posture 5',
                 'Bark/min','Hide','Panting (%)','Sleep','Liplick & lipsmack/min','Blink/min','Yawn/min',
                 'Move (%)','Stand (%)','Active (%)','Lying (%)','Sit (%)']]
y = firework_df['Condition']

# Split the dataset into training and test sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2, random_state=42)

# Train Random Forest model
rf_model = RandomForestClassifier(n_estimators=100, random_state=42)
rf_model.fit(X_train, y_train)
rf_y_pred = rf_model.predict(X_test)
rf_accuracy = accuracy_score(y_test, rf_y_pred)

# Train XGBoost model
xgb_model = XGBClassifier(n_estimators=100, random_state=42)
xgb_model.fit(X_train, y_train)
xgb_y_pred = xgb_model.predict(X_test)
xgb_accuracy = accuracy_score(y_test, xgb_y_pred)

# Train Decision Tree model
dt_model = DecisionTreeClassifier(random_state=42)
dt_model.fit(X_train, y_train)
dt_y_pred = dt_model.predict(X_test)
dt_accuracy = accuracy_score(y_test, dt_y_pred)

# Generate classification reports
rf_report = classification_report(y_test, rf_y_pred, output_dict=True)
xgb_report = classification_report(y_test, xgb_y_pred, output_dict=True)
dt_report = classification_report(y_test, dt_y_pred, output_dict=True)


print("Classification report of Random Forest:")
print(classification_report(y_test, rf_y_pred))
print("Classification report of XGBoost model:")
print(classification_report(y_test, xgb_y_pred))
print("Classification report of Deciscion Tree model:")
print(classification_report(y_test, dt_y_pred))

"""# **Feature Importance of each model**"""

# Extract feature importances from all models
rf_importances = rf_model.feature_importances_
xgb_importances = xgb_model.feature_importances_
dt_importances = dt_model.feature_importances_

# Top 5 features of each model
rf_importances = pd.DataFrame({'Behaviour': X.columns, 'Importance (%)': rf_importances * 100})
rf_importances = rf_importances.sort_values(by='Importance (%)', ascending=False)

xgb_importances = pd.DataFrame({'Behaviour': X.columns, 'Importance (%)': xgb_importances * 100})
xgb_importances = xgb_importances.sort_values(by='Importance (%)', ascending=False)

dt_importances = pd.DataFrame({'Behaviour': X.columns, 'Importance (%)': dt_importances * 100})
dt_importances = dt_importances.sort_values(by='Importance (%)', ascending=False)

importances_top_rf = rf_importances.head(5)
importances_top_xgb = xgb_importances.head(5)
importances_top_dt = dt_importances.head(5)

print("Top 5 features of Random Forest:")
print(importances_top_rf)
print("\nTop 5 features of XGBoost:")
print(importances_top_xgb)
print("\nTop 5 features of Decision Tree:")
print(importances_top_dt)

"""Plot"""

# === Create Figure === #
fig, ax2 = plt.subplots(figsize=(10, 6))  # Only create one subplot

# === Plot 2: Feature Importances === #
# Merge the top 5 features from each model
combined_importances = pd.merge(importances_top_rf, importances_top_xgb, on='Behaviour', how='outer', suffixes=('_RF', '_XGB'))
combined_importances = pd.merge(combined_importances, importances_top_dt, on='Behaviour', how='outer')
combined_importances.columns = ['Behaviour', 'Importance (%)_RF', 'Importance (%)_XGB', 'Importance (%)_DT']
combined_importances.fillna(0, inplace=True)

# Plot the feature importances
combined_importances.set_index('Behaviour').plot(kind='barh', ax=ax2, stacked=False, width=0.7)
ax2.set_xlabel('Feature Importance (%)')
ax2.set_title('Top 5 Features Importance Across Models')
ax2.legend(loc='upper right', title="Model")

# === Final Layout Adjustments === #
plt.tight_layout()
plt.show()

"""# **SHAP Impact of RF**"""

firework_df['Condition'] = firework_df['Condition'].replace({'fireworks': 0, 'control': 1})

# Split the data into features (X) and target (y)
X = firework_df[['Ear 1 (low)','Ear 2','Ear 3','Ear 4','Ear 5','Tail 1','Tail 2','Tail 3','Tail 4','Tail 5', 'Posture 1', 'Posture 2', 'Posture 3', 'Posture 4', 'Posture 5',
                 'Bark/min','Hide','Panting (%)','Sleep','Liplick & lipsmack/min','Blink/min','Yawn/min',
                 'Move (%)','Stand (%)','Active (%)','Lying (%)','Sit (%)']]
y = firework_df['Condition']

# Split the dataset into training and test sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Train Random Forest model
rf_model = RandomForestClassifier(n_estimators=100, random_state=42)
rf_model.fit(X_train, y_train)
rf_y_pred = rf_model.predict(X_test)
rf_accuracy = accuracy_score(y_test, rf_y_pred)

# Print classification reports
print("Random Forest Classification Report:")
print(classification_report(y_test, rf_y_pred))

# === SHAP Analysis for Random Forest === #
import shap

explainer_rf = shap_values = shap.TreeExplainer(rf_model) # SHAP explainer for Random Forest
shap_values_rf = explainer_rf.shap_values(X_test)

# SHAP summary plot for Random Forest
plt.figure(figsize=(10, 6))
shap.summary_plot(shap_values_rf[:,:,0], X_test, feature_names=X.columns, show=False)

"""# **Feature selection for best predictors with K-fold CV**

## RF best: ear positions, blink/min, move (%), lying (%)
"""

from sklearn.model_selection import KFold
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, roc_auc_score

"""K = 10"""

firework_df['Condition'] = firework_df['Condition'].replace({'fireworks': 0, 'control': 1})

# Split the data into features (X) and target (y)
X = firework_df[['Ear_1_low','Ear_2','Ear_3','Ear_4','Ear_5',
                 'Blink_min','Move (%)','Lying (%)']]
y = firework_df['Condition']

# Initialize Stratified K-Fold
k = 10
skf = StratifiedKFold(n_splits=k, shuffle=True, random_state=42)

# Storage for metrics
accuracies = []
precisions = []
recalls = []
f1s = []
roc_aucs = []

# Cross-validation loop
for train_index, test_index in skf.split(X, y):
    X_train, X_test = X.iloc[train_index], X.iloc[test_index]
    y_train, y_test = y.iloc[train_index], y.iloc[test_index]

    model = RandomForestClassifier(n_estimators=100, random_state=42)
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    y_proba = model.predict_proba(X_test)[:, 1]  # For ROC AUC

# Collect metrics
    accuracies.append(accuracy_score(y_test, y_pred))
    precisions.append(precision_score(y_test, y_pred))
    recalls.append(recall_score(y_test, y_pred))
    f1s.append(f1_score(y_test, y_pred))
    roc_aucs.append(roc_auc_score(y_test, y_proba))

# Report average metrics
print(f"Average Accuracy:  {np.mean(accuracies):.3f}")
print(f"Average Precision: {np.mean(precisions):.3f}")
print(f"Average Recall:    {np.mean(recalls):.3f}")
print(f"Average F1 Score:  {np.mean(f1s):.3f}")
print(f"Average ROC AUC:   {np.mean(roc_aucs):.3f}")

"""K = 5"""

firework_df['Condition'] = firework_df['Condition'].replace({'fireworks': 0, 'control': 1})

# Split the data into features (X) and target (y)
X = firework_df[['Ear_1_low','Ear_2','Ear_3','Ear_4','Ear_5',
                 'Blink_min','Move (%)','Lying (%)']]
y = firework_df['Condition']

# Initialize Stratified K-Fold
k = 5
skf = StratifiedKFold(n_splits=k, shuffle=True, random_state=42)

# Storage for metrics
accuracies = []
precisions = []
recalls = []
f1s = []
roc_aucs = []

# Cross-validation loop
for train_index, test_index in skf.split(X, y):
    X_train, X_test = X.iloc[train_index], X.iloc[test_index]
    y_train, y_test = y.iloc[train_index], y.iloc[test_index]

    model = RandomForestClassifier(n_estimators=100, random_state=42)
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    y_proba = model.predict_proba(X_test)[:, 1]  # For ROC AUC

# Collect metrics
    accuracies.append(accuracy_score(y_test, y_pred))
    precisions.append(precision_score(y_test, y_pred))
    recalls.append(recall_score(y_test, y_pred))
    f1s.append(f1_score(y_test, y_pred))
    roc_aucs.append(roc_auc_score(y_test, y_proba))

# Report average metrics
print(f"Average Accuracy:  {np.mean(accuracies):.3f}")
print(f"Average Precision: {np.mean(precisions):.3f}")
print(f"Average Recall:    {np.mean(recalls):.3f}")
print(f"Average F1 Score:  {np.mean(f1s):.3f}")
print(f"Average ROC AUC:   {np.mean(roc_aucs):.3f}")

"""## SHAP best: ear positions, blink/min, move (%), panting (%)

K = 10
"""

firework_df['Condition'] = firework_df['Condition'].replace({'fireworks': 0, 'control': 1})

# Split the data into features (X) and target (y)
X = firework_df[['Ear_1_low','Ear_2','Ear_3','Ear_4','Ear_5',
                 'Blink_min','Move (%)','Panting (%)']]
y = firework_df['Condition']

# Initialize Stratified K-Fold
k = 10
skf = StratifiedKFold(n_splits=k, shuffle=True, random_state=42)

# Storage for metrics
accuracies = []
precisions = []
recalls = []
f1s = []
roc_aucs = []

# Cross-validation loop
for train_index, test_index in skf.split(X, y):
    X_train, X_test = X.iloc[train_index], X.iloc[test_index]
    y_train, y_test = y.iloc[train_index], y.iloc[test_index]

    model = RandomForestClassifier(n_estimators=100, random_state=42)
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    y_proba = model.predict_proba(X_test)[:, 1]  # For ROC AUC

# Collect metrics
    accuracies.append(accuracy_score(y_test, y_pred))
    precisions.append(precision_score(y_test, y_pred))
    recalls.append(recall_score(y_test, y_pred))
    f1s.append(f1_score(y_test, y_pred))
    roc_aucs.append(roc_auc_score(y_test, y_proba))

# Report average metrics
print(f"Average Accuracy:  {np.mean(accuracies):.3f}")
print(f"Average Precision: {np.mean(precisions):.3f}")
print(f"Average Recall:    {np.mean(recalls):.3f}")
print(f"Average F1 Score:  {np.mean(f1s):.3f}")
print(f"Average ROC AUC:   {np.mean(roc_aucs):.3f}")

"""K = 5"""

firework_df['Condition'] = firework_df['Condition'].replace({'fireworks': 0, 'control': 1})

# Split the data into features (X) and target (y)
X = firework_df[['Ear_1_low','Ear_2','Ear_3','Ear_4','Ear_5',
                 'Blink_min','Move (%)','Panting (%)']]
y = firework_df['Condition']

# Initialize Stratified K-Fold
k = 5
skf = StratifiedKFold(n_splits=k, shuffle=True, random_state=42)

# Storage for metrics
accuracies = []
precisions = []
recalls = []
f1s = []
roc_aucs = []

# Cross-validation loop
for train_index, test_index in skf.split(X, y):
    X_train, X_test = X.iloc[train_index], X.iloc[test_index]
    y_train, y_test = y.iloc[train_index], y.iloc[test_index]

    model = RandomForestClassifier(n_estimators=100, random_state=42)
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    y_proba = model.predict_proba(X_test)[:, 1]  # For ROC AUC

# Collect metrics
    accuracies.append(accuracy_score(y_test, y_pred))
    precisions.append(precision_score(y_test, y_pred))
    recalls.append(recall_score(y_test, y_pred))
    f1s.append(f1_score(y_test, y_pred))
    roc_aucs.append(roc_auc_score(y_test, y_proba))

# Report average metrics
print(f"Average Accuracy:  {np.mean(accuracies):.3f}")
print(f"Average Precision: {np.mean(precisions):.3f}")
print(f"Average Recall:    {np.mean(recalls):.3f}")
print(f"Average F1 Score:  {np.mean(f1s):.3f}")
print(f"Average ROC AUC:   {np.mean(roc_aucs):.3f}")

"""## SHAP best & extra: ear positions, blink/min, move (%), panting (%), active (%), sleep

K = 10
"""

firework_df['Condition'] = firework_df['Condition'].replace({'fireworks': 0, 'control': 1})

# Split the data into features (X) and target (y)
X = firework_df[['Ear_1_low','Ear_2','Ear_3','Ear_4','Ear_5',
                 'Blink_min','Move (%)','Panting (%)','Active (%)','Sleep']]
y = firework_df['Condition']

# Initialize Stratified K-Fold
k = 10
skf = StratifiedKFold(n_splits=k, shuffle=True, random_state=42)

# Storage for metrics
accuracies = []
precisions = []
recalls = []
f1s = []
roc_aucs = []

# Cross-validation loop
for train_index, test_index in skf.split(X, y):
    X_train, X_test = X.iloc[train_index], X.iloc[test_index]
    y_train, y_test = y.iloc[train_index], y.iloc[test_index]

    model = RandomForestClassifier(n_estimators=100, random_state=42)
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    y_proba = model.predict_proba(X_test)[:, 1]  # For ROC AUC

# Collect metrics
    accuracies.append(accuracy_score(y_test, y_pred))
    precisions.append(precision_score(y_test, y_pred))
    recalls.append(recall_score(y_test, y_pred))
    f1s.append(f1_score(y_test, y_pred))
    roc_aucs.append(roc_auc_score(y_test, y_proba))

# Report average metrics
print(f"Average Accuracy:  {np.mean(accuracies):.3f}")
print(f"Average Precision: {np.mean(precisions):.3f}")
print(f"Average Recall:    {np.mean(recalls):.3f}")
print(f"Average F1 Score:  {np.mean(f1s):.3f}")
print(f"Average ROC AUC:   {np.mean(roc_aucs):.3f}")

"""K = 5"""

firework_df['Condition'] = firework_df['Condition'].replace({'fireworks': 0, 'control': 1})

# Split the data into features (X) and target (y)
X = firework_df[['Ear_1_low','Ear_2','Ear_3','Ear_4','Ear_5',
                 'Blink_min','Move (%)','Panting (%)','Active (%)','Sleep']]
y = firework_df['Condition']

# Initialize Stratified K-Fold
k = 5
skf = StratifiedKFold(n_splits=k, shuffle=True, random_state=42)

# Storage for metrics
accuracies = []
precisions = []
recalls = []
f1s = []
roc_aucs = []

# Cross-validation loop
for train_index, test_index in skf.split(X, y):
    X_train, X_test = X.iloc[train_index], X.iloc[test_index]
    y_train, y_test = y.iloc[train_index], y.iloc[test_index]

    model = RandomForestClassifier(n_estimators=100, random_state=42)
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    y_proba = model.predict_proba(X_test)[:, 1]  # For ROC AUC

# Collect metrics
    accuracies.append(accuracy_score(y_test, y_pred))
    precisions.append(precision_score(y_test, y_pred))
    recalls.append(recall_score(y_test, y_pred))
    f1s.append(f1_score(y_test, y_pred))
    roc_aucs.append(roc_auc_score(y_test, y_proba))

# Report average metrics
print(f"Average Accuracy:  {np.mean(accuracies):.3f}")
print(f"Average Precision: {np.mean(precisions):.3f}")
print(f"Average Recall:    {np.mean(recalls):.3f}")
print(f"Average F1 Score:  {np.mean(f1s):.3f}")
print(f"Average ROC AUC:   {np.mean(roc_aucs):.3f}")

"""## SHAP extra + posture: ear positions, blink/min, move (%), panting (%), active (%), posture 1, posture 2, posture 3, posture 4, posture 5

K = 10
"""

firework_df['Condition'] = firework_df['Condition'].replace({'fireworks': 0, 'control': 1})

# Split the data into features (X) and target (y)
X = firework_df[['Ear_1_low','Ear_2','Ear_3','Ear_4','Ear_5',
                 'Blink_min','Move (%)','Panting (%)','Active (%)','Posture 1', 'Posture 2', 'Posture 3', 'Posture 4', 'Posture 5']]
y = firework_df['Condition']

# Initialize Stratified K-Fold
k = 10
skf = StratifiedKFold(n_splits=k, shuffle=True, random_state=42)

# Storage for metrics
accuracies = []
precisions = []
recalls = []
f1s = []
roc_aucs = []

# Cross-validation loop
for train_index, test_index in skf.split(X, y):
    X_train, X_test = X.iloc[train_index], X.iloc[test_index]
    y_train, y_test = y.iloc[train_index], y.iloc[test_index]

    model = RandomForestClassifier(n_estimators=100, random_state=42)
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    y_proba = model.predict_proba(X_test)[:, 1]  # For ROC AUC

# Collect metrics
    accuracies.append(accuracy_score(y_test, y_pred))
    precisions.append(precision_score(y_test, y_pred))
    recalls.append(recall_score(y_test, y_pred))
    f1s.append(f1_score(y_test, y_pred))
    roc_aucs.append(roc_auc_score(y_test, y_proba))

# Report average metrics
print(f"Average Accuracy:  {np.mean(accuracies):.3f}")
print(f"Average Precision: {np.mean(precisions):.3f}")
print(f"Average Recall:    {np.mean(recalls):.3f}")
print(f"Average F1 Score:  {np.mean(f1s):.3f}")
print(f"Average ROC AUC:   {np.mean(roc_aucs):.3f}")

"""K = 5"""

firework_df['Condition'] = firework_df['Condition'].replace({'fireworks': 0, 'control': 1})

# Split the data into features (X) and target (y)
X = firework_df[['Ear_1_low','Ear_2','Ear_3','Ear_4','Ear_5',
                 'Blink_min','Move (%)','Panting (%)','Active (%)','Posture 1', 'Posture 2', 'Posture 3', 'Posture 4', 'Posture 5']]
y = firework_df['Condition']

# Initialize Stratified K-Fold
k = 5
skf = StratifiedKFold(n_splits=k, shuffle=True, random_state=42)

# Storage for metrics
accuracies = []
precisions = []
recalls = []
f1s = []
roc_aucs = []

# Cross-validation loop
for train_index, test_index in skf.split(X, y):
    X_train, X_test = X.iloc[train_index], X.iloc[test_index]
    y_train, y_test = y.iloc[train_index], y.iloc[test_index]

    model = RandomForestClassifier(n_estimators=100, random_state=42)
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    y_proba = model.predict_proba(X_test)[:, 1]  # For ROC AUC

# Collect metrics
    accuracies.append(accuracy_score(y_test, y_pred))
    precisions.append(precision_score(y_test, y_pred))
    recalls.append(recall_score(y_test, y_pred))
    f1s.append(f1_score(y_test, y_pred))
    roc_aucs.append(roc_auc_score(y_test, y_proba))

# Report average metrics
print(f"Average Accuracy:  {np.mean(accuracies):.3f}")
print(f"Average Precision: {np.mean(precisions):.3f}")
print(f"Average Recall:    {np.mean(recalls):.3f}")
print(f"Average F1 Score:  {np.mean(f1s):.3f}")
print(f"Average ROC AUC:   {np.mean(roc_aucs):.3f}")

"""## RF feature top 10: ear position, blink/min, move (%), Pant (%), Lying (%), liplick_and_lipsmack, active (%), sit (%), sleep

K = 10
"""

firework_df['Condition'] = firework_df['Condition'].replace({'fireworks': 0, 'control': 1})

# Split the data into features (X) and target (y)
X = firework_df[['Ear_1_low','Ear_2','Ear_3','Ear_4','Ear_5',
                 'Blink_min','Move (%)','Panting (%)','Lying (%)','Liplick & lipsmack/min','Active (%)','Sit (%)','Sleep']]
y = firework_df['Condition']

# Initialize Stratified K-Fold
k = 10
skf = StratifiedKFold(n_splits=k, shuffle=True, random_state=42)

# Storage for metrics
accuracies = []
precisions = []
recalls = []
f1s = []
roc_aucs = []

# Cross-validation loop
for train_index, test_index in skf.split(X, y):
    X_train, X_test = X.iloc[train_index], X.iloc[test_index]
    y_train, y_test = y.iloc[train_index], y.iloc[test_index]

    model = RandomForestClassifier(n_estimators=100, random_state=42)
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    y_proba = model.predict_proba(X_test)[:, 1]  # For ROC AUC

# Collect metrics
    accuracies.append(accuracy_score(y_test, y_pred))
    precisions.append(precision_score(y_test, y_pred))
    recalls.append(recall_score(y_test, y_pred))
    f1s.append(f1_score(y_test, y_pred))
    roc_aucs.append(roc_auc_score(y_test, y_proba))

# Report average metrics
print(f"Average Accuracy:  {np.mean(accuracies):.3f}")
print(f"Average Precision: {np.mean(precisions):.3f}")
print(f"Average Recall:    {np.mean(recalls):.3f}")
print(f"Average F1 Score:  {np.mean(f1s):.3f}")
print(f"Average ROC AUC:   {np.mean(roc_aucs):.3f}")

"""K = 5"""

firework_df['Condition'] = firework_df['Condition'].replace({'fireworks': 0, 'control': 1})

# Split the data into features (X) and target (y)
X = firework_df[['Ear_1_low','Ear_2','Ear_3','Ear_4','Ear_5',
                 'Blink_min','Move (%)','Panting (%)','Lying (%)','Liplick & lipsmack/min','Active (%)','Sit (%)','Sleep']]
y = firework_df['Condition']

# Initialize Stratified K-Fold
k = 5
skf = StratifiedKFold(n_splits=k, shuffle=True, random_state=42)

# Storage for metrics
accuracies = []
precisions = []
recalls = []
f1s = []
roc_aucs = []

# Cross-validation loop
for train_index, test_index in skf.split(X, y):
    X_train, X_test = X.iloc[train_index], X.iloc[test_index]
    y_train, y_test = y.iloc[train_index], y.iloc[test_index]

    model = RandomForestClassifier(n_estimators=100, random_state=42)
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    y_proba = model.predict_proba(X_test)[:, 1]  # For ROC AUC

# Collect metrics
    accuracies.append(accuracy_score(y_test, y_pred))
    precisions.append(precision_score(y_test, y_pred))
    recalls.append(recall_score(y_test, y_pred))
    f1s.append(f1_score(y_test, y_pred))
    roc_aucs.append(roc_auc_score(y_test, y_proba))

# Report average metrics
print(f"Average Accuracy:  {np.mean(accuracies):.3f}")
print(f"Average Precision: {np.mean(precisions):.3f}")
print(f"Average Recall:    {np.mean(recalls):.3f}")
print(f"Average F1 Score:  {np.mean(f1s):.3f}")
print(f"Average ROC AUC:   {np.mean(roc_aucs):.3f}")

"""# **Validation & generalisation: LOSO & k=n_subjects/10**

## Leave-One-Subject-Out: SHAP best features
"""

from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, roc_auc_score
import numpy as np

# Ensure binary labels
firework_df['Condition'] = firework_df['Condition'].replace({'fireworks': 0, 'control': 1})

# Define features, target, and subjects
X = firework_df[['Ear 1 (low)','Ear 2','Ear 3','Ear 4','Ear 5','Blink/min','Move (%)','Panting (%)']]
y = firework_df['Condition']
subjects = firework_df['Name']

# Define specific test subjects
test_subjects = ['Brunner Lucy', 'Zehnder Yaro', 'Gefing Lilly']

# Metrics to store results
accuracies, precisions, recalls, f1s, roc_aucs = [], [], [], [], []

for test_subject in test_subjects:
    test_mask = (subjects == test_subject)

    # Split the data
    X_train, X_test = X[~test_mask], X[test_mask]
    y_train, y_test = y[~test_mask], y[test_mask]

    # Train the Random Forest model
    model = RandomForestClassifier(n_estimators=100, random_state=42)
    model.fit(X_train, y_train)

    # Predict and evaluate
    y_pred = model.predict(X_test)
    y_proba = model.predict_proba(X_test)[:, 1]

    # Store metrics
    accuracies.append(accuracy_score(y_test, y_pred))
    precisions.append(precision_score(y_test, y_pred))
    recalls.append(recall_score(y_test, y_pred))
    f1s.append(f1_score(y_test, y_pred))
    roc_aucs.append(roc_auc_score(y_test, y_proba))

    print(f"Results for {test_subject}:")
    print(f"Accuracy:  {accuracies[-1]:.3f}")
    print(f"Precision: {precisions[-1]:.3f}")
    print(f"Recall:    {recalls[-1]:.3f}")
    print(f"F1 Score:  {f1s[-1]:.3f}")
    print(f"ROC AUC:   {roc_aucs[-1]:.3f}")
    print("-" * 30)

# Print average performance
print(f"\nAverage Results Across {len(test_subjects)} Subjects:")
print(f"Accuracy:  {np.mean(accuracies):.3f}")
print(f"Precision: {np.mean(precisions):.3f}")
print(f"Recall:    {np.mean(recalls):.3f}")
print(f"F1 Score:  {np.mean(f1s):.3f}")
print(f"ROC AUC:   {np.mean(roc_aucs):.3f}")

"""## Leave-One-Subject-Out: All behaviours"""

from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, roc_auc_score
import numpy as np

# Ensure binary labels
firework_df['Condition'] = firework_df['Condition'].replace({'fireworks': 0, 'control': 1})

# Define features, target, and subjects
X = firework_df[['Ear 1 (low)','Ear 2','Ear 3','Ear 4','Ear 5','Tail 1','Tail 2','Tail 3','Tail 4','Tail 5', 'Posture 1', 'Posture 2', 'Posture 3', 'Posture 4', 'Posture 5',
                 'Bark/min','Hide','Panting (%)','Sleep','Liplick & lipsmack/min','Blink/min','Yawn/min',
                 'Move (%)','Stand (%)','Active (%)','Lying (%)','Sit (%)']]
y = firework_df['Condition']
subjects = firework_df['Name']

# Define specific test subjects
test_subjects = ['Brunner Lucy', 'Zehnder Yaro', 'Gefing Lilly']

# Metrics to store results
accuracies, precisions, recalls, f1s, roc_aucs = [], [], [], [], []

for test_subject in test_subjects:
    test_mask = (subjects == test_subject)

    # Split the data
    X_train, X_test = X[~test_mask], X[test_mask]
    y_train, y_test = y[~test_mask], y[test_mask]

    # Train the Random Forest model
    model = RandomForestClassifier(n_estimators=100, random_state=42)
    model.fit(X_train, y_train)

    # Predict and evaluate
    y_pred = model.predict(X_test)
    y_proba = model.predict_proba(X_test)[:, 1]

    # Store metrics
    accuracies.append(accuracy_score(y_test, y_pred))
    precisions.append(precision_score(y_test, y_pred))
    recalls.append(recall_score(y_test, y_pred))
    f1s.append(f1_score(y_test, y_pred))
    roc_aucs.append(roc_auc_score(y_test, y_proba))

    print(f"Results for {test_subject}:")
    print(f"Accuracy:  {accuracies[-1]:.3f}")
    print(f"Precision: {precisions[-1]:.3f}")
    print(f"Recall:    {recalls[-1]:.3f}")
    print(f"F1 Score:  {f1s[-1]:.3f}")
    print(f"ROC AUC:   {roc_aucs[-1]:.3f}")
    print("-" * 30)

# Print average performance
print(f"\nAverage Results Across {len(test_subjects)} Subjects:")
print(f"Accuracy:  {np.mean(accuracies):.3f}")
print(f"Precision: {np.mean(precisions):.3f}")
print(f"Recall:    {np.mean(recalls):.3f}")
print(f"F1 Score:  {np.mean(f1s):.3f}")
print(f"ROC AUC:   {np.mean(roc_aucs):.3f}")

"""## K = 4-fold CV: SHAP best features"""

from sklearn.model_selection import KFold

# Define features, labels, and subject IDs
feature_columns = ['Ear 1 (low)','Ear 2','Ear 3','Ear 4','Ear 5','Blink/min','Move (%)','Panting (%)']
X = firework_df[feature_columns]
y = firework_df['Condition']
subjects = firework_df['Name']

# Define fixed test subjects
fixed_test_subjects = ['Brunner Lucy', 'Zehnder Yaro', 'Gefing Lilly']

# Extract test set (fixed)
test_mask = subjects.isin(fixed_test_subjects)
X_test, y_test = X[test_mask], y[test_mask]

# Remaining subjects for training (will be cross-validated)
train_subjects = np.unique(subjects[~test_mask])
k = 4  # Choose based on dataset size

kf = KFold(n_splits=k, shuffle=True, random_state=42)

# Metrics storage
accuracies, precisions, recalls, f1s, roc_aucs = [], [], [], [], []

# Loop over training folds
for train_idx, _ in kf.split(train_subjects):
    current_train_subjects = train_subjects[train_idx]
    train_mask = subjects.isin(current_train_subjects)

    X_train, y_train = X[train_mask], y[train_mask]

    # Train model
    model = RandomForestClassifier(n_estimators=100, random_state=42)
    model.fit(X_train, y_train)

    # Evaluate only on fixed test set
    y_pred = model.predict(X_test)
    y_proba = model.predict_proba(X_test)[:, 1]

    # Store metrics
    accuracies.append(accuracy_score(y_test, y_pred))
    precisions.append(precision_score(y_test, y_pred, zero_division=0))
    recalls.append(recall_score(y_test, y_pred, zero_division=0))
    f1s.append(f1_score(y_test, y_pred, zero_division=0))

    if len(np.unique(y_test)) == 2:
        roc_aucs.append(roc_auc_score(y_test, y_proba))

    print(f"Fold done (trained on {len(current_train_subjects)} subjects, tested on fixed 3).")

# Final averaged results
print("\nFinal Results (tested on Lilly, Lucy, Yaro):")
print(f"Average Accuracy:  {np.mean(accuracies):.3f}")
print(f"Average Precision: {np.mean(precisions):.3f}")
print(f"Average Recall:    {np.mean(recalls):.3f}")
print(f"Average F1 Score:  {np.mean(f1s):.3f}")
print(f"Average ROC AUC:   {np.nanmean(roc_aucs):.3f}")

"""## K = 4-fold CV: All behaviours"""

# Define features, labels, and subject IDs
feature_columns = ['Ear 1 (low)','Ear 2','Ear 3','Ear 4','Ear 5','Tail 1','Tail 2','Tail 3','Tail 4','Tail 5', 'Posture 1', 'Posture 2', 'Posture 3', 'Posture 4', 'Posture 5',
                 'Bark/min','Hide','Panting (%)','Sleep','Liplick & lipsmack/min','Blink/min','Yawn/min',
                 'Move (%)','Stand (%)','Active (%)','Lying (%)','Sit (%)']
X = firework_df[feature_columns]
y = firework_df['Condition']
subjects = firework_df['Name']

# Define fixed test subjects
fixed_test_subjects = ['Brunner Lucy', 'Zehnder Yaro', 'Gefing Lilly']

# Extract test set (fixed)
test_mask = subjects.isin(fixed_test_subjects)
X_test, y_test = X[test_mask], y[test_mask]

# Remaining subjects for training (will be cross-validated)
train_subjects = np.unique(subjects[~test_mask])
k = 4  # Choose based on dataset size

kf = KFold(n_splits=k, shuffle=True, random_state=42)

# Metrics storage
accuracies, precisions, recalls, f1s, roc_aucs = [], [], [], [], []

# Loop over training folds
for train_idx, _ in kf.split(train_subjects):
    current_train_subjects = train_subjects[train_idx]
    train_mask = subjects.isin(current_train_subjects)

    X_train, y_train = X[train_mask], y[train_mask]

    # Train model
    model = RandomForestClassifier(n_estimators=100, random_state=42)
    model.fit(X_train, y_train)

    # Evaluate only on fixed test set
    y_pred = model.predict(X_test)
    y_proba = model.predict_proba(X_test)[:, 1]

    # Store metrics
    accuracies.append(accuracy_score(y_test, y_pred))
    precisions.append(precision_score(y_test, y_pred, zero_division=0))
    recalls.append(recall_score(y_test, y_pred, zero_division=0))
    f1s.append(f1_score(y_test, y_pred, zero_division=0))

    if len(np.unique(y_test)) == 2:
        roc_aucs.append(roc_auc_score(y_test, y_proba))

    print(f"Fold done (trained on {len(current_train_subjects)} subjects, tested on fixed 3).")

# Final averaged results
print("\nFinal Results (tested on Lilly, Lucy, Yaro):")
print(f"Average Accuracy:  {np.mean(accuracies):.3f}")
print(f"Average Precision: {np.mean(precisions):.3f}")
print(f"Average Recall:    {np.mean(recalls):.3f}")
print(f"Average F1 Score:  {np.mean(f1s):.3f}")
print(f"Average ROC AUC:   {np.nanmean(roc_aucs):.3f}")