Postmortem brain MRI reveals differential associations of subcortical and limbic volumes with cortical thinning and neurodegenerative pathologies.

Code ↔ Paper

11 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 11 matches

[1] § RESULTS › AD, FTLD‐TDP, and FTLD‐tau have lower subcortical and limbic volumes than LBD ↔ subcortical_analyses.ipynb, lines 111–128 · score 0.69 · Lewy body disease, likelihood ratio, ICV normalized postmortem, subcortical volumes, Pairwise, FTLD TDP
[2] § RESULTS › Subcortical and limbic volume loss tracks with cortical thinning differently in each disease ↔ subcortical_cortical_analyses.ipynb, lines 169–176 · score 0.68 · partial Spearman correlations, normalized weighted, ICV normalized, cortical thickness, composite, coefficients
[3] § RESULTS › AD, FTLD‐TDP, and FTLD‐tau have lower subcortical and limbic volumes than LBD ↔ subcortical_cortical_analyses.ipynb, lines 410–432 · score 0.66 · Boxplots compare, Lewy body disease, likelihood ratio, FTLD tau, Pairwise, FTLD TDP
[4] § RESULTS › Polypathological models suggest primary pathology mostly drives subcortical and limbic volume loss ↔ subcortical_analyses.ipynb, lines 2632–2637 · score 0.65 · standardized coefficients, pathology predictors, OLS models, FDR correction, fitting, covariates
[5] § METHODS › Semi‐quantitative neuropathology and histology measures ↔ notebooks/analysis_notebook_10_5_23.ipynb, lines 428–523 · score 0.60 · angular gyrus, cingulate, superior, temporal, field, frontal
[6] § METHODS › Semi‐quantitative neuropathology and histology measures ↔ subcortical_analyses.ipynb, lines 23–54 · score 0.59 · Semi quantitative, antibodies, pons, brainstem, neurons, scores
[7] § METHODS › Statistical analysis › Structural group differences and subcortical–cortical associations ↔ subcortical_cortical_analyses.ipynb, lines 328–407 · score 0.59 · DKT atlas, partial Spearman, cortical region, nucleus, thickness, correlations
[8] § RESULTS › Primary AD is linked to worse subcortical and cortical atrophy than primary LBD in donors with both AD and LBD copathologies ↔ subcortical_analyses.ipynb, lines 1248–1288 · score 0.55 · normalized postmortem subcortical, Pathology burden, Primary pathology, linear, fit, age
[9] § RESULTS › Histologic neuronal loss and gliosis mediate relationships between subcortical and limbic pathology and volume loss ↔ subcortical_analyses.ipynb, lines 545–562 · score 0.53 · regional pathology burden, FDR correction, postmortem subcortical, age, sex, education
[10] § RESULTS › Polypathological models suggest primary pathology mostly drives subcortical and limbic volume loss ↔ subcortical_analyses.ipynb, lines 3070–3218 · score 0.52 · polypathology model, multiple pathologies, FTLD TDP, synuclein, tauopathy, regression
[11] § RESULTS › Primary AD is linked to worse subcortical and cortical atrophy than primary LBD in donors with both AD and LBD copathologies ↔ notebooks/analysis_notebook_10_5_23.ipynb, lines 428–523 · score 0.51 · superior parietal, globus pallidus, cortex, thicknesses, caudate, cortical

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Jupyter notebook · 6,168 lines · 195 KB · no license · 6 matches

