from numpy import var, mean, sqrt
from scipy import stats
from pandas import Series
import pandas as pd
def cohend(d1: Series, d2: Series) -> float:
# calculate the size of samples
n1, n2 = len(d1), len(d2)
# calculate the variance of the samples
s1, s2 = var(d1, ddof=1), var(d2, ddof=1)
# calculate the pooled standard deviation
s = sqrt(((n1 - 1) * s1 + (n2 - 1) * s2) / (n1 + n2 - 2))
# calculate the means of the samples
u1, u2 = mean(d1), mean(d2)
# return the effect size
return (u1 - u2) / s
cpet_data = pd.read_csv('cpet_clinical_data.tsv', sep='\t')
metrics = [
'DBP', # Diastolic Blood Pressure
'HR', # Heart Rate
#'OUES', # Oxygen Uptake Equivalent Slope
'PETCO2', # End-Tidal CO2
'PETO2', # End-Tidal O2
'PP', # Pulse Pressure
'RER', # Respiratory Exchange Ratio
'RPE', # Rating of Perceived Exertion
'RPM', # Pedal RPM
'RPP', # Rate Pressure Product
'RR', # Respiratory Rate
'SBP', # Systolic Blood Pressure
'VCO2', # Carbon Dioxide Production
'VO2', # Oxygen Consumption
'VO2_HR', # Oxygen Consumption per Heart Rate
'VO2_t', # Oxygen consumption (ml/min)
'Ve_BTPS', # Ventilation (BTPS)
'Ve_VCO2', # Ventilatory Equivalent for CO2
'Ve_VO2', # Ventilatory Equivalent for O2
'Vt_BTPS_L', # Tidal Volume (BTPS in Liters)
'wkld', # Workload
'time_sec'
]
# Convert all metric columns to numeric, coercing errors
for metric in metrics:
cpet_data[metric] = pd.to_numeric(cpet_data[metric], errors='coerce')
# Split the data by time points
time_points = ['AT', 'max', 'rest']
# Initialize an empty list to store differences
diff_data = []
for time_point in time_points:
# Filter data for the current time point
df_tp = cpet_data[cpet_data['Time_Point'] == time_point]
# Split the data into Day 1 and Day 2
day1 = df_tp[df_tp['Study_Visit'] == 'D1']
day2 = df_tp[df_tp['Study_Visit'] == 'D2']
# Merge day1 and day2 on ParticipantID
merged = pd.merge(day1, day2, on=['ParticipantID', 'matched_pair', 'sex', 'Time_Point', 'race', 'phenotype'], suffixes=('_D1', '_D2'))
# Calculate the absolute difference for each metric
for metric in metrics:
merged[metric + '_abs_diff'] = merged[metric + '_D2'] - merged[metric + '_D1']
# Calculate the percentage difference for each metric
for metric in metrics:
merged[metric + '_pct_diff'] = ((merged[metric + '_D2'] - merged[metric + '_D1']) / merged[metric + '_D1']) * 100
# Keep only relevant columns and add time point information
diff_columns = ['ParticipantID', 'matched_pair', 'sex', 'Time_Point', 'phenotype'] + [metric + '_pct_diff' for metric in metrics]
merged['Time_Point'] = time_point
diff_data.append(merged)
# Combine the differences for all time points
cpet_diff = pd.concat(diff_data)
# List of absolute difference metrics
abs_diff_metrics = [metric + '_abs_diff' for metric in metrics]
# Dictionary to store DataFrames for each time point
time_point_dfs = {}
# Calculate differences, Cohen's d, and p-values for each time point
for time_point in time_points:
# Filter data for the current time point
tp_data = cpet_diff[cpet_diff['Time_Point'] == time_point]
# Initialize lists to store results
results = []
for metric in abs_diff_metrics:
# Get data for each group
mecfs_data = tp_data[tp_data['phenotype'] == 'MECFS'][metric].dropna()
hc_data = tp_data[tp_data['phenotype'] == 'HC'][metric].dropna()
# Calculate means
mecfs_mean = mecfs_data.mean()
hc_mean = hc_data.mean()
mean_diff = mecfs_mean - hc_mean
cohens_d = cohend(mecfs_data, hc_data)
# Perform t-test
t_stat, p_value = stats.ttest_ind(mecfs_data, hc_data)
# Store results
results.append({
'Metric': metric,
'Mean_Difference': mean_diff,
'Cohens_d': cohens_d,
'p_value': p_value
})
# Convert results to DataFrame
df = pd.DataFrame(results)
# Sort by p-value
df = df.reindex(df['p_value'].sort_values(ascending=True).index)
# Store the DataFrame in the dictionary
time_point_dfs[time_point] = df
at_diff = time_point_dfs['AT']
max_diff = time_point_dfs['max']
rest_diff = time_point_dfs['rest']
