Ethics Section¶

df_wlabs = pd.read_csv("data/combined_features.csv")
df = df_wlabs.drop(columns=['ID'])

X = df.drop(['height_above'], axis=1).select_dtypes(include=[np.number]).values
y = df['height_above'].values
indices = np.arange(len(X))
# np.random.shuffle(indices)
X = X[indices]
y = y[indices]

print(f"\nFinal feature matrix shape: {X.shape}")
print(f"Target variable shape: {y.shape}")

# split datasets in train, validation and test sets
X_train, y_train, X_val, y_val, X_test, y_test  = split_dataset(X, y, test_size=0.3, random_state=RANDOM_STATE)

#z-normalize datasets across splits 
X_train_norm, X_val_norm, X_test_norm = z_normalize(X_train, X_val, X_test)

# Robust Hyperparameter Tuning with Cross-Validation
rf = RandomForestRegressor(random_state=RANDOM_STATE)

# Parameter grid for regression
rf_params = {
    'n_estimators': [50, 100, 200, 300],
    'max_depth': [None, 10, 20, 30],
    'min_samples_split': [2, 5, 10],
    'min_samples_leaf': [1, 2, 4],
    'max_features': ['sqrt', 'log2', None],
    'bootstrap': [True, False]
}

# Cross-validation for regression
rf_cv = RandomizedSearchCV(
    rf, 
    rf_params, 
    n_iter=100,  # More iterations for better search
    cv=5,       # 10-fold CV for more robust evaluation
    scoring='neg_mean_squared_error',  # MSE for regression
    random_state=RANDOM_STATE,
    n_jobs=-1,   # Use all cores
    verbose=1
)

print("Starting hyperparameter optimization...")
rf_optimized = rf_cv.fit(X_train_norm, y_train)

print('Best parameters:', rf_optimized.best_params_)
print('Best cross-validation score (neg_MSE):', rf_optimized.best_score_)

# Model Evaluation on All Sets
best_rf = rf_optimized.best_estimator_
# Predictions on all sets
y_train_pred = best_rf.predict(X_train_norm)
y_val_pred = best_rf.predict(X_val_norm)
y_test_pred = best_rf.predict(X_test_norm)
print("="*60)
print("MODEL PERFORMANCE METRICS")
print("="*60)

# Calculate metrics for all sets
train_metrics = calculate_regression_metrics(y_train, y_train_pred, "Training")
val_metrics = calculate_regression_metrics(y_val, y_val_pred, "Validation")
test_metrics = calculate_regression_metrics(y_test, y_test_pred, "Test")

fig, axes = plt.subplots(1, 3, figsize=(18, 5))
plot_actual_vs_predicted(y_train, y_train_pred, axes[0], 'Training', train_metrics)
plot_actual_vs_predicted(y_val, y_val_pred, axes[1], 'Validation', val_metrics)
plot_actual_vs_predicted(y_test, y_test_pred, axes[2], 'Test', test_metrics)

plt.tight_layout()
plt.show()

fig, axes = plt.subplots(1, 3, figsize=(18, 5))
plot_residuals(y_train, y_train_pred, axes[0], 'Training')
plot_residuals(y_val, y_val_pred, axes[1], 'Validation')
plot_residuals(y_test, y_test_pred, axes[2], 'Test')

plt.tight_layout()
plt.show()

feature_importance, importance_df = feature_importance_analysis(best_rf, df)
importance_df

modal_features_df = pd.read_csv("data/rf_feature_importance_wcat.csv")

# group by 'category' and sum up 'importance'
importance_by_category = (
    modal_features_df
    .groupby("category")["importance"]
    .sum()
    .reset_index()
    .sort_values(by="importance", ascending=False)
)
importance_by_category

plot_importances(importance_by_category)

plot_cumulative_importances(importance_df)

cumulative_importance = np.cumsum(importance_df['importance'].values)
# Find number of features needed for 80% and 90% of importance
idx_80 = np.argmax(cumulative_importance >= 0.8) + 1
idx_90 = np.argmax(cumulative_importance >= 0.9) + 1

