Add mint_stats module for statistical analyses in post-community detection workflows

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

import pandas as pd

import scipy.stats as stats

import seaborn as sns

from scipy.stats import chi2_contingency, shapiro, levene, ttest_ind, mannwhitneyu

from scipy.stats import levene, ttest_ind, chi2_contingency

from scipy.stats import mannwhitneyu

from statsmodels.stats.multitest import multipletests

custom_colors = ["#38363A", "#BAAC97", "#902C44"]

# Load Excel file

file_path = "/Users/samansarraf/Desktop/Stanford/Data_for_statistics.xlsx"

data = pd.read_excel(file_path, sheet_name="SheetName")

# Table with diagnosis

diagnosis_mapping = {1: "CN", 2: "Dementia due to AD", 3: "MCI due to AD", 5: "MCI due to other etiology",

                     -8: "MCI due to AD"}

data["final_DX"] = data["final_DX"].map(diagnosis_mapping)

def load_data(file_path, sheet_name):

    return pd.read_excel(file_path, sheet_name=sheet_name)

# Contingency table

overlap_table = pd.crosstab(data['final_DX'], data['GL Optimized CV'])

print(overlap_table)

# Table with percentages

overlap_percent = overlap_table.div(overlap_table.sum(axis=1), axis=0) * 100

print(overlap_percent)

def preprocess_data(data):

    final_dx_labels = {1: "CN", 2: "Dementia due to AD", 3: "MCI due to AD", 5: "MCI due to other etiology",

                       -8: "MCI due to AD"}

    data['final_DX'] = data['final_DX'].map(final_dx_labels)

    data['AAB'] = pd.to_numeric(data['AAB'], errors='coerce')

    return data

# Sex count

sex_count = data["SEX"].value_counts().reset_index()

sex_count.columns = ["SEX", "Count"]

print(sex_count)

# Chi-square tests

for col in ["APOE", "SEX"]:

    contingency_table = pd.crosstab(data[col], data['GL Optimized CV'])

    chi2, p, _, _ = chi2_contingency(contingency_table)

    print(f"Chi-square test for {col}: chi2 = {chi2:.3f}, p-value = {p:.4f}")

def perform_chi_square_tests(data, variables):

    for var in variables:

        contingency_table = pd.crosstab(data[var], data['GL Optimized CV'])

# Chi-square tests for community 1 and 2 only

filtered_data = data[data['GL Optimized CV'].isin([1, 2])]

for col in ["APOE", "SEX"]:

    contingency_table = pd.crosstab(filtered_data[col], filtered_data['GL Optimized CV'])

    chi2, p, _, _ = chi2_contingency(contingency_table)

        print(f"Chi-square test for {var}: chi2={chi2}, p-value={p}")

    print(f"Chi-square test for {col} (communities 1 & 2): chi2 = {chi2:.3f}, p-value = {p:.4f}")

# Age at baseline summary statistics

data['AAB'] = pd.to_numeric(data['AAB'], errors='coerce')

age_summary = data['AAB'].agg(["count", "max", "min", "mean", "std"]).reset_index()

age_summary.columns = ["Statistic", "Value"]

print(age_summary)

# Statistical comparisons

variables = ["AAB", "SEX", "MMSE", "MOCA", "ADAS13", "CDRSB", "AV45_cs", "ABETA", "TAU", "PTAU", "Hippocampus",

             "Entorhinal"]

data[variables] = data[variables].apply(pd.to_numeric, errors='coerce')

def calculate_statistics(data, variables):

results = []