# %%
"""
Imports core Python libraries for data wrangling, statistical analysis, and visualization, including pandas, NumPy, seaborn, and matplotlib.
Also imports tools for t-tests, Tukey post-hoc comparisons, formatted summary tables, custom plot legends.
"""
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import warnings
import numpy as np
from matplotlib.patches import Patch
from scipy.stats import ttest_rel, ttest_ind
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import seaborn as sns
import pandas as pd
import pingouin as pg
from prettytable import PrettyTable
import warnings
import warnings
from statsmodels.tools.sm_exceptions import ConvergenceWarning
# %%
"""
Converts semi-quantitative regional pathology scores into numeric values for analysis.
The function identifies region-specific pathology columns, replaces ordinal/string scores such as 'Rare', '1+', '2+', and '3+' with numeric values, converts unavailable entries to NaN, and returns a cleaned copy of the dataframe.
"""
def semiq_scores_to_numeric(df):
# All the columns for a given region
# Already done
measures = [ 'Tau', 'ThioPlaques', 'AntibodyPlaques', 'aSyn', 'Ubiquitin',
'Gliosis', 'NeuronLoss', 'TDP43', 'Other', 'Update', 'Angiopathy']
# Find all the columns with regional pathology scores
regions = [ 'Amyg', 'DG', 'CS', 'EC', 'MF', 'Ang', 'SMT', 'Cing', 'OC',
'Neocortical', 'CP', 'GP', 'TS', 'Subcortical', 'MB', 'SN',
'Pons', 'LC', 'Med', 'CB', 'SC', 'Brainstem', 'MC', 'OFC']
# Cross product of these lists
cols_semiq = [x + y for x in regions for y in measures ]
# Conversion rules
repl_semiq = { 'Rare': 0.5, '2+': 2.0, '3+': 3.0, '1+': 1.0, '0': 0.0,
'Presumed 0': 0.0, 'Not Avail': np.nan, 'Not Done': np.nan,
'Not Available': np.nan, '.': np.nan }
# Apply the replacements
df_copy = df.copy()
for m in cols_semiq:
if m in df.columns:
df_copy[m] = pd.to_numeric(df[m].replace(repl_semiq), errors='coerce')
return df_copy
# %%
"""
Defines a lookup dictionary mapping SynthSeg integer label IDs to their corresponding anatomical structure names.
This mapping is used to translate segmentation label values into readable brain-region names for downstream volume extraction, summaries, and plots.
"""
label_to_structure = {
0: "Background",
2: "Left cerebral white matter",
3: "Left cerebral cortex",
4: "Left lateral ventricle",
5: "Left inferior lateral ventricle",
7: "Left cerebellum white matter",
8: "Left cerebellum cortex",
10: "Left thalamus",
11: "Left caudate",
12: "Left putamen",
13: "Left pallidum",
14: "3rd ventricle",
15: "4th ventricle",
16: "Brain-stem",
17: "Left hippocampus",
18: "Left amygdala",
24: "CSF (SynthSeg 2.0 only)",
26: "Left accumbens area",
28: "Left ventral DC",
41: "Right cerebral white matter",
42: "Right cerebral cortex",
43: "Right lateral ventricle",
44: "Right inferior lateral ventricle",
46: "Right cerebellum white matter",
47: "Right cerebellum cortex",
49: "Right thalamus",
50: "Right caudate",
51: "Right putamen",
52: "Right pallidum",
53: "Right hippocampus",
54: "Right amygdala",
58: "Right accumbens area",
60: "Right ventral DC",
}
# %%
"""
Loads the main analysis dataframe from a CSV file.
The resulting dataframe, df, is used as the input table for subsequent cleaning, statistical analysis, and visualization steps.
"""
df = pd.read_csv("data_frame.csv")
# %%
# %%
"""
Computes ICV-normalized postmortem subcortical volumes and performs covariate-adjusted pairwise likelihood-ratio tests across neuropathological diagnostic groups.
The script applies FDR correction to pairwise p-values and generates publication-style box/strip plots with significance annotations for limbic and subcortical structures.
"""
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from itertools import combinations
from scipy.stats import chi2
from statsmodels.formula.api import ols
import pingouin as pg
df_use = df.copy()
order = ["alzheimer's disease", "lewy body disease", "ftld-tdp", "tauopathies"]
# Pretty x-labels
x_labels = [
"Alzheimer’s\ndisease",
"Lewy body\ndisease",
"FTLD-TDP",
"Tauopathies"
]
# Panel groups
panelA = ["hippocampus", "amygdala", "accumbens_area"]
panelB = ["thalamus", "caudate", "putamen", "pallidum"]
all_structures = panelA + panelB
covars = ["AgeatDeath", "Sex", "Education", "PMI"]
palette = ["#3366CC", "#DC3912", "#109618", "#FF9900"]
sns.set(style="whitegrid", context="talk", font_scale=1.2)
# ============================================================
# NORMALIZE
# ============================================================
for s in all_structures:
pm = f"postmortem_{s}"
if pm in df_use.columns:
df_use[f"{s}_norm"] = df_use[pm] / df_use["antemortem_icv"]
needed = [f"{s}_norm" for s in all_structures] + ["NPDx1"] + covars
df_use = df_use[needed].dropna()
# ============================================================
# PAIRWISE LRTs
# ============================================================
pairwise_results = []
for s in all_structures:
ycol = f"{s}_norm"
for g1, g2 in combinations(order, 2):
d = df_use[df_use["NPDx1"].isin([g1, g2])]
if len(d) < 10:
continue
reduced = ols(f"{ycol} ~ " + " + ".join(covars), data=d).fit()
full = ols(f"{ycol} ~ C(NPDx1) + " + " + ".join(covars), data=d).fit()
lr = 2*(full.llf - reduced.llf)
df_diff = full.df_model - reduced.df_model
p = chi2.sf(lr, df_diff)
pairwise_results.append({
"Structure": s, "Group1": g1, "Group2": g2, "p_raw": p
})
pairwise_df = pd.DataFrame(pairwise_results)
reject, p_corr = pg.multicomp(pairwise_df["p_raw"], method="fdr_bh")
pairwise_df["p_FDR"] = p_corr
pairwise_df["Sig"] = pairwise_df["p_FDR"].apply(
lambda p: "***" if p < 0.001 else "**" if p < 0.01 else "*" if p < 0.05 else ""
)
# ============================================================
# PLOTTING
# ============================================================
fig = plt.figure(figsize=(26, 18))
# TOP PANEL (3)
gsA = fig.add_gridspec(
1, 3, left=0.05, right=0.97,
top=0.92, bottom=0.56, wspace=0.33
)
axesA = [fig.add_subplot(gsA[0, k]) for k in range(3)]
# BOTTOM PANEL (4)
gsB = fig.add_gridspec(
1, 4, left=0.05, right=0.97,
top=0.50, bottom=0.12, wspace=0.30
)
axesB = [fig.add_subplot(gsB[0, k]) for k in range(4)]
axes = axesA + axesB
def plot_struct(ax, s):
ycol = f"{s}_norm"
d = df_use[["NPDx1", ycol]].dropna()
# -------------- PLOT GROUPS --------------
for idx, g in enumerate(order):
vals = d.loc[d["NPDx1"] == g, ycol]
tmp = pd.DataFrame({"group": [g]*len(vals), "y": vals})
sns.boxplot(
data=tmp, x="group", y="y",
color=palette[idx], ax=ax,
width=0.55, fliersize=0,
linewidth=1.3, boxprops=dict(alpha=0.72)
)
sns.stripplot(
data=tmp, x="group", y="y",
color="black", size=4, alpha=0.55,
ax=ax, jitter=0.15
)
# REMOVE DEFAULT LABELS
ax.set_xlabel("")
ax.set_ylabel("")
# Set pretty labels
ax.set_xticklabels(x_labels, fontsize=11)
ymin, ymax = d[ycol].min(), d[ycol].max()
yr = ymax - ymin
ax.set_ylim(ymin - 0.06*yr, ymax + 0.45*yr)
ax.set_title(s.replace("_", " ").capitalize(), fontsize=16, fontweight="bold")
ax.grid(axis="y", linestyle=":", alpha=0.45)
# ---------------- SIG BARS -----------------
pairs = pairwise_df[pairwise_df["Structure"] == s]
y_offset = 0.020 * yr
y_pos = ymax + 0.10 * yr
for _, row in pairs.iterrows():
if row["Sig"]:
g1, g2 = row["Group1"], row["Group2"]
x1 = order.index(g1)
x2 = order.index(g2)
ax.plot([x1, x1, x2, x2],
[y_pos, y_pos+y_offset, y_pos+y_offset, y_pos],
lw=1.35, color="black")
ax.text((x1+x2)/2, y_pos + y_offset*0.8,
row["Sig"], ha="center",
fontsize=13, fontweight="bold")
y_pos += y_offset * 1.9
# Draw all structures
for ax, s in zip(axes, all_structures):
plot_struct(ax, s)
# ============================================================
# GLOBAL LABELS + TITLE
# ============================================================
fig.text(0.001, 0.55, "Normalized volume",
va="center", rotation="vertical",
fontsize=20, fontweight="bold")
fig.text(0.50, 0.06, "Diagnostic groups",
ha="center", fontsize=20, fontweight="bold")
plt.subplots_adjust(top=0.93)
fig.suptitle(
"Postmortem subcortical volumes differentiates neuropathological groups",
fontsize=26, fontweight="bold"
)
plt.tight_layout()
plt.savefig("Postmortem subcortical volumes differentiates neuropathological groups.png", dpi=600, bbox_inches="tight")
plt.show()
# %%
# %%
"""
Generates a supplementary DOCX table of pairwise covariate-adjusted likelihood-ratio tests
for ICV-normalized postmortem subcortical volumes across neuropathological diagnostic groups.
The script computes adjusted means, mean differences, standard errors, 95% confidence intervals,
raw p-values, FDR-corrected p-values, and significance labels, then exports publication-ready
tables to a Word document.
"""
import pandas as pd
import numpy as np
from itertools import combinations
from statsmodels.formula.api import ols
from scipy.stats import chi2, shapiro
from statsmodels.stats.diagnostic import het_breuschpagan
from docx import Document
from docx.shared import Pt
from docx.enum.table import WD_TABLE_ALIGNMENT
# ============================================================
# SETTINGS
# ============================================================
order = ["alzheimer's disease", "lewy body disease", "ftld-tdp", "tauopathies"]
covars = ["AgeatDeath", "Sex", "Education", "PMI"]
structures = [
"hippocampus","amygdala","accumbens_area",
"thalamus","caudate","putamen","pallidum"
]
rows = []
# ============================================================
# LOOP OVER STRUCTURES AND GROUP PAIRS
# ============================================================
for s in structures:
ycol = f"{s}_norm"
for g1, g2 in combinations(order, 2):
d = df_use[df_use["NPDx1"].isin([g1, g2])].copy()
if len(d) < 10:
continue
n_g1 = (d["NPDx1"] == g1).sum()
n_g2 = (d["NPDx1"] == g2).sum()
# ----------------- MODELS -----------------
reduced = ols(f"{ycol} ~ " + " + ".join(covars), data=d).fit()
full = ols(f"{ycol} ~ C(NPDx1) + " + " + ".join(covars), data=d).fit()
# ----------------- LRT -----------------
ll_full = full.llf
ll_reduced = reduced.llf
LR = 2 * (ll_full - ll_reduced)
df_diff = full.df_model - reduced.df_model
p_raw = chi2.sf(LR, df_diff)
# ----------------- ADJUSTED MEANS -----------------
d["_pred"] = full.fittedvalues
mean_g1 = d.loc[d["NPDx1"] == g1, "_pred"].mean()
mean_g2 = d.loc[d["NPDx1"] == g2, "_pred"].mean()
diff = mean_g1 - mean_g2
# ----------------- SE + CI -----------------
contrast = np.zeros(len(full.params))
idx1 = list(full.params.index).index(f"C(NPDx1)[T.{g2}]") if f"C(NPDx1)[T.{g2}]" in full.params else None
if idx1 is not None:
contrast[idx1] = -1
se = np.sqrt(contrast @ full.cov_params() @ contrast.T)
ci_low = diff - 1.96 * se
ci_high = diff + 1.96 * se
# ----------------- DIAGNOSTICS -----------------
res = full.resid
shapiro_p = shapiro(res)[1]
lm, lm_pvalue, fval, f_pvalue = het_breuschpagan(res, full.model.exog)
bp_p = lm_pvalue
cooks = full.get_influence().cooks_distance[0].max()
rows.append({
"Structure": s,
"Group1": g1,
"Group2": g2,
"n_g1": n_g1,
"n_g2": n_g2,
"AdjMean_Group1": mean_g1,
"AdjMean_Group2": mean_g2,
"AdjMean_Diff": diff,
"SE_Diff": se,
"CI95_low": ci_low,
"CI95_high": ci_high,
"LR": LR,
"p_raw": p_raw
})
# ============================================================
# DATAFRAME + FDR
# ============================================================
sup_df = pd.DataFrame(rows)
from pingouin import multicomp
_, p_corr = multicomp(sup_df["p_raw"], method="fdr_bh")
sup_df["p_FDR"] = p_corr
sup_df["Sig"] = sup_df["p_FDR"].apply(
lambda p: "***" if p < 0.001 else "**" if p < 0.01 else "*" if p < 0.05 else ""
)
###############################################################################
# DOCX EXPORT — CLEAN PUBLICATION-READY
###############################################################################
doc = Document()
# -------------------------------------------------------------
# Global title + caption
# -------------------------------------------------------------
doc.add_heading("Supplementary Table 1", level=1)
caption_text = (
"Pairwise covariate-adjusted likelihood ratio tests across diagnostic groups "
"for each postmortem subcortical structure. Covariate-adjusted marginal means, "
"adjusted mean differences (Δ), standard errors, and 95% confidence intervals "
"are reported. FDR correction was applied globally across all pairwise tests. "
"Scientific notation is used for all numeric values. Significant FDR-corrected "
"p-values are bolded and marked with *, **, or ***."
)
doc.add_paragraph(caption_text)
doc.add_paragraph("\n")
# -------------------------------------------------------------
# Abbreviations
# -------------------------------------------------------------
abbr = {
"alzheimer's disease": "AD",
"lewy body disease": "LBD",
"ftld-tdp": "FTLD-TDP",
"tauopathies": "Tau"
}
# -------------------------------------------------------------
# Helper functions
# -------------------------------------------------------------
def sci(x):
try:
return f"{float(x):.2e}"
except:
return str(x)
def bold_run(cell, text):
p = cell.paragraphs[0]
run = p.add_run(text)
run.bold = True
def add_text(cell, text):
cell.paragraphs[0].add_run(text)
# -------------------------------------------------------------
# Generate tables
# -------------------------------------------------------------
structures_sorted = [
"hippocampus","amygdala","accumbens_area",
"thalamus","caudate","putamen","pallidum"
]
for s in structures_sorted:
sub = sup_df[sup_df["Structure"] == s].copy()
if len(sub) == 0:
continue
sub = sub.sort_values(by=["p_FDR", "LR"], ascending=[True, False])
doc.add_heading(f"{s.capitalize()}", level=2)
# Build columns
sub["Comparison"] = sub.apply(
lambda r: f"{abbr[r['Group1']]} (n={int(r['n_g1'])}) vs {abbr[r['Group2']]} (n={int(r['n_g2'])})",
axis=1
)
sub["AdjMeans"] = sub.apply(
lambda r: f"{abbr[r['Group1']]} = {sci(r['AdjMean_Group1'])}\n"
f"{abbr[r['Group2']]} = {sci(r['AdjMean_Group2'])}",
axis=1
)
sub["Diff"] = sub.apply(lambda r: sci(r["AdjMean_Diff"]), axis=1)
sub["CI"] = sub.apply(
lambda r: f"[{sci(r['CI95_low'])}, {sci(r['CI95_high'])}]",
axis=1
)
sub["p_FDR_sup"] = sub.apply(lambda r: sci(r["p_FDR"]) + r["Sig"], axis=1)
final_cols = [
"Comparison", "AdjMeans", "Diff",
"SE_Diff", "CI", "LR", "p_raw", "p_FDR_sup"
]
header_labels = {
"Diff": "Δ (Adj Diff)",
"SE_Diff": "SE",
"p_FDR_sup": "p_FDR"
}
table = doc.add_table(rows=1, cols=len(final_cols))
table.alignment = WD_TABLE_ALIGNMENT.CENTER
table.style = "Table Grid"
hdr = table.rows[0].cells
for j, col in enumerate(final_cols):
bold_run(hdr[j], header_labels.get(col, col))
for _, row in sub.iterrows():
cells = table.add_row().cells
for j, col in enumerate(final_cols):
val = row[col]
if col in ["SE_Diff", "LR", "p_raw"]:
val = sci(val)
if col == "p_FDR_sup" and any(x in str(val) for x in ["*", "**", "***"]):
bold_run(cells[j], str(val))
else:
add_text(cells[j], str(val))
doc.add_paragraph("\n")
# -------------------------------------------------------------
# Save output
# -------------------------------------------------------------
doc.save("group_wise_comparisons.docx")
print("Saved polished Word file: group_wise_comparisons.docx")
# %%
"""
Generates a 4-group grid of partial Spearman correlations between regional pathology burden and ICV-normalized postmortem subcortical volumes.
For each diagnostic group and structure, the script selects the relevant pathology marker, adjusts correlations for age at death, sex, PMI, and education, applies FDR correction within each diagnostic group, and saves a publication-style multi-panel figure.
"""
matplotlib.rcParams.update({
"font.family": "DejaVu Sans",
"axes.labelsize": 14,
"axes.titlesize": 14,
"xtick.labelsize": 12,
"ytick.labelsize": 12,
})
sns.set(style="whitegrid", context="talk")
warnings.filterwarnings("ignore", category=RuntimeWarning)
np.seterr(divide="ignore", invalid="ignore")
# ============================================================
# INPUT DATA
# ============================================================
df_use = df.copy()
df_use["NPDx1"] = df_use["NPDx1"].astype(str).str.strip().str.lower()
df_use["Sex"] = df_use["Sex"].astype("category").cat.codes
# ============================================================
# DISEASE → PRIMARY PATHOLOGY MAPPING (4 GROUPS)
# ============================================================
pathology_markers = {
"alzheimer's disease": "Tau",
"lewy body disease": "aSyn",
"ftld-tdp": "TDP43",
"tauopathies": "Tau",
}
# Labels shown on x-axis
display_labels = {
"Tau": "p-tau",
"aSyn": "α-synuclein",
"TDP43": "TDP-43",
}
# Row labels if you want to use them later
row_titles = {
"alzheimer's disease": "Alzheimer’s disease (AD pathology)",
"lewy body disease": "Lewy body disease (α-synuclein)",
"ftld-tdp": "FTLD-TDP (TDP-43)",
"tauopathies": "FTLD-Tau (p-tau)",
}
# ============================================================
# PREFIX MAP
# ============================================================
region_map = {
"hippocampus": "EC_CS_DG",
"amygdala": "Amyg",
"caudate": "CP",
"putamen": "CP",
"thalamus": "TS",
"pallidum": "GP",
}
structures = ["hippocampus", "amygdala", "caudate", "putamen", "thalamus", "pallidum"]
palette = [
"#56B4E9", "#E69F00", "#009E73",
"#CC79A7", "#F0E442", "#0072B2"
]
# ============================================================
# NORMALIZE POSTMORTEM VOLUMES
# ============================================================
for s in structures:
pm_col = f"postmortem_{s}"
if pm_col in df_use.columns and "antemortem_icv" in df_use.columns:
df_use[f"{pm_col}_norm"] = df_use[pm_col] / df_use["antemortem_icv"]
# ============================================================
# FIGURE SETUP
# ============================================================
groups = list(pathology_markers.keys()) # 4 groups
nrows, ncols = len(groups), len(structures)
fig, axes = plt.subplots(
nrows, ncols,
figsize=(5.3 * ncols, 4.4 * nrows),
sharey=False
)
# Make axes always 2D
if nrows == 1:
axes = np.expand_dims(axes, axis=0)
plt.subplots_adjust(hspace=0.60, wspace=0.55, top=0.92, bottom=0.07, left=0.09, right=0.97)
# ============================================================
# MAIN LOOP (PARTIAL SPEARMAN)
# ============================================================
for row_idx, (group_name, marker_suffix) in enumerate(pathology_markers.items()):
gdf = df_use[df_use["NPDx1"] == group_name].copy()
# center row title across the full row of subplots
left_ax = axes[row_idx, 0]
right_ax = axes[row_idx, -1]
left_pos = left_ax.get_position()
right_pos = right_ax.get_position()
x_center = (left_pos.x0 + right_pos.x1) / 2
y_top = left_pos.y1 + 0.015
fig.text(
x_center,
y_top,
row_titles[group_name],
ha="center",
va="bottom",
fontsize=22
)
all_stats = []
all_pvals = []
# ---------- compute stats ----------
for col_idx, s in enumerate(structures):
ax = axes[row_idx, col_idx]
post_col = f"postmortem_{s}_norm"
prefix = region_map[s]
if post_col not in gdf.columns:
ax.text(0.5, 0.5, "No volume", ha="center", va="center", fontsize=12)
ax.set_axis_off()
continue
path_candidates = [c for c in gdf.columns if c.lower().startswith(prefix.lower())]
path_cols = [c for c in path_candidates if marker_suffix.lower() in c.lower()]
if not path_cols:
ax.text(0.5, 0.5, "No pathology", ha="center", va="center", fontsize=12)
ax.set_axis_off()
continue
path_col = path_cols[0]
needed_cols = [path_col, post_col, "AgeatDeath", "Sex", "PMI", "Education"]
missing_cols = [c for c in needed_cols if c not in gdf.columns]
if missing_cols:
ax.text(0.5, 0.5, "Missing covariates", ha="center", va="center", fontsize=12)
ax.set_axis_off()
continue
d = gdf[needed_cols].dropna()
if len(d) < 5:
ax.text(0.5, 0.5, "N too small", ha="center", va="center", fontsize=12)
ax.set_axis_off()
continue
res = pg.partial_corr(
data=d,
x=path_col,
y=post_col,
covar=["AgeatDeath", "Sex", "PMI", "Education"],
method="spearman"
)
r = res["r"].iloc[0]
p = res["p-val"].iloc[0]
# Preserve original subplot column index
all_stats.append({
"col_idx": col_idx,
"structure": s,
"n": len(d),
"r": r,
"p_raw": p,
"path_col": path_col,
"post_col": post_col,
})
all_pvals.append(p)
# ---------- FDR correction ----------
if len(all_pvals) > 0:
reject, p_corr = pg.multicomp(all_pvals, method="fdr_bh")
else:
reject, p_corr = [], []
# attach corrected p-values back to stats
for i in range(len(all_stats)):
all_stats[i]["p_fdr"] = p_corr[i]
all_stats[i]["reject"] = reject[i]
# ---------- plot ----------
for stat in all_stats:
col_idx = stat["col_idx"]
s = stat["structure"]
r = stat["r"]
p_raw = stat["p_raw"]
p_fdr = stat["p_fdr"]
path_col = stat["path_col"]
post_col = stat["post_col"]
ax = axes[row_idx, col_idx]
d = gdf[[path_col, post_col, "AgeatDeath", "Sex", "PMI", "Education"]].dropna()
sns.regplot(
data=d,
x=path_col,
y=post_col,
scatter_kws=dict(alpha=0.7, s=45),
line_kws=dict(color=palette[col_idx % len(palette)], lw=2),
color=palette[col_idx % len(palette)],
ax=ax
)
sig = "***" if p_fdr < 0.001 else "**" if p_fdr < 0.01 else "*" if p_fdr < 0.05 else ""
ax.set_title(
f"{s.capitalize()} (ρ={r:.2f}, FDR p={p_fdr:.3f}{sig})",
fontsize=16,
fontweight="bold"
)
ax.set_xlabel(display_labels[marker_suffix], fontsize=16)
ax.set_ylabel("")
ax.grid(True, linestyle=":", alpha=0.5)
ax.margins(x=0.12)
ax.tick_params(axis="y", pad=10)
# ============================================================
# GLOBAL Y-AXIS LABEL
# ============================================================
fig.text(
0.055, 0.5,
"Normalized volume",
va="center", ha="center",
rotation=90,
fontsize=18, fontweight="bold"
)
# ============================================================
# GLOBAL TITLE
# ============================================================
fig.suptitle(
"Partial Spearman correlation between regional pathology burden and normalized postmortem volume",
fontsize=22,
y=0.99
)
plt.savefig(
"Partial_Spearman_pathology_vs_volume_4groups_corrected.png",
dpi=600,
bbox_inches="tight"
)
plt.show()
# %%
"""
Generates an LBD-only 3×3 grid of partial Spearman correlations between thalamic α-synuclein burden and ICV-normalized postmortem limbic/subcortical volumes.
The script adjusts correlations for age at death, sex, PMI, and education, applies FDR correction across tested structures, and saves a publication-style multi-panel regression figure.
"""
# ============================================================
# GLOBAL AESTHETICS
# ============================================================
matplotlib.rcParams.update({
"font.family": "DejaVu Sans",
"axes.labelsize": 14,
"axes.titlesize": 14,
"xtick.labelsize": 12,
"ytick.labelsize": 12,
})
sns.set(style="whitegrid", context="talk")
warnings.filterwarnings("ignore", category=RuntimeWarning)
np.seterr(divide="ignore", invalid="ignore")
# ============================================================
# INPUT DATA
# ============================================================
df_use = df.copy()
df_use["NPDx1"] = df_use["NPDx1"].astype(str).str.strip().str.lower()
df_use["Sex"] = df_use["Sex"].astype("category").cat.codes
# ============================================================
# KEEP ONLY LBD
# ============================================================
gdf = df_use[df_use["NPDx1"] == "lewy body disease"].copy()
# ============================================================
# STRUCTURES
# ============================================================
structures = [
"hippocampus", "amygdala", "caudate",
"putamen", "thalamus", "pallidum",
"accumbens_area"
]
palette = [
"#56B4E9", "#E69F00", "#009E73",
"#CC79A7", "#F0E442", "#0072B2",
"#D55E00"
]
# ============================================================
# NORMALIZE POSTMORTEM VOLUMES
# ============================================================
for s in structures:
pm_col = f"postmortem_{s}"
if pm_col in gdf.columns and "antemortem_icv" in gdf.columns:
gdf[f"{pm_col}_norm"] = gdf[pm_col] / gdf["antemortem_icv"]
# ============================================================
# FIND THALAMIC aSyn COLUMN
# ============================================================
thal_asyn_candidates = [c for c in gdf.columns if c.lower().startswith("ts")]
thal_asyn_cols = [c for c in thal_asyn_candidates if "asyn" in c.lower()]
if len(thal_asyn_cols) == 0:
raise ValueError("No thalamic aSyn column found.")
thal_asyn_col = thal_asyn_cols[0]
print("Using thalamic aSyn column:", thal_asyn_col)
# ============================================================
# FIGURE SETUP: 3 x 3
# ============================================================
fig, axes = plt.subplots(3, 3, figsize=(16, 14), sharey=False)
axes = axes.flatten()
plt.subplots_adjust(hspace=0.5, wspace=0.4, top=0.88, bottom=0.08, left=0.08, right=0.98)
all_stats = []
all_pvals = []
# ============================================================
# COMPUTE STATS
# ============================================================
for idx, s in enumerate(structures):
ax = axes[idx]
post_col = f"postmortem_{s}_norm"
if post_col not in gdf.columns:
ax.text(0.5, 0.5, f"No volume\n{post_col}", ha="center", va="center", fontsize=12)
ax.set_axis_off()
continue
needed_cols = [thal_asyn_col, post_col, "AgeatDeath", "Sex", "PMI", "Education"]
missing_cols = [c for c in needed_cols if c not in gdf.columns]
if missing_cols:
ax.text(0.5, 0.5, "Missing covariates", ha="center", va="center", fontsize=12)
ax.set_axis_off()
continue
d = gdf[needed_cols].dropna()
if len(d) < 5:
ax.text(0.5, 0.5, "N too small", ha="center", va="center", fontsize=12)
ax.set_axis_off()
continue
res = pg.partial_corr(
data=d,
x=thal_asyn_col,
y=post_col,
covar=["AgeatDeath", "Sex", "PMI", "Education"],
method="spearman"
)
r = res["r"].iloc[0]
p = res["p-val"].iloc[0]
all_stats.append({
"idx": idx,
"structure": s,
"n": len(d),
"r": r,
"p_raw": p,
"post_col": post_col,
})
all_pvals.append(p)
# ============================================================
# FDR CORRECTION
# ============================================================
if len(all_pvals) > 0:
reject, p_corr = pg.multicomp(all_pvals, method="fdr_bh")
else:
reject, p_corr = [], []
for i in range(len(all_stats)):
all_stats[i]["p_fdr"] = p_corr[i]
all_stats[i]["reject"] = reject[i]
# ============================================================
# PLOT
# ============================================================
for stat in all_stats:
idx = stat["idx"]
s = stat["structure"]
r = stat["r"]
p_fdr = stat["p_fdr"]
post_col = stat["post_col"]
ax = axes[idx]
d = gdf[[thal_asyn_col, post_col, "AgeatDeath", "Sex", "PMI", "Education"]].dropna()
sns.regplot(
data=d,
x=thal_asyn_col,
y=post_col,
scatter_kws=dict(alpha=0.7, s=45),
line_kws=dict(color=palette[idx % len(palette)], lw=2),
color=palette[idx % len(palette)],
ax=ax
)
sig = "***" if p_fdr < 0.001 else "**" if p_fdr < 0.01 else "*" if p_fdr < 0.05 else ""
title_name = s.replace("_", " ").title()
ax.set_title(f"{title_name} (ρ={r:.2f}, p={stat['p_raw']:.3f}{sig})", fontsize=14, fontweight="bold")
ax.set_xlabel("Thalamic α-synuclein", fontsize=13)
ax.set_ylabel("Normalized volume", fontsize=13)
ax.grid(True, linestyle=":", alpha=0.5)
ax.margins(x=0.12)
# ============================================================
# TURN OFF UNUSED PANELS
# ============================================================
for j in range(len(structures), 9):
axes[j].axis("off")
# ============================================================
# GLOBAL TITLE
# ============================================================
fig.suptitle(
"LBD only: thalamic α-synuclein vs subcortical and limbic volumes",
fontsize=20,
y=0.97
)
plt.savefig(
"LBD_thalamic_aSyn_vs_all_structures_3x3_with_accumbens.png",
dpi=600,
bbox_inches="tight"
)
plt.show()
# %%
"""
Fits a linear mixed-effects model testing associations between regional pathology markers and normalized postmortem subcortical volume across caudate, putamen, thalamus, and pallidum.
The script builds a long-format dataframe with one row per subject–structure pair, z-scores the normalized volume outcome, includes Tau, TDP-43, and α-synuclein pathology scores plus covariates as fixed effects, and models participant-level random intercepts.
"""
import pandas as pd
import numpy as np
import warnings
import statsmodels.formula.api as smf
warnings.filterwarnings("ignore")
# ============================================================
# INPUT DATA
# ============================================================
df_use = df.copy()
# Required columns check
SUBJECT_COL = "INDDID"
needed_base = [SUBJECT_COL, "AgeatDeath", "Sex", "PMI", "Education"]
missing = [c for c in needed_base if c not in df_use.columns]
if missing:
raise ValueError(f"Missing required columns: {missing}")
# Clean / encode
df_use["Sex"] = df_use["Sex"].astype("category").cat.codes
# ============================================================
# ONLY THESE 4 STRUCTURES
# ============================================================
structures = ["caudate", "putamen", "thalamus", "pallidum"]
# region/prefix map (your convention)
region_map = {
"caudate": "CP",
"putamen": "CP",
"thalamus": "TS",
"pallidum": "GP",
}
# ============================================================
# PRINT PATHOLOGY COLUMNS ENDING WITH "Tau"
# (case-insensitive; also strips whitespace)
# ============================================================
tau_cols = [c for c in df_use.columns if str(c).strip().lower().endswith("tau")]
print("\n===== Columns ending with 'Tau' (case-insensitive) =====")
for c in tau_cols:
print(c)
# ============================================================
# NORMALIZE POSTMORTEM VOLUMES (if not already)
# expects: postmortem_<structure> and antemortem_icv
# ============================================================
if "antemortem_icv" not in df_use.columns:
raise ValueError("Missing 'antemortem_icv' needed for volume normalization.")
for s in structures:
pm_col = f"postmortem_{s}"
if pm_col not in df_use.columns:
raise ValueError(f"Missing volume column: {pm_col}")
df_use[f"{pm_col}_norm"] = df_use[pm_col] / df_use["antemortem_icv"]
# ============================================================
# Helper: pick pathology column for a structure + marker
# Markers: "Tau", "TDP43", "aSyn"
# Enforce CP split: caudate vs putamen by preferring colnames containing "caud" or "put"
# ============================================================
def pick_path_col(columns, prefix, marker, structure=None):
"""
columns: list-like of df column names
prefix: e.g. "CP", "TS", "GP"
marker: "Tau" / "TDP43" / "aSyn"
structure: optional; used only to enforce CP caudate vs putamen separation
"""
cols = list(columns)
prefix_l = prefix.lower()
marker_l = marker.lower()
# 1) candidate columns: start with prefix AND contain marker anywhere
cand = [c for c in cols
if str(c).lower().startswith(prefix_l) and marker_l in str(c).lower()]
if not cand:
return None
# 2) enforce CP split for caudate vs putamen
if prefix_l == "cp" and structure is not None:
s_l = structure.lower()
if s_l == "caudate":
better = [c for c in cand if ("caud" in str(c).lower()) or ("caudate" in str(c).lower())]
if better:
cand = better
if s_l == "putamen":
better = [c for c in cand if ("put" in str(c).lower()) or ("putamen" in str(c).lower())]
if better:
cand = better
# 3) prefer columns that END with marker (e.g., "...Tau") if available
end_match = [c for c in cand if str(c).strip().lower().endswith(marker_l)]
if end_match:
cand = end_match
# 4) deterministic pick (alphabetical)
cand = sorted(cand)
return cand[0]
# ============================================================
# BUILD LONG DATAFRAME
# Each row = (INDDID, structure)
# Columns: VolNorm, Tau, TDP43, aSyn + covariates
# ============================================================
rows = []
chosen_cols = [] # debug table
for s in structures:
prefix = region_map[s]
vol_col = f"postmortem_{s}_norm"
tau_col = pick_path_col(df_use.columns, prefix, "Tau", structure=s)
tdp_col = pick_path_col(df_use.columns, prefix, "TDP43", structure=s)
asyn_col = pick_path_col(df_use.columns, prefix, "aSyn", structure=s)
chosen_cols.append({
"Structure": s,
"Prefix": prefix,
"VolCol": vol_col,
"TauCol": tau_col,
"TDP43Col": tdp_col,
"aSynCol": asyn_col
})
# If any pathology marker is missing for this structure, we still build rows
# but those rows will be dropped later if predictors are NaN.
use_cols = [SUBJECT_COL, vol_col, "AgeatDeath", "Sex", "PMI", "Education"]
rename_map = {SUBJECT_COL: "Subject", vol_col: "VolNorm"}
if tau_col is not None:
use_cols.append(tau_col)
rename_map[tau_col] = "Tau"
else:
# create placeholder column later
pass
if tdp_col is not None:
use_cols.append(tdp_col)
rename_map[tdp_col] = "TDP43"
if asyn_col is not None:
use_cols.append(asyn_col)
rename_map[asyn_col] = "aSyn"
d = df_use[use_cols].copy().rename(columns=rename_map)
d["Structure"] = s # kept for debugging only (NOT in model)
# Ensure missing predictors exist as columns
for pred in ["Tau", "TDP43", "aSyn"]:
if pred not in d.columns:
d[pred] = np.nan
rows.append(d)
long_df = pd.concat(rows, ignore_index=True)
print("\n===== Chosen pathology columns per structure =====")
chosen_cols_df = pd.DataFrame(chosen_cols)
print(chosen_cols_df)
# ============================================================
# DROP MISSING + Z-SCORE OUTCOME ONLY
# ============================================================
# keep only complete cases for model variables
model_cols = ["Subject", "VolNorm", "Tau", "TDP43", "aSyn", "Education", "AgeatDeath", "Sex", "PMI"]
long_df = long_df.dropna(subset=model_cols).copy()
# z-score outcome (global)
long_df["VolNorm_z"] = (long_df["VolNorm"] - long_df["VolNorm"].mean()) / long_df["VolNorm"].std(ddof=0)
# categorical subject
long_df["Subject"] = long_df["Subject"].astype("category")
print("\nN rows:", len(long_df), "| N participants:", long_df["Subject"].nunique())
print(long_df[["Subject","Structure","VolNorm","VolNorm_z","Tau","TDP43","aSyn"]].head())
# ============================================================
# FIT MIXED EFFECTS MODEL
# - No Structure term in formula (as you requested)
# - Random intercept per participant
# ============================================================
formula = "VolNorm_z ~ Tau + TDP43 + aSyn + Education + AgeatDeath + Sex + PMI"
print("\n===== Fitting MixedLM =====")
print("Formula:", formula)
m = smf.mixedlm(
formula,
data=long_df,
groups=long_df["Subject"]
).fit(reml=False, method="lbfgs")
print(m.summary())
# ============================================================
# OPTIONAL: show coefficients with more precision (helps when values look like 0.000)
# ============================================================
coefs = pd.DataFrame({
"coef": m.params,
"se": m.bse,
"z": m.tvalues,
"p": m.pvalues
})
print("\n===== Coefs with more precision =====")
print(coefs.to_string(float_format=lambda x: f"{x:.6g}"))
# %%
"""
Fits disease-wise linear mixed-effects models testing whether the primary pathology burden for each diagnostic group is associated with ICV-normalized postmortem subcortical volume.
The script builds a long-format dataset across caudate, putamen, thalamus, and pallidum, selects the disease-specific primary pathology marker, z-scores the volume outcome only, adjusts for education, age at death, sex, and PMI, and includes a random intercept for each participant.
"""
import pandas as pd
import numpy as np
import warnings
import statsmodels.formula.api as smf
warnings.filterwarnings("ignore")
# ============================================================
# 0) LOAD
# ============================================================
df = pd.read_csv("")df_use = df.copy()
# ============================================================
# 1) FIND DISEASE GROUP COLUMN (you had KeyError for NPDx1)
# We'll auto-detect from common options.
# ============================================================
CAND_GROUP_COLS = ["NPDx1", "Group", "group", "NPDx", "Dx", "Diagnosis", "diagnosis"]
GROUP_COL = next((c for c in CAND_GROUP_COLS if c in df_use.columns), None)
if GROUP_COL is None:
raise ValueError(
"❌ Could not find a disease group column. "
"Tried: " + ", ".join(CAND_GROUP_COLS) + "\n"
"Available columns (first 50):\n" + str(list(df_use.columns)[:50])
)
# standardize group strings
df_use[GROUP_COL] = df_use[GROUP_COL].astype(str).str.strip().str.lower()
# ============================================================
# 2) REQUIRED COLUMNS
# ============================================================
SUBJECT_COL = "INDDID"
if SUBJECT_COL not in df_use.columns:
raise ValueError(f"❌ Subject column '{SUBJECT_COL}' not found in df columns.")
# sex -> numeric codes
if "Sex" not in df_use.columns:
raise ValueError("❌ Missing required covariate column: Sex")
df_use["Sex"] = df_use["Sex"].astype("category").cat.codes
for cov in ["AgeatDeath", "PMI", "Education"]:
if cov not in df_use.columns:
raise ValueError(f"❌ Missing required covariate column: {cov}")
# ICV needed for normalization
if "antemortem_icv" not in df_use.columns:
raise ValueError("❌ Missing 'antemortem_icv' needed for volume normalization.")
# ============================================================
# 3) STRUCTURES (ONLY these 4)
# ============================================================
structures = ["caudate", "putamen", "thalamus", "pallidum"]
# disease -> primary marker
# (same as your earlier mapping)
primary_marker_by_group = {
"alzheimer's disease": "Tau",
"lewy body disease": "aSyn",
"ftld-tdp": "TDP43",
"tauopathies": "Tau",
}
# region prefixes (your convention)
# NOTE: CP will be split caudate vs putamen using extra logic below
region_prefix = {
"caudate": "CP",
"putamen": "CP",
"thalamus": "TS",
"pallidum": "GP",
}
# ============================================================
# 4) NORMALIZE VOLUMES (RIGHT)
# ============================================================
for s in structures:
pm_col = f"postmortem_{s}"
if pm_col not in df_use.columns:
raise ValueError(f"❌ Missing volume column: {pm_col}")
df_use[f"{pm_col}_norm"] = df_use[pm_col] / df_use["antemortem_icv"]
# ============================================================
# 5) HELPERS: pick pathology column for a given structure + marker
# CP split: enforce caudate vs putamen when possible
# ============================================================
def pick_pathology_col(columns, structure, marker):
"""
Returns the best-matching pathology column name or None.
marker in {"Tau","TDP43","aSyn"} (case-insensitive match).
"""
cols = list(columns)
marker_l = marker.lower()
struct_l = structure.lower()
# candidate pool by prefix
pref = region_prefix[structure].lower()
pref_candidates = [c for c in cols if c.lower().startswith(pref)]
# keep only those with marker in name
marker_candidates = [c for c in pref_candidates if marker_l in c.lower()]
if len(marker_candidates) == 0:
return None
# --- CP split enforcement ---
# try to choose a caudate-specific vs putamen-specific CP column
if structure == "caudate":
# prefer explicit caudate hints
caud_pref = [c for c in marker_candidates if ("caud" in c.lower() or "caudate" in c.lower())]
if len(caud_pref) > 0:
return caud_pref[0]
if structure == "putamen":
put_pref = [c for c in marker_candidates if ("put" in c.lower() or "putamen" in c.lower())]
if len(put_pref) > 0:
return put_pref[0]
# fallback: take the first marker match within prefix
return marker_candidates[0]
# ============================================================
# 6) BUILD LONG DF: one row per (Subject, Structure)
# PrimaryPath depends on disease group
# ============================================================
rows = []
chosen_cols_log = [] # track what got picked
for group_name, marker in primary_marker_by_group.items():
gdf = df_use[df_use[GROUP_COL] == group_name].copy()
if gdf.empty:
continue
for s in structures:
vol_col = f"postmortem_{s}_norm"
path_col = pick_pathology_col(gdf.columns, structure=s, marker=marker)
chosen_cols_log.append((group_name, s, marker, path_col, vol_col))
if path_col is None:
continue
d = gdf[[SUBJECT_COL, path_col, vol_col, "AgeatDeath", "Sex", "PMI", "Education"]].copy()
d = d.rename(columns={
SUBJECT_COL: "Subject",
path_col: "PrimaryPath",
vol_col: "VolNorm",
})
d["Disease"] = group_name # keep disease label
d["Structure"] = s # kept for debugging (NOT in formula)
rows.append(d)
if len(rows) == 0:
raise ValueError("❌ No long-format rows created. Check pathology column naming / prefixes / markers.")
mB_df = pd.concat(rows, ignore_index=True)
# drop NAs & reset index to avoid statsmodels IndexError
need = ["Subject", "PrimaryPath", "VolNorm", "AgeatDeath", "Sex", "PMI", "Education", "Disease"]
mB_df = mB_df.dropna(subset=need).reset_index(drop=True)
# Z-score OUTCOME only (across the whole dataset used for each disease model)
# (You asked: don't z-score predictors)
mB_df["VolNorm_z"] = (mB_df["VolNorm"] - mB_df["VolNorm"].mean()) / mB_df["VolNorm"].std(ddof=0)
# cast subject categorical
mB_df["Subject"] = mB_df["Subject"].astype("category")
# ============================================================
# 7) DEBUG OUTPUTS YOU ASKED FOR
# ============================================================
print("\n===== GROUP COLUMN USED =====")
print("GROUP_COL =", GROUP_COL)
print("\n===== Chosen pathology columns (first 30 rows) =====")
log_df = pd.DataFrame(chosen_cols_log, columns=["Disease", "Structure", "Marker", "ChosenPathCol", "VolCol"])
print(log_df.head(30))
print("\n===== PrimaryPath value counts (including NA) =====")
print(mB_df["PrimaryPath"].value_counts(dropna=False))
print("\n===== Preview mB_df =====")
print(mB_df.head())
# ============================================================
# 8) FIT 4 DISEASE-WISE LME MODELS (PRIMARY PATH ONLY)
# No Structure term in formula (as requested)
# ============================================================
formula = "VolNorm_z ~ PrimaryPath + Education + AgeatDeath + Sex + PMI"
print("\n===== Formula =====")
print(formula)
models = {}
for disease in primary_marker_by_group.keys():
ddf = mB_df[mB_df["Disease"] == disease].copy()
# sanity checks
n_sub = ddf["Subject"].nunique()
n_rows = len(ddf)
print(f"\n------------------------------\nDISEASE: {disease}\nrows={n_rows} | subjects={n_sub}")
if n_rows < 20 or n_sub < 5:
print("⚠️ Skipping (too few rows or subjects).")
continue
# Important: reset index again per subset (prevents out-of-bounds indexing)
ddf = ddf.reset_index(drop=True)
# Fit random-intercept model
fit = smf.mixedlm(
formula=formula,
data=ddf,
groups=ddf["Subject"]
).fit(reml=False, method="lbfgs")
models[disease] = fit
print(fit.summary())
# ============================================================
# 9) (OPTIONAL) QUICK: show two random subjects' rows (debug)
# ============================================================
rng = np.random.default_rng(0)
subj_two = rng.choice(mB_df["Subject"].astype(str).unique(), size=2, replace=False)
print("\n===== Two random subjects =====")
print("Subjects:", subj_two.tolist())
print(mB_df[mB_df["Subject"].astype(str).isin(subj_two)].sort_values(["Subject", "Structure"]))
# %%
"""
Fits disease-wise polypathology linear mixed-effects models testing whether Tau, TDP-43, and α-synuclein pathology are associated with ICV-normalized postmortem subcortical volume.
The script builds a long-format dataset across caudate, putamen, thalamus, and pallidum, z-scores the normalized volume outcome only, adjusts for age at death, sex, PMI, and education, and includes participant-level random intercepts while fitting separate models within each diagnostic group.
"""
import pandas as pd
import numpy as np
import re
import warnings
import statsmodels.formula.api as smf
warnings.filterwarnings("ignore")
# ============================================================
# 0) LOAD DATA
# ============================================================
df = pd.read_csv("")print("Loaded:", df.shape)
print(df.head())
# ============================================================
# 1) CONFIG
# ============================================================
SUBJECT_COL = "INDDID" # confirmed by you
DX_COL = "NPDx1" # if your file doesn't have NPDx1, see fallback below
# 4 subcortical structures you want
structures = ["caudate", "putamen", "thalamus", "pallidum"]
# region prefixes you were using
# NOTE: caudate/putamen share CP prefix in your data
region_prefix = {
"caudate": "CP",
"putamen": "CP",
"thalamus": "TS",
"pallidum": "GP",
}
# marker suffixes in your sheet
markers = ["Tau", "TDP43", "aSyn"]
# outcome columns (normalized volumes)
VOL_BASE = {s: f"postmortem_{s}" for s in structures}
ICV_COL = "antemortem_icv"
# covariates
COVARS = ["AgeatDeath", "Sex", "PMI", "Education"]
# diseases of interest (edit these to match your values)
disease_list = [
"alzheimer's disease",
"lewy body disease",
"ftld-tdp",
"tauopathies",
]
# ============================================================
# 2) SAFETY / CLEANING
# ============================================================
df_use = df.copy()
# --- Subject ID must exist ---
if SUBJECT_COL not in df_use.columns:
raise ValueError(f"❌ Missing subject column '{SUBJECT_COL}' in dataframe.")
# --- Dx column fallback (if NPDx1 not present) ---
if DX_COL not in df_use.columns:
# If your file uses a different column name, set DX_COL above.
# Otherwise, try a few common alternatives:
for alt in ["Group", "NPDx", "Dx", "Diagnosis", "NPDx1_clean"]:
if alt in df_use.columns:
DX_COL = alt
print(f"⚠️ Using '{DX_COL}' as diagnosis column (fallback).")
break
else:
raise ValueError("❌ Could not find a diagnosis column. Set DX_COL to the correct column name.")
# Normalize dx strings
df_use[DX_COL] = df_use[DX_COL].astype(str).str.strip().str.lower()
# Sex: make numeric scalar codes (avoid patsy “>1-dimensional” issues)
# If Sex is already numeric 0/1, this will keep it numeric.
if df_use["Sex"].dtype.name == "category" or df_use["Sex"].dtype == object:
df_use["Sex"] = df_use["Sex"].astype("category").cat.codes
# Force numeric covariates (this prevents PatsyError if any are weird objects/arrays)
for c in ["AgeatDeath", "PMI", "Education"]:
if c in df_use.columns:
df_use[c] = pd.to_numeric(df_use[c], errors="coerce")
# ============================================================
# 3) PRINT “Tau columns” (names ending with Tau)
# ============================================================
tau_cols = [c for c in df_use.columns if re.search(r"tau\s*$", str(c), flags=re.IGNORECASE)]
print("\n===== Columns ending with 'Tau' =====")
print(tau_cols[:200])
print("Count:", len(tau_cols))
# ============================================================
# 4) NORMALIZE VOLUMES and BUILD LONG DF
# ============================================================
# create normalized volume columns
for s in structures:
vol_col = VOL_BASE[s]
if vol_col in df_use.columns and ICV_COL in df_use.columns:
df_use[f"{vol_col}_norm"] = df_use[vol_col] / df_use[ICV_COL]
else:
print(f"⚠️ Missing {vol_col} or {ICV_COL}; cannot normalize for {s}.")
def choose_path_col(columns, prefix, marker, structure):
"""
Choose pathology column for a given (prefix, marker, structure).
- primary match: startswith(prefix) AND endswith(marker)
- CP special case: prefer columns containing 'caud' for caudate and 'put' for putamen
"""
cols = list(columns)
# candidate columns by prefix
pref = [c for c in cols if str(c).lower().startswith(prefix.lower())]
if not pref:
return None
# keep those ending with marker (Tau / TDP43 / aSyn)
# (ending-with constraint as you requested)
mk = [c for c in pref if re.search(rf"{re.escape(marker)}\s*$", str(c), flags=re.IGNORECASE)]
if not mk:
# if none strictly end with marker, loosen slightly (marker anywhere)
mk = [c for c in pref if marker.lower() in str(c).lower()]
if not mk:
return None
# enforce CP split caudate vs putamen
if prefix.lower() == "cp":
if structure.lower() == "caudate":
# prefer columns mentioning caudate-ish strings
for key in ["caud", "caudate", "head"]:
hit = [c for c in mk if key in str(c).lower()]
if hit:
return hit[0]
if structure.lower() == "putamen":
for key in ["put", "putamen"]:
hit = [c for c in mk if key in str(c).lower()]
if hit:
return hit[0]
# fallback: if no structure-specific keyword, still return first match
return mk[0]
# non-CP: just first best match
return mk[0]
rows = []
debug_choices = []
for s in structures:
vol_norm = f"{VOL_BASE[s]}_norm"
if vol_norm not in df_use.columns:
continue
prefix = region_prefix[s]
# choose columns for each marker
col_tau = choose_path_col(df_use.columns, prefix, "Tau", s)
col_tdp = choose_path_col(df_use.columns, prefix, "TDP43", s)
col_asyn = choose_path_col(df_use.columns, prefix, "aSyn", s)
debug_choices.append((s, prefix, col_tau, col_tdp, col_asyn, vol_norm))
keep = [SUBJECT_COL, DX_COL, vol_norm] + COVARS
# add pathology cols if exist
for cc in [col_tau, col_tdp, col_asyn]:
if cc is not None:
keep.append(cc)
d = df_use[keep].copy()
# rename to standard names
rename_map = {SUBJECT_COL: "Subject", DX_COL: "Dx", vol_norm: "VolNorm"}
if col_tau is not None: rename_map[col_tau] = "Tau"
if col_tdp is not None: rename_map[col_tdp] = "TDP43"
if col_asyn is not None: rename_map[col_asyn] = "aSyn"
d = d.rename(columns=rename_map)
d["Structure"] = s
# ensure pathology columns exist (if missing, create with NaN so formula is stable)
for m in ["Tau", "TDP43", "aSyn"]:
if m not in d.columns:
d[m] = np.nan
rows.append(d)
long_df = pd.concat(rows, ignore_index=True)
print("\n===== Column choices per structure =====")
debug_df = pd.DataFrame(debug_choices, columns=["Structure", "Prefix", "Tau_col", "TDP43_col", "aSyn_col", "VolNorm_col"])
print(debug_df)
# ============================================================
# 5) CLEAN LONG DF + Z-SCORE OUTCOME ONLY (as you requested)
# ============================================================
# keep only diseases we care about (optional but recommended)
long_df["Dx"] = long_df["Dx"].astype(str).str.strip().str.lower()
long_df = long_df[long_df["Dx"].isin(disease_list)].copy()
# force numeric predictors (avoid PatsyError)
for m in ["Tau", "TDP43", "aSyn"]:
long_df[m] = pd.to_numeric(long_df[m], errors="coerce")
for c in ["AgeatDeath", "PMI", "Education", "VolNorm"]:
long_df[c] = pd.to_numeric(long_df[c], errors="coerce")
# drop rows missing essential fields (but allow missing pathology in some rows; you can tighten later)
essential = ["Subject", "Dx", "VolNorm"] + COVARS
long_df = long_df.dropna(subset=essential).copy()
# z-score OUTCOME only (within entire stacked dataset)
long_df["VolNorm_z"] = (long_df["VolNorm"] - long_df["VolNorm"].mean()) / long_df["VolNorm"].std(ddof=0)
# IMPORTANT for statsmodels MixedLM: reset index AFTER drops to avoid your IndexError
long_df = long_df.reset_index(drop=True)
print("\n===== LONG DF summary =====")
print("Rows:", len(long_df), "| Subjects:", long_df["Subject"].nunique())
print(long_df[["Subject", "Dx", "Structure", "VolNorm", "VolNorm_z", "Tau", "TDP43", "aSyn"]].head())
# ============================================================
# 6) QUICK DEBUG HELPERS (you asked for these)
# ============================================================
# Show any two random subjects (raw rows)
two_subj = np.random.choice(long_df["Subject"].unique(), size=min(2, long_df["Subject"].nunique()), replace=False)
print("\n===== Two random subjects =====")
print(two_subj)
print(long_df[long_df["Subject"].isin(two_subj)][
["Subject", "Dx", "Structure", "VolNorm", "VolNorm_z", "Tau", "TDP43", "aSyn", "AgeatDeath", "Sex", "PMI", "Education"]
].sort_values(["Subject", "Structure"]))
# Check PrimaryPath-style column completeness (here: Tau/TDP43/aSyn)
print("\n===== Value counts (including NaN) for Tau / TDP43 / aSyn =====")
print("Tau:")
print(long_df["Tau"].isna().value_counts(dropna=False))
print("TDP43:")
print(long_df["TDP43"].isna().value_counts(dropna=False))
print("aSyn:")
print(long_df["aSyn"].isna().value_counts(dropna=False))
# ============================================================
# 7) FIT POLYPATHOLOGY LME PER DISEASE (4 models)
# - No Structure term in formula (you explicitly asked)
# - Random intercept per Subject
# ============================================================
# NOTE: If some diseases have pathology columns totally missing (all NaN),
# MixedLM will fail. We handle this by dropping predictors that are all-NaN within each disease.
base_cov = "AgeatDeath + Sex + PMI + Education"
base_paths = ["Tau", "TDP43", "aSyn"]
results = {}
for dx in disease_list:
dxd = long_df[long_df["Dx"] == dx].copy()
# pick which pathology predictors are usable in this disease subset
usable = [p for p in base_paths if dxd[p].notna().sum() > 5] # require at least 5 non-NaN points
if len(usable) == 0:
print(f"\n⚠️ Skipping {dx}: no pathology predictors have enough non-missing values.")
continue
# drop rows with missing in *usable* predictors (model needs complete cases)
model_df = dxd.dropna(subset=usable + COVARS + ["VolNorm_z", "Subject"]).copy()
# also reset index to avoid IndexError in statsmodels
model_df = model_df.reset_index(drop=True)
if model_df["Subject"].nunique() < 3 or len(model_df) < 20:
print(f"\n⚠️ Skipping {dx}: too few subjects/rows after filtering. "
f"Rows={len(model_df)}, Subjects={model_df['Subject'].nunique()}")
continue
# build formula (NO Structure)
path_term = " + ".join(usable)
formula = f"VolNorm_z ~ {path_term} + {base_cov}"
print("\n" + "="*90)
print(f"Fitting polypathology LME for: {dx}")
print("Rows:", len(model_df), "| Subjects:", model_df["Subject"].nunique())
print("Formula:", formula)
# Fit
try:
m = smf.mixedlm(
formula=formula,
data=model_df,
groups=model_df["Subject"]
).fit(reml=False, method="lbfgs")
print(m.summary())
results[dx] = m
except Exception as e:
print(f"❌ Model failed for {dx}: {e}")
print("\nDone. Models fit:", list(results.keys()))
##########################################################################################
# %%
##########################################################################################
########## Addtional markers
##########################################################################################
# %%
"""
Generates separate publication-style figures for partial Spearman correlations between global pathology severity markers and ICV-normalized postmortem limbic/subcortical volumes.
For each pathology marker, diagnostic group, and structure, the script adjusts for age at death, sex, PMI, and education, applies FDR correction across structures within each group, annotates significance with superscript asterisks, and saves both PNG and PDF outputs.
"""
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib
import pingouin as pg
from matplotlib.backends.backend_pdf import PdfPages
matplotlib.rcParams['font.family'] = 'DejaVu Sans'
sns.set(style="whitegrid")
np.seterr(divide='ignore', invalid='ignore')
##########################################################################################
# DATA PREP
##########################################################################################
df_use = df.copy()
df_use["NPDx1"] = df_use["NPDx1"].astype(str).str.strip().str.lower()
df_use["Sex"] = df_use["Sex"].astype("category").cat.codes
##########################################################################################
# STRUCTURES + NORMALIZATION
##########################################################################################
structures = ["hippocampus", "amygdala", "caudate",
"putamen", "thalamus", "pallidum", "accumbens_area"]
pretty_structure = {
"hippocampus": "Hippocampus",
"amygdala": "Amygdala",
"caudate": "Caudate",
"putamen": "Putamen",
"thalamus": "Thalamus",
"pallidum": "Pallidum",
"accumbens_area": "Accumbens area"
}
for s in structures:
pm = f"postmortem_{s}"
if pm in df_use.columns:
df_use[f"{pm}_norm"] = df_use[pm] / df_use["antemortem_icv"]
##########################################################################################
# PATHOLOGY MARKERS
##########################################################################################
markers = ["ABeta", "CERAD", "Braak06", "CAA",
"Arteriolosclerosis", "Atherosclerosis"]
pretty_marker = {
"ABeta": "Aβ (Amyloid-β)",
"CERAD": "CERAD plaques",
"Braak06": "Braak stage",
"CAA": "CAA",
"Arteriolosclerosis": "Arteriolosclerosis",
"Atherosclerosis": "Atherosclerosis"
}
severity_map = {
"0": 0, "None": 0, "None/Normal": 0,
"1": 1, "Mild": 1,
"2": 2, "Moderate": 2,
"3": 3, "Severe": 3
}
for m in markers:
if m in df_use.columns:
df_use[m] = (
df_use[m]
.astype(str)
.replace("Unknown", np.nan)
.replace(severity_map)
.replace("nan", np.nan)
)
df_use[m] = pd.to_numeric(df_use[m], errors="coerce")
##########################################################################################
# GROUPS + COLORS
##########################################################################################
groups = ["alzheimer's disease", "lewy body disease", "ftld-tdp", "tauopathies"]
pretty_group = {
"alzheimer's disease": "AD",
"lewy body disease": "LBD",
"ftld-tdp": "FTLD-TDP",
"tauopathies": "Tauopathies"
}
group_color = {
"alzheimer's disease": "#4C72B0",
"lewy body disease": "#DD8452",
"ftld-tdp": "#55A868",
"tauopathies": "#CBAF00"
}
##########################################################################################
# SIGNIFICANCE → SUPERSCRIPT ASTERISKS
##########################################################################################
def p_to_superscript(p_fdr):
if p_fdr < 0.001:
return "$^{***}$"
elif p_fdr < 0.01:
return "$^{**}$"
elif p_fdr < 0.05:
return "$^{*}$"
else:
return ""
##########################################################################################
# MAIN LOOP — ONE FIGURE PER MARKER
##########################################################################################
for marker in markers:
fig, axes = plt.subplots(
len(structures), len(groups),
figsize=(4.2 * len(groups), 3.0 * len(structures)),
sharex=False, sharey=False
)
# Compute for each disease group
for col, g in enumerate(groups):
gdf = df_use[df_use["NPDx1"] == g]
# FIRST PASS: collect all raw p-values for FDR correction
raw_pvals = []
for s in structures:
post_col = f"postmortem_{s}_norm"
if marker not in gdf.columns or post_col not in gdf.columns:
raw_pvals.append(np.nan)
continue
d = gdf[[marker, post_col, "AgeatDeath", "Sex", "PMI", "Education"]]
d = d.apply(pd.to_numeric, errors="coerce").dropna()
if len(d) < 5:
raw_pvals.append(np.nan)
continue
res = pg.partial_corr(
data=d, x=marker, y=post_col,
covar=["AgeatDeath", "Sex", "PMI", "Education"],
method="spearman"
)
raw_pvals.append(res["p-val"].iloc[0])
# FDR
valid_mask = ~pd.isna(raw_pvals)
valid_p = np.array(raw_pvals)[valid_mask]
if len(valid_p) > 0:
_, p_fdr_valid = pg.multicomp(valid_p, method="fdr_bh")
else:
p_fdr_valid = []
# Expand back
p_fdr_full = []
idx = 0
for v in valid_mask:
if v:
p_fdr_full.append(p_fdr_valid[idx])
idx += 1
else:
p_fdr_full.append(np.nan)
# SECOND PASS: plotting
for row, s in enumerate(structures):
ax = axes[row, col]
post_col = f"postmortem_{s}_norm"
# Remove default labels
ax.set_ylabel("")
# Set clean y-label only for left-most column
if col == 0:
ax.set_ylabel(pretty_structure[s], fontsize=13)
# Prepare data
if marker not in gdf.columns or post_col not in gdf.columns:
ax.text(0.5, 0.5, "No data", ha="center", va="center")
ax.set_xticks([])
ax.set_yticks([])
continue
d = gdf[[marker, post_col, "AgeatDeath", "Sex", "PMI", "Education"]]
d = d.apply(pd.to_numeric, errors="coerce").dropna()
if len(d) < 5:
ax.text(0.5, 0.5, "N too small", ha="center", va="center")
ax.set_xticks([])
ax.set_yticks([])
continue
# Partial Spearman
res = pg.partial_corr(
data=d, x=marker, y=post_col,
covar=["AgeatDeath", "Sex", "PMI", "Education"],
method="spearman"
)
r = res["r"].iloc[0]
p_raw = res["p-val"].iloc[0]
p_fdr = p_fdr_full[row]
sup = p_to_superscript(p_fdr)
# Scatter + line
sns.regplot(
data=d,
x=marker,
y=post_col,
scatter_kws=dict(color=group_color[g], s=50, alpha=0.75),
line_kws=dict(color=group_color[g], lw=2.3, alpha=0.9),
ax=ax
)
# Reset y-label after regplot overwrites it
ax.set_ylabel("")
if col == 0:
ax.set_ylabel(pretty_structure[s], fontsize=13)
# Pretty x-axis label
if row == len(structures) - 1:
ax.set_xlabel(pretty_marker[marker], fontsize=12)
else:
ax.set_xlabel("")
# Title with superscript asterisks
ax.set_title(
f"{pretty_group[g]} (ρ = {r:.2f}; p = {p_raw:.3f}{sup})",
fontsize=12, pad=6
)
ax.grid(True, linestyle=":", alpha=0.5)
# Figure title
plt.suptitle(
f"Partial Spearman correlations — {pretty_marker[marker]}",
fontsize=18, weight="bold"
)
plt.tight_layout(rect=[0, 0, 1, 0.95])
# SAVE PNG
png_name = f"Postmortem_Partial_Spearman_{marker}.png"
plt.savefig(png_name, dpi=300, bbox_inches="tight")
print(f"SAVED PNG: {png_name}")
# SAVE PDF
pdf_name = f"Postmortem_Partial_Spearman_{marker}.pdf"
with PdfPages(pdf_name) as pdf:
pdf.savefig(fig, bbox_inches="tight")
print(f"SAVED PDF: {pdf_name}")
plt.show()
print("\n✔✔✔ ALL FIGURES SAVED (PNG 300dpi + PDF)")
# %%
"""
Computes partial Spearman correlations between global pathology/vascular severity markers and ICV-normalized postmortem limbic/subcortical volumes across diagnostic groups.
The script adjusts for age at death, sex, PMI, and education, applies FDR correction across structures within each marker–group comparison, and visualizes the resulting correlation coefficients as clean 3×3 heatmaps with significance annotations and a separate colorbar legend.
"""
import pandas as pd, numpy as np, seaborn as sns, matplotlib.pyplot as plt, pingouin as pg, warnings
from matplotlib.transforms import Bbox
import matplotlib
matplotlib.rcParams['font.family'] = 'DejaVu Sans'
warnings.filterwarnings("ignore", category=RuntimeWarning)
np.seterr(divide='ignore', invalid='ignore')
# ---------------------------------------------------------------------
# USE df_use EXACTLY AS IS (NO RELABELING OF GROUPS)
# ---------------------------------------------------------------------
df_use = df.copy()
df_use["NPDx1"] = df_use["NPDx1"].astype(str).str.strip().str.lower()
df_use["Sex"] = df_use["Sex"].astype("category").cat.codes
# ---------------------------------------------------------------------
# Structures
# ---------------------------------------------------------------------
structures = ["hippocampus", "amygdala", "caudate",
"putamen", "thalamus", "pallidum", "accumbens_area"]
# Normalize volumes
for s in structures:
pm = f"postmortem_{s}"
if pm in df_use.columns:
df_use[f"{pm}_norm"] = df_use[pm] / df_use["antemortem_icv"]
# ---------------------------------------------------------------------
# Pathology markers
# ---------------------------------------------------------------------
pathology_markers = [
"ABeta", "Braak06", "CERAD", "CAA",
"Arteriolosclerosis", "Atherosclerosis"
]
severity_map = {
"0": 0, "None": 0, "None/Normal": 0,
"1": 1, "Mild": 1,
"2": 2, "Moderate": 2,
"3": 3, "Severe": 3
}
for m in pathology_markers:
if m in df_use.columns:
df_use[m] = (
df_use[m].astype(str)
.replace("Unknown", np.nan)
.replace(severity_map)
.replace("nan", np.nan)
)
df_use[m] = pd.to_numeric(df_use[m], errors="coerce")
# ---------------------------------------------------------------------
# Settings
# ---------------------------------------------------------------------
groups = ["alzheimer's disease", "lewy body disease", "ftld-tdp", "tauopathies"]
group_labels = ["AD", "LBD", "FTLD-TDP", "Tauopathies"]
covars = ["AgeatDeath", "Sex", "PMI", "Education"]
sns.set(style="white", context="talk")
# ---------------------------------------------------------------------
# Compute partial correlations
# ---------------------------------------------------------------------
results_all = []
for marker in pathology_markers:
for g in groups:
sub = df_use[df_use["NPDx1"] == g].copy()
if sub.empty:
continue
all_rho, all_p = [], []
for s in structures:
post_col = f"postmortem_{s}_norm"
if post_col not in sub.columns:
all_rho.append(np.nan)
all_p.append(np.nan)
continue
cols = [marker, post_col] + covars
d = sub[cols].apply(pd.to_numeric, errors="coerce").dropna()
if len(d) < 5:
all_rho.append(np.nan)
all_p.append(np.nan)
continue
res = pg.partial_corr(
data=d, x=marker, y=post_col,
covar=covars, method="spearman"
)
rho = res["r"].iloc[0]
p = res["p-val"].iloc[0]
all_rho.append(rho)
all_p.append(p)
valid = ~pd.isna(all_p)
pvals = np.array(all_p)[valid]
if len(pvals) > 0:
_, p_corr = pg.multicomp(pvals, method="fdr_bh")
else:
p_corr = []
idx = 0
for s, rho, p in zip(structures, all_rho, all_p):
if not pd.isna(p):
pFDR = p_corr[idx]
idx += 1
else:
pFDR = np.nan
results_all.append({
"Marker": marker,
"Group": g,
"Structure": s,
"rho": rho,
"p_val": p,
"p_FDR": pFDR
})
results_master_df = pd.DataFrame(results_all)
print("\n✅ Partial Spearman correlations computed + FDR corrected.\n")
# ---------------------------------------------------------------------
# Plotting: 3×3 heatmaps
# ---------------------------------------------------------------------
sns.set(style="white", context="talk")
cmap = sns.diverging_palette(280, 150, s=90, l=60, as_cmap=True)
vmin, vmax = -0.6, 0.6
def fdr_star(p):
if p < 0.001: return "***"
elif p < 0.01: return "**"
elif p < 0.05: return "*"
return ""
titles = {
"ABeta": "Aβ (Amyloid-β)",
"Braak06": "Braak stage",
"CERAD": "CERAD plaques",
"CAA": "CAA",
"Arteriolosclerosis": "Arteriolosclerosis",
"Atherosclerosis": "Atherosclerosis"
}
structure_order = structures
fig, axes = plt.subplots(3, 3, figsize=(18, 18))
axes = axes.flatten()
for i, marker in enumerate(pathology_markers):
ax = axes[i]
sub = results_master_df[results_master_df["Marker"] == marker]
if sub.empty:
ax.axis("off")
continue
rho_mat = sub.pivot(index="Structure", columns="Group", values="rho") \
.reindex(index=structure_order, columns=groups)
p_mat = sub.pivot(index="Structure", columns="Group", values="p_val") \
.reindex(index=structure_order, columns=groups)
pFDR_mat = sub.pivot(index="Structure", columns="Group", values="p_FDR") \
.reindex(index=structure_order, columns=groups)
annot = rho_mat.copy().astype(str)
for r in rho_mat.index:
for c in rho_mat.columns:
rv = rho_mat.loc[r, c]
pv = p_mat.loc[r, c]
pf = pFDR_mat.loc[r, c]
if not np.isnan(rv):
annot.loc[r, c] = f"{rv:.2f}{fdr_star(pf)}\n(p={pv:.3f})"
else:
annot.loc[r, c] = ""
sns.heatmap(
rho_mat, annot=annot, fmt="", cmap=cmap,
vmin=vmin, vmax=vmax, cbar=False,
linewidths=0.5, linecolor="white",
ax=ax, annot_kws={"fontsize": 9, "color": "black"}
)
ax.set_title(titles[marker], fontsize=13, color="black", pad=6)
ax.set_xticklabels(group_labels, fontsize=10, rotation=0, color="black")
ax.set_xlabel("")
ax.set_ylabel("")
if i % 3 == 0:
ax.set_yticklabels(
[s.capitalize().replace("_area", "") for s in structure_order],
fontsize=10, color="black"
)
else:
ax.set_yticklabels([])
# Turn off unused plots
for j in range(len(pathology_markers), 9):
axes[j].axis("off")
# ---------------------------------------------------------------------
# Figure title
# ---------------------------------------------------------------------
plt.suptitle(
"Partial Spearman correlation between postmortem MRI volumes\n"
"and global markers of pathology, degeneration and vascular burden",
fontsize=16, color="black", y=0.96
)
plt.tight_layout(rect=[0, 0, 1, 0.94])
plt.savefig("heatmap_addiotnal_markers.png", dpi=600, bbox_inches="tight")
plt.show()
##########################################################################################
# Standalone colorbar (exact matching)
##########################################################################################
fig, ax = plt.subplots(figsize=(1.2, 5))
sm = plt.cm.ScalarMappable(
cmap=cmap,
norm=plt.Normalize(vmin=vmin, vmax=vmax)
)
sm.set_array([])
cbar = plt.colorbar(sm, cax=ax)
cbar.set_label("Spearman ρ", fontsize=12, color="black", labelpad=6)
cbar.ax.tick_params(labelsize=10, colors="black", length=3)
for spine in ax.spines.values():
spine.set_visible(False)
plt.tight_layout()
plt.savefig("heatmap_addiotnal_markers_legend.png", dpi=600, bbox_inches="tight")
plt.show()
# %%
"""
Generates standardized OLS β heatmaps for vascular pathology markers across four diagnostic groups and postmortem limbic/subcortical volumes.
For each group, marker, and structure, the script z-scores the ICV-normalized volume outcome, pathology marker, and covariates, fits covariate-adjusted OLS models, applies BH-FDR correction, and visualizes standardized β values with raw p-values and FDR-based significance stars.
"""
# ---------------------------------------------------------------------
# GLOBAL SETTINGS
# ---------------------------------------------------------------------
matplotlib.rcParams["font.family"] = "DejaVu Sans"
warnings.filterwarnings("ignore", category=RuntimeWarning)
np.seterr(divide="ignore", invalid="ignore")
sns.set(style="white", context="talk")
# ---------------------------------------------------------------------
# USE df AS INPUT
# ---------------------------------------------------------------------
df_use = df.copy()
df_use["NPDx1"] = df_use["NPDx1"].astype(str).str.strip().str.lower()
df_use["Sex"] = df_use["Sex"].astype("category").cat.codes
# ---------------------------------------------------------------------
# Groups (4)
# ---------------------------------------------------------------------
groups = [
"alzheimer's disease",
"lewy body disease",
"ftld-tdp",
"tauopathies"
]
pretty_group = {
"alzheimer's disease": "AD",
"lewy body disease": "LBD",
"ftld-tdp": "FTLD-TDP",
"tauopathies": "FTLD-Tau"
}
group_labels = [pretty_group[g] for g in groups]
# ---------------------------------------------------------------------
# Structures + normalization
# ---------------------------------------------------------------------
structures = [
"hippocampus", "amygdala", "caudate",
"putamen", "thalamus", "pallidum", "accumbens_area"
]
pretty_y = {
"hippocampus": "Hippocampus",
"amygdala": "Amygdala",
"caudate": "Caudate",
"putamen": "Putamen",
"thalamus": "Thalamus",
"pallidum": "Pallidum",
"accumbens_area": "Accumbens"
}
for s in structures:
pm = f"postmortem_{s}"
if pm in df_use.columns and "antemortem_icv" in df_use.columns:
df_use[f"{pm}_norm"] = df_use[pm] / df_use["antemortem_icv"]
# ---------------------------------------------------------------------
# Markers (3)
# ---------------------------------------------------------------------
pathology_markers = ["CAA", "Arteriolosclerosis", "Atherosclerosis"]
titles = {
"CAA": "CAA",
"Arteriolosclerosis": "Arteriolosclerosis",
"Atherosclerosis": "Atherosclerosis"
}
# ---------------------------------------------------------------------
# Convert vascular severity text to numeric
# ---------------------------------------------------------------------
severity_map = {
"0": 0, "None": 0, "None/Normal": 0,
"1": 1, "Mild": 1,
"2": 2, "Moderate": 2,
"3": 3, "Severe": 3
}
for m in pathology_markers:
if m in df_use.columns:
df_use[m] = (
df_use[m].astype(str)
.replace("Unknown", np.nan)
.replace(severity_map)
.replace("nan", np.nan)
)
df_use[m] = pd.to_numeric(df_use[m], errors="coerce")
# ---------------------------------------------------------------------
# Covariates
# ---------------------------------------------------------------------
covars = ["AgeatDeath", "Sex", "PMI", "Education"]
# ---------------------------------------------------------------------
# Helpers
# ---------------------------------------------------------------------
def fdr_star(p):
if pd.isna(p):
return ""
if p < 0.001:
return "***"
if p < 0.01:
return "**"
if p < 0.05:
return "*"
return ""
def zscore_inplace(df_in, cols):
out = df_in.copy()
for c in cols:
sd = out[c].std(ddof=0)
if (sd == 0) or np.isnan(sd):
out[c] = np.nan
else:
out[c] = (out[c] - out[c].mean()) / sd
return out
# ---------------------------------------------------------------------
# Compute OLS standardized β for each (marker, group, structure)
# Standardization happens WITHIN EACH GROUP for:
# y, marker, and covariates
# ---------------------------------------------------------------------
rows = []
MIN_N = 8 # can change if needed
for marker in pathology_markers:
for g in groups:
sub = df_use[df_use["NPDx1"] == g].copy()
if sub.empty:
for s in structures:
rows.append({
"Marker": marker,
"Group": g,
"Structure": s,
"beta_std": np.nan,
"p_raw": np.nan,
"N": 0
})
continue
for s in structures:
y = f"postmortem_{s}_norm"
if (marker not in sub.columns) or (y not in sub.columns):
rows.append({
"Marker": marker,
"Group": g,
"Structure": s,
"beta_std": np.nan,
"p_raw": np.nan,
"N": 0
})
continue
d = sub[[marker, y] + covars].apply(pd.to_numeric, errors="coerce").dropna()
if len(d) < MIN_N:
rows.append({
"Marker": marker,
"Group": g,
"Structure": s,
"beta_std": np.nan,
"p_raw": np.nan,
"N": int(len(d))
})
continue
# z-score within group
dz = zscore_inplace(d, [marker, y] + covars).dropna()
if len(dz) < MIN_N:
rows.append({
"Marker": marker,
"Group": g,
"Structure": s,
"beta_std": np.nan,
"p_raw": np.nan,
"N": int(len(dz))
})
continue
# OLS: z(y) ~ z(marker) + z(covariates)
X = sm.add_constant(dz[[marker] + covars], has_constant="add")
fit = sm.OLS(dz[y], X).fit()
rows.append({
"Marker": marker,
"Group": g,
"Structure": s,
"beta_std": float(fit.params.get(marker, np.nan)),
"p_raw": float(fit.pvalues.get(marker, np.nan)),
"N": int(len(dz))
})
ols_master = pd.DataFrame(rows)
# ---------------------------------------------------------------------
# FDR correction:
# For each marker, correct across ALL (structure × group) tests
# Stars come from p_FDR, while raw p is displayed
# ---------------------------------------------------------------------
ols_master["p_FDR"] = np.nan
for marker in pathology_markers:
idx = ols_master["Marker"] == marker
pvals = ols_master.loc[idx, "p_raw"].values
ok = ~np.isnan(pvals)
if ok.sum() > 0:
_, p_corr = pg.multicomp(pvals[ok], method="fdr_bh")
out = np.full_like(pvals, np.nan, dtype=float)
out[ok] = p_corr
ols_master.loc[idx, "p_FDR"] = out
# ---------------------------------------------------------------------
# Plot heatmaps (3 panels, one per marker) with 4 groups (columns)
# ---------------------------------------------------------------------
cmap = sns.diverging_palette(280, 150, s=90, l=60, as_cmap=True)
absmax = np.nanmax(np.abs(ols_master["beta_std"].values))
if not np.isfinite(absmax) or absmax == 0:
absmax = 0.5
vmin, vmax = -absmax, absmax
fig, axes = plt.subplots(1, 3, figsize=(24, 10))
axes = np.array(axes).flatten()
for i, marker in enumerate(pathology_markers):
ax = axes[i]
sub = ols_master[ols_master["Marker"] == marker].copy()
beta_mat = (
sub.pivot(index="Structure", columns="Group", values="beta_std")
.reindex(index=structures, columns=groups)
)
p_mat = (
sub.pivot(index="Structure", columns="Group", values="p_raw")
.reindex(index=structures, columns=groups)
)
pf_mat = (
sub.pivot(index="Structure", columns="Group", values="p_FDR")
.reindex(index=structures, columns=groups)
)
annot = beta_mat.copy().astype(object)
for r in beta_mat.index:
for c in beta_mat.columns:
b = beta_mat.loc[r, c]
pr = p_mat.loc[r, c]
pf = pf_mat.loc[r, c]
if pd.notna(b) and pd.notna(pr):
annot.loc[r, c] = f"β={b:.2f}{fdr_star(pf)}\n(p={pr:.3f})"
else:
annot.loc[r, c] = ""
sns.heatmap(
beta_mat,
annot=annot,
fmt="",
cmap=cmap,
vmin=vmin,
vmax=vmax,
cbar=(i == 2),
linewidths=0.6,
linecolor="white",
ax=ax,
annot_kws={"fontsize": 22, "color": "black"}
)
ax.set_title(titles[marker], fontsize=22, pad=10, color="black")
ax.set_xticklabels(group_labels, fontsize=16, rotation=0, color="black")
ax.set_xlabel("")
ax.set_ylabel("")
if i == 0:
ax.set_yticklabels([pretty_y[s] for s in structures], fontsize=18, color="black")
else:
ax.set_yticklabels([])
# ---------------------------------------------------------------------
# Format colorbar
# ---------------------------------------------------------------------
cbar = axes[-1].collections[0].colorbar
cbar.set_label("Standardized β (marker effect)", fontsize=18, labelpad=10)
cbar.ax.tick_params(labelsize=16)
plt.tight_layout(rect=[0, 0, 1, 0.95])
plt.savefig("heatmap_beta_rawp_FDRstars_vascular_4groups.jpg", dpi=600, bbox_inches="tight")
plt.show()
# %%
##########################################################################################
########## Linear Mixed Effects Model: Polypathology
##########################################################################################
# %%
"""
Generates 2×2 disease-specific panels showing standardized pathology–volume associations across postmortem limbic/subcortical structures.
For each diagnostic group and structure, the script fits covariate-adjusted standardized OLS models using Tau, α-synuclein, and TDP-43 pathology predictors, applies BH-FDR correction within each disease group, and visualizes standardized β coefficients as grouped bar plots with significance annotations.
"""
# ---------------------------------------------------------------------
# Use your cleaned dataframe exactly as-is
# ---------------------------------------------------------------------
df_use = df.copy() # ← very important
df_use["NPDx1"] = df_use["NPDx1"].astype(str).str.strip().str.lower()
df_use["Sex"] = pd.to_numeric(df_use["Sex"], errors="coerce")
disease_groups = [
"alzheimer's disease",
"lewy body disease",
"ftld-tdp",
"tauopathies"
]
df_use = df_use[df_use["NPDx1"].isin(disease_groups)].copy()
# ---------------------------------------------------------------------
# Structures + region prefixes
# ---------------------------------------------------------------------
structures = ["hippocampus", "amygdala", "caudate", "putamen", "thalamus", "pallidum"]
region_map = {
"hippocampus": "EC_CS_DG",
"amygdala": "Amyg",
"caudate": "CP",
"putamen": "CP",
"thalamus": "TS",
"pallidum": "GP",
}
# Normalize volumes
if "antemortem_icv" in df_use.columns:
for s in structures:
pm = f"postmortem_{s}"
if pm in df_use.columns:
df_use[f"{pm}_norm"] = df_use[pm] / df_use["antemortem_icv"]
# Display names
pretty_group = {
"alzheimer's disease": "Alzheimer’s Disease",
"lewy body disease": "Lewy Body Disease",
"ftld-tdp": "FTLD-TDP",
"tauopathies": "Tauopathies"
}
# Colors per structure
structure_palette = {
"hippocampus": "#6a51a3",
"amygdala": "#9e9ac8",
"caudate": "#807dba",
"putamen": "#bcbddc",
"thalamus": "#cbc9e2",
"pallidum": "#dadaeb"
}
# ---------------------------------------------------------------------
# Fit LME for a single structure
# ---------------------------------------------------------------------
def fit_structure(df_sub, s):
prefix = region_map[s]
# outcome (normalized preferred)
candidates = [f"postmortem_{s}_norm", f"postmortem_{s}"]
outcome = next((x for x in candidates if x in df_sub.columns), None)
if outcome is None:
return pd.DataFrame()
# pathology predictors: *only Tau, aSyn, TDP43*
path_cols = [
c for c in df_sub.columns
if c.startswith(prefix)
and any(x in c for x in ["Tau", "aSyn", "TDP43"])
]
if not path_cols:
return pd.DataFrame()
covars = [c for c in ["AgeatDeath", "Sex", "Education", "PMI"] if c in df_sub.columns]
keep = ["INDDID", outcome] + path_cols + covars
d = df_sub[keep].dropna()
if len(d) < 10:
return pd.DataFrame()
# Standardize numeric columns
for col in d.select_dtypes(include=[np.number]).columns:
if d[col].std(ddof=0) > 0:
d[col] = (d[col] - d[col].mean()) / d[col].std(ddof=0)
formula = outcome + " ~ " + " + ".join(path_cols + covars)
try:
model = smf.ols(formula, d).fit()
except:
return pd.DataFrame()
out = []
for param, coef, pval in zip(model.params.index, model.params.values, model.pvalues.values):
if param == "Intercept":
continue
if any(x in param for x in ["Tau", "aSyn", "TDP43"]):
out.append({
"Structure": s,
"Marker": param.replace(prefix, ""),
"StdBeta": coef,
"p": pval,
"N": len(d)
})
return pd.DataFrame(out)
# ---------------------------------------------------------------------
# 2×2 PANEL FIGURE (one panel per disease)
# ---------------------------------------------------------------------
fig, axes = plt.subplots(2, 2, figsize=(16, 12), sharey=True)
axes = axes.flatten()
for idx, dx in enumerate(disease_groups):
ax = axes[idx]
gdf = df_use[df_use["NPDx1"] == dx].copy()
if gdf.empty:
ax.text(0.5, 0.5, "No data", ha="center", va="center")
ax.axis("off")
continue
all_out = []
for s in structures:
r = fit_structure(gdf, s)
if not r.empty:
all_out.append(r)
if not all_out:
ax.text(0.5, 0.5, "No valid models", ha="center", va="center")
ax.axis("off")
continue
df_dx = pd.concat(all_out, ignore_index=True)
# FDR per disease
_, p_corr = pg.multicomp(df_dx["p"], method="fdr_bh")
df_dx["p_FDR"] = p_corr
df_dx["Sig"] = df_dx["p_FDR"].apply(
lambda p: "***" if p < 0.001 else "**"
if p < 0.01 else "*" if p < 0.05 else ""
)
# Short clean marker names
df_dx["MarkerClean"] = df_dx["Marker"].replace({
"Tau": "Tau",
"aSyn": "aSyn",
"TDP43": "TDP43"
})
# Barplot
sns.barplot(
data=df_dx,
x="MarkerClean", y="StdBeta",
hue="Structure",
palette=structure_palette,
errorbar=None, ax=ax
)
# Add significance stars
for _, r in df_dx.iterrows():
x = ["Tau", "aSyn", "TDP43"].index(r["MarkerClean"])
y = r["StdBeta"]
ax.text(
x + (structures.index(r["Structure"]) * 0.015),
y + np.sign(y) * 0.03,
r["Sig"],
ha="center", fontsize=11, weight="bold"
)
ax.axhline(0, color="black", lw=1)
ax.set_title(pretty_group[dx], fontsize=15, weight="bold")
ax.set_xlabel("")
if idx % 2 == 0:
ax.set_ylabel("Std β")
else:
ax.set_ylabel("")
# Legend
handles, labels = axes[0].get_legend_handles_labels()
fig.legend(
handles, labels, title="Structure",
bbox_to_anchor=(0.5, -0.02), loc="lower center",
ncol=6, frameon=False
)
plt.suptitle("Postmortem Polypathology–Volume Effects (LME, standardized β)", fontsize=17)
plt.tight_layout(rect=[0.05, 0.05, 1, 0.95])
plt.show()
# %%
# %%
"""
Runs disease- and structure-specific standardized postmortem pathology–volume models for limbic and subcortical regions.
For each diagnostic group and structure, the script fits covariate-adjusted standardized OLS models using Tau, TDP-43, and α-synuclein pathology predictors, extracts β coefficients, standard errors, 95% confidence intervals, raw and FDR-corrected p-values, and visualizes pathology-specific standardized effects in 2×2 disease panels.
"""
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
import seaborn as sns
import matplotlib.pyplot as plt
import pingouin as pg
import warnings, matplotlib
warnings.filterwarnings("ignore")
matplotlib.rcParams['font.family'] = 'DejaVu Sans'
# ===============================================================
# 🔧 1. Load YOUR dataset
# ===============================================================
df_use = df.copy()
df_use["NPDx1"] = df_use["NPDx1"].astype(str).str.strip().str.lower()
if "Sex" in df_use.columns:
df_use["Sex"] = pd.to_numeric(df_use["Sex"], errors="coerce")
groups = ["alzheimer's disease", "lewy body disease", "ftld-tdp", "tauopathies"]
df_use = df_use[df_use["NPDx1"].isin(groups)].copy()
# ===============================================================
# 🔧 2. Structures + pathology prefixes
# ===============================================================
structures = ["hippocampus", "amygdala", "caudate", "putamen", "thalamus", "pallidum"]
region_map = {
"hippocampus": "EC_CS_DG",
"amygdala": "Amyg",
"caudate": "CP",
"putamen": "CP",
"thalamus": "TS",
"pallidum": "GP",
}
# Normalize volumes by ICV
if "antemortem_icv" in df_use.columns:
for s in structures:
pm = f"postmortem_{s}"
if pm in df_use.columns:
df_use[f"{pm}_norm"] = df_use[pm] / df_use["antemortem_icv"]
# ===============================================================
# 🔧 3. Run LME for each disease group × structure
# ===============================================================
def run_lme_by_group(df_sub, group_name):
gdf = df_sub[df_sub["NPDx1"] == group_name].copy()
results = []
if gdf.empty:
print(f"[skip] {group_name}: no subjects")
return pd.DataFrame()
for s in structures:
# pick normalized → raw
out_candidates = [f"postmortem_{s}_norm", f"postmortem_{s}"]
outcome = next((x for x in out_candidates if x in gdf.columns), None)
if outcome is None:
continue
prefix = region_map[s]
# pathology predictors: Tau / aSyn / TDP43
path_cols = [
c for c in gdf.columns
if c.startswith(prefix) and any(k in c for k in ["Tau", "aSyn", "TDP43"])
]
if not path_cols:
print(f"[skip] {group_name} {s}: no matching pathology variables ({prefix})")
continue
covars = [c for c in ["AgeatDeath", "Sex", "Education", "PMI"] if c in gdf.columns]
keep = ["INDDID", outcome] + path_cols + covars
d = gdf[keep].dropna()
if len(d) < 10:
print(f"[skip] {group_name} {s}: insufficient N={len(d)}")
continue
# --------------------
# Standardize
# --------------------
for col in d.select_dtypes(include=[np.number]).columns:
if d[col].std(ddof=0) > 0:
d[col] = (d[col] - d[col].mean()) / d[col].std(ddof=0)
rhs = path_cols + covars
formula = f"{outcome} ~ " + " + ".join(rhs)
print("\n--------------------------------------------------------")
print(f"GROUP: {group_name.upper()} | STRUCTURE: {s.upper()}")
print("MODEL FORMULA:")
print(formula)
print("--------------------------------------------------------")
try:
model = smf.ols(formula, d).fit()
except Exception as e:
print("ERROR:", e)
continue
"""
for p, coef, pval in zip(model.params.index, model.params.values, model.pvalues.values):
if p == "Intercept":
continue
results.append({
"Group": group_name,
"Structure": s,
"Parameter": p,
"StdBeta": coef,
"p-value": pval,
"N": len(d)
})
"""
# inside the loop: for p, coef, pval in zip(...)
for p, coef, pval in zip(model.params.index, model.params.values, model.pvalues.values):
if p == "Intercept":
continue
se = model.bse.get(p, np.nan)
ci_low = coef - 1.96 * se
ci_high = coef + 1.96 * se
results.append({
"Group": group_name,
"Structure": s,
"Parameter": p,
"StdBeta": coef,
"SE": se,
"CI_low": ci_low,
"CI_high": ci_high,
"p-value": pval,
"N": len(d)
})
res_df = pd.DataFrame(results)
if not res_df.empty:
_, p_corr = pg.multicomp(res_df["p-value"], method="fdr_bh")
res_df["p-FDR"] = p_corr
res_df["Sig(FDR)"] = res_df["p-FDR"].apply(
lambda p: "***" if p < 0.001 else "**" if p < 0.01 else "*" if p < 0.05 else ""
)
return res_df
# ===============================================================
# 🔧 4. Run all models
# ===============================================================
all_results = []
for g in groups:
print(f"\n\n====================== {g.upper()} ======================")
out = run_lme_by_group(df_use, g)
if not out.empty:
all_results.append(out)
print(out.to_string(index=False))
else:
print("NO MODELS FIT.")
# ===============================================================
# 🔧 5. Combine results
# ===============================================================
if not all_results:
raise SystemExit("❌ No results to plot")
results_df = pd.concat(all_results, ignore_index=True)
# only pathology rows
subset = results_df[
results_df["Parameter"].str.contains("Tau|TDP43|aSyn", case=False)
]
# ===============================================================
# 🔧 6. Plot — 2×2 GRID (one panel per disease)
# ===============================================================
sns.set(style="whitegrid", context="talk")
fig, axes = plt.subplots(2, 2, figsize=(12, 8), sharey=True)
axes = axes.flatten()
for i, g in enumerate(groups):
ax = axes[i]
sub = subset[subset["Group"] == g]
if sub.empty:
ax.text(0.5, 0.5, "No data", ha="center", va="center")
ax.axis("off")
continue
sns.barplot(
data=sub,
x="Parameter", y="StdBeta",
hue="Structure",
errorbar=None, ax=ax
)
ax.axhline(0, color="black", lw=1)
ax.tick_params(axis="x", rotation=90, labelsize=8)
ax.set_title(g.title(), fontsize=13)
ax.set_xlabel("")
ax.set_ylabel("Std β")
handles, labels = axes[0].get_legend_handles_labels()
fig.legend(handles, labels, title="Structure",
loc="lower center", bbox_to_anchor=(0.5, -0.05),
ncol=6, frameon=False)
plt.suptitle("Postmortem Polypathology–Volume Associations (Standardized β)", fontsize=16)
plt.tight_layout(rect=[0,0.05,1,0.95])
plt.show()
# %%
"""
Exports postmortem polypathology model results to a publication-ready Word summary table.
The script formats disease-specific tables reporting structure, predictor, standardized β, SE, 95% CI, raw p-values, and FDR-corrected p-values with superscript significance markers, then saves the results as a DOCX file.
"""
import pandas as pd
from docx import Document
from docx.enum.table import WD_TABLE_ALIGNMENT
from docx.oxml import OxmlElement
from docx.oxml.ns import qn
df_in = results_df.copy()
# ======================================================================================
# Pretty names
# ======================================================================================
pretty_group = {
"alzheimer's disease": "Alzheimer’s disease",
"lewy body disease": "Lewy body disease",
"ftld-tdp": "FTLD-TDP",
"tauopathies": "Tauopathies",
}
pretty_struct = {
"hippocampus": "Hippocampus",
"amygdala": "Amygdala",
"caudate": "Caudate",
"putamen": "Putamen",
"thalamus": "Thalamus",
"pallidum": "Pallidum",
}
pretty_param = {
"AgeatDeath": "Age at death",
"Sex": "Sex",
"Education": "Education",
"PMI": "PMI",
}
def clean_param(p):
if "Tau" in p: return "p-tau"
if "aSyn" in p: return "α-synuclein"
if "TDP43" in p: return "TDP-43"
return pretty_param.get(p, p)
# ======================================================================================
# Superscript helper
# ======================================================================================
def add_superscript(run, text):
for ch in text:
r = OxmlElement('w:r')
t = OxmlElement('w:t')
t.text = ch
r.append(t)
rpr = OxmlElement('w:rPr')
vert = OxmlElement('w:vertAlign')
vert.set(qn('w:val'), 'superscript')
rpr.append(vert)
r.append(rpr)
run._r.addnext(r)
# ======================================================================================
# DOCX init
# ======================================================================================
doc = Document()
title = doc.add_heading("Polypathology LME – Postmortem MRI", level=1)
title.alignment = 1
cap = (
"Linear regression models predicting postmortem MRI volumes from multiple pathology "
"burdens (p-tau, α-synuclein, TDP-43) and covariates (age at death, sex, education, PMI). "
"Reported: standardized β, standard error (SE), 95% confidence interval (CI), raw p-value "
"and FDR-corrected p-value with superscript significance (*, **, ***). "
"One table per disease group."
)
doc.add_paragraph(cap)
doc.add_page_break()
# ======================================================================================
# Build one table per disease group
# ======================================================================================
for g in pretty_group.keys():
sub = df_in[df_in["Group"] == g].copy()
if sub.empty:
continue
doc.add_heading(pretty_group[g], level=2)
sub["ParamPrint"] = sub["Parameter"].apply(clean_param)
sub.sort_values(["Structure", "ParamPrint"], inplace=True)
structs = sub["Structure"].unique()
total_rows = sum(len(sub[sub["Structure"] == s]) for s in structs)
# ⭐ NOW 7 columns instead of 6
table = doc.add_table(rows=total_rows + 1, cols=7)
table.style = "Table Grid"
table.alignment = WD_TABLE_ALIGNMENT.CENTER
# Header
hdr = table.rows[0].cells
hdr[0].text = "Structure"
hdr[1].text = "Predictor"
hdr[2].text = "Std β"
hdr[3].text = "SE"
hdr[4].text = "CI (95%)"
hdr[5].text = "p(raw)"
hdr[6].text = "p(FDR)"
row_idx = 1
for s in structs:
block = sub[sub["Structure"] == s]
span = len(block)
base = table.cell(row_idx, 0)
for _ in range(span - 1):
base.merge(table.cell(row_idx + 1, 0))
base.text = pretty_struct[s]
for _, r in block.iterrows():
cells = table.rows[row_idx].cells
cells[1].text = r["ParamPrint"]
cells[2].text = f"{r['StdBeta']:.3f}"
cells[3].text = f"{r['SE']:.3f}"
cells[4].text = f"[{r['CI_low']:.3f}, {r['CI_high']:.3f}]"
# p(raw)
cells[5].text = f"{r['p-value']:.2e}"
# p(FDR) with superscripts
run = cells[6].paragraphs[0].add_run(f"{r['p-FDR']:.2e}")
if r["Sig(FDR)"]:
add_superscript(run, r["Sig(FDR)"])
row_idx += 1
doc.add_page_break()
# ======================================================================================
# SAVE
# ======================================================================================
doc.save("polypathology_LME_results.docx")
print("Saved: polypathology_LME_results.docx")
# %%
# %%
##########################################################################################
# Polypathology LME Heatmaps (FINAL VERSION)
# Updated with:
# - β + raw p in cell
# - FDR stars (***) based on p-FDR
# - "Not applicable" explicitly shown for missing cells
# - Missing cells are masked so they do NOT get misleading colors
##########################################################################################
import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import matplotlib
from matplotlib import colors
matplotlib.rcParams["font.family"] = "DejaVu Sans"
sns.set(style="white", context="talk")
def plot_polypathology_heatmaps_dfuse(
results_df,
groups=(
"alzheimer's disease",
"lewy body disease",
"ftld-tdp",
"tauopathies"
),
out_png="lme_polypathology_heatmaps_clean_final.png"
):
# ----------------------------------------------------------
# 1) Ensure Group column exists
# ----------------------------------------------------------
if "Group" not in results_df.columns:
if "NPDx1" in results_df.columns:
results_df = results_df.rename(columns={"NPDx1": "Group"})
else:
raise ValueError("results_df must contain Group or NPDx1 column.")
# ----------------------------------------------------------
# 2) Filter predictors (Tau, TDP43, aSyn)
# ----------------------------------------------------------
mask = results_df["Parameter"].str.contains("Tau|TDP43|aSyn", case=False, na=False)
gdf = results_df[mask].copy()
def to_pred(p):
p = str(p)
if "TDP43" in p:
return "TDP-43"
if "Tau" in p:
return "p-tau"
if "aSyn" in p:
return "α-synuclein"
return None
gdf["Predictor"] = gdf["Parameter"].apply(to_pred)
gdf = gdf[~gdf["Predictor"].isna()].copy()
# ----------------------------------------------------------
# 3) Structure ordering
# ----------------------------------------------------------
struct_order = ["hippocampus", "amygdala", "caudate",
"putamen", "thalamus", "pallidum"]
struct_label_map = {s: s.capitalize() for s in struct_order}
gdf = gdf[gdf["Structure"].isin(struct_order)].copy()
gdf["StructureLabel"] = gdf["Structure"].map(struct_label_map)
predictor_order = ["TDP-43", "p-tau", "α-synuclein"]
# ----------------------------------------------------------
# 4) Color map (soft PRGn)
# ----------------------------------------------------------
base_cmap = sns.color_palette("PRGn", as_cmap=True)
cmap = colors.LinearSegmentedColormap.from_list(
"PRGn_soft",
base_cmap(np.linspace(0.1, 0.9, 256))
)
# ----------------------------------------------------------
# 5) Pretty titles
# ----------------------------------------------------------
title_map = {
"alzheimer's disease": "Alzheimer’s disease",
"lewy body disease": "Lewy body disease",
"ftld-tdp": "FTLD-TDP",
"tauopathies": "FTLD-Tau"
}
# ----------------------------------------------------------
# 6) FDR star helper
# ----------------------------------------------------------
def fdr_star(pf):
if pd.isna(pf):
return ""
if pf < 0.001:
return "***"
if pf < 0.01:
return "**"
if pf < 0.05:
return "*"
return ""
# ----------------------------------------------------------
# 7) Create figure (2×2)
# ----------------------------------------------------------
fig, axes = plt.subplots(2, 2, figsize=(12.5, 9))
axes = axes.flatten()
for i, g in enumerate(groups):
ax = axes[i]
sub = gdf[gdf["Group"] == g].copy()
if sub.empty:
ax.text(0.5, 0.5, "No data", ha="center", va="center")
ax.axis("off")
continue
# ----------------------------------------------------------
# Pivot tables
# ----------------------------------------------------------
beta_pivot = sub.pivot_table(
index="StructureLabel", columns="Predictor", values="StdBeta"
).reindex(
index=[struct_label_map[s] for s in struct_order],
columns=predictor_order
)
pval_pivot = sub.pivot_table(
index="StructureLabel", columns="Predictor", values="p-value"
).reindex(index=beta_pivot.index, columns=beta_pivot.columns)
pfdr_pivot = sub.pivot_table(
index="StructureLabel", columns="Predictor", values="p-FDR"
).reindex(index=beta_pivot.index, columns=beta_pivot.columns)
# ----------------------------------------------------------
# Mask missing β (these cells become "Not applicable")
# ----------------------------------------------------------
na_mask = beta_pivot.isna()
# ----------------------------------------------------------
# Heatmap background (masked so NA cells are blank/white)
# ----------------------------------------------------------
sns.heatmap(
beta_pivot,
ax=ax,
cmap=cmap,
center=0,
vmin=-1, vmax=1,
mask=na_mask,
linewidths=0.6,
linecolor="white",
cbar=False,
annot=False
)
# ----------------------------------------------------------
# Axis formatting
# ----------------------------------------------------------
ax.set_yticks(np.arange(len(struct_order)) + 0.5)
ax.set_yticklabels([struct_label_map[s] for s in struct_order], fontsize=14)
# x-axis only on bottom row panels
if i in [2, 3]:
ax.set_xticks(np.arange(len(predictor_order)) + 0.5)
ax.set_xticklabels(predictor_order, fontsize=14)
else:
ax.set_xticks([])
ax.set_xticklabels([])
ax.set_xlabel("")
ax.set_ylabel("")
ax.tick_params(axis="both", length=0)
# ----------------------------------------------------------
# Annotate text: valid -> β + raw p + FDR stars
# missing -> "Not applicable"
# ----------------------------------------------------------
for y in range(beta_pivot.shape[0]):
for x in range(beta_pivot.shape[1]):
beta = beta_pivot.iloc[y, x]
pval = pval_pivot.iloc[y, x]
pfdr = pfdr_pivot.iloc[y, x]
# Not applicable if no beta OR no p-value
if pd.isna(beta) or pd.isna(pval) or np.isclose(beta, 0, atol=1e-6):
ax.text(
x + 0.5, y + 0.5,
"Not\napplicable",
ha="center", va="center",
fontsize=13,
color="black"
)
continue
star = fdr_star(pfdr)
text = f"β={beta:.2f}{star}\n(p={pval:.3f})"
ax.text(
x + 0.5, y + 0.5,
text,
ha="center", va="center",
fontsize=15,
color="black"
)
ax.set_title(title_map.get(g, g), fontsize=15, pad=12)
# Turn off any unused axes (if fewer than 4 groups)
for j in range(len(groups), 4):
axes[j].axis("off")
# ----------------------------------------------------------
# Global labels
# ----------------------------------------------------------
fig.text(0.04, 0.5, "Structure", va="center", rotation=90, fontsize=14)
fig.text(0.52, 0.06, "Predictor", ha="center", fontsize=14)
# ----------------------------------------------------------
# Shared colorbar
# ----------------------------------------------------------
cbar_ax = fig.add_axes([0.93, 0.28, 0.02, 0.45])
sm = plt.cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=-1, vmax=1))
sm.set_array([])
cbar = plt.colorbar(sm, cax=cbar_ax)
cbar.set_label("Std β", fontsize=12)
plt.tight_layout(rect=[0.08, 0.10, 0.9, 0.95])
plt.savefig(out_png, dpi=600, bbox_inches="tight")
plt.show()
# ---- RUN ----
plot_polypathology_heatmaps_dfuse(results_df)
# %%
##########################################################################################
########## Linear Mixed Effects Model: Primary pathology
##########################################################################################
# %%
"""
Fits standardized primary-pathology models linking disease-specific pathology burden to postmortem limbic/subcortical volumes while adjusting for key covariates.
The script runs one model set per diagnostic group, extracts standardized β estimates with SE, 95% CI, raw and FDR-corrected p-values, and exports the results as publication-ready DOCX tables.
"""
##########################################################################################
# 1) Load df_use
##########################################################################################
df_use = df.copy()
df_use["NPDx1"] = df_use["NPDx1"].astype(str).str.strip().str.lower()
df_use["Sex"] = pd.to_numeric(df_use["Sex"], errors="coerce")
disease_groups = [
"alzheimer's disease",
"lewy body disease",
"ftld-tdp",
"tauopathies"
]
df_use = df_use[df_use["NPDx1"].isin(disease_groups)].copy()
##########################################################################################
# 2) Structures + prefixes
##########################################################################################
structures = ["hippocampus", "amygdala", "caudate", "putamen",
"thalamus", "pallidum"]
region_prefix_map = {
"hippocampus": "EC_CS_DG",
"amygdala": "Amyg",
"caudate": "CP",
"putamen": "CP",
"thalamus": "TS",
"pallidum": "GP",
}
primary_marker_suffix = {
"alzheimer's disease": "Tau",
"lewy body disease": "aSyn",
"ftld-tdp": "TDP43",
"tauopathies": "Tau",
}
covars = ["AgeatDeath", "Sex", "Education", "PMI"]
##########################################################################################
# 3) Normalize volumes
##########################################################################################
if "antemortem_icv" in df_use.columns:
for s in structures:
pm = f"postmortem_{s}"
if pm in df_use.columns:
df_use[f"{pm}_norm"] = df_use[pm] / df_use["antemortem_icv"]
##########################################################################################
# 4) LME WITH CI
##########################################################################################
def run_primary_lme_by_group(df_use, group_name, structures):
results_all = []
marker_suffix = primary_marker_suffix[group_name]
gdf = df_use[df_use["NPDx1"] == group_name].copy()
if gdf.empty:
return pd.DataFrame()
for s in structures:
region_prefix = region_prefix_map[s]
# find pathology for this structure
path_candidates = [c for c in gdf.columns if c.lower().startswith(region_prefix.lower())]
path_cols = [c for c in path_candidates if marker_suffix.lower() in c.lower()]
if not path_cols:
continue
path_col = path_cols[0]
outcome = next(
(c for c in (f"postmortem_{s}_norm", f"postmortem_{s}")
if c in gdf.columns),
None
)
if outcome is None:
continue
covars_here = [c for c in covars if c in gdf.columns]
d = gdf[["INDDID", "NPDx1", outcome, path_col] + covars_here].dropna()
if len(d) < 8:
continue
# Z-score numeric vars
for col in d.select_dtypes(include=[np.number]).columns:
if d[col].std(ddof=0) > 0:
d[col] = (d[col] - d[col].mean()) / d[col].std(ddof=0)
rhs = [path_col] + covars_here
formula = f"{outcome} ~ " + " + ".join(rhs)
res = smf.ols(formula, d).fit()
# -------------------------------
# PRIMARY predictor + CI
# -------------------------------
coef = res.params.get(path_col, np.nan)
se = res.bse.get(path_col, np.nan)
ci_low = coef - 1.96 * se
ci_high = coef + 1.96 * se
results_all.append({
"Group": group_name,
"Structure": s,
"Parameter": path_col,
"Type": "Primary",
"StdBeta": coef,
"CI_low": ci_low,
"CI_high": ci_high,
"p-value": res.pvalues.get(path_col, np.nan),
"N": len(d["INDDID"].unique())
})
# -------------------------------
# Covariates + CI
# -------------------------------
for cov in covars_here:
coef = res.params.get(cov, np.nan)
se = res.bse.get(cov, np.nan)
ci_low = coef - 1.96 * se
ci_high = coef + 1.96 * se
results_all.append({
"Group": group_name,
"Structure": s,
"Parameter": cov,
"Type": "Covariate",
"StdBeta": coef,
"CI_low": ci_low,
"CI_high": ci_high,
"p-value": res.pvalues.get(cov, np.nan),
"N": len(d["INDDID"].unique())
})
res_df = pd.DataFrame(results_all)
# FDR correction
if not res_df.empty:
_, pf = pg.multicomp(res_df["p-value"], method="fdr_bh")
res_df["p-FDR"] = pf
res_df["Sig(FDR)"] = res_df["p-FDR"].apply(
lambda p: "***" if p < 0.001 else
"**" if p < 0.01 else
"*" if p < 0.05 else ""
)
return res_df
##########################################################################################
# 5) Run all groups
##########################################################################################
all_res = []
for g in disease_groups:
r = run_primary_lme_by_group(df_use, g, structures)
if not r.empty:
all_res.append(r)
primary_results_df = pd.concat(all_res, ignore_index=True)
##########################################################################################
# 6) Pretty printing
##########################################################################################
marker_clean = {
"EC_CS_DGTau": "p-tau (EC/CS/DG)",
"AmygTau": "p-tau (Amygdala)",
"CPTau": "p-tau (Caudate/Putamen)",
"TSTau": "p-tau (Thalamus)",
"GPTau": "p-tau (Pallidum)",
"EC_CS_DGaSyn": "α-syn (EC/CS/DG)",
"EC_CS_DGTDP43": "TDP-43 (EC/CS/DG)"
}
pretty_cov = {
"AgeatDeath": "Age at death",
"Sex": "Sex",
"Education": "Education",
"PMI": "PMI"
}
def pretty_param(p):
if p in marker_clean: return marker_clean[p]
if p in pretty_cov: return pretty_cov[p]
return p
struct_pretty = {
"hippocampus": "Hippocampus",
"amygdala": "Amygdala",
"caudate": "Caudate",
"putamen": "Putamen",
"thalamus": "Thalamus",
"pallidum": "Pallidum"
}
pretty_group_short = {
"alzheimer's disease": "AD",
"lewy body disease": "LBD",
"ftld-tdp": "FTLD-TDP",
"tauopathies": "Tauopathies"
}
##########################################################################################
# 7) Build final table with CI
##########################################################################################
df2 = primary_results_df.copy()
df2["GroupPrint"] = df2["Group"].map(pretty_group_short)
df2["StructurePrint"] = df2["Structure"].map(struct_pretty)
df2["ParameterPrint"] = df2["Parameter"].apply(pretty_param)
df2["StdBeta_fmt"] = df2["StdBeta"].map(lambda x: f"{x:.3f}")
df2["CI_fmt"] = df2.apply(lambda r: f"[{r.CI_low:.3f}, {r.CI_high:.3f}]", axis=1)
df2["p_fmt"] = df2["p-value"].map(lambda x: f"{x:.4f}")
df2["pFDR_fmt"] = df2["p-FDR"].map(lambda x: f"{x:.4f}")
df2 = df2.sort_values(["GroupPrint", "StructurePrint", "Type"])
final_table = df2[[
"GroupPrint", "StructurePrint", "ParameterPrint",
"StdBeta_fmt", "CI_fmt", "p_fmt", "pFDR_fmt",
"Sig(FDR)", "N"
]]
##########################################################################################
# 8) DOCX EXPORT
##########################################################################################
doc = Document()
title = doc.add_heading("Primary Pathology → Postmortem Volume (Standardized LME)", level=1)
title.alignment = 1
doc.add_paragraph(
"Each table reports standardized β, 95% confidence interval (CI), raw p-value, "
"FDR-corrected p-value, significance code, and sample size."
)
doc.add_page_break()
groups_in_order = ["alzheimer's disease", "lewy body disease", "ftld-tdp", "tauopathies"]
for g in groups_in_order:
sub = final_table[final_table["GroupPrint"] == pretty_group_short[g]]
if sub.empty:
continue
h = doc.add_heading(pretty_group[g], level=2)
tbl = doc.add_table(rows=1, cols=len(final_table.columns))
tbl.style = "Table Grid"
hdr = tbl.rows[0].cells
for j, colname in enumerate(final_table.columns):
hdr[j].text = colname
for _, row in sub.iterrows():
rw = tbl.add_row().cells
for j, colname in enumerate(final_table.columns):
rw[j].text = str(row[colname])
doc.add_page_break()
output_file = "primary_pathology_LME_results_by_group_with_CI.docx"
doc.save(output_file)
print("Saved DOCX:", output_file)
# %%
"""
Heatmaps
"""
def plot_primary_pathology_heatmaps_clean(primary_results_df):
# -------------------------------------------------------
# Filter *only PRIMARY PATHOLOGY rows*
# -------------------------------------------------------
dfp = primary_results_df[primary_results_df["Type"] == "Primary"].copy()
# -------------------------------------------------------
# Structures (ordered)
# -------------------------------------------------------
struct_order = ["hippocampus", "amygdala", "caudate",
"putamen", "thalamus", "pallidum"]
struct_labels = {s: s.capitalize() for s in struct_order}
dfp["StructureLabel"] = dfp["Structure"].map(struct_labels)
# -------------------------------------------------------
# Pretty pathology names for prefix-based parameters
# -------------------------------------------------------
pretty_path = {
"EC_CS_DGTau": "p-tau",
"AmygTau": "p-tau",
"CPTau": "p-tau",
"TSTau": "p-tau",
"GPTau": "p-tau",
"EC_CS_DGaSyn": "α-syn",
"EC_CS_DGTDP43": "TDP-43"
}
# Parameter → pretty label
dfp["PathLabel"] = dfp["Parameter"].map(pretty_path).fillna("")
# -------------------------------------------------------
# Correct group order and labels
# -------------------------------------------------------
group_order = [
"alzheimer's disease",
"lewy body disease",
"ftld-tdp",
"tauopathies"
]
column_labels = [
"AD\n(p-tau)",
"LBD\n(α-syn)",
"FTLD-TDP\n(TDP-43)",
"Tauopathies\n(p-tau)"
]
# -------------------------------------------------------
# Create matrices (StdBeta, p, stars)
# -------------------------------------------------------
M = []
P = []
Stars = []
for s in struct_order:
beta_row = []
p_row = []
star_row = []
for grp in group_order:
sub = dfp[(dfp["Group"] == grp) & (dfp["Structure"] == s)]
if sub.empty:
beta_row.append(np.nan)
p_row.append(np.nan)
star_row.append("")
continue
beta = sub["StdBeta"].values[0]
pval = sub["p-value"].values[0]
pfdr = sub["p-FDR"].values[0]
# stars based on FDR
if pfdr < 0.001:
sig = "***"
elif pfdr < 0.01:
sig = "**"
elif pfdr < 0.05:
sig = "*"
else:
sig = ""
beta_row.append(beta)
p_row.append(pval)
star_row.append(sig)
M.append(beta_row)
P.append(p_row)
Stars.append(star_row)
M = np.array(M)
P = np.array(P)
Stars = np.array(Stars)
# -------------------------------------------------------
# Colormap
# -------------------------------------------------------
base_cmap = sns.color_palette("PRGn", as_cmap=True)
cmap = colors.LinearSegmentedColormap.from_list(
"PRGn_soft", base_cmap(np.linspace(0.15, 0.85, 256))
)
# -------------------------------------------------------
# Plot
# -------------------------------------------------------
fig, ax = plt.subplots(figsize=(13, 6))
sns.heatmap(
M,
ax=ax,
cmap=cmap,
vmin=-1, vmax=1,
linewidths=0.6,
linecolor="white",
cbar=True,
cbar_kws={"label": "Std β"},
annot=False
)
ax.set_yticks(np.arange(len(struct_order)) + 0.5)
ax.set_yticklabels([struct_labels[s] for s in struct_order], rotation=0)
ax.set_xticks(np.arange(len(group_order)) + 0.5)
ax.set_xticklabels(column_labels, rotation=0, ha="center")
# -------------------------------------------------------
# Annotate cells
# -------------------------------------------------------
for y in range(M.shape[0]):
for x in range(M.shape[1]):
beta = M[y, x]
pval = P[y, x]
sig = Stars[y, x]
if np.isnan(beta):
continue
txt = f"β={beta:.2f}{sig}\np={pval:.3f}"
ax.text(
x + 0.5,
y + 0.5,
txt,
ha="center",
va="center",
fontsize=8.5
)
ax.set_title("Primary Pathology → Postmortem Volume (LME)", fontsize=16, pad=18)
ax.set_ylabel("Brain Structure", fontsize=13)
plt.tight_layout()
plt.savefig("primary_pathology_heatmap_clean.png", dpi=600, bbox_inches="tight")
plt.show()
# ===================== RUN =====================
plot_primary_pathology_heatmaps_clean(primary_results_df)
# %%
##########################################################################################
##########################################################################################
####### Gliosis and NeuronLoss: Mediation Analyses
##########################################################################################
##########################################################################################
# %%
"""
Generates disease-specific partial Spearman correlation plots between gliosis/neuron-loss pathology scores and ICV-normalized postmortem limbic/subcortical volumes.
For each diagnostic group, pathology marker, and structure, the script adjusts for age at death, sex, PMI, and education, applies FDR correction within each marker–group set, and displays regression plots annotated with ρ, raw p-value, significance stars, and sample size.
"""
# ---------------------------------------------------------------------
# USE YOUR CLEANED DF EXACTLY AS IS
# ---------------------------------------------------------------------
df_use = df.copy()
df_use["NPDx1"] = df_use["NPDx1"].astype(str).str.strip().str.lower()
df_use["Sex"] = df_use["Sex"].astype("category").cat.codes
if "antemortem_icv" not in df_use.columns:
raise ValueError("Missing 'antemortem_icv' column!")
# ---------------------------------------------------------------------
# Structures and mappings
# ---------------------------------------------------------------------
structures = ["hippocampus", "amygdala", "caudate",
"putamen", "thalamus", "pallidum"]
region_map = {
"caudate": "CP",
"putamen": "CP",
"thalamus": "TS",
"pallidum": "GP",
"hippocampus": "EC_CS_DG",
"amygdala": "Amyg",
}
# Normalize postmortem volumes by ICV
for s in structures:
pm_col = f"postmortem_{s}"
if pm_col in df_use.columns:
df_use[f"{pm_col}_norm"] = df_use[pm_col] / df_use["antemortem_icv"]
# ---------------------------------------------------------------------
# Pathology features
# ---------------------------------------------------------------------
pathology_features = ["Gliosis", "NeuronLoss"]
# Your 4 disease groups already harmonized in df_use
groups = ["alzheimer's disease", "lewy body disease",
"ftld-tdp", "tauopathies"]
sns.set(style="whitegrid", context="talk")
color = "tab:purple"
# ---------------------------------------------------------------------
# Run and plot per pathology marker × group
# ---------------------------------------------------------------------
for marker_suffix in pathology_features:
print(f"\n===== {marker_suffix.upper()} =====")
for g in groups:
subdf = df_use[df_use["NPDx1"] == g]
if subdf.empty:
print(f"[skip] No subjects for group: {g}")
continue
print(f"\n--- {g.upper()} ---")
fig, axes = plt.subplots(
1, len(structures), figsize=(6 * len(structures), 5)
)
if len(structures) == 1:
axes = [axes]
all_stats, all_p = [], []
# -----------------------------
# Compute partial Spearman stats
# -----------------------------
for i, s in enumerate(structures):
ax = axes[i]
post_col = f"postmortem_{s}_norm"
if post_col not in subdf.columns:
ax.text(0.5, 0.5, "Missing volume", ha="center", va="center")
ax.axis("off")
continue
# Find pathology column
prefix = region_map[s]
path_candidates = [c for c in subdf.columns
if c.lower().startswith(prefix.lower())]
path_cols = [
c for c in path_candidates
if marker_suffix.lower() in c.lower()
]
if not path_cols:
ax.text(0.5, 0.5, f"No {prefix}{marker_suffix}",
ha="center", va="center")
ax.axis("off")
continue
path_col = path_cols[0]
cols = [path_col, post_col, "AgeatDeath", "Sex", "PMI", "Education"]
d = subdf[cols].copy().apply(pd.to_numeric, errors="coerce").dropna()
if len(d) < 5:
ax.text(0.5, 0.5, "Insufficient data", ha="center", va="center")
ax.axis("off")
continue
res = pg.partial_corr(
data=d,
x=path_col, y=post_col,
covar=["AgeatDeath", "Sex", "PMI", "Education"],
method="spearman"
)
r, p = res["r"].iloc[0], res["p-val"].iloc[0]
all_stats.append((s, r, p, len(d)))
all_p.append(p)
# -----------------------------
# FDR correction
# -----------------------------
if len(all_p) > 0:
reject, p_corr = pg.multicomp(all_p, method="fdr_bh")
else:
reject, p_corr = [], []
# -----------------------------
# Final plotting with p_FDR stars
# -----------------------------
for (s, r, p, n), pc, rej, ax in zip(all_stats, p_corr, reject, axes):
post_col = f"postmortem_{s}_norm"
prefix = region_map[s]
path_candidates = [c for c in subdf.columns
if c.lower().startswith(prefix.lower())]
path_cols = [
c for c in path_candidates
if marker_suffix.lower() in c.lower()
]
path_col = path_cols[0]
d = subdf[
[path_col, post_col, "AgeatDeath", "Sex", "PMI", "Education"]
].copy().apply(pd.to_numeric, errors="coerce").dropna()
sig = "***" if pc < 0.001 else "**" if pc < 0.01 else "*" if pc < 0.05 else ""
sns.regplot(
data=d, x=path_col, y=post_col,
scatter_kws=dict(alpha=0.7, s=55),
line_kws=dict(color=color, lw=2),
color=color, ax=ax
)
ax.set_title(
f"{s.capitalize()} ({marker_suffix})\n"
f"ρ = {r:.2f}, p = {p:.3f}{sig}, n = {n}",
fontsize=11, fontweight="bold", pad=10
)
ax.set_xlabel(marker_suffix)
ax.set_ylabel("Postmortem Volume / ICV")
ax.grid(True, linestyle=":", alpha=0.5)
plt.suptitle(
f"{marker_suffix} vs Postmortem Volume — {g.title()} (Partial Spearman, FDR-corrected)",
fontsize=16, fontweight="bold"
)
plt.tight_layout(rect=[0, 0, 1, 0.94])
plt.show()
# %%
##########################################################################################
# Final Postmortem Partial Spearman Heatmaps (Gliosis & Neuronal loss)
##########################################################################################
# ---------------------------------------------------------------------
# Use your cleaned df exactly as-is
# ---------------------------------------------------------------------
df_use = df.copy()
df_use["NPDx1"] = df_use["NPDx1"].astype(str).str.strip().str.lower()
df_use["Sex"] = df_use["Sex"].astype("category").cat.codes
# ---------------------------------------------------------------------
# Structures, groups, pathology features
# ---------------------------------------------------------------------
structures = ["hippocampus", "amygdala", "caudate", "putamen", "thalamus", "pallidum"]
region_map = {
"caudate": "CP",
"putamen": "CP",
"thalamus": "TS",
"pallidum": "GP",
"hippocampus": "EC_CS_DG",
"amygdala": "Amyg",
}
groups = ["alzheimer's disease", "lewy body disease",
"ftld-tdp", "tauopathies"]
group_labels = [
"Alzheimer's\ndisease", "Lewy body\ndisease",
"FTLD-TDP", "Tauopathies"
]
pathology_features = [
("Gliosis", "Gliosis"),
("NeuronLoss", "Neuronal loss")
]
covars = ["AgeatDeath", "Sex", "PMI", "Education"]
# ---------------------------------------------------------------------
# Normalize postmortem volumes by ICV
# ---------------------------------------------------------------------
if "antemortem_icv" not in df_use.columns:
raise ValueError("Missing 'antemortem_icv' column!")
for s in structures:
pm_col = f"postmortem_{s}"
if pm_col in df_use.columns:
df_use[f"{pm_col}_norm"] = df_use[pm_col] / df_use["antemortem_icv"]
sns.set(style="white")
# ---------------------------------------------------------------------
# Compute partial correlations
# ---------------------------------------------------------------------
results = []
for g in groups:
subdf = df_use[df_use["NPDx1"] == g].copy()
for marker_col, marker_label in pathology_features:
for s in structures:
post_col = f"postmortem_{s}_norm"
prefix = region_map[s]
# pathology columns beginning with prefix (EC_CS_DG, Amyg, CP, TS, GP)
path_candidates = [
c for c in subdf.columns
if c.lower().startswith(prefix.lower())
]
# match Gliosis / NeuronLoss case-insensitive
path_cols = [
c for c in path_candidates
if marker_col.lower() in c.lower()
]
if not path_cols or post_col not in subdf.columns:
continue
path_col = path_cols[0]
cols = [path_col, post_col] + covars
d = subdf[cols].apply(pd.to_numeric, errors="coerce").dropna()
if len(d) < 5:
continue
res = pg.partial_corr(
data=d, x=path_col, y=post_col,
covar=covars, method="spearman"
)
results.append({
"Group": g,
"Pathology": marker_label,
"Structure": s,
"r": res["r"].iloc[0],
"p": res["p-val"].iloc[0]
})
res_df = pd.DataFrame(results)
# ---------------------------------------------------------------------
# FDR correction within each group × pathology
# ---------------------------------------------------------------------
res_df["p_FDR"] = np.nan
for g in groups:
for marker_label in [m[1] for m in pathology_features]:
sub = res_df[(res_df["Group"] == g) &
(res_df["Pathology"] == marker_label)]
if len(sub) > 0:
_, p_corr = pg.multicomp(sub["p"], method="fdr_bh")
res_df.loc[sub.index, "p_FDR"] = p_corr
res_df["sig"] = res_df["p_FDR"].apply(
lambda p: "***" if p < 0.001 else
"**" if p < 0.01 else
"*" if p < 0.05 else ""
)
# ---------------------------------------------------------------------
# Plot two side-by-side heatmaps
# ---------------------------------------------------------------------
fig, axes = plt.subplots(1, 2, figsize=(13.2, 6.8), sharey=True)
cmap = sns.color_palette("PRGn", as_cmap=True)
vmin, vmax = -1, 1
for ax, marker_label in zip(axes, [m[1] for m in pathology_features]):
sub = res_df[res_df["Pathology"] == marker_label]
# Pivot into structure x group matrix
pivot_r = sub.pivot(index="Structure", columns="Group", values="r")
pivot_r = pivot_r.reindex(index=structures, columns=groups)
# Build annotation text matrix
annot_df = sub.set_index(["Structure", "Group"])
annot_text = []
for s in structures:
row = []
for g in groups:
try:
r = annot_df.loc[(s, g), "r"]
p = annot_df.loc[(s, g), "p"]
sig = annot_df.loc[(s, g), "sig"]
row.append(f"{r:.2f}\n(p={p:.3f}){sig}" if not pd.isna(r) else "")
except KeyError:
row.append("")
annot_text.append(row)
sns.heatmap(
pivot_r,
annot=np.array(annot_text), fmt="",
cmap=cmap, center=0, vmin=vmin, vmax=vmax,
linewidths=1, linecolor="white",
cbar=(ax == axes[-1]),
cbar_kws={"label": "ρ (partial Spearman)"},
ax=ax,
annot_kws={"fontsize": 9, "color": "black"}
)
ax.set_title(marker_label, fontsize=12, pad=10, color="black")
ax.set_xticklabels(group_labels, rotation=0, fontsize=9)
ax.set_yticklabels([s.capitalize() for s in structures],
rotation=0, fontsize=9)
ax.set_xlabel("")
ax.set_ylabel("")
# ---------------------------------------------------------------------
# Unified axis labels
# ---------------------------------------------------------------------
fig.text(0.5, 0.035, "Disease group", ha="center",
fontsize=11, color="black")
fig.text(0.06, 0.5, "Structure", va="center",
rotation=90, fontsize=11, color="black")
# ---------------------------------------------------------------------
# Main title
# ---------------------------------------------------------------------
plt.suptitle(
"Partial Spearman correlation between postmortem MRI volumes\n"
"and regional neurodegeneration markers (Gliosis & Neuronal loss)",
fontsize=14, y=0.98
)
plt.tight_layout(rect=[0.06, 0.05, 0.95, 0.93], w_pad=1.2)
plt.savefig("heatmpa_gliosis_nl.png", dpi=600, bbox_inches="tight")
plt.show()
# %%
##########################################################################################
# FINAL — Scatter Plots for Gliosis and Neuronal Loss
# Same layout as ABeta/CERAD/Braak scatter plots
# 6 rows (structures) × 4 columns (disease groups)
# 2 figures total (Gliosis, Neuronal Loss)
##########################################################################################
##########################################################################################
# DATA PREP
##########################################################################################
df_use = df.copy()
df_use["NPDx1"] = df_use["NPDx1"].astype(str).str.strip().str.lower()
df_use["Sex"] = df_use["Sex"].astype("category").cat.codes
##########################################################################################
# STRUCTURES + NORMALIZATION
##########################################################################################
structures = ["hippocampus", "amygdala", "caudate",
"putamen", "thalamus", "pallidum"]
pretty_structure = {
"hippocampus": "Hippocampus",
"amygdala": "Amygdala",
"caudate": "Caudate",
"putamen": "Putamen",
"thalamus": "Thalamus",
"pallidum": "Pallidum"
}
for s in structures:
pm = f"postmortem_{s}"
if pm in df_use.columns:
df_use[f"{pm}_norm"] = df_use[pm] / df_use["antemortem_icv"]
##########################################################################################
# PATHOLOGY REGIONS (same map as heatmap)
##########################################################################################
region_map = {
"hippocampus": "EC_CS_DG",
"amygdala": "Amyg",
"caudate": "CP",
"putamen": "CP",
"thalamus": "TS",
"pallidum": "GP"
}
pathology_features = [
("Gliosis", "Gliosis"),
("NeuronLoss", "Neuronal loss")
]
##########################################################################################
# DISEASE GROUPS + COLORS
##########################################################################################
groups = ["alzheimer's disease", "lewy body disease",
"ftld-tdp", "tauopathies"]
pretty_group = {
"alzheimer's disease": "AD",
"lewy body disease": "LBD",
"ftld-tdp": "FTLD-TDP",
"tauopathies": "Tauopathies"
}
group_color = {
"alzheimer's disease": "#4C72B0",
"lewy body disease": "#DD8452",
"ftld-tdp": "#55A868",
"tauopathies": "#CBAF00"
}
##########################################################################################
# SUPERSCRIPT ASTERISKS FOR SIGNIFICANCE
##########################################################################################
def p_to_superscript(p_fdr):
if p_fdr < 0.001:
return "$^{***}$"
elif p_fdr < 0.01:
return "$^{**}$"
elif p_fdr < 0.05:
return "$^{*}$"
else:
return ""
##########################################################################################
# MAIN LOOP — TWO FIGURES (Gliosis, Neuronal loss)
##########################################################################################
for marker_col, marker_label in pathology_features:
fig, axes = plt.subplots(
len(structures), len(groups),
figsize=(4.2 * len(groups), 3.0 * len(structures)),
sharex=False, sharey=False
)
# For each disease group
for col, g in enumerate(groups):
gdf = df_use[df_use["NPDx1"] == g].copy()
# FIRST PASS: gather raw p-values for FDR
raw_pvals = []
for s in structures:
post_col = f"postmortem_{s}_norm"
prefix = region_map[s]
# find pathology column
path_candidates = [c for c in gdf.columns if c.lower().startswith(prefix.lower())]
path_cols = [c for c in path_candidates if marker_col.lower() in c.lower()]
if not path_cols or post_col not in gdf.columns:
raw_pvals.append(np.nan)
continue
path_col = path_cols[0]
cols = [path_col, post_col, "AgeatDeath", "Sex", "PMI", "Education"]
d = gdf[cols].apply(pd.to_numeric, errors="coerce").dropna()
if len(d) < 5:
raw_pvals.append(np.nan)
continue
res = pg.partial_corr(
data=d, x=path_col, y=post_col,
covar=["AgeatDeath", "Sex", "PMI", "Education"],
method="spearman"
)
raw_pvals.append(res["p-val"].iloc[0])
# FDR — within group × pathology
valid_mask = ~pd.isna(raw_pvals)
valid_p = np.array(raw_pvals)[valid_mask]
if len(valid_p) > 0:
_, p_fdr_valid = pg.multicomp(valid_p, method="fdr_bh")
else:
p_fdr_valid = []
# Expand back
p_fdr_full = []
idx = 0
for v in valid_mask:
if v:
p_fdr_full.append(p_fdr_valid[idx])
idx += 1
else:
p_fdr_full.append(np.nan)
# SECOND PASS: plotting
for row, s in enumerate(structures):
ax = axes[row, col]
# Remove seaborn y-label
ax.set_ylabel("")
if col == 0:
ax.set_ylabel(pretty_structure[s], fontsize=13)
post_col = f"postmortem_{s}_norm"
prefix = region_map[s]
path_candidates = [c for c in gdf.columns if c.lower().startswith(prefix.lower())]
path_cols = [c for c in path_candidates if marker_col.lower() in c.lower()]
if not path_cols or post_col not in gdf.columns:
ax.text(0.5, 0.5, "No data", ha="center", va="center")
ax.set_xticks([]); ax.set_yticks([])
continue
path_col = path_cols[0]
cols = [path_col, post_col, "AgeatDeath", "Sex", "PMI", "Education"]
d = gdf[cols].apply(pd.to_numeric, errors="coerce").dropna()
if len(d) < 5:
ax.text(0.5, 0.5, "N too small", ha="center", va="center")
ax.set_xticks([]); ax.set_yticks([])
continue
# Stats
res = pg.partial_corr(
data=d, x=path_col, y=post_col,
covar=["AgeatDeath", "Sex", "PMI", "Education"],
method="spearman"
)
r = res["r"].iloc[0]
p_raw = res["p-val"].iloc[0]
p_fdr = p_fdr_full[row]
sup = p_to_superscript(p_fdr)
# Plot
sns.regplot(
data=d,
x=path_col,
y=post_col,
scatter_kws=dict(color=group_color[g], s=50, alpha=0.75),
line_kws=dict(color=group_color[g], lw=2.3, alpha=0.9),
ax=ax
)
# Clean axis after regplot overwrites
ax.set_ylabel("")
if col == 0:
ax.set_ylabel(pretty_structure[s], fontsize=13)
# X-label only bottom row
if row == len(structures)-1:
ax.set_xlabel(marker_label, fontsize=12)
else:
ax.set_xlabel("")
ax.set_title(
f"{pretty_group[g]} (ρ = {r:.2f}; p = {p_raw:.3f}{sup})",
fontsize=12, pad=6
)
ax.grid(True, linestyle=":", alpha=0.5)
# FIGURE TITLE
plt.suptitle(
f"Partial Spearman correlations — {marker_label}",
fontsize=18, weight="bold"
)
plt.tight_layout(rect=[0, 0, 1, 0.95])
# SAVE PNG
png_name = f"Scatter_{marker_col}.png"
plt.savefig(png_name, dpi=300, bbox_inches="tight")
print(f"SAVED PNG: {png_name}")
# SAVE PDF
pdf_name = f"Scatter_{marker_col}.pdf"
with PdfPages(pdf_name) as pdf:
pdf.savefig(fig, bbox_inches="tight")
print(f"SAVED PDF: {pdf_name}")
plt.show()
print("\n✔✔✔ ALL SCATTER FIGURES DONE")
# %%
# %%
##########################################################################################
# OLS heatmaps (Gliosis & Neuronal loss): STANDARDIZED β as color + annotation
# Annotation: β_std + FDR stars, raw p shown
# df_use version — NO regrouping, NO relabeling, PLOT ONLY
##########################################################################################
import pandas as pd, numpy as np, seaborn as sns, matplotlib.pyplot as plt
import statsmodels.api as sm
import pingouin as pg, warnings, matplotlib
matplotlib.rcParams['font.family'] = 'DejaVu Sans'
warnings.filterwarnings("ignore", category=RuntimeWarning)
np.seterr(divide='ignore', invalid='ignore')
# ---------------------------------------------------------------------
# Use your cleaned df exactly as-is
# ---------------------------------------------------------------------
df_use = df.copy()
df_use["NPDx1"] = df_use["NPDx1"].astype(str).str.strip().str.lower()
df_use["Sex"] = df_use["Sex"].astype("category").cat.codes
# ---------------------------------------------------------------------
# Structures, groups, pathology features
# ---------------------------------------------------------------------
structures = ["hippocampus", "amygdala", "caudate", "putamen", "thalamus", "pallidum"]
region_map = {
"caudate": "CP",
"putamen": "CP",
"thalamus": "TS",
"pallidum": "GP",
"hippocampus": "EC_CS_DG",
"amygdala": "Amyg",
}
groups = ["alzheimer's disease", "lewy body disease", "ftld-tdp", "tauopathies"]
group_labels = [
"AD",
"LBD",
"FTLD-TDP",
"FTLD-Tau"
]
pathology_features = [
("Gliosis", "Gliosis"),
("NeuronLoss", "Neuronal loss")
]
covars = ["AgeatDeath", "Sex", "PMI", "Education"]
# ---------------------------------------------------------------------
# Normalize postmortem volumes by ICV
# ---------------------------------------------------------------------
if "antemortem_icv" not in df_use.columns:
raise ValueError("Missing 'antemortem_icv' column!")
for s in structures:
pm_col = f"postmortem_{s}"
if pm_col in df_use.columns:
df_use[f"{pm_col}_norm"] = df_use[pm_col] / df_use["antemortem_icv"]
sns.set(style="white")
# ---------------------------------------------------------------------
# Fit OLS per group × pathology × structure
# RAW model is optional; we compute it but PLOT standardized β
# Standardized β computed by z-scoring y and x (pathology) within subset d
# ---------------------------------------------------------------------
results = []
for g in groups:
subdf = df_use[df_use["NPDx1"] == g].copy()
for marker_col, marker_label in pathology_features:
for s in structures:
y_col = f"postmortem_{s}_norm"
prefix = region_map[s]
# candidates: columns beginning with prefix
path_candidates = [c for c in subdf.columns if c.lower().startswith(prefix.lower())]
# match marker substring case-insensitive
path_cols = [c for c in path_candidates if marker_col.lower() in c.lower()]
if (not path_cols) or (y_col not in subdf.columns):
continue
x_col = path_cols[0]
cols = [x_col, y_col] + covars
d = subdf[cols].apply(pd.to_numeric, errors="coerce").dropna()
if len(d) < 6:
continue
# ---- RAW fit (kept for reference/debug) ----
X = sm.add_constant(d[[x_col] + covars], has_constant="add")
y = d[y_col].astype(float)
# ---- Standardized fit (what we plot) ----
d_std = d.copy()
x_sd = d_std[x_col].std(ddof=0)
y_sd = d_std[y_col].std(ddof=0)
beta_raw, p_raw, beta_std = np.nan, np.nan, np.nan
try:
fit = sm.OLS(y, X).fit()
beta_raw = float(fit.params.get(x_col, np.nan))
p_raw = float(fit.pvalues.get(x_col, np.nan))
except Exception:
pass
# standardize only if variability exists
if (x_sd is not None) and (y_sd is not None) and (x_sd > 0) and (y_sd > 0):
d_std[x_col] = (d_std[x_col] - d_std[x_col].mean()) / x_sd
d_std[y_col] = (d_std[y_col] - d_std[y_col].mean()) / y_sd
X_std = sm.add_constant(d_std[[x_col] + covars], has_constant="add")
y_std = d_std[y_col].astype(float)
try:
fit_std = sm.OLS(y_std, X_std).fit()
beta_std = float(fit_std.params.get(x_col, np.nan))
except Exception:
beta_std = np.nan
results.append({
"Group": g,
"Pathology": marker_label,
"Structure": s,
"PathCol": x_col,
"beta": beta_raw,
"beta_std": beta_std, # <-- PLOTTED
"p": p_raw, # <-- RAW p (for display + FDR)
"n": int(len(d))
})
res_df = pd.DataFrame(results)
# ---------------------------------------------------------------------
# FDR correction within each group × pathology (across structures)
# Stars from FDR; raw p printed
# ---------------------------------------------------------------------
res_df["p_FDR"] = np.nan
for g in groups:
for marker_label in [m[1] for m in pathology_features]:
sub = res_df[(res_df["Group"] == g) &
(res_df["Pathology"] == marker_label) &
res_df["p"].notna()]
if len(sub) > 0:
_, p_corr = pg.multicomp(sub["p"].values, method="fdr_bh")
res_df.loc[sub.index, "p_FDR"] = p_corr
res_df["sig"] = res_df["p_FDR"].apply(
lambda p: "***" if pd.notna(p) and p < 0.001 else
"**" if pd.notna(p) and p < 0.01 else
"*" if pd.notna(p) and p < 0.05 else ""
)
# ---------------------------------------------------------------------
# Plot two side-by-side heatmaps (color = STANDARDIZED β)
# Annotation: β_std + FDR stars, raw p shown
# ---------------------------------------------------------------------
fig, axes = plt.subplots(1, 2, figsize=(13.2, 6.8), sharey=True)
cmap = sns.color_palette("PRGn", as_cmap=True)
# symmetric limits based on standardized betas
betas_std = res_df["beta_std"].replace([np.inf, -np.inf], np.nan).dropna()
if len(betas_std) == 0:
raise ValueError("No standardized betas to plot (beta_std is empty). Check model fits / variance.")
v = float(np.nanmax(np.abs(betas_std.values)))
vmin, vmax = -v, v
for ax, marker_label in zip(axes, [m[1] for m in pathology_features]):
sub = res_df[res_df["Pathology"] == marker_label].copy()
# heatmap values = standardized β
pivot_b = sub.pivot(index="Structure", columns="Group", values="beta_std")
pivot_b = pivot_b.reindex(index=structures, columns=groups)
annot_df = sub.set_index(["Structure", "Group"])
annot_text = []
for s in structures:
row = []
for g in groups:
try:
bs = annot_df.loc[(s, g), "beta_std"]
p = annot_df.loc[(s, g), "p"] # RAW p
sig = annot_df.loc[(s, g), "sig"] # from FDR
if pd.isna(bs) or pd.isna(p):
row.append("")
else:
row.append(f"β={bs:.2f}{sig}\n(p={p:.3f})")
except KeyError:
row.append("")
annot_text.append(row)
sns.heatmap(
pivot_b,
annot=np.array(annot_text), fmt="",
cmap=cmap, center=0, vmin=-1.0, vmax=1.0,
linewidths=1, linecolor="white",
cbar=(ax == axes[-1]),
cbar_kws={
"label": "Std β"
},
ax=ax,
annot_kws={"fontsize": 16, "color": "black"}
)
ax.set_title(marker_label, fontsize=12, pad=10, color="black")
ax.set_xticklabels(group_labels, rotation=0, fontsize=9)
ax.set_yticklabels([s.capitalize() for s in structures], rotation=0, fontsize=9)
ax.set_xlabel("")
ax.set_ylabel("")
# unified labels
#fig.text(0.5, 0.035, "Disease group", ha="center", fontsize=12, color="black")
#fig.text(0.06, 0.5, "Structure", va="center", rotation=90, fontsize=11, color="black")
plt.suptitle(
"OLS association between postmortem MRI volumes (ICV-normalized)",
fontsize=13, y=0.98
)
plt.tight_layout(rect=[0.06, 0.05, 0.95, 0.93], w_pad=1.2)
plt.savefig("heatmap_OLS_stdBeta_gliosis_nl.png", dpi=600, bbox_inches="tight")
plt.show()
# %%
##########################################################################################
# DOCX TABLE EXPORT — GLIOSIS & NEURONAL LOSS
##########################################################################################
# ---------------------------------------------------------------------
# GLOBAL SETTINGS
# ---------------------------------------------------------------------
matplotlib.rcParams["font.family"] = "DejaVu Sans"
warnings.filterwarnings("ignore", category=RuntimeWarning)
np.seterr(divide="ignore", invalid="ignore")
# ---------------------------------------------------------------------
# CLEAN INPUT DATA
# ---------------------------------------------------------------------
df_use = df.copy()
df_use["NPDx1"] = df_use["NPDx1"].astype(str).str.strip().str.lower()
df_use["Sex"] = df_use["Sex"].astype("category").cat.codes
# ---------------------------------------------------------------------
# STRUCTURES / GROUPS / PATHOLOGY FEATURES
# ---------------------------------------------------------------------
structures = [
"hippocampus",
"amygdala",
"caudate",
"putamen",
"thalamus",
"pallidum"
]
region_map = {
"caudate": "CP",
"putamen": "CP",
"thalamus": "TS",
"pallidum": "GP",
"hippocampus": "EC_CS_DG",
"amygdala": "Amyg",
}
groups = [
"alzheimer's disease",
"lewy body disease",
"ftld-tdp",
"tauopathies"
]
group_labels = {
"alzheimer's disease": "AD",
"lewy body disease": "LBD",
"ftld-tdp": "FTLD-TDP",
"tauopathies": "FTLD-Tau"
}
pathology_features = [
("Gliosis", "Gliosis"),
("NeuronLoss", "Neuronal loss")
]
covars = [
"AgeatDeath",
"Sex",
"PMI",
"Education"
]
# ---------------------------------------------------------------------
# NORMALIZE POSTMORTEM VOLUMES BY ICV
# ---------------------------------------------------------------------
if "antemortem_icv" not in df_use.columns:
raise ValueError("Missing 'antemortem_icv' column!")
for s in structures:
pm_col = f"postmortem_{s}"
if pm_col in df_use.columns:
df_use[f"{pm_col}_norm"] = df_use[pm_col] / df_use["antemortem_icv"]
# ---------------------------------------------------------------------
# PRETTY NAMES
# ---------------------------------------------------------------------
struct_pretty = {
"hippocampus": "Hippocampus",
"amygdala": "Amygdala",
"caudate": "Caudate",
"putamen": "Putamen",
"thalamus": "Thalamus",
"pallidum": "Pallidum"
}
pretty_cov = {
"AgeatDeath": "Age at death",
"Sex": "Sex",
"PMI": "PMI",
"Education": "Education"
}
def pretty_param_from_pathcol(pathcol):
x = str(pathcol)
region_pretty = {
"EC_CS_DG": "EC/CS/DG",
"Amyg": "Amygdala",
"CP": "Caudate/Putamen",
"TS": "Thalamus",
"GP": "Pallidum",
}
marker_pretty = {
"Gliosis": "Gliosis",
"NeuronLoss": "Neuronal loss"
}
region_name = None
marker_name = None
for rp in region_pretty:
if x.lower().startswith(rp.lower()):
region_name = region_pretty[rp]
break
for mk in marker_pretty:
if mk.lower() in x.lower():
marker_name = marker_pretty[mk]
break
if region_name is not None and marker_name is not None:
return f"{marker_name} ({region_name})"
return x
def fdr_star(p):
if pd.isna(p):
return ""
if p < 0.001:
return "***"
if p < 0.01:
return "**"
if p < 0.05:
return "*"
return ""
# ---------------------------------------------------------------------
# Z-SCORE HELPER
# ---------------------------------------------------------------------
def zscore_subset(df_in, cols):
out = df_in.copy()
for c in cols:
sd = out[c].std(ddof=0)
if pd.isna(sd) or sd == 0:
out[c] = np.nan
else:
out[c] = (out[c] - out[c].mean()) / sd
return out
# ---------------------------------------------------------------------
# FIT OLS FOR EACH GROUP × PATHOLOGY × STRUCTURE
# ---------------------------------------------------------------------
results_all = []
MIN_N = 6
for g in groups:
subdf = df_use[df_use["NPDx1"] == g].copy()
for marker_col, marker_label in pathology_features:
for s in structures:
y_col = f"postmortem_{s}_norm"
prefix = region_map[s]
path_candidates = [
c for c in subdf.columns
if c.lower().startswith(prefix.lower())
]
path_cols = [
c for c in path_candidates
if marker_col.lower() in c.lower()
]
if (not path_cols) or (y_col not in subdf.columns):
continue
x_col = path_cols[0]
cols = [x_col, y_col] + covars
d = subdf[cols].apply(pd.to_numeric, errors="coerce").dropna()
if len(d) < MIN_N:
continue
dz = zscore_subset(d, [x_col, y_col] + covars).dropna()
if len(dz) < MIN_N:
continue
X = sm.add_constant(dz[[x_col] + covars], has_constant="add")
y = dz[y_col].astype(float)
try:
fit = sm.OLS(y, X).fit()
except Exception:
continue
# ---------------------------------------------------------
# Primary pathology row
# ---------------------------------------------------------
coef = fit.params.get(x_col, np.nan)
se = fit.bse.get(x_col, np.nan)
ci_low = coef - 1.96 * se
ci_high = coef + 1.96 * se
p_raw = fit.pvalues.get(x_col, np.nan)
results_all.append({
"Group": g,
"Pathology": marker_label,
"Structure": s,
"Parameter": x_col,
"Type": "Primary",
"StdBeta": float(coef),
"SE": float(se),
"CI_low": float(ci_low),
"CI_high": float(ci_high),
"p_raw": float(p_raw),
"N": int(len(dz))
})
# ---------------------------------------------------------
# Covariate rows
# ---------------------------------------------------------
for cov in covars:
coef = fit.params.get(cov, np.nan)
se = fit.bse.get(cov, np.nan)
ci_low = coef - 1.96 * se
ci_high = coef + 1.96 * se
p_raw = fit.pvalues.get(cov, np.nan)
results_all.append({
"Group": g,
"Pathology": marker_label,
"Structure": s,
"Parameter": cov,
"Type": "Covariate",
"StdBeta": float(coef),
"SE": float(se),
"CI_low": float(ci_low),
"CI_high": float(ci_high),
"p_raw": float(p_raw),
"N": int(len(dz))
})
res_df = pd.DataFrame(results_all)
if res_df.empty:
raise ValueError("No results were generated. Check input columns, group labels, and minimum sample size.")
# ---------------------------------------------------------------------
# FDR CORRECTION FOR ALL ROWS
#
# Primary rows corrected together within group × pathology.
# Covariate rows corrected together within group × pathology.
#
# This gives p(FDR) for every row.
# ---------------------------------------------------------------------
res_df["p_FDR"] = np.nan
for g in groups:
for marker_label in [m[1] for m in pathology_features]:
for row_type in ["Primary", "Covariate"]:
idx = (
(res_df["Group"] == g) &
(res_df["Pathology"] == marker_label) &
(res_df["Type"] == row_type) &
res_df["p_raw"].notna()
)
pvals = res_df.loc[idx, "p_raw"].values
if len(pvals) > 0:
_, p_corr = pg.multicomp(pvals, method="fdr_bh")
res_df.loc[idx, "p_FDR"] = p_corr
res_df["Sig(FDR)"] = res_df["p_FDR"].apply(fdr_star)
# ---------------------------------------------------------------------
# PRETTY TABLE
# ---------------------------------------------------------------------
table_df = res_df.copy()
table_df["GroupPrint"] = table_df["Group"].map(group_labels)
table_df["StructurePrint"] = table_df["Structure"].map(struct_pretty)
table_df["ParameterPrint"] = table_df["Parameter"].apply(
lambda x: pretty_cov[x] if x in pretty_cov else pretty_param_from_pathcol(x)
)
table_df["StdBeta_fmt"] = table_df["StdBeta"].map(lambda x: f"{x:.3f}")
table_df["SE_fmt"] = table_df["SE"].map(lambda x: f"{x:.3f}")
table_df["CI_fmt"] = table_df.apply(
lambda r: f"[{r.CI_low:.3f}, {r.CI_high:.3f}]",
axis=1
)
table_df["p_fmt"] = table_df["p_raw"].map(lambda x: f"{x:.4f}")
table_df["pFDR_fmt"] = table_df["p_FDR"].map(lambda x: f"{x:.4f}" if pd.notna(x) else "")
table_df["N_fmt"] = table_df["N"].astype(int).astype(str)
final_table = table_df[[
"GroupPrint",
"Pathology",
"StructurePrint",
"ParameterPrint",
"Type",
"StdBeta_fmt",
"SE_fmt",
"CI_fmt",
"p_fmt",
"pFDR_fmt",
"Sig(FDR)",
"N_fmt"
]].sort_values([
"GroupPrint",
"Pathology",
"StructurePrint",
"Type",
"ParameterPrint"
])
final_table = final_table.rename(columns={
"GroupPrint": "Group",
"StructurePrint": "Structure",
"ParameterPrint": "Predictor",
"StdBeta_fmt": "β",
"SE_fmt": "SE",
"CI_fmt": "CI (95%)",
"p_fmt": "p(raw)",
"pFDR_fmt": "p(FDR)",
"N_fmt": "N"
})
# ---------------------------------------------------------------------
# DOCX EXPORT — ONE TABLE PER GROUP
# ---------------------------------------------------------------------
doc = Document()
title = doc.add_heading(
"Gliosis and Neuronal Loss → Postmortem Volume (Standardized OLS)",
level=1
)
title.alignment = WD_ALIGN_PARAGRAPH.CENTER
doc.add_paragraph(
"Each table reports standardized β, standard error (SE), 95% confidence interval (CI), "
"raw p-value, FDR-corrected p-value, significance code, and sample size. "
"FDR correction is applied to all rows, including covariates, within each disease group "
"× pathology × row-type block."
)
doc.add_page_break()
group_order = [
"alzheimer's disease",
"lewy body disease",
"ftld-tdp",
"tauopathies"
]
for g in group_order:
gp = group_labels[g]
sub = final_table[final_table["Group"] == gp].copy()
if sub.empty:
continue
h = doc.add_heading(gp, level=2)
h.alignment = WD_ALIGN_PARAGRAPH.LEFT
tbl = doc.add_table(rows=1, cols=len(sub.columns))
tbl.style = "Table Grid"
# Header row
hdr = tbl.rows[0].cells
for j, colname in enumerate(sub.columns):
hdr[j].text = str(colname)
# Body rows
for _, row in sub.iterrows():
rw = tbl.add_row().cells
for j, colname in enumerate(sub.columns):
rw[j].text = str(row[colname])
doc.add_page_break()
docx_output = "gliosis_neuronloss_OLS_results_by_group_WITH_SE_AND_ALL_FDR.docx"
doc.save(docx_output)
print("Saved DOCX:", docx_output)
# %%
##########################################################################################
# Postmortem-only Monte Carlo Mediation Analysis
##########################################################################################
import pandas as pd, numpy as np, seaborn as sns, matplotlib.pyplot as plt, pingouin as pg
from statsmodels.formula.api import ols
import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)
np.seterr(divide='ignore', invalid='ignore')
# ---------------------------------------------------------------------
# Use your cleaned dataframe exactly as-is
# ---------------------------------------------------------------------
df_use = df.copy()
# Ensure clean NPDx1 and numeric postmortem volumes
df_use["NPDx1"] = df_use["NPDx1"].astype(str).str.strip().str.lower()
for c in df_use.columns:
if c.startswith("postmortem_"):
df_use[c] = pd.to_numeric(df_use[c], errors="coerce")
# ---------------------------------------------------------------------
# Settings
# ---------------------------------------------------------------------
disease_groups = ["alzheimer's disease", "ftld-tdp", "lewy body disease", "tauopathies"]
region_map = {
"hippocampus": "EC_CS_DG",
"amygdala": "Amyg",
"caudate": "CP",
"putamen": "CP",
"thalamus": "TS",
"pallidum": "GP"
}
covars = ["AgeatDeath", "Sex", "Education", "PMI"]
mediators = ["NeuronLoss", "Gliosis"]
# ---------------------------------------------------------------------
# Monte Carlo bootstrap function
# ---------------------------------------------------------------------
def montecarlo_indirect(df_in, path_col, med_col, vol_col, covars, n_iter=5000, seed=42):
np.random.seed(seed)
try:
# a-path
a_mod = ols(f"{med_col} ~ {path_col} + {' + '.join(covars)}", data=df_in).fit()
# b-path
b_mod = ols(f"{vol_col} ~ {path_col} + {med_col} + {' + '.join(covars)}", data=df_in).fit()
# total effect
c_mod = ols(f"{vol_col} ~ {path_col} + {' + '.join(covars)}", data=df_in).fit()
a = a_mod.params.get(path_col, np.nan)
b = b_mod.params.get(med_col, np.nan)
c_prime = b_mod.params.get(path_col, np.nan)
c_total = c_mod.params.get(path_col, np.nan)
a_draws = np.random.normal(a, a_mod.bse.get(path_col, np.nan), n_iter)
b_draws = np.random.normal(b, b_mod.bse.get(med_col, np.nan), n_iter)
ab_samples = a_draws * b_draws
indirect = np.mean(ab_samples)
ci_low, ci_high = np.percentile(ab_samples, [2.5, 97.5])
# two-sided Monte Carlo p-value
p_val = 2 * min(np.mean(ab_samples < 0), np.mean(ab_samples > 0))
prop_med = (indirect / c_total) if c_total not in [0, np.nan] else np.nan
return c_prime, indirect, prop_med, ci_low, ci_high, p_val
except Exception:
return [np.nan]*6
# ---------------------------------------------------------------------
# Run postmortem mediation analysis
# ---------------------------------------------------------------------
results = []
for dx in disease_groups:
gdf = df_use[df_use["NPDx1"] == dx].copy()
if gdf.empty:
continue
# pathology suffix selection
if dx == "alzheimer's disease":
suffix = "Tau"
elif dx == "ftld-tdp":
suffix = "TDP43"
elif dx == "lewy body disease":
suffix = "aSyn"
elif dx == "tauopathies":
suffix = "Tau"
# loop structures
for s, prefix in region_map.items():
path_col = f"{prefix}{suffix}"
for med in mediators:
med_col = f"{prefix}{med}"
vol_col = f"postmortem_{s}"
# must exist
if not (path_col in gdf.columns and med_col in gdf.columns and vol_col in gdf.columns):
continue
cols = [path_col, med_col, vol_col] + covars
d = gdf[cols].apply(pd.to_numeric, errors="coerce").dropna()
if len(d) < 15:
continue
# z-score all continuous variables
for c in cols:
if d[c].std(ddof=0) > 0:
d[c] = (d[c] - d[c].mean()) / d[c].std(ddof=0)
cprime, indirect, prop, cil, cih, pv = montecarlo_indirect(
d, path_col, med_col, vol_col, covars
)
if np.isnan(indirect):
continue
results.append({
"Disease": dx,
"Structure": s,
"Mediator": med,
"Direct(c')": cprime,
"Indirect(a*b)": indirect,
"PropMediated": prop,
"CI_low": cil,
"CI_high": cih,
"p_val": pv,
"N": len(d)
})
# ---------------------------------------------------------------------
# Assemble dataframe + FDR by disease group
# ---------------------------------------------------------------------
results_df = pd.DataFrame(results)
if results_df.empty:
print("❌ No valid postmortem mediation results found.")
else:
results_df["p_FDR"] = np.nan
for dx in disease_groups:
sub = results_df[results_df["Disease"] == dx]
if sub.empty:
continue
reject, p_corr = pg.multicomp(sub["p_val"], method="fdr_bh")
results_df.loc[sub.index, "p_FDR"] = p_corr
# significance annotation
results_df["Sig(FDR)"] = results_df["p_FDR"].apply(
lambda p: "***" if p < 0.001 else
"**" if p < 0.01 else
"*" if p < 0.05 else ""
)
print("\n✅ Postmortem mediation results (FDR corrected within each group):\n")
print(results_df.round(4))
# %%
################################################################################
# Postmortem Mediation Analysis + DOCX Supplementary Table
################################################################################
import numpy as np
import pandas as pd
import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)
from statsmodels.formula.api import ols
import pingouin as pg
from docx import Document
from docx.enum.table import WD_TABLE_ALIGNMENT
from docx.oxml import OxmlElement
from docx.oxml.ns import qn
np.seterr(divide="ignore", invalid="ignore")
# ============================================================================
# 1) MEDIATION ANALYSIS
# ============================================================================
df_use = df.copy()
# Clean group label and postmortem columns
df_use["NPDx1"] = df_use["NPDx1"].astype(str).str.strip().str.lower()
for c in df_use.columns:
if c.startswith("postmortem_"):
df_use[c] = pd.to_numeric(df_use[c], errors="coerce")
disease_groups = ["alzheimer's disease", "ftld-tdp",
"lewy body disease", "tauopathies"]
region_map = {
"hippocampus": "EC_CS_DG",
"amygdala": "Amyg",
"caudate": "CP",
"putamen": "CP",
"thalamus": "TS",
"pallidum": "GP"
}
covars = ["AgeatDeath", "Sex", "Education", "PMI"]
mediators = ["NeuronLoss", "Gliosis"]
# --------------------------- Monte Carlo helper -----------------------------
def montecarlo_indirect(df_in, path_col, med_col, vol_col, covars,
n_iter=5000, seed=42):
"""
Return:
c_prime, p_direct,
indirect, p_MC,
prop_med, ci_low, ci_high
"""
np.random.seed(seed)
# a-path
a_mod = ols(f"{med_col} ~ {path_col} + {' + '.join(covars)}", data=df_in).fit()
# b-path (direct + mediator)
b_mod = ols(f"{vol_col} ~ {path_col} + {med_col} + {' + '.join(covars)}",
data=df_in).fit()
# total effect
c_mod = ols(f"{vol_col} ~ {path_col} + {' + '.join(covars)}", data=df_in).fit()
a = a_mod.params.get(path_col, np.nan)
b = b_mod.params.get(med_col, np.nan)
c_prime = b_mod.params.get(path_col, np.nan)
p_direct = b_mod.pvalues.get(path_col, np.nan)
c_total = c_mod.params.get(path_col, np.nan)
# MC samples for a*b
a_se = a_mod.bse.get(path_col, np.nan)
b_se = b_mod.bse.get(med_col, np.nan)
if np.isnan(a) or np.isnan(b) or np.isnan(a_se) or np.isnan(b_se):
return [np.nan]*7
a_draws = np.random.normal(a, a_se, n_iter)
b_draws = np.random.normal(b, b_se, n_iter)
ab_samples = a_draws * b_draws
indirect = np.mean(ab_samples)
ci_low, ci_high = np.percentile(ab_samples, [2.5, 97.5])
# two-sided MC p-value
p_MC = 2 * min(np.mean(ab_samples < 0), np.mean(ab_samples > 0))
prop_med = (indirect / c_total) if (c_total not in [0, np.nan]) else np.nan
return c_prime, p_direct, indirect, p_MC, prop_med, ci_low, ci_high
# --------------------------- Run mediation -----------------------------
results = []
for dx in disease_groups:
gdf = df_use[df_use["NPDx1"] == dx].copy()
if gdf.empty:
continue
if dx == "alzheimer's disease":
suffix = "Tau"
elif dx == "ftld-tdp":
suffix = "TDP43"
elif dx == "lewy body disease":
suffix = "aSyn"
elif dx == "tauopathies":
suffix = "Tau"
for s, prefix in region_map.items():
path_col = f"{prefix}{suffix}"
vol_col = f"postmortem_{s}"
for med in mediators:
med_col = f"{prefix}{med}"
if not (path_col in gdf.columns and
med_col in gdf.columns and
vol_col in gdf.columns):
continue
cols = [path_col, med_col, vol_col] + covars
d = gdf[cols].apply(pd.to_numeric, errors="coerce").dropna()
if len(d) < 15:
continue
# z-score continuous variables
for c in cols:
if d[c].std(ddof=0) > 0:
d[c] = (d[c] - d[c].mean()) / d[c].std(ddof=0)
(c_prime, p_direct,
indirect, p_MC,
prop, cil, cih) = montecarlo_indirect(
d, path_col, med_col, vol_col, covars
)
if np.isnan(indirect):
continue
results.append({
"Disease": dx,
"Structure": s,
"Mediator": med,
"Direct(c')": c_prime,
"p_direct": p_direct,
"Indirect(a*b)": indirect,
"p_MC": p_MC,
"PropMediated": prop,
"CI_low": cil,
"CI_high": cih,
"N": len(d)
})
results_df = pd.DataFrame(results)
if results_df.empty:
raise RuntimeError("No valid mediation results computed.")
# --------------------------- FDR corrections -----------------------------
results_df["p_FDR_direct"] = np.nan
results_df["p_FDR_MC"] = np.nan
for dx in disease_groups:
sub = results_df[results_df["Disease"] == dx]
if sub.empty:
continue
# FDR for direct effects
rej_d, p_corr_d = pg.multicomp(sub["p_direct"], method="fdr_bh")
results_df.loc[sub.index, "p_FDR_direct"] = p_corr_d
# FDR for indirect (MC) effects
rej_i, p_corr_i = pg.multicomp(sub["p_MC"], method="fdr_bh")
results_df.loc[sub.index, "p_FDR_MC"] = p_corr_i
# significance codes
def sig_star(p):
if p < 0.001:
return "***"
elif p < 0.01:
return "**"
elif p < 0.05:
return "*"
else:
return ""
results_df["Sig_direct"] = results_df["p_FDR_direct"].apply(sig_star)
results_df["Sig_MC"] = results_df["p_FDR_MC"].apply(sig_star)
print("Mediation results (head):")
print(results_df.head())
# ============================================================================
# 2) BUILD DOCX — mediation_analyses.docx
# ============================================================================
group_display = {
"alzheimer's disease": "Alzheimer’s disease",
"lewy body disease": "Lewy body disease",
"ftld-tdp": "FTLD-TDP",
"tauopathies": "Tauopathies"
}
mediator_display = {
"Gliosis": "Gliosis",
"NeuronLoss": "Neuronal loss"
}
structures_order = [
"hippocampus", "amygdala", "caudate",
"putamen", "thalamus", "pallidum"
]
# ---------- small helpers ----------
def sci(x):
try:
return f"{float(x):.2e}"
except Exception:
return ""
def add_bold(cell, text):
run = cell.paragraphs[0].add_run(text)
run.bold = True
def add_superscript(cell, text):
run = cell.paragraphs[0].add_run(text)
run.font.superscript = True
def remove_top_border(cell):
"""
Remove ONLY the top border of a cell (to hide line
between Gliosis and Neuronal loss rows).
"""
tc = cell._tc
tcPr = tc.get_or_add_tcPr()
tcBorders = tcPr.find(qn('w:tcBorders'))
if tcBorders is None:
tcBorders = OxmlElement('w:tcBorders')
tcPr.append(tcBorders)
top = tcBorders.find(qn('w:top'))
if top is None:
top = OxmlElement('w:top')
tcBorders.append(top)
top.set(qn('w:val'), 'nil') # no border
# ---------- create doc ----------
doc = Document()
doc.add_heading("Supplementary Table — Postmortem Mediation Analyses", level=1)
intro = (
"For each diagnostic group, Monte Carlo mediation analyses were performed "
"with regional pathology as the predictor, gliosis or neuronal loss as "
"mediators, and postmortem MRI volume as the outcome. Models adjust for "
"age at death, sex, education, and PMI. FDR correction was applied within "
"each diagnostic group separately for direct effects and indirect "
"mediation effects. Asterisks (*, **, ***) indicate FDR-corrected "
"significance levels."
)
doc.add_paragraph(intro)
doc.add_paragraph("\n")
for dx in results_df["Disease"].unique():
sub = results_df[results_df["Disease"] == dx].copy()
if sub.empty:
continue
doc.add_heading(group_display.get(dx, dx), level=2)
sub["Structure"] = pd.Categorical(sub["Structure"], structures_order)
sub = sub.sort_values(["Structure", "Mediator"])
# table columns
cols = [
"Structure", "Mediator",
"Direct(c')", "p_direct", "p_FDR_direct",
"Indirect(a*b)", "p_MC", "p_FDR_MC",
"PropMediated", "CI"
]
headers = {
"Structure": "Structure",
"Mediator": "Mediator",
"Direct(c')": "Direct (c′)",
"p_direct": "p(c′)",
"p_FDR_direct": "p_FDR(c′)",
"Indirect(a*b)": "Indirect (a×b)",
"p_MC": "p(MC)",
"p_FDR_MC": "p_FDR(MC)",
"PropMediated": "Proportion mediated",
"CI": "95% CI [low, high]"
}
table = doc.add_table(rows=1, cols=len(cols))
table.style = "Table Grid"
table.alignment = WD_TABLE_ALIGNMENT.CENTER
hdr = table.rows[0].cells
for j, c in enumerate(cols):
add_bold(hdr[j], headers[c])
# fill rows: structures × mediators (Gliosis, NeuronLoss)
for s in structures_order:
s_sub = sub[sub["Structure"] == s]
if s_sub.empty:
continue
for i, med in enumerate(["Gliosis", "NeuronLoss"]):
rowdata = s_sub[s_sub["Mediator"] == med]
if rowdata.empty:
continue
r = rowdata.iloc[0]
cells = table.add_row().cells
# structure label only on first mediator row
if i == 0:
cells[0].text = s.capitalize()
else:
cells[0].text = ""
# remove top border for all cells in second mediator row
for c in cells:
remove_top_border(c)
cells[1].text = mediator_display[med]
# direct effect + p
cells[2].text = sci(r["Direct(c')"])
cells[3].text = sci(r["p_direct"])
cells[4].text = sci(r["p_FDR_direct"])
if r["Sig_direct"]:
add_superscript(cells[4], r["Sig_direct"])
# indirect effect + p(MC)
cells[5].text = sci(r["Indirect(a*b)"])
cells[6].text = sci(r["p_MC"])
cells[7].text = sci(r["p_FDR_MC"])
if r["Sig_MC"]:
add_superscript(cells[7], r["Sig_MC"])
# proportion mediated + CI
cells[8].text = sci(r["PropMediated"])
cells[9].text = f"[{sci(r['CI_low'])}, {sci(r['CI_high'])}]"
doc.add_paragraph("\n")
out_name = "mediation_analyses.docx"
doc.save(out_name)
print(f"\nSaved Word file: {out_name}")
# %%
##########################################################################################
# Postmortem Mediation
##########################################################################################
import pandas as pd, numpy as np, seaborn as sns, matplotlib.pyplot as plt, pingouin as pg
from statsmodels.formula.api import ols
import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)
np.seterr(divide='ignore', invalid='ignore')
# ---------------------------------------------------------------------
# Use your already-cleaned dataframe
# ---------------------------------------------------------------------
df_use = df.copy() # <-- ADAPTED (NO re-cleaning)
# Ensure numeric postmortem volumes (safe cast only)
for c in df_use.columns:
if c.startswith("postmortem_"):
df_use[c] = pd.to_numeric(df_use[c], errors="coerce")
# ---------------------------------------------------------------------
# Settings
# ---------------------------------------------------------------------
disease_groups = ["alzheimer's disease", "ftld-tdp", "lewy body disease", "tauopathies"]
region_map = {
"Hippocampus": "EC_CS_DG",
"Amygdala": "Amyg",
"Caudate": "CP",
"Putamen": "CP",
"Thalamus": "TS",
"Pallidum": "GP"
}
covars = ["AgeatDeath", "Sex", "Education", "PMI"]
mediators = ["Gliosis", "NeuronLoss"]
# Mapping from disease → pathology suffix
path_suffix_map = {
"alzheimer's disease": "Tau",
"ftld-tdp": "TDP43",
"lewy body disease": "aSyn",
"tauopathies": "Tau"
}
# ---------------------------------------------------------------------
# Monte Carlo bootstrap indirect effect
# ---------------------------------------------------------------------
def montecarlo_indirect(df_mc, path_col, med_col, vol_col, covars, n_iter=5000, seed=42):
np.random.seed(seed)
try:
a_model = ols(f"{med_col} ~ {path_col} + {' + '.join(covars)}", data=df_mc).fit()
b_model = ols(f"{vol_col} ~ {path_col} + {med_col} + {' + '.join(covars)}", data=df_mc).fit()
c_model = ols(f"{vol_col} ~ {path_col} + {' + '.join(covars)}", data=df_mc).fit()
a = a_model.params.get(path_col, np.nan)
b = b_model.params.get(med_col, np.nan)
c_prime = b_model.params.get(path_col, np.nan)
c_total = c_model.params.get(path_col, np.nan)
a_draws = np.random.normal(a, a_model.bse.get(path_col, np.nan), n_iter)
b_draws = np.random.normal(b, b_model.bse.get(med_col, np.nan), n_iter)
ab_samples = a_draws * b_draws
indirect = np.mean(ab_samples)
ci_low, ci_high = np.percentile(ab_samples, [2.5, 97.5])
p_val = 2 * min(np.mean(ab_samples < 0), np.mean(ab_samples > 0))
prop = (indirect / c_total) if c_total not in [0, np.nan] else np.nan
return c_prime, indirect, prop, ci_low, ci_high, p_val
except:
return [np.nan] * 6
# ---------------------------------------------------------------------
# Run mediation analysis
# ---------------------------------------------------------------------
results = []
for dx in disease_groups:
gdf = df_use[df_use["NPDx1"] == dx].copy()
if gdf.empty:
continue
suffix = path_suffix_map[dx]
for s, prefix in region_map.items():
for med in mediators:
path_col = f"{prefix}{suffix}"
med_col = f"{prefix}{med}"
vol_col = f"postmortem_{s.lower()}"
if not all(c in gdf.columns for c in [path_col, med_col, vol_col]):
continue
cols = [path_col, med_col, vol_col] + covars
d = gdf[cols].apply(pd.to_numeric, errors="coerce").dropna()
if len(d) < 15:
continue
# Z-score numeric values
for c in cols:
if d[c].std(ddof=0) > 0:
d[c] = (d[c] - d[c].mean()) / d[c].std(ddof=0)
cprime, indirect, prop, cil, cih, pv = montecarlo_indirect(
d, path_col, med_col, vol_col, covars
)
if np.isnan(indirect):
continue
results.append({
"Disease": dx, "Structure": s, "Mediator": med,
"Direct(c')": cprime, "Indirect(a*b)": indirect,
"PropMediated": prop, "CI_low": cil, "CI_high": cih,
"p_val": pv, "N": len(d)
})
results_df = pd.DataFrame(results)
if results_df.empty:
raise SystemExit("❌ No valid postmortem mediation results found.")
# ---------------------------------------------------------------------
# FDR correction within each disease group
# ---------------------------------------------------------------------
results_df["p_FDR"] = np.nan
for dx in disease_groups:
sub = results_df[results_df["Disease"] == dx]
if sub.empty:
continue
_, p_corr = pg.multicomp(sub["p_val"], method="fdr_bh")
results_df.loc[sub.index, "p_FDR"] = p_corr
results_df["Sig(FDR)"] = results_df["p_FDR"].apply(
lambda p: "***" if p < 0.001 else "**" if p < 0.01 else "*" if p < 0.05 else ""
)
# ---------------------------------------------------------------------
# Publication-Ready Figure
# ---------------------------------------------------------------------
sns.set(style="white", context="talk")
effect_palette = {"Direct(c')": "#bca0dc", "Indirect(a*b)": "#a3d9a5"}
pretty_names = {
"alzheimer's disease": "Alzheimer’s Disease",
"ftld-tdp": "FTLD-TDP",
"lewy body disease": "Lewy Body Disease",
"tauopathies": "Tauopathies"
}
fig, axes = plt.subplots(2, 4, figsize=(20, 7.5),
sharey=True,
gridspec_kw={'hspace': 0.35, 'wspace': 0.10})
for row_i, med in enumerate(mediators):
for col_i, dx in enumerate(disease_groups):
ax = axes[row_i, col_i]
sub = results_df[(results_df["Mediator"] == med) &
(results_df["Disease"] == dx)]
if sub.empty:
ax.axis("off")
continue
plot_df = sub.melt(
id_vars=["Structure", "p_FDR", "Sig(FDR)"],
value_vars=["Direct(c')", "Indirect(a*b)"],
var_name="EffectType",
value_name="EffectSize"
)
# shorter x-labels
plot_df["Structure"] = plot_df["Structure"].replace({"Hippocampus": "Hippo"})
sns.barplot(
data=plot_df, x="Structure", y="EffectSize",
hue="EffectType", hue_order=["Direct(c')", "Indirect(a*b)"],
palette=effect_palette, dodge=True, edgecolor=None,
errorbar=None, ax=ax
)
ax.axhline(0, color="black", lw=0.6)
ax.tick_params(axis="x", rotation=0, labelsize=9)
ax.set_xlabel("")
ax.set_ylim(-0.65, 0.35)
if row_i == 0:
ax.set_title(pretty_names[dx], fontsize=14, weight="semibold")
# significance stars
ylim = ax.get_ylim()
for patch, (_, row) in zip(ax.patches, plot_df.iterrows()):
star = row["Sig(FDR)"]
if star:
height = patch.get_height()
y_pos = np.clip(height + (0.015 if height >= 0 else -0.015),
ylim[0] + 0.01, ylim[1] - 0.01)
ax.text(patch.get_x() + patch.get_width()/2, y_pos,
star, ha="center",
va="bottom" if height >= 0 else "top",
fontsize=11, weight="bold")
if ax.get_legend():
ax.get_legend().remove()
# Shared x/y labels
fig.text(0.5, 0.03, "Structure", ha="center", fontsize=15, weight="semibold")
fig.text(0.08, 0.5, "β Effect", va="center", ha="center",
rotation=90, fontsize=15, weight="semibold")
# Row titles
fig.text(0.5, 0.94, "Gliosis mediated effects", ha="center", fontsize=13)
fig.text(0.5, 0.46, "Neuron loss mediated effects", ha="center", fontsize=13)
# Legend
handles, labels = axes[0, 0].get_legend_handles_labels()
label_map = {"Direct(c')": "Direct", "Indirect(a*b)": "Indirect"}
labels = [label_map.get(l, l) for l in labels]
fig.legend(handles, labels, title="Effect", loc="lower center",
bbox_to_anchor=(0.5, -0.05), ncol=2,
frameon=False, fontsize=11, title_fontsize=12)
plt.suptitle("Postmortem Mediation: direct vs indirect effects", fontsize=16, y=0.99)
plt.tight_layout(rect=[0.05, 0.10, 1, 0.92])
print("✅ Figure saved as 'postmortem_mediation_final.png'")
plt.savefig("postmortem_mediation_final.png", dpi=600, bbox_inches="tight")
plt.show()
# %%
# %%
######### Antemortem box plots
# ============================================================
# PREP
# ============================================================
df_use = merged_df.copy()
df_use["NPDx1"] = df_use["NPDx1"].astype(str).str.strip().str.lower()
order = ["alzheimer's disease", "lewy body disease", "ftld-tdp", "tauopathies"]
# Pretty x-labels
x_labels = [
"AD",
"LBD",
"FTLD-TDP",
"FTLD-Tau"
]
covars = ["AgeatDeath", "Sex", "Education", "AMI"]
palette = ["#3366CC", "#DC3912", "#109618", "#FF9900"]
sns.set(style="whitegrid", context="talk", font_scale=1.2)
# ============================================================
# 7 REGIONS ONLY
# ============================================================
regions = [
"hippocampus",
"amygdala",
"accumbens_area",
"thalamus",
"caudate",
"putamen",
"pallidum"
]
regions = [r for r in regions if f"antemortem_avg_{r}" in df_use.columns]
if "antemortem_icv" not in df_use.columns:
raise ValueError("antemortem_icv column not found")
print("Using regions:", regions)
# ============================================================
# NORMALIZE BY ANTEMORTEM ICV
# ============================================================
for r in regions:
avg_col = f"antemortem_avg_{r}"
norm_col = f"{r}_norm"
df_use[norm_col] = df_use[avg_col] / df_use["antemortem_icv"]
# ============================================================
# PAIRWISE LRTs
# ============================================================
pairwise_results = []
for r in regions:
ycol = f"{r}_norm"
for g1, g2 in combinations(order, 2):
d = df_use[df_use["NPDx1"].isin([g1, g2])].copy()
covars_here = [c for c in covars if c in d.columns]
d = d[["NPDx1", ycol] + covars_here].dropna()
if len(d) < 10:
continue
reduced = ols(f"{ycol} ~ " + " + ".join(covars_here), data=d).fit()
full = ols(f"{ycol} ~ C(NPDx1) + " + " + ".join(covars_here), data=d).fit()
lr = 2 * (full.llf - reduced.llf)
df_diff = full.df_model - reduced.df_model
p = chi2.sf(lr, df_diff)
pairwise_results.append({
"Region": r,
"Group1": g1,
"Group2": g2,
"p_raw": p
})
pairwise_df = pd.DataFrame(pairwise_results)
if not pairwise_df.empty:
reject, p_corr = pg.multicomp(pairwise_df["p_raw"], method="fdr_bh")
pairwise_df["p_FDR"] = p_corr
pairwise_df["Sig"] = pairwise_df["p_FDR"].apply(
lambda p: "***" if p < 0.001 else "**" if p < 0.01 else "*" if p < 0.05 else ""
)
else:
pairwise_df = pd.DataFrame(columns=["Region", "Group1", "Group2", "p_raw", "p_FDR", "Sig"])
print(pairwise_df)
# ============================================================
# PLOTTING
# ============================================================
fig = plt.figure(figsize=(26, 18))
# TOP PANEL (3)
gsA = fig.add_gridspec(
1, 3, left=0.05, right=0.97,
top=0.92, bottom=0.56, wspace=0.33
)
axesA = [fig.add_subplot(gsA[0, k]) for k in range(3)]
# BOTTOM PANEL (4)
gsB = fig.add_gridspec(
1, 4, left=0.05, right=0.97,
top=0.50, bottom=0.12, wspace=0.30
)
axesB = [fig.add_subplot(gsB[0, k]) for k in range(4)]
axes = axesA + axesB
def plot_struct(ax, r):
ycol = f"{r}_norm"
d = df_use[["NPDx1", ycol]].dropna()
for idx, g in enumerate(order):
vals = d.loc[d["NPDx1"] == g, ycol]
tmp = pd.DataFrame({"group": [g] * len(vals), "y": vals})
sns.boxplot(
data=tmp, x="group", y="y",
color=palette[idx], ax=ax,
width=0.55, fliersize=0,
linewidth=1.3, boxprops=dict(alpha=0.72)
)
sns.stripplot(
data=tmp, x="group", y="y",
color="black", size=4, alpha=0.55,
ax=ax, jitter=0.15
)
ax.set_xlabel("")
ax.set_ylabel("")
ax.set_xticklabels(x_labels, fontsize=16, fontweight="bold")
ymin, ymax = d[ycol].min(), d[ycol].max()
yr = ymax - ymin if ymax > ymin else 1.0
ax.set_ylim(ymin - 0.06 * yr, ymax + 0.45 * yr)
ax.set_title(r.replace("_", " ").capitalize(), fontsize=18, fontweight="bold")
ax.grid(axis="y", linestyle=":", alpha=0.45)
pairs = pairwise_df[pairwise_df["Region"] == r]
y_offset = 0.020 * yr
y_pos = ymax + 0.10 * yr
for _, row in pairs.iterrows():
if row["Sig"]:
g1, g2 = row["Group1"], row["Group2"]
x1 = order.index(g1)
x2 = order.index(g2)
ax.plot([x1, x1, x2, x2],
[y_pos, y_pos + y_offset, y_pos + y_offset, y_pos],
lw=1.35, color="black")
ax.text((x1 + x2) / 2, y_pos + y_offset * 0.8,
row["Sig"], ha="center",
fontsize=14, fontweight="bold")
y_pos += y_offset * 1.9
for ax, r in zip(axes, regions):
plot_struct(ax, r)
fig.text(0.0001, 0.55, "Normalized volume",
va="center", rotation="vertical",
fontsize=25, fontweight="bold")
fig.text(0.50, 0.06, "Disease groups",
ha="center", fontsize=25, fontweight="bold")
plt.subplots_adjust(top=0.93)
fig.suptitle(
"Antemortem subcortical and limbic volumes differentiates neuropathological groups",
fontsize=26, fontweight="bold"
)
plt.tight_layout()
plt.savefig("antemortem_avg_subcortical_volumes_normalized_by_icv.png",
dpi=600, bbox_inches="tight")
plt.show()
# %%
############################
############################
"""
For the 5 groups, the groups to be useed are:
order = [
"alzheimer's disease",
"lewy body disease",
"ftld-tdp",
"3R-tau",
"4R-tau"
]
"""

