Evaluation of equations to test data

[gm89uk]

I found have found it v.time consuming to find the sweet spot of complexity that has the best fit for test data with the higher complexities on the training data overfitting and performing poorly on the test data.

In addition I also needed to apply the equations to different subgroups of the data and see how each equation applied to those.

I have put together some python code to make my life easier, and someone else might find it useful.
It assumes that in the main database (xlsx), has multiple otherwise identical worksheets, that only differ with the data. For example you can have "train", "test_group1", "test_group2" etc.

The output would look something like this:
This table is for all training and test data combined

The line graphs separate train and test, and present combined data. Separate line graphs are produced for each subcategory

Apologies if code is highly unoptimised!

Example:

import pandas as pd
import numpy as np
import math
import warnings
import matplotlib.pyplot as plt
from sklearn.metrics import mean_absolute_error, mean_squared_error
from tabulate import tabulate

# Suppress the performance warning
warnings.simplefilter(action='ignore', category=pd.errors.PerformanceWarning)

# Define the file paths and sheet names
data_file_path = "filepath_database.xlsx"
sheet_names = ["Train","Test"] #works with as many sheet names as you want tested. 
new_label = sheet_names + ["All Cases"]
subgroup_categories = ['All', 'Short', 'Medium', 'Long', 'Eq_Wt'] #Eq_Wt means all groups are given equal weighting even if unequal n. 
 
# Define the target variable and subcategory criteria
target_variable = 'Postoperative_BCVA' #Y
subcategory_variable = 'Preoperative_BCVA' #column in database that subgroup is based on
short_threshold = 0.3 #<0.3
medium_threshold_lower = short_threshold #>=0.3
medium_threshold_upper = 0.6 #<0.6
long_threshold = medium_threshold_upper #>=0.6

equations_df = pd.read_csv("filepath_to_hall_of_fame_csv.csv")
equation_strs = equations_df["Equation"].str.lower().tolist() #lower case to prevent frustrating errors, remove if not relevant. 

all_dfs = {i: pd.read_excel(data_file_path, sheet_name=sheet_name) for i, sheet_name in enumerate(sheet_names)}

# Define mapping between variable names and column names
variable_mapping = { #left is variable names in equations, x0,x1,x2,x3 etc. Right are column names in database (xlsx)
    'phaco_ppv': 'Phaco_PPV',
    'eifl_present': 'EIFL_Present',
    'eifl_t': 'EIFL_T',
    'eifl_v': 'EIFL_V',
    'preoperative_bcva': 'Preoperative_BCVA',
    'crt': 'CRT',
    'r_mv': 'R_MV',
    'r_fv': 'R_FV',
    'onl_t': 'ONL_T',
    'onl_mv': 'ONL_MV',
    'onl_fv': 'ONL_FV',
    'inl_t': 'INL_T',
    'inl_mv': 'INL_MV',
    'inl_fv': 'INL_FV'
}

relevant_columns = list(variable_mapping.values())

def get_sample_values(row):
    sample_values = {key: row[col] for key, col in variable_mapping.items()}
    sample_values.update({ #Add unary variables used in pysr
        'sin': math.sin,
        'relu': lambda x: max(0, x),
        'cond': lambda x, y: y if x > 0 else 0,
        'cos': math.cos,
        'tan': math.tan,
        'exp': math.exp,
        'log': math.log,
        'ln': math.log,
        'sqrt': math.sqrt,
        'square': lambda x: x ** 2,
        'inv': lambda x: 1 / x,
    })
    return sample_values

#------------------no further customisation required below this line for it to work

# Eval the equations
def evaluate_equation(row, equation_str):
    try:
        sample_values = get_sample_values(row)
        return eval(equation_str, sample_values)
    except Exception as e:
        #print(f"Error evaluating equation on row {row.name}: {e}")
        return np.nan

for i, equation_str in enumerate(equation_strs):
    for x, df in all_dfs.items():
        col_e = f'Eval_{i+1}'
        df[col_e] = df.apply(lambda row: evaluate_equation(row, equation_str), axis=1)

# Combine all individual DF into a single combined DF
combined_df = pd.concat(all_dfs.values(), ignore_index=True)
all_dfs['combined'] = combined_df

# Custom error functions that handle NaN values
def custom_mean_absolute_error(y_true, y_pred):
    mask = ~np.isnan(y_true) & ~np.isnan(y_pred)
    return np.mean(np.abs(y_true[mask] - y_pred[mask]))

def custom_mean_squared_error(y_true, y_pred):
    mask = ~np.isnan(y_true) & ~np.isnan(y_pred)
    return np.sqrt(np.mean((y_true[mask] - y_pred[mask]) ** 2))

def get_target_variable(df, category):
    if category == 'All':
        return df[target_variable]
    elif category == 'Short':
        return df[df[subcategory_variable] < short_threshold][target_variable]
    elif category == 'Medium':
        return df[(df[subcategory_variable] >= medium_threshold_lower) & (df[subcategory_variable] < medium_threshold_upper)][target_variable]
    elif category == 'Long':
        return df[df[subcategory_variable] >= long_threshold][target_variable]
    elif category == 'Eq_Wt':
        return df[target_variable]