for var in variables:

        if var in data.columns:

            for g1 in range(1, 3):

                for g2 in range(g1 + 1, 4):

    for g1 in [1, 2]:

        for g2 in [2, 3]:

            group1 = data[data['GL Optimized CV'] == g1][var].dropna()

            group2 = data[data['GL Optimized CV'] == g2][var].dropna()

            if len(group1) < 15 or len(group2) < 15:

                test_type = "Mann-Whitney U test"

                stat, p_value = mannwhitneyu(group1, group2, alternative='two-sided')

                        test_type = "Mann-Whitney U"

            else:

                        normal_g1 = shapiro(group1)[1] > 0.05

                        normal_g2 = shapiro(group2)[1] > 0.05

                        equal_var = levene(group1, group2)[1] > 0.05

                norm_g1 = stats.shapiro(group1).pvalue > 0.05

                norm_g2 = stats.shapiro(group2).pvalue > 0.05

                var_equal = levene(group1, group2).pvalue > 0.05

                        if normal_g1 and normal_g2:

                            stat, p_value = ttest_ind(group1, group2, equal_var=equal_var)

                            test_type = "Student's t-test" if equal_var else "Welch's t-test"

                if norm_g1 and norm_g2:

                    test_type = "Student's t-test" if var_equal else "Welch's t-test"

                    stat, p_value = ttest_ind(group1, group2, equal_var=var_equal)

                else:

                    test_type = "Mann-Whitney U test"

                    stat, p_value = mannwhitneyu(group1, group2, alternative='two-sided')

                            test_type = "Mann-Whitney U"

                    results.append([var, f"Group {g1} vs {g2}", test_type, round(stat, 3), round(p_value, 4)])

            results.append([var, f"Group {g1} vs {g2}", test_type, round(stat, 3), p_value])

    results_df = pd.DataFrame(results, columns=["Variable", "Comparison", "Test Type", "Statistic", "p-value"])

    results_df["Corrected p-value"] = np.minimum(results_df["p-value"] * len(results_df), 1.0)

# Convert results to DataFrame

results_df = pd.DataFrame(results, columns=["Variable", "Comparison", "Test_Type", "t_or_w_value", "p_value"])

# Bonferroni correction

results_df["corrected_p_value"] = multipletests(results_df["p_value"], method='bonferroni')[1]

results_df = results_df.sort_values(by=["Variable", "Comparison"])

print(results_df)

# See amyloid positivity numbers

cross_table = pd.crosstab(data['AV45_cs_pos'], data['GL Optimized CV'])

print(cross_table)

# Load Excel file

file_path = "/Users/samansarraf/Desktop/Stanford/Data_for_statistics.xlsx"

data = pd.read_excel(file_path, sheet_name="SheetName")

def boxplot_variable(data, var, ylabel):

    plt.figure(figsize=(8, 5))

    sns.boxplot(x='GL Optimized CV', y=var, hue='final_DX', data=data, palette=custom_colors)

# Custom colors

custom_colors = ["#38363A", "#BAAC97", "#902C44"]

# Ensure correct data format

data['GL Optimized CV'] = data['GL Optimized CV'].astype(str)

data['AV45_cs'] = pd.to_numeric(data['AV45_cs'], errors='coerce')

data['TAU'] = pd.to_numeric(data['TAU'], errors='coerce')

data['PTAU'] = pd.to_numeric(data['PTAU'], errors='coerce')

data['Hippocampus'] = pd.to_numeric(data['Hippocampus'], errors='coerce')

# Filter data to include only final_DX = 1, 3, and 2

filtered_data = data[data['final_DX'].isin([1, 3, 2])]

# Order final_DX

final_dx_mapping = {1: "CN", 3: "MCI", 2: "AD"}

filtered_data['final_DX'] = filtered_data['final_DX'].map(final_dx_mapping)

filtered_data['final_DX'] = pd.Categorical(filtered_data['final_DX'], categories=["CN", "MCI", "AD"], ordered=True)

# Boxplot function

def plot_boxplot(variable, y_label, title):

    plt.figure(figsize=(8, 6))

    sns.boxplot(x='GL Optimized CV', y=variable, data=filtered_data, hue='final_DX', palette=custom_colors)

    plt.xlabel("GL Optimized CV")

    plt.ylabel(ylabel)

    plt.title(f"Boxplot for {ylabel} by community")

    plt.legend(title="Diagnosis")

    # plt.savefig(f"boxplot_{var}.png", dpi=300)

    plt.ylabel(y_label)

    plt.title(title)

    plt.legend(title="Diagnosis", loc='best')

    plt.show()

if __name__ == "__main__":

    file_path = "/Users/samansarraf/Desktop/Stanford/Data_for_statistics.xlsx"

    sheet_name = "SheetName"

    data = load_data(file_path, sheet_name)

    data = preprocess_data(data)