subcortical_analyses.ipynb at commit 6abfd1d, no license · at the source

Overview

Authors: Pulkit Khandelwal^1,2, Michael Tran Duong¹, Lisa M Levorse¹, Winifred Trotman¹, Alejandra Bahena¹, Sydney A Lim¹, Amanda E Denning¹, Eunice Chung¹, Christopher A Olm¹, Hamsanandini Radhakrishnan¹, Ranjit Ittyerah¹, Karthik Prabhakaran¹, Gabor Mizsei¹, Theresa Schuck¹, Sheina Emrani¹, Joaquin A Vizcarra¹, John Robinson¹, Daniel T Ohm¹, Jeffrey S Phillips¹, Jesse Cohen³

and 11 other authors

Laura E M Wisse⁴, John A Detre¹, Ilya M Nasrallah¹, Christopher A Brown¹, Sandhitsu R Das¹, Edward B Lee¹, M Dylan Tisdall¹, David J Irwin¹, Corey T McMillan¹, David A Wolk¹, Paul A Yushkevich¹

ORCID iDs: Michael Tran Duong

Institutions: Harvard University (United States); Massachusetts General Hospital (United States); Athinoula A. Martinos Center for Biomedical Imaging (United States); University of Pennsylvania (United States); University of Florida (United States); Lund University (Sweden)