print(f"\nFeature Selection Insights:")
print(f"Number of features needed for 80% of importance: {idx_80}")
print(f"Number of features needed for 90% of importance: {idx_90}")
print(f"Total number of features: {len(feature_importance)}")

# Save feature importance to CSV
importance_path = r'data/rf_feature_importance.csv'
importance_df.to_csv(importance_path, index=False)
print(f"\nFeature importance saved to: {importance_path}")

# Final Model Summary
print("="*80)
print("FINAL MODEL SUMMARY")
print("="*80)
print(f"Best Random Forest Parameters: {rf_optimized.best_params_}")
print(f"Best Cross-Validation Score (neg_MSE): {rf_optimized.best_score_:.4f}")
print(f"\nFinal Test Set Performance:")
print(f"  R² Score: {test_metrics['r2']:.4f}")
print(f"  RMSE: {test_metrics['rmse']:.4f}")
print(f"  MAE: {test_metrics['mae']:.4f}")
print(f"  MSE: {test_metrics['mse']:.4f}")

# Check for overfitting
print(f"\nOverfitting Check:")
print(f"  Training R²: {train_metrics['r2']:.4f}")
print(f"  Validation R²: {val_metrics['r2']:.4f}")
print(f"  Test R²: {test_metrics['r2']:.4f}")
print(f"  Train-Val R² Gap: {train_metrics['r2'] - val_metrics['r2']:.4f}")
print(f"  Train-Test R² Gap: {train_metrics['r2'] - test_metrics['r2']:.4f}")

if train_metrics['r2'] - test_metrics['r2'] > 0.2:
    print("  WARNING: Potential overfitting detected (>20% R² gap between train and test)")
else:
    print("  Model shows good generalization")

print("\n" + "="*80)

# Extract residuals with IDs for all data splits
print("="*60)
print("EXTRACTING RESIDUALS WITH IDs")
print("="*60)

# We need to reconstruct the full dataset with predictions to get residuals with IDs
# First, let's get predictions for the entire dataset
full_X = df.drop(['height_above'], axis=1).select_dtypes(include=[np.number]).values
full_y = df_wlabs['height_above'].values
original_ids = df_wlabs['ID'].values

# Standardize the full dataset using training statistics (same as used for model training)
X_mean = np.nanmean(X_train, axis=0)
X_std = np.nanstd(X_train, axis=0)
full_X_norm = (full_X - X_mean) / X_std

# Get predictions for all samples
full_predictions = best_rf.predict(full_X_norm)

# Calculate residuals (actual - predicted)
residuals = full_y - full_predictions

# Create residuals dataframe with IDs
residuals_df = pd.DataFrame({
    'ID': original_ids,
    'Actual_height_above': full_y,
    'Predicted_height_above': full_predictions,
    'Residual': residuals,
    'Absolute_Residual': np.abs(residuals)
})

# Sort by ID for easier analysis
residuals_df = residuals_df.sort_values('ID')

# Display basic statistics
print(f"Total samples: {len(residuals_df)}")
print(f"Mean residual: {residuals_df['Residual'].mean():.4f}")
print(f"Std residual: {residuals_df['Residual'].std():.4f}")
print(f"Mean absolute residual: {residuals_df['Absolute_Residual'].mean():.4f}")

# Show first few rows
print(f"\nFirst 10 rows of residuals:")
print(residuals_df.head(10))

# Save to CSV
output_path = 'data/rf_residuals.csv'
residuals_df.to_csv(output_path, index=False)
print(f"\nResiduals saved to: {output_path}")

# Additional analysis: identify largest residuals
print(f"\nTop 10 largest positive residuals (model under-predicted):")
print(residuals_df.nlargest(10, 'Residual')[['ID', 'Actual_height_above', 'Predicted_height_above', 'Residual']])

print(f"\nTop 10 largest negative residuals (model over-predicted):")
print(residuals_df.nsmallest(10, 'Residual')[['ID', 'Actual_height_above', 'Predicted_height_above', 'Residual']])

print(f"\nTop 10 largest absolute residuals:")
print(residuals_df.nlargest(10, 'Absolute_Residual')[['ID', 'Actual_height_above', 'Predicted_height_above', 'Residual']])