# Create boxplots

plot_boxplot("AV45_cs", "Aβ PET", "Amyloid Levels")

plot_boxplot("TAU", "CSF TAU", "Tau Levels")

plot_boxplot("PTAU", "CSF PTAU", "pTau Levels")

plot_boxplot("Hippocampus", "Hippocampal Volume", "Hippocampal Levels")

    print(pd.crosstab(data['final_DX'], data['GL Optimized CV']))

    print(pd.crosstab(data['final_DX'], data['GL Optimized CV'], normalize='index') * 100)

    print(data['SEX'].value_counts())

# Convert 'GL Optimized CV' to categorical

data['GL Optimized CV'] = data['GL Optimized CV'].astype(str)

    perform_chi_square_tests(data, ['APOE', 'SEX'])

    variables = ["AAB", "SEX", "MMSE", "MOCA", "ADAS13", "CDRSB", "AV45_cs", "ABETA", "TAU", "PTAU", "Hippocampus",

                 "Entorhinal"]

    calculate_statistics(data, variables)

def plot_boxplot_extra(data, variable, y_label, title):

    plt.figure(figsize=(8, 6))

    sns.boxplot(x=data['GL Optimized CV'].astype(str), y=data[variable], palette=["lightblue"])

    plt.title(title)

    plt.xlabel("Communities")

    plt.ylabel(y_label)

    plt.show()

# Generate boxplots for specified variables

variables_to_plot = [

    ("MMSE", "MMSE"),

    ("MOCA", "MOCA"),

    ("CDRSB", "CDR"),

    ("ADAS13", "ADAS13"),

    ("PTAU", "PTAU"),

    ("TAU", "TAU"),

    ("AV45_cs", "AV45"),

    ("ICV", "ICV"),

    ("Hippocampus", "Hippocampus")

]

for var, label in variables_to_plot:

    plot_boxplot_extra(data, var, label, f"Boxplot for {label} by community")

    for var in ["MMSE", "MOCA", "CDRSB", "ADAS13", "PTAU", "TAU", "AV45_cs", "Hippocampus"]:

        boxplot_variable(data, var, var)

# Filter data for final_DX == 3

filtered_data = data[data['final_DX'] == 3]

# Define variables to test and ensure they are numeric

variables = ["AV45_cs", "TAU", "PTAU", "ABETA", "Hippocampus", "Entorhinal"]

filtered_data[variables] = filtered_data[variables].apply(pd.to_numeric, errors='coerce')

# Perform Mann-Whitney U tests for each variable

results = {}

cv_groups = [1, 2, 3]

for var in variables:

    var_results = {}

    for i in range(len(cv_groups) - 1):

        for j in range(i + 1, len(cv_groups)):

            group1 = filtered_data[filtered_data['GL Optimized CV'] == cv_groups[i]][var].dropna()

            group2 = filtered_data[filtered_data['GL Optimized CV'] == cv_groups[j]][var].dropna()

            if len(group1) > 0 and len(group2) > 0:

                stat, p = mannwhitneyu(group1, group2, alternative='two-sided')

                var_results[f"CV {cv_groups[i]} vs CV {cv_groups[j]}"] = {"U-value": stat, "p-value": p}

    results[var] = var_results

# Print results with U-value and p-value

for var, comparisons in results.items():

    print(f"\n{var}")

    for comp, res in comparisons.items():

        print(f"{comp}\n  U-value: {res['U-value']}\n  p-value: {res['p-value']}")

import matplotlib.pyplot as plt

import pandas as pd

import scipy.stats as stats

import seaborn as sns

from scipy.stats import levene, ttest_ind, chi2_contingency

from scipy.stats import mannwhitneyu

from statsmodels.stats.multitest import multipletests