Journal: Alzheimer's & dementia : the journal of the Alzheimer's Association, volume 22, issue 7, article e71649

Dates: received 7 January 2026; accepted 6 June 2026; published online 8 July 2026; in print July 2026

Type: Research article · Language: English

License: CC BY-NC

Identifiers: DOI 10.1002/alz.71649 · PMID 42418450 · PMCID PMC13344898 · OpenAlex W7167693549

Open access: hybrid, a free copy (OpenAlex)

Status: code verified

Categories: structural MRI / diffusion (modality), histology / microscopy (modality), human (organism), other condition (population), Alzheimer's / dementia (population), Parkinson's (population), cellular / molecular (subfield)

Methods: Statistics, Smoothing, state filtering, decompositions, Connectivity, Preprocessing

Keywords: ADRD, alzheimer's disease, alzheimer's disease and related dementias, copathologies, entorhinal cortex, ex vivo magnetic resonance imaging (MRI), frontotemporal lobar degeneration with transactive response DNA binding protein 43 (FTLD‐TDP), FTLD‐tau, lewy body dementia, neuroimaging, neuropathologies, postmortem imaging, subcortical structures, synuclein, tau, TDP‐43

MeSH: Alzheimer Disease*, Brain*, Cerebral Cortical Thinning*, Frontotemporal Lobar Degeneration*, Lewy Body Disease*, Limbic System*, Magnetic Resonance Imaging*, Neurodegenerative Diseases*, Aged, Aged, 80 and over, alpha-Synuclein, Atrophy, DNA-Binding Proteins, Female, Humans, Male, Middle Aged, tau Proteins (* major topic)