num_equations = len(equation_strs)
num_dfs = len(all_dfs)
num_categories = len(subgroup_categories)
mae_lists = np.zeros((num_equations, num_dfs, num_categories))
rmse_lists = np.zeros((num_equations, num_dfs, num_categories))

for i, equation_str in enumerate(equation_strs):
    for x, (key, df) in enumerate(all_dfs.items()):
        for j, category in enumerate(subgroup_categories):
            true_data = get_target_variable(df, category)
            estimate_data = df[f'Eval_{i+1}'].reindex(true_data.index)  

            if j == len(subgroup_categories) - 1:  # If Eq_Wt category
                mae_lists[i][x][j] = np.mean(mae_lists[i][x][j - 3:j])
                rmse_lists[i][x][j] = np.mean(rmse_lists[i][x][j - 3:j])
            else:
                mae_lists[i][x][j] = custom_mean_absolute_error(true_data, estimate_data)
                rmse_lists[i][x][j] = custom_mean_squared_error(true_data, estimate_data)

min_mae = np.nanmin(mae_lists, axis=0)  
min_rmse = np.nanmin(rmse_lists, axis=0) 

min_mae_eq_num = np.argmin(mae_lists, axis=0)  # Equation number for min MAE
min_rmse_eq_num = np.argmin(rmse_lists, axis=0)  # Equation number for min RMSE

round_precision = 5

# Define the data for the table
data = [["Subgroup Category", "Min RMSE", "Best Equation Number (RMSE)"] + [category + " (RMSE)" for category in subgroup_categories] + [category + " (MAE)" for category in subgroup_categories]]

#RMSE table
for category in subgroup_categories:
    category_index = subgroup_categories.index(category)
    min_rmse_value = round(min_rmse[len(all_dfs) - 1, category_index], round_precision)
    equation_number_rmse = min_rmse_eq_num[len(all_dfs) - 1, category_index] + 1
    row_data = [category, min_rmse_value, equation_number_rmse]
    for i in subgroup_categories:
        row_data.append(round(rmse_lists[equation_number_rmse-1, len(all_dfs) - 1, subgroup_categories.index(i)], round_precision))
    for i in subgroup_categories:
        row_data.append(round(mae_lists[equation_number_rmse-1, len(all_dfs) - 1, subgroup_categories.index(i)], round_precision))
    data.append(row_data)

print(tabulate(data, tablefmt="grid"))

data = [["Subgroup Category", "Min MAE", "Best Equation Number (MAE)"] + [category + " (RMSE)" for category in subgroup_categories] + [category + " (MAE)" for category in subgroup_categories]]

#MAE table
for category in subgroup_categories:
    category_index = subgroup_categories.index(category)
    min_mae_value = round(min_mae[len(all_dfs) - 1, category_index], round_precision)
    equation_number_mae = min_mae_eq_num[len(all_dfs) - 1, category_index] + 1
    row_data = [category, min_mae_value, equation_number_mae]
    for i in subgroup_categories:
        row_data.append(round(rmse_lists[equation_number_mae-1, len(all_dfs) - 1, subgroup_categories.index(i)], round_precision))
    for i in subgroup_categories:
        row_data.append(round(mae_lists[equation_number_mae-1, len(all_dfs) - 1, subgroup_categories.index(i)], round_precision))
    data.append(row_data)

print(tabulate(data, tablefmt="grid"))

# Define x-axis (equation numbers)
x_axis = range(1, len(equation_strs) + 1)
# Calculate the width of the figure based on the number of equations
min_width = 10  # Minimum width of the figure
additional_width_per_equation = 0.5  # Additional width per equation
num_equations = len(equation_strs)
total_width = min_width + additional_width_per_equation * num_equations

# Plot RMSE and MAE for each eye length category
for j, subgroup_category in enumerate(subgroup_categories):
    plt.figure(figsize=(total_width, 10)) #10 is height of figure, can increase if lots of equations
    for x, df in enumerate(all_dfs):
        label = new_label[x]
        color = plt.cm.tab10(x)
        plt.plot(x_axis, rmse_lists[:, x, j], marker='o', label=label, color=color)
        if label != 'Min Test RMSE':
            plt.axhline(y=min_rmse[x, j], linestyle='--', color=color)
    plt.xlabel('Equation Number')
    plt.ylabel('RMSE')
    plt.title(f'RMSE for Each Equation ({subgroup_category} Eyes)')
    plt.xticks(x_axis)
    plt.legend()
    plt.grid(True)
    plt.show()
    
    plt.figure(figsize=(total_width, 10))
    for x, df in enumerate(all_dfs):
        label = new_label[x]
        color = plt.cm.tab10(x)
        plt.plot(x_axis, mae_lists[:, x, j], marker='o', label=label, color=color)
        if label != 'Min Test MAE':
            plt.axhline(y=min_mae[x, j], linestyle='--', color=color)
    plt.xlabel('Equation Number')
    plt.ylabel('MAE')
    plt.title(f'MAE for Each Equation ({subgroup_category} Eyes)')
    plt.xticks(x_axis)
    plt.legend()
    plt.grid(True)
    plt.show()