def stats_analysis(file_path, sheet_name):

    # Load Excel file

    data_stat = pd.read_excel(file_path, sheet_name)

    data_stat_copy = data_stat.copy()

    # Table with diagnosis

    diagnosis_mapping = {1: "CN", 2: "Dementia due to AD", 3: "MCI due to AD", 5: "MCI due to other etiology",

                         -8: "MCI due to AD"}

    data_stat["final_DX"] = data_stat["final_DX"].map(diagnosis_mapping)

    # Contingency table

    overlap_table = pd.crosstab(data_stat['final_DX'], data_stat['GL Optimized CV'])

    print(overlap_table)

    # Table with percentages

    overlap_percent = overlap_table.div(overlap_table.sum(axis=1), axis=0) * 100

    print(overlap_percent)

    # Sex count

    sex_count = data_stat["SEX"].value_counts().reset_index()

    sex_count.columns = ["SEX", "Count"]

    print(sex_count)

    # Chi-square tests

    for col in ["APOE", "SEX"]:

        contingency_table = pd.crosstab(data_stat[col], data_stat['GL Optimized CV'])

        chi2, p, _, _ = chi2_contingency(contingency_table)

        print(f"Chi-square test for {col}: chi2 = {chi2:.3f}, p-value = {p:.4f}")

    # Chi-square tests for community 1 and 2 only

    filtered_data = data_stat[data_stat['GL Optimized CV'].isin([1, 2])]

    for col in ["APOE", "SEX"]:

        contingency_table = pd.crosstab(filtered_data[col], filtered_data['GL Optimized CV'])

        chi2, p, _, _ = chi2_contingency(contingency_table)

        print(f"Chi-square test for {col} (communities 1 & 2): chi2 = {chi2:.3f}, p-value = {p:.4f}")

    # Age at baseline summary statistics

    data_stat['AAB'] = pd.to_numeric(data_stat['AAB'], errors='coerce')

    age_summary = data_stat['AAB'].agg(["count", "max", "min", "mean", "std"]).reset_index()

    age_summary.columns = ["Statistic", "Value"]

    print(age_summary)

    # Statistical comparisons

    variables = ["AAB", "SEX", "MMSE", "MOCA", "ADAS13", "CDRSB", "AV45_cs", "ABETA", "TAU", "PTAU", "Hippocampus",

                 "Entorhinal"]

    data_stat[variables] = data_stat[variables].apply(pd.to_numeric, errors='coerce')

    results = []

    for var in variables:

        for g1 in [1, 2]:

            for g2 in [2, 3]:

                group1 = data_stat[data_stat['GL Optimized CV'] == g1][var].dropna()

                group2 = data_stat[data_stat['GL Optimized CV'] == g2][var].dropna()

                if len(group1) < 15 or len(group2) < 15:

                    test_type = "Mann-Whitney U test"

                    stat, p_value = mannwhitneyu(group1, group2, alternative='two-sided')

                else:

                    norm_g1 = stats.shapiro(group1).pvalue > 0.05

                    norm_g2 = stats.shapiro(group2).pvalue > 0.05

                    var_equal = levene(group1, group2).pvalue > 0.05

                    if norm_g1 and norm_g2:

                        test_type = "Student's t-test" if var_equal else "Welch's t-test"

                        stat, p_value = ttest_ind(group1, group2, equal_var=var_equal)

                    else:

                        test_type = "Mann-Whitney U test"

                        stat, p_value = mannwhitneyu(group1, group2, alternative='two-sided')

                results.append([var, f"Group {g1} vs {g2}", test_type, round(stat, 3), p_value])

    # Convert results to DataFrame

    results_df = pd.DataFrame(results, columns=["Variable", "Comparison", "Test_Type", "t_or_w_value", "p_value"])

    # Bonferroni correction

    results_df["corrected_p_value"] = multipletests(results_df["p_value"], method='bonferroni')[1]

    results_df = results_df.sort_values(by=["Variable", "Comparison"])

    print(results_df)

    # See amyloid positivity numbers

    cross_table = pd.crosstab(data_stat['AV45_cs_pos'], data_stat['GL Optimized CV'])

    print(cross_table)

    # Custom colors

    custom_colors = ["#38363A", "#BAAC97", "#902C44"]

    data_stat = data_stat_copy

    # Ensure correct data_stat format

    data_stat['GL Optimized CV'] = data_stat['GL Optimized CV'].astype(str)

    data_stat['AV45_cs'] = pd.to_numeric(data_stat['AV45_cs'], errors='coerce')

    data_stat['TAU'] = pd.to_numeric(data_stat['TAU'], errors='coerce')

    data_stat['PTAU'] = pd.to_numeric(data_stat['PTAU'], errors='coerce')

    data_stat['Hippocampus'] = pd.to_numeric(data_stat['Hippocampus'], errors='coerce')

    # Filter data_stat to include only final_DX = 1, 3, and 2

    filtered_data = data_stat[data_stat['final_DX'].isin([1, 3, 2])]

    # Order final_DX

    final_dx_mapping = {1: "CN",  2: "AD", 3: "MCI"}

    filtered_data['final_DX'] = filtered_data['final_DX'].map(final_dx_mapping)

    filtered_data['final_DX'] = pd.Categorical(filtered_data['final_DX'], categories=["CN", "MCI", "AD"], ordered=True)

    # Create box plots

    sf = True

    plot_boxplot(filtered_data, "AV45_cs", "Aβ PET", "Amyloid Levels", custom_colors, save_fig=sf)

    plot_boxplot(filtered_data, "TAU", "CSF TAU", "Tau Levels", custom_colors, save_fig=sf)

    plot_boxplot(filtered_data, "PTAU", "CSF PTAU", "pTau Levels", custom_colors, save_fig=sf)

    plot_boxplot(filtered_data, "Hippocampus", "Hippocampal Volume", "Hippocampal Levels", custom_colors, save_fig=sf)

    data_stat['GL Optimized CV'] = data_stat['GL Optimized CV'].astype(str)

    # Generate boxplots for specified variables

    variables_to_plot = [

        ("MMSE", "MMSE"),

        ("MOCA", "MOCA"),

        ("CDRSB", "CDR"),

        ("ADAS13", "ADAS13"),

        ("PTAU", "PTAU"),

        ("TAU", "TAU"),

        ("AV45_cs", "AV45"),

        ("ICV", "ICV"),

        ("Hippocampus", "Hippocampus")

    ]

    sfe = True

    for var, label in variables_to_plot:

        plot_boxplot_extra(data_stat, var, label, f"Boxplot for {label} by community", save_fig=sfe)

    # Filter data_stat for final_DX == 3

    filtered_data = data_stat[data_stat['final_DX'] == 3]

    # Define variables to test and ensure they are numeric

    variables = ["AV45_cs", "TAU", "PTAU", "ABETA", "Hippocampus", "Entorhinal"]

    filtered_data[variables] = filtered_data[variables].apply(pd.to_numeric, errors='coerce')

    # Perform Mann-Whitney U tests for each variable

    results = {}

    cv_groups = [1, 2, 3]

    for var in variables:

        var_results = {}

        for i in range(len(cv_groups) - 1):

            for j in range(i + 1, len(cv_groups)):

                group1 = filtered_data[filtered_data['GL Optimized CV'] == cv_groups[i]][var].dropna()

                group2 = filtered_data[filtered_data['GL Optimized CV'] == cv_groups[j]][var].dropna()

                if len(group1) > 0 and len(group2) > 0:

                    stat, p = mannwhitneyu(group1, group2, alternative='two-sided')

                    var_results[f"CV {cv_groups[i]} vs CV {cv_groups[j]}"] = {"U-value": stat, "p-value": p}

        results[var] = var_results

    # Print results with U-value and p-value

    for var, comparisons in results.items():

        print(f"\n{var}")

        for comp, res in comparisons.items():

            print(f"{comp}\n  U-value: {res['U-value']}\n  p-value: {res['p-value']}")

    return cross_table

def plot_boxplot(data_plot, variable, y_label, title, colors, save_fig=False):

    plt.figure(figsize=(8, 6))

    sns.boxplot(x='GL Optimized CV', y=variable, data=data_plot,

                hue='final_DX', palette=colors, order=sorted(data_plot['GL Optimized CV'].unique()))

    plt.xlabel("GL Optimized CV")

    plt.ylabel(y_label)

    plt.title(title)

    plt.legend(title="Diagnosis", loc='best')

    if save_fig:

        plt.savefig(f"{variable}_diagnosis_boxplot.png")