Topic: Dementia and Cognitive Impairment Research (Psychiatry and Mental health, Medicine), according to OpenAlex

Funding: NIH HHS (P01 AG017586, F30 AG074524, P30 AG072979, U19 AG062418, P01 AG066597, R01 NS109260, RF1 AG056014, P01 AG084497, R01 AG054519, R01 AG069474); National Institutes of Health (R01 AG054519, P01 AG084497, U19 AG062418, P01 AG017586, P01 AG066597, R01 NS109260, F30 AG074524, P30 AG072979, R01 AG069474, RF1 AG056014)

Citations: not cited yet (Europe PMC); 74 references in the paper

Abstract

INTRODUCTION: The impact of different neuropathologies on deep brain structures remains to be understood. We examine subcortical and limbic volumetry in neurodegenerative diseases involving phosphorylated tau (p‐tau), α‐synuclein, and transactive response DNA binding protein 43 (TDP‐43).

METHODS: We acquired neuropathological measures and brain segmentations from postmortem analysis of 132 donors with Alzheimer's disease (AD), Lewy body disease (LBD), frontotemporal lobar degeneration with TDP‐43 (FTLD‐TDP), and FTLD‐tau.

RESULTS: LBD had the least subcortical, limbic, and cortical atrophy compared to AD, FTLD‐TDP, and FTLD‐tau. In donors with both AD and LBD pathologies, primary LBD was associated with less atrophy than primary AD. While AD had cortico‐subcortical and cortico‐limbic morphometric associations, LBD had more limited parieto‐occipital cortico‐limbic associations. FTLD‐TDP had cortico‐subcortical while FTLD‐tau had cortico‐subcortical and cortico‐limbic associations. In AD and FTLD‐tau, hippocampal volumes correlated with p‐tau burden, neuron loss, and gliosis. In LBD, thalamic α‐synuclein severity was associated with subcortical/limbic volumes.

DISCUSSION: Postmortem neuroimaging reveals disease‐ and region‐specific structure–pathology relationships.

Reproduced under the paper's license (CC BY-NC), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 11 matches between paragraphs and lines of code.

Pulkit-Khandelwal/purple-mri

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 671f35b6d8eb41e3f75c4ff3b16cedc1986a9b24, 10 September 2026

Languages: Shell (39), Python (29), C++ (1), Jupyter (1)

Size: 178 files, 70 scripts

Software Heritage: not archived

Found in: the text, “Image analysis”

Holds: README, environment (pyproject.toml, docs/requirements.txt, pkg_src/pyproject.toml, docker/nighres_docker/Dockerfile, docker/segmentation_docker/Dockerfile, docker/segmentation_docker/requirements.txt, docs/source/_static/singularity.png), documentation, 1 notebook

Not found: license file, CITATION.cff, tests, continuous integration

Tools: NumPy (24 files), FreeSurfer (21 files), PyTorch (14 files), SciPy (13 files), SimpleITK (12 files), nnU-Net (11 files), OpenCV (9 files), pandas (7 files), NiBabel (6 files), ANTs (2 files), Matplotlib (2 files), Pillow (2 files), FSL (1 file), Nighres (1 file), Pingouin (1 file), Plotly (1 file), scikit-image (1 file), scikit-learn (1 file), seaborn (1 file), statannotations (1 file), statsmodels (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

71 files

docker/nighres_docker/clean_labels_final.sh, Shell, 19 lines
docker/nighres_docker/nighres_topology_corection.py, Python, 55 lines
docker/nighres_docker/prepare_cruise_files.sh, Shell, 37 lines
docker/segmentation_docker/commands_nnunet_inference.sh, Shell, 43 lines
docker/segmentation_docker/sitecustomize.py, Python, 45 lines
docs/source/conf.py, Python, 36 lines
glm/get_roi_statistics.sh, Shell, 21 lines
glm/glm.sh, Shell, 32 lines
glm/glm_dbm.sh, Shell, 76 lines
glm/prepare_vtk_files_for_glm.py, Python, 171 lines
glm/warp_to_template_space.sh, Shell, 78 lines
misc_scripts/c3d_commands.sh, Shell, 5 lines
misc_scripts/differences_between_lists.sh, Shell, 40 lines
misc_scripts/dots_to_template/archive_create_multi_comp_vtk_files.sh, Shell, 38 lines
misc_scripts/dots_to_template/dilation_split_dots.sh, Shell, 22 lines
misc_scripts/dots_to_template/dots_to_surface_template.sh, Shell, 87 lines
misc_scripts/dots_to_template/prepare_vtk_file_for_merge.py, Python, 158 lines
misc_scripts/dots_to_template/split_labels.py, Python, 42 lines
misc_scripts/flip_conform.sh, Shell, 23 lines
misc_scripts/flywheel_download.py, Python, 27 lines
misc_scripts/get_flash_reslice_and_flash_dots.sh, Shell, 22 lines
misc_scripts/get_patches.py, Python, 55 lines
misc_scripts/mris_place_surface.cpp, C++, 2,174 lines
misc_scripts/my_functions.py, Python, 70 lines
misc_scripts/nifti_to_dicom.py, Python, 99 lines
misc_scripts/perform_bias_correction.sh, Shell, 44 lines
misc_scripts/reorient_secondary.py, Python, 48 lines
nnU-Net-archs/generic_AnatomyCENet.py, Python, 702 lines
nnU-Net-archs/generic_AnatomyProjectNet.py, Python, 792 lines
nnU-Net-archs/generic_AnatomySECENet.py, Python, 702 lines
nnU-Net-archs/generic_AnatomySENet.py, Python, 702 lines
nnU-Net-archs/generic_AttentionNet.py, Python, 522 lines
nnU-Net-archs/generic_UNet.py, Python, 471 lines
nnU-Net-archs/generic_UNet_DP.py, Python, 124 lines
nnU-Net-archs/generic_UPKNet.py, Python, 640 lines
nnU-Net-archs/generic_VNet.py, Python, 615 lines
nnU-Net-archs/generic_VoxResNet.py, Python, 420 lines
notebooks/analysis_notebook_10_5_23.ipynb, Jupyter, 615 lines, 2 matches
pkg_src/purple_mri/__init__.py, Python, 1 line
pkg_src/purple_mri/__main__.py, Python, 13 lines
pkg_src/purple_mri/scripts/create_mesh.sh, Shell, 25 lines
pkg_src/purple_mri/scripts/decimate_mesh.py, Python, 296 lines
pkg_src/purple_mri/scripts/do_mni_txs.sh, Shell, 32 lines
pkg_src/purple_mri/scripts/make_fs_directories.sh, Shell, 43 lines
pkg_src/purple_mri/scripts/parcellation.sh, Shell, 236 lines
pkg_src/purple_mri/scripts/prepare_segm_files.sh, Shell, 37 lines
pkg_src/purple_mri/scripts/run_surface_pipeline.sh, Shell, 43 lines
pkg_src/purple_mri/scripts/topology_correction.sh, Shell, 39 lines
postsynth/generate.py, Python, 359 lines
postsynth/helper.py, Python, 822 lines
purple_mri/create_mesh.sh, Shell, 26 lines
purple_mri/decimate_mesh.py, Python, 310 lines
purple_mri/do_mni_txs.sh, Shell, 31 lines
purple_mri/infant_scripts/do_mni_txs_mcribs.sh, Shell, 31 lines
purple_mri/infant_scripts/parcellation_mcribs.sh, Shell, 121 lines
purple_mri/make_fs_directories.sh, Shell, 43 lines
purple_mri/parcellation.sh, Shell, 289 lines
purple_mri/prepare_segm_files.sh, Shell, 57 lines
purple_mri/run_surface_pipeline.sh, Shell, 45 lines
purple_mri/topology_correction.sh, Shell, 39 lines
scripts/exvivo_invivo_greedy_registration.sh, Shell, 135 lines
scripts/intensity_template/greedy_build_template.sh, Shell, 588 lines
scripts/intensity_template/pybatch.sh, Shell, 121 lines
scripts/intensity_template/warp_and_mri_init_template.sh, Shell, 59 lines
scripts/register_flash_t2w_ciss.sh, Shell, 46 lines
scripts/registration_exvivo_invivo_greedy_v2.sh, Shell, 217 lines
visualization/download_for_screenshots_fs.sh, Shell, 33 lines
visualization/ggseg_sample.py, Python, 26 lines
visualization/screenshots_fs_viz.py, Python, 18 lines
visualization/screenshots_fs_viz.sh, Shell, 19 lines
README.md, Text, 111 lines

At the source: github.com/Pulkit-Khandelwal/purple-mri

Pulkit-Khandelwal/postmortem-subcortical-limbic-pathologies

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 6abfd1d98fa06009578e0e292a72996ba3e46a58, 5 May 2026

Languages: Jupyter (2)

Size: 3 files, 2 scripts

Software Heritage: not archived

Found in: the text, “Statistical analysis”

Holds: README, 2 notebooks

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: Matplotlib (2 files), NumPy (2 files), pandas (2 files), Pingouin (2 files), SciPy (2 files), seaborn (2 files), statsmodels (2 files), NiBabel (1 file), Pillow (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

3 files

subcortical_analyses.ipynb, Jupyter, 6,168 lines, 6 matches
subcortical_cortical_analyses.ipynb, Jupyter, 1,182 lines, 3 matches
README.md, Text, 2 lines

At the source: github.com/Pulkit-Khandelwal/postmortem-subcortical-limbic-pathologies

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
72 scripts, each with its path and the digest of its content;
11 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

purple-mri.readthedocs.io, at purple-mri.readthedocs.io; found in the text, “Image analysis”

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 31 authors, 16 keywords, 18 MeSH terms, 2 funders, 71 references.

Cite

This paper

Khandelwal, P., Duong, M. T., Levorse, L. M., Trotman, W., Bahena, A., Lim, S. A., Denning, A. E., Chung, E., Olm, C. A., Radhakrishnan, H., Ittyerah, R., Prabhakaran, K., Mizsei, G., Schuck, T., Emrani, S., Vizcarra, J. A., Robinson, J., Ohm, D. T., Phillips, J. S., . . . Yushkevich, P. A. (2026). Postmortem brain MRI reveals differential associations of subcortical and limbic volumes with cortical thinning and neurodegenerative pathologies. Alzheimer's & dementia : the journal of the Alzheimer's Association, 22(7), e71649. https://doi.org/10.1002/alz.71649

BibTeX

@article{khandelwal2026postmortem, author = {Khandelwal, Pulkit and Duong, Michael Tran and Levorse, Lisa M and Trotman, Winifred and Bahena, Alejandra and Lim, Sydney A and Denning, Amanda E and Chung, Eunice and Olm, Christopher A and Radhakrishnan, Hamsanandini and Ittyerah, Ranjit and Prabhakaran, Karthik and Mizsei, Gabor and Schuck, Theresa and Emrani, Sheina and Vizcarra, Joaquin A and Robinson, John and Ohm, Daniel T and Phillips, Jeffrey S and Cohen, Jesse and Wisse, Laura E M and Detre, John A and Nasrallah, Ilya M and Brown, Christopher A and Das, Sandhitsu R and Lee, Edward B and Tisdall, M Dylan and Irwin, David J and McMillan, Corey T and Wolk, David A and Yushkevich, Paul A}, title = {{Postmortem brain MRI reveals differential associations of subcortical and limbic volumes with cortical thinning and neurodegenerative pathologies}}, journal = {Alzheimer's \& dementia : the journal of the Alzheimer's Association}, year = {2026}, month = jul, volume = {22}, number = {7}, pages = {e71649}, publisher = {Wiley}, issn = {1552-5260}, doi = {10.1002/alz.71649}, url = {https://doi.org/10.1002/alz.71649}, pmid = {42418450}, pmcid = {PMC13344898} }

RIS

TY - JOUR AU - Khandelwal, Pulkit AU - Duong, Michael Tran AU - Levorse, Lisa M AU - Trotman, Winifred AU - Bahena, Alejandra AU - Lim, Sydney A AU - Denning, Amanda E AU - Chung, Eunice AU - Olm, Christopher A AU - Radhakrishnan, Hamsanandini AU - Ittyerah, Ranjit AU - Prabhakaran, Karthik AU - Mizsei, Gabor AU - Schuck, Theresa AU - Emrani, Sheina AU - Vizcarra, Joaquin A AU - Robinson, John AU - Ohm, Daniel T AU - Phillips, Jeffrey S AU - Cohen, Jesse AU - Wisse, Laura E M AU - Detre, John A AU - Nasrallah, Ilya M AU - Brown, Christopher A AU - Das, Sandhitsu R AU - Lee, Edward B AU - Tisdall, M Dylan AU - Irwin, David J AU - McMillan, Corey T AU - Wolk, David A AU - Yushkevich, Paul A TI - Postmortem brain MRI reveals differential associations of subcortical and limbic volumes with cortical thinning and neurodegenerative pathologies T2 - Alzheimer's & dementia : the journal of the Alzheimer's Association J2 - Alzheimers Dement PY - 2026 DA - 2026/07/01 VL - 22 IS - 7 SP - e71649 SN - 1552-5260 PB - Wiley DO - 10.1002/alz.71649 UR - https://doi.org/10.1002/alz.71649 LA - en ER -

CSL-JSON

{ "id": "10.1002/alz.71649", "type": "article-journal", "title": "Postmortem brain MRI reveals differential associations of subcortical and limbic volumes with cortical thinning and neurodegenerative pathologies", "container-title": "Alzheimer's & dementia : the journal of the Alzheimer's Association", "author": [ { "family": "Khandelwal", "given": "Pulkit" }, { "family": "Duong", "given": "Michael Tran" }, { "family": "Levorse", "given": "Lisa M" }, { "family": "Trotman", "given": "Winifred" }, { "family": "Bahena", "given": "Alejandra" }, { "family": "Lim", "given": "Sydney A" }, { "family": "Denning", "given": "Amanda E" }, { "family": "Chung", "given": "Eunice" }, { "family": "Olm", "given": "Christopher A" }, { "family": "Radhakrishnan", "given": "Hamsanandini" }, { "family": "Ittyerah", "given": "Ranjit" }, { "family": "Prabhakaran", "given": "Karthik" }, { "family": "Mizsei", "given": "Gabor" }, { "family": "Schuck", "given": "Theresa" }, { "family": "Emrani", "given": "Sheina" }, { "family": "Vizcarra", "given": "Joaquin A" }, { "family": "Robinson", "given": "John" }, { "family": "Ohm", "given": "Daniel T" }, { "family": "Phillips", "given": "Jeffrey S" }, { "family": "Cohen", "given": "Jesse" }, { "family": "Wisse", "given": "Laura E M" }, { "family": "Detre", "given": "John A" }, { "family": "Nasrallah", "given": "Ilya M" }, { "family": "Brown", "given": "Christopher A" }, { "family": "Das", "given": "Sandhitsu R" }, { "family": "Lee", "given": "Edward B" }, { "family": "Tisdall", "given": "M Dylan" }, { "family": "Irwin", "given": "David J" }, { "family": "McMillan", "given": "Corey T" }, { "family": "Wolk", "given": "David A" }, { "family": "Yushkevich", "given": "Paul A" } ], "container-title-short": "Alzheimers Dement", "volume": "22", "issue": "7", "page": "e71649", "DOI": "10.1002/alz.71649", "PMID": "42418450", "PMCID": "PMC13344898", "ISSN": "1552-5260", "publisher": "Wiley", "URL": "https://doi.org/10.1002/alz.71649", "language": "en", "issued": { "date-parts": [ [ 2026, 7, 1 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.21203/rs.3.rs-9326213/v1 [code]: Multi-task fMRI outperforms resting-state fMRI for revealing task-invariant organization of the human brain
Journal: Research Square (preprint)
In common: statannotations, ANTs, FreeSurfer, 11 other tools, 1 reference
[2] doi:10.64898/2026.03.09.710558 [code]: Multi-task fMRI outperforms resting-state fMRI for revealing task-invariant organization of the human brain
Journal: bioRxiv (preprint)
In common: statannotations, ANTs, FreeSurfer, 11 other tools, 1 reference
[3] doi:10.1038/s41598-026-55397-w [code]: Fast surface reconstruction of human brain MRI: benchmarking deep-learning based morphometry tools.
Journal: Scientific reports
In common: Nighres, SimpleITK, ANTs, 8 other tools, structural MRI / diffusion, 2 references
[4] doi:10.1002/hipo.70124 [code]: Association Between Anterior Hippocampal Gyrification and Episodic Memory Performance in Neurotypical Young Adults.
Journal: Hippocampus
In common: Nighres, nnU-Net, SimpleITK, 9 other tools, structural MRI / diffusion, 1 reference
[5] doi:10.1162/imag.a.1164 [code]: Bias and generalizability of brain age prediction models: A multi-cohort evaluation with anatomical and interpretability insights.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: ANTs, FreeSurfer, FSL, 11 other tools, Alzheimer's / dementia, structural MRI / diffusion
[6] doi:10.3389/frai.2026.1771088 [code]: Few-shot deployment of pretrained MRI transformers in brain imaging tasks.
Journal: Frontiers in artificial intelligence
In common: nnU-Net, SimpleITK, OpenCV, 10 other tools, structural MRI / diffusion, 1 reference
[7] doi:10.2463/mrms.mp.2024-0149 [code]: Image Distortion Correction for Diffusion MR Imaging Using a Transformer-based U-Net.
Journal: Magnetic resonance in medical sciences : MRMS : an official journal of Japan Society of Magnetic Resonance in Medicine
In common: nnU-Net, ANTs, FreeSurfer, 10 other tools, structural MRI / diffusion
[8] doi:10.1371/journal.pcbi.1014555 [code]: Body surface potential driven personalisation of electrophysiological digital twins in hypertrophic cardiomyopathy.
Journal: PLoS computational biology
In common: nnU-Net, SimpleITK, ANTs, 10 other tools, structural MRI / diffusion
[9] doi:10.3389/fnins.2026.1870124 [code]: An end-to-end pipeline for automated fetal brain segmentation and biometry from 3D SSFP MRI.
Journal: Frontiers in neuroscience
In common: nnU-Net, SimpleITK, FSL, 10 other tools, structural MRI / diffusion
[10] doi:10.1016/j.crmeth.2026.101473 [code]: AmygdalaGo-BOLT for boundary-aware segmentation of the human amygdala.
Journal: Cell reports methods
In common: SimpleITK, ANTs, FreeSurfer, 9 other tools, structural MRI / diffusion

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.