OSCR

Quantifying generalization error in machine learning prediction of cognitive decline.

Code ↔ Paper

4 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 4 matches
  1. [1] § Methods › Structural MRI data ↔ 01_Analysis.py, lines 512–600 · score 0.99 · right cerebellar white, mid anterior, mid posterior, cerebral white matter, cortical thickness, right lateral
  2. [2] § Results › Most predictive features ↔ 01_Analysis.py, lines 512–600 · score 0.86 · word fluency, right thalamus, right amygdala, right hippocampus, corpus callosum, accumbens
  3. [3] § Methods › Predictive analysis ↔ 01_Analysis.py, lines 120–163 · score 0.61 · Scikit learn, random forest, regression, MMSE, SOB, predictive
  4. [4] § Methods › Model evaluation ↔ 01_Analysis.py, lines 760–792 · score 0.54 · squared error, absolute error, metrics, MSE, MAE, splits

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 authors' code

Python · 1,371 lines · 65 KB · MIT · 4 matches

  1. # %%
  2. # # 01_Analysis
  4. # %%
  5. # In summary, 01_Analysis.py contains code for predictive modeling and statistical analyses comparing cognitive decline prediction across OASIS-3 and ADNI datasets. The code is structured as follows:
  6. #
  7. # 0. Setup and Helper Functions
  8. # 0.1 Package imports
  9. # 0.2 Helper functions for data processing, model fitting, and evaluation
  10. #
  11. # 1. Data Loading and Within-Dataset Predictions
  12. # 1.1 Data loading and preprocessing
  13. # 1.2 Within-dataset predictions (OASIS-3 and ADNI)
  14. # 1.3 Results aggregation and summary
  15. #
  16. # 2. Permutation Importance Analysis
  17. # 2.1 Permutation importance for individual features
  18. # 2.2 Feature importance visualization and analysis
  19. # 2.3 Coalition-based feature importance
  20. #
  21. # 3. Across-Dataset Predictions
  22. # 3.1 Cross-dataset prediction (clinical, structural, combined, top-15)
  23. # 3.2 Cross-dataset results aggregation
  24. #
  25. # 4. Statistical Model Comparisons
  26. # 4.1 Statistical Testing Framework and Model Comparisons
  27. # 4.2 Absolute Error Comparisons Between Models
  28. #
  29. # 5. Subgroup Analysis
  30. # 5.1 Absolute error analysis preparation
  31. # 5.2 Subgroup comparisons by diagnosis, sex, APOE E4, and age
  32. # 5.3 Results formatting and export
  33. #
  34. # 6. Post-Hoc Matched Samples Analysis
  35. # 6.1 Quantile-based sample matching
  36. # 6.2 Matched sample predictions and evaluation
  37. # 6.3 Results for both combined and top-15 models
  38. #
  39. # Note: Running the predictive models is computationally expensive with n_splits=1000. Consider reducing splits for testing.
  40. # Full results tables are available in the 02_Supplementary_Material folder at https://osf.io/up65f/files/osfstorage.
  43. # %%
  44. # ## 0. Setup and Helper Functions
  45. # ### 0.1 Package Imports
  47. import joblib
  48. import pandas as pd
  49. import numpy as np
  50. import seaborn as sns
  51. import seaborn.objects as so
  53. import sklearn
  54. from sklearn import metrics
  55. from sklearn.pipeline import make_pipeline
  56. from sklearn.ensemble import RandomForestRegressor
  57. from sklearn.multioutput import MultiOutputRegressor
  58. from sklearn.experimental import enable_iterative_imputer
  59. from sklearn.impute import IterativeImputer
  60. from sklearn.metrics import r2_score
  61. from sklearn.model_selection import ShuffleSplit
  62. from sklearn.metrics import r2_score, mean_absolute_error, mean_squared_error
  63. from sklearn.inspection import permutation_importance
  65. import matplotlib
  66. import matplotlib.pyplot as plt
  67. matplotlib.rcParams['font.family'] = ['Arial']
  69. from scipy import stats
  70. from scipy.stats import wilcoxon
  71. from scipy.stats import mannwhitneyu
  73. from tqdm.notebook import tqdm
  75. from pathlib import Path
  76. from joblib import Parallel, delayed
  78. np.random.seed(21) # by setting a seed, you can ensure that every time you run the code, the sequence of random numbers generated will be the same
  80. # %%
  81. # ### 0.2 Helper Functions for Data Processing, Model Fitting, and Evaluation
  83. def get_data(clinical_features_path, data_slopes_path, structural_data_path, missing_cols=[]):
  84. clinical_features = pd.read_pickle(clinical_features_path)
  85. data_slopes = pd.read_pickle(data_slopes_path)
  86. structural_data = pd.read_pickle(structural_data_path)
  88. # clinical_features.columns.str.strip() str --> columns as strings not as indexes
  89. clinical_features.columns = clinical_features.columns.str.removeprefix("clin__npsy__").str.removeprefix("clin__risk__").str.removeprefix("clin__assess__")
  91. clinical_features["demo_sex"] = clinical_features["demo_sex"].map({'F': 0, 'M': 1, 0: 0, 1: 1})
  92. clinical_features["diag"] = clinical_features["diag"].map({"hc": 0, "mci": 1, "dem": 2.0, "NL": 0.0, "MCI": 1.0, "Dementia": 2.0})
  94. # defining the dependent and independent variables and removing the variable "subject" from the dataframes
  95. y = data_slopes.copy().set_index("subject")
  96. X_clin = clinical_features.copy().set_index("subject")
  97. X_fs = structural_data.copy().set_index("subject")
  98. X_clin_fs = pd.concat([X_clin, X_fs], axis=1)
  100. # remove missing columns if specified
  101. if missing_cols:
  102. X_clin = X_clin.drop(columns=missing_cols, errors='ignore')
  103. X_fs = X_fs.drop(columns=missing_cols, errors='ignore')
  104. X_clin_fs = X_clin_fs.drop(columns=missing_cols, errors='ignore')
  106. return y, X_clin, X_fs, X_clin_fs
  109. def get_iterative_imputer(X_train, X_test, random_state):
  110. imputer = IterativeImputer(add_indicator=True, random_state=random_state, n_nearest_features=50)
  111. missing_cols_train = [c + "_missing" for c in X_train.columns[X_train.isna().any()]]
  112. missing_cols_test = [c + "_missing" for c in X_test.columns[X_test.isna().any()]]
  113. feature_names_train = X_train.columns.to_list() + missing_cols_train
  114. feature_names_test = X_test.columns.to_list() + missing_cols_test
  115. return imputer, feature_names_train, feature_names_test
  117. # %%
  120. def fit_pipeline(X, y, train_idx, test_idx, root_dir):
  121. X_train, X_test = X.iloc[train_idx], X.iloc[test_idx]
  122. y_train, y_test = y.iloc[train_idx], y.iloc[test_idx]
  124. y_train = pd.DataFrame(y_train)
  125. y_test = pd.DataFrame(y_test)
  127. # iterative imputation
  129. # define imputer
  130. imputer, feature_names_train, feature_names_test = get_iterative_imputer(X_train, X_test, random_state=0)
  132. # define pipeline
  133. pipeline = make_pipeline(imputer, RandomForestRegressor(random_state=0))
  135. # fit pipeline
  136. pipeline.fit(X_train, y_train)
  138. # predict
  139. y_pred = pd.DataFrame(pipeline.predict(X_test)).rename(columns={0: 'mmse_pred', 1: 'sob_pred'})
  140. y_pred.index = y_test.index.values
  142. # concatenate
  143. df = pd.concat([y_pred, y_test], axis = 1)
  145. # get MAE, MSE, R2 for the predictions with scikit-learn
  147. metrics = {
  148. 'r2': r2_score,
  149. 'mae': mean_absolute_error,
  150. 'mse': mean_squared_error,
  151. }
  152. results_df = pd.DataFrame({
  153. var: {name: func(df[f'{var}_slope'], df[f'{var}_pred']) for name, func in metrics.items()}
  154. for var in ['sob', 'mmse']
  155. }).T
  159. # save pipeline and dataframe
  160. if root_dir:
  161. save_pipeline_and_df(pipeline, df, results_df, root_dir)
  163. return df, pipeline, results_df
  165. def save_pipeline_and_df(pipeline, df, results_df, root_dir):
  166. root_dir = Path(root_dir)
  167. root_dir.mkdir(parents=True, exist_ok=True)
  168. joblib.dump(pipeline, f"{root_dir}/pipeline.joblib", compress=('xz', 3))
  169. df.to_csv(f"{root_dir}/predictions.csv", index=True, header=True)
  170. results_df.to_csv(f"{root_dir}/results.csv", index=True, header=True)
  172. def compress_pipeline_file(root_dir):
  173. pipeline_path = Path(root_dir) / "pipeline.joblib"
  174. if pipeline_path.exists():
  175. with open(pipeline_path, 'rb') as f:
  176. pipeline = joblib.load(f)
  177. with open(pipeline_path, 'wb') as f:
  178. joblib.dump(pipeline, f, compress=('xz', 3))
  181. def split_and_fit(X, y, root_dir, n_splits=1000, test_size=0.2):
  182. ss = ShuffleSplit(n_splits=n_splits, test_size=test_size, random_state=0)
  184. # split data
  185. for iter, (train_idx, test_idx) in tqdm(enumerate(ss.split(X, y)), total=n_splits):
  186. path = Path(f"{root_dir}/iteration_{iter+1}")
  187. path.mkdir(parents=True, exist_ok=True)
  188. # check if both pipeline and df already exist
  189. if (path / "pipeline.joblib").exists() and (path / "predictions.csv").exists():
  190. continue
  192. fit_pipeline(X, y, train_idx, test_idx, path)
  194. def get_results(root, iters=range(1000)):
  195. dfs = [pd.read_csv(f"{root}/iteration_{i+1}/results.csv", index_col=0) for i in iters]
  196. return pd.concat(dfs, keys=iters, names=["iteration", "target"])
  198. def get_results_summary(root, iters=range(1000)):
  199. results = get_results(root, iters)
  200. summary = results.groupby(level=0).agg(['mean', 'median'])
  201. summary.columns = ['_'.join(col).strip() for col in summary.columns.values]
  202. return summary
  204. # %%
  205. # ## 1. Data Loading and Within-Dataset Predictions
  206. # ### 1.1 Data Loading and Preprocessing
  208. # Load data
  209. oasis_y, oasis_X_clin, oasis_X_fs, oasis_X_clin_fs = get_data("OASIS_clinical_features_X.pkl", "OASIS_data_slopes_y.pkl", "OASIS_structGlobScort.pkl")
  210. adni_y, adni_X_clin, adni_X_fs, adni_X_clin_fs = get_data("ADNI_1237_data_clinical_features_X_clinsessionn_2025.pkl", "ADNI_1237_data_slopes_y_2025.pkl", "ADNI_1237_data_structGlobScort_2025.pkl")
  212. missing_adni = ["smoke_PACKSPER", "TMTA_TRAILALI", "smoke_TOBAC100", "smoke_SMOKYRS", "TMTB_TRAILBLI", "smoke_TOBAC30", "WMSds_DIGIF", "WMSds_DIGIFLEN", "WF_VEG", "WMSds_DIGIB", "WMSds_DIGIBLEN", "WAIS_WAIS", "familyHist_DADDEM", "familyHist_MOMDEM"]
  213. missing_oasis = ["familyHist_sumSibDem","familyHist_ratioSibDem","TMTB_TRAILBRR","TMTB_TRAILBLI","TMTA_TRAILALI","TMTA_TRAILARR"]
  215. missing_cols = missing_adni + missing_oasis
  217. oasis_y, oasis_X_clin_filtered, oasis_X_fs_filtered, oasis_X_clin_fs_filtered = get_data("OASIS_clinical_features_X.pkl", "OASIS_data_slopes_y.pkl", "OASIS_structGlobScort.pkl", missing_cols=missing_cols)
  218. adni_y, adni_X_clin_filtered, adni_X_fs_filtered, adni_X_clin_fs_filtered = get_data("ADNI_1237_data_clinical_features_X_clinsessionn_2025.pkl", "ADNI_1237_data_slopes_y_2025.pkl", "ADNI_1237_data_structGlobScort_2025.pkl", missing_cols=missing_cols)
  220. data_combinations = [
  221. (oasis_y, oasis_X_clin, "OASIS_clin_predictions"),
  222. (oasis_y, oasis_X_fs, "OASIS_struct_predictions"),
  223. (oasis_y, oasis_X_clin_fs, "OASIS_clin_struct_predictions"),
  224. (adni_y, adni_X_clin, "ADNI_clin_predictions"),
  225. (adni_y, adni_X_fs, "ADNI_struct_predictions"),
  226. (adni_y, adni_X_clin_fs, "ADNI_clin_struct_predictions"),
  227. (oasis_y, oasis_X_clin_filtered, "OASIS_clin_predictions_filtered"),
  228. (oasis_y, oasis_X_fs_filtered, "OASIS_struct_predictions_filtered"),
  229. (oasis_y, oasis_X_clin_fs_filtered, "OASIS_clin_struct_predictions_filtered"),
  230. (adni_y, adni_X_clin_filtered, "ADNI_clin_predictions_filtered"),
  231. (adni_y, adni_X_fs_filtered, "ADNI_struct_predictions_filtered"),
  232. (adni_y, adni_X_clin_fs_filtered, "ADNI_clin_struct_predictions_filtered"),
  233. ]
  235. oasis_X_clin_filtered.to_csv("OASIS_clin_filtered.csv", index=True, header=True)
  236. adni_X_clin_filtered.to_csv("ADNI_clin_filtered.csv", index=True, header=True)
  238. # ### 1.2 Within-Dataset Predictions (OASIS-3 and ADNI)
  240. for y, X, root_dir in data_combinations:
  241. print(f"Processing {root_dir}...")
  242. split_and_fit(X, y, root_dir, n_splits=1000)
  244. # %%
  245. # ### 1.3 Results Aggregation and Summary
  247. def get_dataset_results(clin_file, struct_file, clin_struct_file, iters=range(1000)):
  248. clin_results = get_results(clin_file, iters)
  249. struct_results = get_results(struct_file, iters)
  250. clin_struct_results = get_results(clin_struct_file, iters)
  252. results = pd.concat([clin_results, struct_results, clin_struct_results], keys=["clin", "struct", "clin_struct"], names=["modality", "iteration", "target"])
  254. return results
  256. oasis_results = get_dataset_results("OASIS_clin_predictions", "OASIS_struct_predictions", "OASIS_clin_struct_predictions")
  257. adni_results = get_dataset_results("ADNI_clin_predictions", "ADNI_struct_predictions", "ADNI_clin_struct_predictions")
  258. oasis_filtered_results = get_dataset_results("OASIS_clin_predictions_filtered", "OASIS_struct_predictions_filtered", "OASIS_clin_struct_predictions_filtered")
  259. adni_filtered_results = get_dataset_results("ADNI_clin_predictions_filtered", "ADNI_struct_predictions_filtered", "ADNI_clin_struct_predictions_filtered")
  261. results = pd.concat([oasis_results, adni_results, oasis_filtered_results, adni_filtered_results], keys=["OASIS", "ADNI", "OASIS_filtered", "ADNI_filtered"], names=["dataset", "modality", "iteration", "target"])
  262. results.to_csv("results_within.csv")
  264. results.groupby(level=["dataset", "modality", "target"]).agg(['mean', 'median'])
  267. # %%
  268. # ## 2. Permutation Importance Analysis
  269. # ### 2.1 Permutation Importance for Individual Features
  271. def load_pipeline(iter, root_dir):
  272. path = Path(f"{root_dir}/iteration_{iter+1}/pipeline.joblib")
  273. if path.exists():
  274. return joblib.load(path)
  275. else:
  276. raise FileNotFoundError(f"Pipeline for iteration {iter+1} not found in {root_dir}.")
  279. def get_permutation_importance_iter(pipeline, X_test, y_test, n_repeats=100, random_state=0):
  280. pi = permutation_importance(pipeline, X_test, y_test, n_repeats=n_repeats, scoring='r2', random_state=random_state)
  281. return pd.DataFrame({"feature": X_test.columns, "permutation_importance": pi.importances_mean})
  283. def get_permutation_importance(X, y, root_dir, n_splits=1000, n_repeats=100):
  285. # quickly check that every path exists
  286. root_dir = Path(root_dir)
  287. for iter in range(n_splits):
  288. path = root_dir / f"iteration_{iter+1}"
  289. if not path.exists():
  290. raise FileNotFoundError(f"Path {path} does not exist. Make sure to run split_and_fit first.")
  292. def process_iteration(iter):
  293. path = root_dir / f"iteration_{iter+1}"
  294. pi_path = path / f"permutation_importance_repeats-{n_repeats}.csv"
  295. if pi_path.exists():
  296. print(f"Skipping iteration {iter+1}, already processed.")
  297. return pd.read_csv(pi_path, index_col=0)
  298. else:
  299. ss = ShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=0)
  300. train_indices = []
  301. test_indices = []
  302. for i, (train_idx, test_idx) in enumerate(ss.split(X, y)):
  303. if i == iter:
  304. train_indices = train_idx
  305. test_indices = test_idx
  306. break
  308. pipeline = load_pipeline(iter, root_dir)
  309. X_test = X.iloc[test_indices]
  310. y_test = y.iloc[test_indices]
  312. pi_df = get_permutation_importance_iter(pipeline, X_test, y_test, n_repeats=n_repeats)
  313. pi_df.to_csv(pi_path)
  314. return pi_df
  316. # Run the processing in parallel with 4 jobs
  317. results = Parallel(n_jobs=4)(
  318. delayed(process_iteration)(iter) for iter in tqdm(range(n_splits), desc="Processing iterations")
  319. )
  321. results = pd.concat(results, keys=range(n_splits), names=["iteration", "feature"])
  322. results = results.droplevel("feature")
  324. return results
  326. oasis_pi = get_permutation_importance(oasis_X_clin_fs, oasis_y,"OASIS_clin_struct_predictions", n_splits=1000, n_repeats=5)
  327. adni_pi = get_permutation_importance(adni_X_clin_fs, adni_y, "ADNI_clin_struct_predictions", n_splits=1000, n_repeats=5)
  328. oasis_filtered_pi = get_permutation_importance(oasis_X_clin_fs_filtered, oasis_y, "OASIS_clin_struct_predictions_filtered", n_splits=1000, n_repeats=5)
  329. adni_filtered_pi = get_permutation_importance(adni_X_clin_fs_filtered, adni_y, "ADNI_clin_struct_predictions_filtered", n_splits=1000, n_repeats=5)
  331. oasis_filtered_pi.set_index("feature", append=True).unstack("feature").droplevel(0, axis=1).to_csv("OASIS_clin_struct_permutation_importance.csv")
  332. adni_filtered_pi.set_index("feature", append=True).unstack("feature").droplevel(0, axis=1).to_csv("ADNI_clin_struct_permutation_importance.csv")
  334. # ### 2.2 Feature Importance Visualization and Analysis
  336. # %%
  337. def get_across_permutation_importance(X, y, root_dir, n_repeats=100, label="permutation_importance"):
  338. # quickly check that every path exists
  339. path = Path(root_dir)
  341. pi_path = path / f"{label}_repeats-{n_repeats}.csv"
  342. if pi_path.exists():
  343. print(f"Skipping {label}, already processed.")
  344. pi_df = pd.read_csv(pi_path, index_col=0)
  345. else:
  346. pipeline = joblib.load(f"{root_dir}/pipeline.joblib")
  347. pi = permutation_importance(pipeline, X, y, n_repeats=n_repeats, scoring='r2', random_state=0)
  348. pi_df = pd.DataFrame({"feature": X.columns, "permutation_importance": pi.importances_mean, "permutation_importance_std": pi.importances_std})
  349. pi_df.to_csv(pi_path)
  350. return pi_df
  352. oasis2adni_pi = get_across_permutation_importance(adni_X_clin_fs_filtered, adni_y, "OASIS_to_ADNI_clin_struct", n_repeats=100)
  353. adni2oasis_pi = get_across_permutation_importance(oasis_X_clin_fs_filtered, oasis_y, "ADNI_to_OASIS_clin_struct", n_repeats=100)
  356. # ### 2.3 Coalition-Based Feature Importance
  358. # %%
  359. feature_coalitions = pd.read_csv("features.csv", header=0)
  361. # group by coalition and feature type
  362. coalition_dict = feature_coalitions.groupby("coalition")["feature"].apply(list).to_dict()
  363. feature_type_dict = feature_coalitions.groupby("feature-type")["feature"].apply(list).to_dict()
  366. # get importance by feature coalition
  367. def importance_by_coalition(pipeline, X_test, y_test, coalition_dict, random_state=0):
  368. # manual reimplementation of permutation_importance
  369. result = {}
  371. for coalition, features_in in coalition_dict.items():
  372. # permute the rows within the coalition
  373. X_test_coalition = X_test.copy()
  374. X_test_coalition[features_in] = X_test_coalition[features_in].sample(frac=1, random_state=random_state).values
  375. pred = pipeline.predict(X_test_coalition)
  376. result[coalition] = metrics.r2_score(y_test, pred)
  378. pred = pipeline.predict(X_test)
  380. all_score = metrics.r2_score(y_test, pred)
  382. result = pd.DataFrame(result, index=['permutation_importance'])
  383. result = all_score - result.T
  385. return result
  388. def get_importance_by_coalition(X, y, root_dir, coalition_dict, n_splits=1000, n_repeats=5, label="coalition_importance"):
  389. # quickly check that every path exists
  390. root_dir = Path(root_dir)
  391. for i in range(n_splits):
  392. path = root_dir / f"iteration_{i+1}"
  393. if not path.exists():
  394. raise FileNotFoundError(f"Path {path} does not exist. Make sure to run split_and_fit first.")
  396. def process_iteration(iter):
  397. path = root_dir / f"iteration_{iter+1}"
  398. pi_path = path / f"{label}_repeats-{n_repeats}.csv"
  399. if pi_path.exists():
  400. return pd.read_csv(pi_path, index_col=0)
  401. else:
  402. ss = ShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=0)
  403. test_indices = []
  404. for i, (train_idx, test_idx) in enumerate(ss.split(X, y)):
  405. if i == iter:
  406. test_indices = test_idx
  407. break
  409. pipeline = load_pipeline(iter, root_dir)
  410. X_test = X.iloc[test_indices]
  411. y_test = y.iloc[test_indices]
  413. pi_df = [importance_by_coalition(pipeline, X_test, y_test, coalition_dict, random_state=i) for i in range(n_repeats)]
  414. pi_df = pd.concat(pi_df).groupby(level=0).mean()
  415. pi_df.to_csv(pi_path)
  416. return pi_df
  418. # Run the processing in parallel with 4 jobs
  419. results = Parallel(n_jobs=4)(
  420. delayed(process_iteration)(iter) for iter in tqdm(range(n_splits), desc="Processing iterations")
  421. )
  422. results = pd.concat(results, keys=range(n_splits), names=["iteration", "coalition"])
  424. return results
  426. oasis_filtered_ci = get_importance_by_coalition(oasis_X_clin_fs_filtered, oasis_y, "OASIS_clin_struct_predictions_filtered", coalition_dict, n_splits=1000, n_repeats=5)
  427. adni_filtered_ci = get_importance_by_coalition(adni_X_clin_fs_filtered, adni_y, "ADNI_clin_struct_predictions_filtered", coalition_dict, n_splits=1000, n_repeats=5)
  428. oasis_filtered_type_ci = get_importance_by_coalition(oasis_X_clin_fs_filtered, oasis_y, "OASIS_clin_struct_predictions_filtered", feature_type_dict, n_splits=1000, n_repeats=5, label="feature_type_importance")
  429. adni_filtered_type_ci = get_importance_by_coalition(adni_X_clin_fs_filtered, adni_y, "ADNI_clin_struct_predictions_filtered", feature_type_dict, n_splits=1000, n_repeats=5, label="feature_type_importance")
  431. # %%
  432. def get_mean_ci(df, alpha=0.05, df_dof=4):
  433. mean = df.mean()
  434. ci = stats.t.interval(alpha, df_dof, loc=mean, scale=stats.sem(df))
  435. ci_lower, ci_upper = ci
  436. return pd.Series({"mean": mean, "ci_lower": ci_lower, "ci_upper": ci_upper})
  438. def get_mean_ci_nb(df, alpha=0.05, VIF=(1 + 0.2/0.8)):
  439. R = len(df)
  440. mean = df.mean()
  441. var = df.var(ddof=1) # unbiased variance
  442. inflation = 1 + VIF
  443. se = np.sqrt(var * inflation / R)
  444. dof = R - 1
  445. t_val = stats.t.ppf(1 - alpha/2, dof)
  446. ci_lower = mean - t_val * se
  447. ci_upper = mean + t_val * se
  449. return pd.Series({"mean": mean, "ci_lower": ci_lower, "ci_upper": ci_upper})
  451. def get_ci_nb(df, alpha=0.05, VIF=(1 + 0.2/0.8)):
  452. res = get_mean_ci_nb(df, alpha, VIF)
  453. ci_lower = res['ci_lower']
  454. ci_upper = res['ci_upper']
  455. return (ci_lower, ci_upper)
  457. # %%
  459. # plot these as point plots with error bars
  460. def plot_permutation_importance(pi_df, title):
  461. plt.figure(figsize=(10, 6))
  462. sns.pointplot(
  463. data=pi_df.reset_index(),
  464. x='feature',
  465. y='permutation_importance',
  466. capsize=0.1,
  467. errorbar=get_ci_nb,
  468. markers='o',
  469. color='blue',
  470. markersize=2,
  471. linestyle='none',
  472. linewidth=1
  473. )
  474. plt.xticks(rotation=90)
  475. plt.title(title)
  476. plt.xlabel('Feature')
  477. plt.ylabel('Permutation Importance (Mean ± CI)')
  478. plt.tight_layout()
  479. plt.show()
  481. top_15_features_oasis = oasis_pi.groupby("feature").mean()["permutation_importance"].sort_values(ascending=False).head(15).index.tolist()
  482. top_15_features_adni = adni_pi.groupby("feature").mean()["permutation_importance"].sort_values(ascending=False).head(15).index.tolist()
  483. top_15_features = list(set(top_15_features_oasis + top_15_features_adni))
  485. plot_permutation_importance(oasis_pi.loc[oasis_pi["feature"].isin(top_15_features)], "OASIS-3: Top 15 Features Permutation Importance")
  486. plot_permutation_importance(adni_pi.loc[adni_pi["feature"].isin(top_15_features)], "ADNI: Top 15 Features Permutation Importance")
  488. # for across datasets:
  490. # %%
  492. def plot_across_permutation_importance(pi_df, title, ax=None):
  493. pi_df = pi_df.sort_values("permutation_importance", ascending=False).head(15).copy()
  494. pi_df['feature'] = pd.Categorical(pi_df['feature'], categories=pi_df['feature'], ordered=True)
  495. pi_df["ymin"] = pi_df["permutation_importance"] - pi_df["permutation_importance_std"]
  496. pi_df["ymax"] = pi_df["permutation_importance"] + pi_df["permutation_importance_std"]
  498. p = so.Plot(data=pi_df, y='feature', x='permutation_importance', xmin='ymin', xmax='ymax')
  499. p = p.add(so.Dot())
  500. p = p.add(so.Range())
  501. p = p.label(title=title, y="Feature", x="Permutation Importance ($R^2$ decrease)")
  503. # Set font to Arial for all text elements
  504. p = p.theme({
  505. "axes.labelweight": "normal",
  506. "axes.labelsize": 10,
  507. "axes.titlesize": 12,
  508. })
  510. p.on(ax).plot()
  512. # Create a figure with two subplots arranged vertically
  513. fig, axes = plt.subplots(2, 1, figsize=(8, 12), constrained_layout=True, sharex=True)
  516. feature_dict = {
  517. 'session_n': "Number of previous visits",
  518. 'diag': "Diagnosis",
  519. 'demo_sex': "Sex",
  520. 'demo_age': "Age",
  521. 'demo_education': "Education",
  522. 'diabetes_DIABETES': "Diabetes",
  523. 'hypercho_HYPERCHO': "Hypercholesterolemia",
  524. 'cvasc_HYPERTEN': "Hypertension",
  525. 'cvasc_CBSTROKE': "Stroke",
  526. 'cvasc_CBTIA': "Transient Ischemic Attack",
  527. 'cvasc_CVHATT': "Heart Attack",
  528. 'cvasc_CVAFIB': "Atrial Fibrillation",
  529. 'cvasc_CVANGIO': "Angiography",
  530. 'cvasc_CVOTHR': "Other Cardiovascular",
  531. 'cdr_commun': "CDR (Communication)",
  532. 'cdr_homehobb': "CDR (Home and Hobbies)",
  533. 'cdr_judgment': "CDR (Judgment)",
  534. 'cdr_memory': "CDR (Memory)",
  535. 'cdr_orient': "CDR (Orientation)",
  536. 'cdr_perscare': "CDR (Personal Care)",
  537. 'cdr_sob': "CDR-SOB",
  538. 'cdr_cdrGlobal': "CDR (Global)",
  539. 'mmse_mmse': "MMSE (Total)",
  540. 'gds_gdsSum': "GDS (Sum)",
  541. 'faq_BILLS': "FAQ (Bills)",
  542. 'faq_TAXES': "FAQ (Taxes)",
  543. 'faq_SHOPPING': "FAQ (Shopping)",
  544. 'faq_GAMES': "FAQ (Games)",
  545. 'faq_STOVE': "FAQ (Stove)",
  546. 'faq_MEALPREP': "FAQ (Meal Preparation)",
  547. 'faq_EVENTS': "FAQ (Events)",
  548. 'faq_PAYATTN': "FAQ (Pay Attention)",
  549. 'faq_REMDATES': "FAQ (Remember Dates)",
  550. 'faq_TRAVEL': "FAQ (Travel)",
  551. 'faq_faqSum': "FAQ (Total)",
  552. 'npiq_npiqPresSum': "NPI-Q (Presence Sum)",
  553. 'npiq_npiqSevSum': "NPI-Q (Severity Sum)",
  554. 'apoe_e2count': "APOE ε2 Count",
  555. 'apoe_e3count': "APOE ε3 Count",
  556. 'apoe_e4count': "APOE ε4 Count",
  557. 'WMSlm_LOGIMEM': "WMS (Logical Memory)",
  558. 'WMSlm_MEMUNITS': "WMS (Memory Units)",
  559. 'WMSlm_MEMTIME': "WMS (Memory Time)",
  560. 'WF_ANIMALS': "Word Fluency (Animals)",
  561. 'TMTA_TRAILA': "Trail Making Test A",
  562. 'TMTB_TRAILB': "Trail Making Test B",
  563. 'TMTB_TRAILBnorm': "Trail Making Test B (normalized)",
  564. 'BOSTON_BOSTON': "Boston Naming Test (Total)",
  565. 'fs__globalVolume__3rd-Ventricle': "Volume of 3rd Ventricle",
  566. 'fs__globalVolume__4th-Ventricle': "Volume of 4th Ventricle",
  567. 'fs__globalVolume__CC_Anterior': "Volume of Corpus Callosum (Anterior)",
  568. 'fs__globalVolume__CC_Central': "Volume of Corpus Callosum (Central)",
  569. 'fs__globalVolume__CC_Mid_Anterior': "Volume of Corpus Callosum (Mid Anterior)",
  570. 'fs__globalVolume__CC_Mid_Posterior': "Volume of Corpus Callosum (Mid Posterior)",
  571. 'fs__globalVolume__CC_Posterior': "Volume of Corpus Callosum (Posterior)",
  572. 'fs__globalVolume__Left-Cerebellum-Cortex': "Volume of Cerebellum Cortex (Left)",
  573. 'fs__globalVolume__Left-Cerebellum-White-Matter': "Volume of Cerebellum White Matter (Left)",
  574. 'fs__globalVolume__Left-Lateral-Ventricle': "Volume of Lateral Ventricle (Left)",
  575. 'fs__globalVolume__Right-Cerebellum-Cortex': "Volume of Cerebellum Cortex (Right)",
  576. 'fs__globalVolume__Right-Cerebellum-White-Matter': "Volume of Cerebellum White Matter (Right)",
  577. 'fs__globalVolume__Right-Lateral-Ventricle': "Volume of Lateral Ventricle (Right)",
  578. 'fs__globalVolume__SubCortGrayVol': "Volume of Subcortical Gray Matter",
  579. 'fs__globalVolume__TotalGrayVol': "Volume of Total Gray Matter",
  580. 'fs__globalVolume__lhCerebralWhiteMatterVol': "Volume of Cerebral White Matter (Left)",
  581. 'fs__globalVolume__lhCortexVol': "Volume of Cortex (Left)",
  582. 'fs__globalVolume__lh_MeanThickness_thickness': "Mean Cortical Thickness (Left)",
  583. 'fs__globalVolume__rhCerebralWhiteMatterVol': "Volume of Cerebral White Matter (Right)",
  584. 'fs__globalVolume__rhCortexVol': "Volume of Cortex (Right)",
  585. 'fs__globalVolume__rh_MeanThickness_thickness': "Mean Cortical Thickness (Right)",
  586. 'fs__subcortVolume__Left-Accumbens-area': "Volume of Accumbens (Left)",
  587. 'fs__subcortVolume__Left-Amygdala': "Volume of Amygdala (Left)",
  588. 'fs__subcortVolume__Left-Caudate': "Volume of Caudate (Left)",
  589. 'fs__subcortVolume__Left-Hippocampus': "Volume of Hippocampus (Left)",
  590. 'fs__subcortVolume__Left-Pallidum': "Volume of Pallidum (Left)",
  591. 'fs__subcortVolume__Left-Putamen': "Volume of Putamen (Left)",
  592. 'fs__subcortVolume__Left-Thalamus-Proper': "Volume of Thalamus (Left)",
  593. 'fs__subcortVolume__Right-Accumbens-area': "Volume of Accumbens (Right)",
  594. 'fs__subcortVolume__Right-Amygdala': "Volume of Amygdala (Right)",
  595. 'fs__subcortVolume__Right-Caudate': "Volume of Caudate (Right)",
  596. 'fs__subcortVolume__Right-Hippocampus': "Volume of Hippocampus (Right)",
  597. 'fs__subcortVolume__Right-Pallidum': "Volume of Pallidum (Right)",
  598. 'fs__subcortVolume__Right-Putamen': "Volume of Putamen (Right)",
  599. 'fs__subcortVolume__Right-Thalamus-Proper': "Volume of Thalamus (Right)",
  600. }
  602. oasis2adni_pi['feature_names'] = oasis2adni_pi['feature'].map(feature_dict)
  603. adni2oasis_pi['feature_names'] = adni2oasis_pi['feature'].map(feature_dict)
  605. # Plot for OASIS to ADNI
  606. plot_across_permutation_importance(oasis2adni_pi.rename(columns={"feature_names": "feature", "feature": "feature_i"}), "A) OASIS-3 → ADNI: Feature Importance", ax=axes[0])
  608. # Plot for ADNI to OASIS
  609. plot_across_permutation_importance(adni2oasis_pi.rename(columns={"feature_names": "feature", "feature": "feature_i"}), "B) ADNI → OASIS-3: Feature Importance", ax=axes[1])
  612. plt.savefig("cross_dataset_feature_importance.png", dpi=300, bbox_inches="tight")
  614. top15_features_oasis2adni = oasis2adni_pi.sort_values("permutation_importance", ascending=False).head(15)["feature"].to_list()
  615. top15_features_adni2oasis = adni2oasis_pi.sort_values("permutation_importance", ascending=False).head(15)["feature"].to_list()
  617. # %%
  619. top_15_features_oasis_filtered = oasis_filtered_pi.groupby("feature").mean()["permutation_importance"].sort_values(ascending=False).head(15).index.tolist()
  620. top_15_features_adni_filtered = adni_filtered_pi.groupby("feature").mean()["permutation_importance"].sort_values(ascending=False).head(15).index.tolist()
  621. top_15_features_filtered = list(set(top_15_features_oasis_filtered + top_15_features_adni_filtered))
  622. top_15_features_intersection = list(set(top_15_features_oasis_filtered) & set(top_15_features_adni_filtered))
  624. adni_overlap = set(top15_features_adni2oasis) & set(top_15_features_adni_filtered)
  625. adni_difference1 = set(top_15_features_adni_filtered) - set(top15_features_adni2oasis)
  626. adni_difference2 = set(top15_features_adni2oasis) - set(top_15_features_adni_filtered)
  628. oasis_overlap = set(top15_features_oasis2adni) & set(top_15_features_oasis_filtered)
  629. oasis_difference1 = set(top_15_features_oasis_filtered) - set(top15_features_oasis2adni)
  630. oasis_difference2 = set(top15_features_oasis2adni) - set(top_15_features_oasis_filtered)
  632. plot_permutation_importance(oasis_filtered_pi.loc[oasis_filtered_pi["feature"].isin(top_15_features)], "OASIS-3 (Filtered): Top 15 Features Permutation Importance")
  633. plot_permutation_importance(adni_filtered_pi.loc[adni_filtered_pi["feature"].isin(top_15_features)], "ADNI (Filtered): Top 15 Features Permutation Importance")
  635. plot_permutation_importance(oasis_filtered_ci.reset_index("coalition").rename(columns={"coalition": "feature"}), "OASIS-3 (Filtered): Permutation Importance by Coalition")
  636. plot_permutation_importance(adni_filtered_ci.reset_index("coalition").rename(columns={"coalition": "feature"}), "ADNI (Filtered): Permutation Importance by Coalition")
  638. # %%
  639. def round(x, n_figs):
  640. x = float(x)
  641. power = 10 ** np.floor(np.log10(abs(x)))
  642. rounded = np.round(x / power, n_figs - 1) * power
  643. rounded = float(rounded)
  645. # Use 'f' format to preserve trailing zeros if needed
  646. digits_after_decimal = max(n_figs - int(np.floor(np.log10(abs(rounded)))) - 1, 0)
  647. format_str = f"{{:.{digits_after_decimal}f}}"
  648. result = format_str.format(rounded)
  649. return result
  651. # oasis_filtered_pi_top15 = oasis_filtered_pi.loc[oasis_filtered_pi["feature"].isin(top_15_features_oasis_filtered)].copy()
  652. oasis_filtered_pi_top15 = oasis_filtered_pi.copy()
  653. oasis_filtered_pi_top15 = oasis_filtered_pi_top15.groupby("feature")["permutation_importance"].apply(get_mean_ci_nb).unstack()
  654. oasis_filtered_pi_top15 = oasis_filtered_pi_top15.sort_values("mean", ascending=False)
  656. oasis_filtered_pi_top15.loc[top_15_features_oasis_filtered].applymap(round, n_figs=2).apply(lambda x: f"{x['mean']} ({x['ci_lower']}, {x['ci_upper']})", axis=1)
  657. oasis_filtered_pi_top15.applymap(round, n_figs=2).apply(lambda x: f"{x['mean']} ({x['ci_lower']}, {x['ci_upper']})", axis=1).to_csv("OASIS_importance.csv")
  660. adni_filtered_pi_top15 = adni_filtered_pi.copy()
  661. adni_filtered_pi_top15 = adni_filtered_pi_top15.groupby("feature")["permutation_importance"].apply(get_mean_ci_nb).unstack()
  662. adni_filtered_pi_top15 = adni_filtered_pi_top15.sort_values("mean", ascending=False)
  665. adni_filtered_pi_top15.loc[top_15_features_adni_filtered].applymap(round, n_figs=2).apply(lambda x: f"{x['mean']} ({x['ci_lower']}, {x['ci_upper']})", axis=1)
  666. adni_filtered_pi_top15.applymap(round, n_figs=2).apply(lambda x: f"{x['mean']} ({x['ci_lower']}, {x['ci_upper']})", axis=1).to_csv("ADNI_importance.csv")
  669. # what is unique to top_15_features_oasis_filtered?
  670. set(top_15_features_oasis_filtered) - set(top_15_features_adni_filtered)
  671. # what is unique to top_15_features_adni_filtered?
  672. set(top_15_features_adni_filtered) - set(top_15_features_oasis_filtered)
  674. oasis_filtered_ci_table = oasis_filtered_type_ci.groupby("coalition")["permutation_importance"].apply(get_mean_ci_nb).unstack()
  675. oasis_filtered_ci_table = oasis_filtered_ci_table.sort_values("mean", ascending=False)
  676. oasis_filtered_ci_table.applymap(round, n_figs=2).apply(lambda x: f"{x['mean']} ({x['ci_lower']}, {x['ci_upper']})", axis=1)
  678. adni_filtered_ci_table = adni_filtered_type_ci.groupby("coalition")["permutation_importance"].apply(get_mean_ci_nb).unstack()
  679. adni_filtered_ci_table = adni_filtered_ci_table.sort_values("mean", ascending=False)
  680. adni_filtered_ci_table.applymap(round, n_figs=2).apply(lambda x: f"{x['mean']} ({x['ci_lower']}, {x['ci_upper']})", axis=1)
  683. # %%
  684. # ## 3. Across-Dataset Predictions
  685. # ### 3.1 Cross-Dataset Prediction Functions and Execution
  687. def predict_across_datasets(X_train, y_train, X_test, y_test, root_dir):
  688. root_dir = Path(root_dir)
  689. root_dir.mkdir(parents=True, exist_ok=True)
  691. if root_dir / "pipeline.joblib" in root_dir.iterdir() and root_dir / "predictions.csv" in root_dir.iterdir():
  692. print(f"Skipping {root_dir}, already processed.")
  693. return
  694. X = pd.concat([X_train, X_test])
  695. y = pd.concat([y_train, y_test])
  696. train_idx = np.arange(len(X_train))
  697. test_idx = np.arange(len(X_train), len(X))
  698. # Train a model on the combined data
  699. fit_pipeline(X, y, train_idx, test_idx, root_dir)
  701. # train on oasis, test on adni
  702. across_data_combinations = [
  703. (oasis_y, oasis_X_clin_filtered, adni_y, adni_X_clin_filtered, "OASIS_to_ADNI_clin"),
  704. (oasis_y, oasis_X_fs_filtered, adni_y, adni_X_fs_filtered, "OASIS_to_ADNI_struct"),
  705. (oasis_y, oasis_X_clin_fs_filtered, adni_y, adni_X_clin_fs_filtered, "OASIS_to_ADNI_clin_struct"),
  706. (adni_y, adni_X_clin_filtered, oasis_y, oasis_X_clin_filtered, "ADNI_to_OASIS_clin"),
  707. (adni_y, adni_X_fs_filtered, oasis_y, oasis_X_fs_filtered, "ADNI_to_OASIS_struct"),
  708. (adni_y, adni_X_clin_fs_filtered, oasis_y, oasis_X_clin_fs_filtered, "ADNI_to_OASIS_clin_struct"),
  709. (oasis_y, oasis_X_clin_fs_filtered[top_15_features_oasis_filtered], adni_y, adni_X_clin_fs_filtered[top_15_features_oasis_filtered], "OASIS_to_ADNI_clin_struct_top15"),
  710. (adni_y, adni_X_clin_fs_filtered[top_15_features_adni_filtered], oasis_y, oasis_X_clin_fs_filtered[top_15_features_adni_filtered], "ADNI_to_OASIS_clin_struct_top15"),
  711. ]
  713. for y_train, X_train, y_test, X_test, root_dir in across_data_combinations:
  714. print(f"Processing {root_dir}...")
  715. predict_across_datasets(X_train, y_train, X_test, y_test, root_dir)
  717. # ### 3.2 Cross-Dataset Results Aggregation
  719. def get_across_data_results(root):
  720. results = pd.read_csv(f"{root}/results.csv", index_col=0)
  721. summary = results.groupby(level=0).agg(['mean', 'median'])
  722. return summary
  725. oasis2adni_results = [get_across_data_results(f"OASIS_to_ADNI_{c}") for c in ["clin", "struct", "clin_struct", "clin_struct_top15"]]
  726. oasis2adni_results = pd.concat(oasis2adni_results, keys=["clin", "struct", "clin_struct", "clin_struct_top15"], names=["modality", "target"])
  727. adni2oasis_results = [get_across_data_results(f"ADNI_to_OASIS_{c}") for c in ["clin", "struct", "clin_struct", "clin_struct_top15"]]
  728. adni2oasis_results = pd.concat(adni2oasis_results, keys=["clin", "struct", "clin_struct", "clin_struct_top15"], names=["modality", "target"])
  730. across_results = pd.concat([oasis2adni_results, adni2oasis_results], axis=0,
  731. keys=["OASIS_to_ADNI", "ADNI_to_OASIS"], names=["dataset", "modality", "target"])
  733. # %%
  734. # ### 3.3 Cross-Dataset Training Size Sensitivity Analysis
  736. def predict_across_datasets_with_train_size(X_train, y_train, X_test, y_test, root_dir, train_size, n_splits=1000):
  737. if train_size > len(X_train):
  738. raise ValueError(f"train_size={train_size} exceeds the available training sample ({len(X_train)}).")
  740. ss = ShuffleSplit(n_splits=n_splits, train_size=train_size, random_state=0)
  742. for iter, (train_idx, _) in tqdm(enumerate(ss.split(X_train, y_train)), total=n_splits):
  743. path = Path(f"{root_dir}/iteration_{iter+1}")
  744. path.mkdir(parents=True, exist_ok=True)
  746. if (path / "pipeline.joblib").exists() and (path / "predictions.csv").exists():
  747. continue
  749. X_train_iter = X_train.iloc[train_idx]
  750. y_train_iter = y_train.iloc[train_idx]
  752. X = pd.concat([X_train_iter, X_test])
  753. y = pd.concat([y_train_iter, y_test])
  754. train_idx_full = np.arange(len(X_train_iter))
  755. test_idx_full = np.arange(len(X_train_iter), len(X))
  757. fit_pipeline(X, y, train_idx_full, test_idx_full, path)
  760. def fit_split_multioutput_pipeline(X_train, y_train, X_test, y_test, root_dir):
  761. root_dir = Path(root_dir)
  762. root_dir.mkdir(parents=True, exist_ok=True)
  764. if (root_dir / "pipeline.joblib").exists() and (root_dir / "predictions.csv").exists():
  765. return
  767. y_train = pd.DataFrame(y_train)
  768. y_test = pd.DataFrame(y_test)
  770. imputer, _, _ = get_iterative_imputer(X_train, X_test, random_state=0)
  771. X_train_imputed = imputer.fit_transform(X_train)
  772. X_test_imputed = imputer.transform(X_test)
  774. estimator = MultiOutputRegressor(RandomForestRegressor(random_state=0, n_estimators=50))
  775. estimator.fit(X_train_imputed, y_train)
  777. y_pred = pd.DataFrame(estimator.predict(X_test_imputed)).rename(columns={0: 'mmse_pred', 1: 'sob_pred'})
  778. y_pred.index = y_test.index.values
  780. df = pd.concat([y_pred, y_test], axis=1)
  782. metrics = {
  783. 'r2': r2_score,
  784. 'mae': mean_absolute_error,
  785. 'mse': mean_squared_error,
  786. }
  787. results_df = pd.DataFrame({
  788. var: {name: func(df[f'{var}_slope'], df[f'{var}_pred']) for name, func in metrics.items()}
  789. for var in ['sob', 'mmse']
  790. }).T
  792. save_pipeline_and_df({"imputer": imputer, "estimator": estimator}, df, results_df, root_dir)
  795. def predict_across_datasets_split_multioutput(X_train, y_train, X_test, y_test, root_dir):
  796. fit_split_multioutput_pipeline(X_train, y_train, X_test, y_test, root_dir)
  799. def predict_across_datasets_with_train_size_split_multioutput(X_train, y_train, X_test, y_test, root_dir, train_size, n_splits=1000):
  800. if train_size > len(X_train):
  801. raise ValueError(f"train_size={train_size} exceeds the available training sample ({len(X_train)}).")
  803. ss = ShuffleSplit(n_splits=n_splits, train_size=train_size, random_state=0)
  805. for iter, (train_idx, _) in tqdm(enumerate(ss.split(X_train, y_train)), total=n_splits):
  806. path = Path(f"{root_dir}/iteration_{iter+1}")
  807. path.mkdir(parents=True, exist_ok=True)
  809. if (path / "pipeline.joblib").exists() and (path / "predictions.csv").exists():
  810. continue
  812. fit_split_multioutput_pipeline(X_train.iloc[train_idx], y_train.iloc[train_idx], X_test, y_test, path)
  814. between_train_sizes = [60, 100, 200, 300, 600]
  815. across_train_size_combinations = [
  816. (oasis_y, oasis_X_clin_fs_filtered, adni_y, adni_X_clin_fs_filtered, "OASIS_to_ADNI_clin_struct_train_size_new"),
  817. (adni_y, adni_X_clin_fs_filtered, oasis_y, oasis_X_clin_fs_filtered, "ADNI_to_OASIS_clin_struct_train_size_new"),
  818. ]
  819. across_train_size_split_multioutput_combinations = [
  820. (oasis_y, oasis_X_clin_fs_filtered, adni_y, adni_X_clin_fs_filtered, "OASIS_to_ADNI_clin_struct_train_size_split_multioutput_new"),
  821. (adni_y, adni_X_clin_fs_filtered, oasis_y, oasis_X_clin_fs_filtered, "ADNI_to_OASIS_clin_struct_train_size_split_multioutput_new"),
  822. ]
  824. for y_train, X_train, y_test, X_test, root_dir in across_train_size_combinations:
  825. for train_size in between_train_sizes:
  826. sized_root_dir = f"{root_dir}_{train_size}"
  827. print(f"Processing {sized_root_dir}...")
  828. predict_across_datasets_with_train_size(X_train, y_train, X_test, y_test, sized_root_dir, train_size=train_size, n_splits=10)
  830. for y_train, X_train, y_test, X_test, root_dir in across_train_size_split_multioutput_combinations:
  831. for train_size in between_train_sizes:
  832. sized_root_dir = f"{root_dir}_{train_size}"
  833. print(f"Processing {sized_root_dir}...")
  834. predict_across_datasets_with_train_size_split_multioutput(X_train, y_train, X_test, y_test, sized_root_dir, train_size=train_size, n_splits=10)
  837. def get_across_train_size_results(root_stub, train_sizes=between_train_sizes, iters=range(10)):
  838. results = [get_results(f"{root_stub}_{train_size}", iters) for train_size in train_sizes]
  839. return pd.concat(results, keys=train_sizes, names=["train_size", "iteration", "target"])
  842. def get_across_train_size_model_results(native_root_stub, split_root_stub, train_sizes=between_train_sizes, iters=range(10)):
  843. native_results = get_across_train_size_results(native_root_stub, train_sizes, iters)
  844. split_results = get_across_train_size_results(split_root_stub, train_sizes, iters)
  845. return pd.concat(
  846. [native_results, split_results],
  847. keys=["native multioutput", "split multioutput"],
  848. names=["model", "train_size", "iteration", "target"]
  849. )
  852. oasis2adni_train_size_results = get_across_train_size_model_results(
  853. "OASIS_to_ADNI_clin_struct_train_size_new",
  854. "OASIS_to_ADNI_clin_struct_train_size_split_multioutput_new"
  855. )
  856. adni2oasis_train_size_results = get_across_train_size_model_results(
  857. "ADNI_to_OASIS_clin_struct_train_size_new",
  858. "ADNI_to_OASIS_clin_struct_train_size_split_multioutput_new"
  859. )
  861. across_train_size_results = pd.concat(
  862. [oasis2adni_train_size_results, adni2oasis_train_size_results],
  863. keys=["OASIS_to_ADNI", "ADNI_to_OASIS"],
  864. names=["dataset", "model", "train_size", "iteration", "target"]
  865. )
  866. across_train_size_results.to_csv("results_between_train_sizes.csv")
  869. def get_train_size_rsq_plot(data, y, x, hue, colors, title="", xlabel="R^2", ax=None, ref_val=None, full_sample_points=None, digits=2, full_sample_label="full sample"):
  870. if ax is None:
  871. ax = plt.gca()
  873. categories = data[y].drop_duplicates().tolist()
  874. hue_order = data[hue].drop_duplicates().tolist()
  875. palette = dict(zip(hue_order, colors))
  876. offsets = np.linspace(-0.2, 0.2, len(hue_order)) if len(hue_order) > 1 else np.array([0.0])
  878. ax.axvline(0, color='black', linestyle='--', linewidth=1.2)
  879. if ref_val is not None:
  880. ax.axvline(ref_val, color='black', linestyle='-', linewidth=1.2)
  881. [ax.axvline(p, color='gray', linestyle='--', linewidth=0.3) for p in np.arange(-0.6, 0.8, 0.2)]
  883. sns.boxplot(data=data, y=y, x=x, hue=hue, palette=palette, dodge=True, ax=ax, showfliers=False, linewidth=0.5, order=categories, hue_order=hue_order)
  885. medians = data.groupby([y, hue], sort=False)[x].median()
  886. for i, category in enumerate(categories):
  887. for offset, model in zip(offsets, hue_order):
  888. key = (category, model)
  889. if key in medians.index:
  890. median = medians.loc[key]
  891. ax.text(median, i + offset, f"{median:.{digits}f}", ha='center', va='center', fontdict={'fontsize': 11, 'fontweight':'bold', 'fontname': 'Arial'})
  893. if full_sample_points is not None:
  894. full_sample_pos = len(categories)
  895. for offset, model in zip(offsets, hue_order):
  896. if model in full_sample_points:
  897. point_val = full_sample_points[model]
  898. ax.scatter(point_val, full_sample_pos + offset, color=palette[model], s=35, zorder=5)
  899. ax.text(point_val, full_sample_pos + offset, f"{point_val:.{digits}f}", ha='left', va='center', fontdict={'fontsize': 11, 'fontweight':'bold', 'fontname': 'Arial'})
  900. ax.set_yticks(range(len(categories) + 1), labels=categories + [full_sample_label])
  901. ax.set_ylim(len(categories) + 0.5, -0.5)
  903. ax.tick_params(axis='both', labelfontfamily='Arial', labelsize=12)
  904. ax.set_xlabel(xlabel, fontdict={'fontsize': 12, 'fontstyle':'italic', 'fontweight':'bold', 'fontname': 'Arial'})
  905. ax.set_ylabel('')
  906. ax.set_title(title, fontdict={'fontsize': 16, 'fontweight':'bold', 'fontname': 'Arial'})
  907. ax.legend(title='')
  910. across_train_size_results_plot = across_train_size_results.reset_index()
  911. across_train_size_results_plot["train_size"] = across_train_size_results_plot["train_size"].map(lambda train_size: f"n = {train_size}")
  913. fig, axs = plt.subplots(2, 2, figsize=(10, 6), sharex=True, sharey=False, gridspec_kw={'hspace': 0.45, 'wspace': 0.4})
  914. plot_specs = [
  915. ("OASIS_to_ADNI", "sob", "OASIS-3 → ADNI: CDR-SOB", axs[0, 0]),
  916. ("OASIS_to_ADNI", "mmse", "OASIS-3 → ADNI: MMSE", axs[0, 1]),
  917. ("ADNI_to_OASIS", "sob", "ADNI → OASIS-3: CDR-SOB", axs[1, 0]),
  918. ("ADNI_to_OASIS", "mmse", "ADNI → OASIS-3: MMSE", axs[1, 1]),
  919. ]
  921. for dataset, target, title, ax in plot_specs:
  922. plot_df = across_train_size_results_plot.query("dataset == @dataset and target == @target")
  923. full_sample_points = {
  924. "native multioutput": pd.read_csv(f"{dataset}_clin_struct/results.csv", index_col=0).loc[target, "r2"],
  925. "split multioutput": pd.read_csv(f"{dataset}_clin_struct_split_multioutput/results.csv", index_col=0).loc[target, "r2"],
  926. }
  927. ref_val = full_sample_points["native multioutput"]
  928. get_train_size_rsq_plot(plot_df, "train_size", "r2", "model", ["#E89D9D", "#7DB7E8"], title, "R\u00B2", ax=ax, ref_val=ref_val, full_sample_points=full_sample_points)
  930. # remove legend from all and add a single legend for the entire figure
  931. handles, labels = axs[0, 0].get_legend_handles_labels()
  932. fig.legend(handles, labels, title='Model', loc='upper center', ncol=2, fontsize=12, title_fontsize=12, frameon=False)
  933. # add some margin to the top of the figure to accommodate the legend
  934. fig.subplots_adjust(top=0.85)
  936. # edit all axes to have same xlim
  937. for ax in axs.flatten():
  938. ax.set_xlim(-0.05, 0.55)
  939. ax.get_legend().remove()
  942. plt.tight_layout()
  943. # plt.savefig("r2_perf_between_train_sizes.pdf", dpi=300)
  945. # %%
  946. # ## 4. Statistical Model Comparisons
  947. # ### 4.1 Statistical Testing Framework and Model Comparisons
  949. def r2_test(df, subset1, subset2=None, metric='r2', null_value=0):
  950. dataset1, modality1, target1 = subset1
  952. test_func = wilcoxon
  953. df1 = df.loc[(dataset1, modality1, slice(None), target1)]
  954. if subset2 is not None:
  955. dataset2, modality2, target2 = subset2
  956. df2 = df.loc[(dataset2, modality2, slice(None), target2)]
  957. res = test_func(df1[metric], df2[metric], alternative='two-sided')
  958. proportion = (df1[metric] > df2[metric]).mean()
  959. else:
  960. res = test_func(df1[metric] - null_value, alternative='two-sided')
  961. proportion = (df1[metric] > null_value).mean()
  962. statistic, p_value = res.statistic, res.pvalue
  963. return [statistic, p_value, proportion, metric]
  965. model_comparisons = []
  967. # Within dataset comparisons
  968. for metric in ["r2", "mse", "mae"]:
  969. for target in ["sob", "mmse"]:
  970. for modality1, modality2 in [("clin", "struct"), ("struct", "clin_struct"), ("clin_struct", "clin")]:
  971. for dataset in ["OASIS_filtered", "ADNI_filtered"]:
  972. comparison_label = f"{dataset}_{modality1}_vs_{modality2}_{target}"
  973. comparison = r2_test(results, (dataset, modality1, target), (dataset, modality2, target), metric=metric)
  974. comparison = pd.Series(comparison, index=["statistic", "p_value", "proportion", "metric"], name=comparison_label)
  975. model_comparisons.append(comparison)
  977. # Between dataset comparisons
  978. for target in ["sob", "mmse"]:
  979. for modality in ["clin", "struct", "clin_struct"]:
  980. comparison_label = f"OASIS_filtered_{modality}_vs_ADNI_filtered_{modality}_{target}"
  981. comparison = r2_test(results, ("OASIS_filtered", modality, target), ("ADNI_filtered", modality, target), metric=metric)
  982. comparison = pd.Series(comparison, index=["statistic", "p_value", "proportion", "metric"], name=comparison_label)
  983. model_comparisons.append(comparison)
  985. # Within vs across comparisons
  986. for target in ["sob", "mmse"]:
  987. for modality in ["clin", "struct", "clin_struct"]:
  988. for dataset_across in ["OASIS_to_ADNI", "ADNI_to_OASIS"]:
  989. for dataset_within in ["OASIS_filtered", "ADNI_filtered"]:
  990. comparison_label = f"{dataset_across}_{modality}_vs_{dataset_within}_{modality}_{target}"
  991. null = across_results.loc[dataset_across, modality, target].loc["r2", "mean"]
  992. comparison = r2_test(results, (dataset_within, modality, target), None, null_value=null, metric=metric)
  993. comparison = pd.Series(comparison, index=["statistic", "p_value", "proportion", "metric"], name=comparison_label)
  994. model_comparisons.append(comparison)
  996. model_comparisons = pd.DataFrame(model_comparisons).rename_axis("comparison")
  998. model_comparisons.round(3).loc["OASIS_filtered_clin_vs_struct_sob"]
  999. model_comparisons.round(3).loc["OASIS_filtered_clin_vs_struct_mmse"]
  1000. model_comparisons.loc["OASIS_filtered_clin_struct_vs_clin_sob"]
  1001. model_comparisons.loc["OASIS_filtered_clin_struct_vs_clin_mmse"]
  1002. model_comparisons.loc["OASIS_filtered_struct_vs_clin_struct_sob"]
  1003. model_comparisons.loc["OASIS_filtered_struct_vs_clin_struct_mmse"]
  1004. model_comparisons.loc["OASIS_filtered_clin_vs_struct_sob"]
  1005. model_comparisons.loc["OASIS_filtered_clin_vs_struct_mmse"]
  1007. model_comparisons.loc["ADNI_filtered_clin_struct_vs_clin_sob"]
  1008. model_comparisons.loc["ADNI_filtered_clin_struct_vs_clin_mmse"]
  1009. model_comparisons.loc["ADNI_filtered_struct_vs_clin_struct_sob"]
  1010. model_comparisons.loc["ADNI_filtered_struct_vs_clin_struct_mmse"]
  1011. model_comparisons.loc["ADNI_filtered_clin_vs_struct_sob"]
  1012. model_comparisons.loc["ADNI_filtered_clin_vs_struct_mmse"]
  1014. model_comparisons.loc["OASIS_filtered_clin_struct_vs_ADNI_filtered_clin_struct_sob"]
  1015. model_comparisons.loc["OASIS_filtered_clin_struct_vs_ADNI_filtered_clin_struct_mmse"]
  1017. model_comparisons.loc["OASIS_filtered_clin_vs_ADNI_filtered_clin_sob"]
  1018. model_comparisons.loc["OASIS_filtered_clin_vs_ADNI_filtered_clin_mmse"]
  1020. model_comparisons.loc["OASIS_filtered_struct_vs_ADNI_filtered_struct_sob"]
  1021. model_comparisons.loc["OASIS_filtered_struct_vs_ADNI_filtered_struct_mmse"]
  1023. model_comparisons.loc[
  1024. ["OASIS_to_ADNI_clin_struct_vs_OASIS_filtered_clin_struct_sob",
  1025. "OASIS_to_ADNI_clin_struct_vs_OASIS_filtered_clin_struct_mmse",
  1026. "OASIS_to_ADNI_clin_vs_OASIS_filtered_clin_sob",
  1027. "OASIS_to_ADNI_clin_vs_OASIS_filtered_clin_mmse",
  1028. "OASIS_to_ADNI_struct_vs_OASIS_filtered_struct_sob",
  1029. "OASIS_to_ADNI_struct_vs_OASIS_filtered_struct_mmse",
  1030. "ADNI_to_OASIS_clin_struct_vs_OASIS_filtered_clin_struct_sob",
  1031. "ADNI_to_OASIS_clin_struct_vs_OASIS_filtered_clin_struct_mmse",
  1032. "ADNI_to_OASIS_clin_vs_OASIS_filtered_clin_sob",
  1033. "ADNI_to_OASIS_clin_vs_OASIS_filtered_clin_mmse",
  1034. "ADNI_to_OASIS_struct_vs_OASIS_filtered_struct_sob",
  1035. "ADNI_to_OASIS_struct_vs_OASIS_filtered_struct_mmse"]
  1036. ].query("metric == 'r2'")
  1038. model_comparisons.loc[
  1039. ["OASIS_to_ADNI_clin_struct_vs_ADNI_filtered_clin_struct_sob",
  1040. "OASIS_to_ADNI_clin_struct_vs_ADNI_filtered_clin_struct_mmse",
  1041. "OASIS_to_ADNI_clin_vs_ADNI_filtered_clin_sob",
  1042. "OASIS_to_ADNI_clin_vs_ADNI_filtered_clin_mmse",
  1043. "OASIS_to_ADNI_struct_vs_ADNI_filtered_struct_sob",
  1044. "OASIS_to_ADNI_struct_vs_ADNI_filtered_struct_mmse",
  1045. "ADNI_to_OASIS_clin_struct_vs_ADNI_filtered_clin_struct_sob",
  1046. "ADNI_to_OASIS_clin_struct_vs_ADNI_filtered_clin_struct_mmse",
  1047. "ADNI_to_OASIS_clin_vs_ADNI_filtered_clin_sob",
  1048. "ADNI_to_OASIS_clin_vs_ADNI_filtered_clin_mmse",
  1049. "ADNI_to_OASIS_struct_vs_ADNI_filtered_struct_sob",
  1050. "ADNI_to_OASIS_struct_vs_ADNI_filtered_struct_mmse"]
  1051. ].query("metric == 'r2'")
  1053. target="mmse"
  1054. (results
  1055. .groupby(level=["dataset", "modality", "target"])
  1056. .agg(['mean', 'median'])
  1057. .xs("median", level=1, axis=1)
  1058. .xs(target, level="target", axis=0)
  1059. .loc[["OASIS_filtered", "ADNI_filtered"], ["clin", "struct", "clin_struct"], :]
  1060. .round(2)
  1061. [["r2", "mse", "mae"]])
  1063. across_results.round(2).xs(target, level="target", axis=0).xs("median", level=1, axis=1)[["r2", "mse", "mae"]]
  1065. r2_comparisons = model_comparisons.set_index("metric", append=True)
  1066. r2_comparisons = r2_comparisons.xs("r2", level="metric", axis=0)
  1068. ##
  1069. r2_comparisons.loc[f"OASIS_filtered_clin_vs_struct_{target}"]
  1070. r2_comparisons.loc[f"OASIS_filtered_clin_struct_vs_clin_{target}"]
  1071. r2_comparisons.loc[f"OASIS_filtered_struct_vs_clin_struct_{target}"]
  1072. ##
  1073. r2_comparisons.loc[f"ADNI_filtered_clin_vs_struct_{target}"]
  1074. r2_comparisons.loc[f"ADNI_filtered_clin_struct_vs_clin_{target}"]
  1075. r2_comparisons.loc[f"ADNI_filtered_struct_vs_clin_struct_{target}"]
  1076. ##
  1077. r2_comparisons.loc[f"OASIS_filtered_clin_vs_ADNI_filtered_clin_{target}"]
  1078. r2_comparisons.loc[f"OASIS_filtered_struct_vs_ADNI_filtered_struct_{target}"]
  1079. r2_comparisons.loc[f"OASIS_filtered_clin_struct_vs_ADNI_filtered_clin_struct_{target}"]
  1080. ##
  1081. r2_comparisons.loc[f"OASIS_to_ADNI_clin_vs_OASIS_filtered_clin_{target}"]
  1082. r2_comparisons.loc[f"OASIS_to_ADNI_struct_vs_OASIS_filtered_struct_{target}"]
  1083. r2_comparisons.loc[f"OASIS_to_ADNI_clin_struct_vs_OASIS_filtered_clin_struct_{target}"]
  1084. ##
  1085. r2_comparisons.loc[f"OASIS_to_ADNI_clin_vs_ADNI_filtered_clin_{target}"]
  1086. r2_comparisons.loc[f"OASIS_to_ADNI_struct_vs_ADNI_filtered_struct_{target}"]
  1087. r2_comparisons.loc[f"OASIS_to_ADNI_clin_struct_vs_ADNI_filtered_clin_struct_{target}"]
  1088. ##
  1089. r2_comparisons.loc[f"ADNI_to_OASIS_clin_vs_OASIS_filtered_clin_{target}"]
  1090. r2_comparisons.loc[f"ADNI_to_OASIS_struct_vs_OASIS_filtered_struct_{target}"]
  1091. r2_comparisons.loc[f"ADNI_to_OASIS_clin_struct_vs_OASIS_filtered_clin_struct_{target}"]
  1092. ##
  1093. r2_comparisons.loc[f"ADNI_to_OASIS_clin_vs_ADNI_filtered_clin_{target}"]
  1094. r2_comparisons.loc[f"ADNI_to_OASIS_struct_vs_ADNI_filtered_struct_{target}"]
  1095. r2_comparisons.loc[f"ADNI_to_OASIS_clin_struct_vs_ADNI_filtered_clin_struct_{target}"]
  1096. ##
  1099. # %%
  1100. # ### 4.2 Absolute Error Comparisons Between Models
  1102. def get_abserr(root_dir):
  1103. abserr = pd.read_csv(f"{root_dir}/predictions.csv", index_col=0)
  1104. abserr['sob_abs_error'] = (abserr['sob_slope'] - abserr['sob_pred']).abs()
  1105. abserr['mmse_abs_error'] = (abserr['mmse_slope'] - abserr['mmse_pred']).abs()
  1106. return abserr
  1108. def get_abserr_comparison(root_dir1, root_dir2, paired=False):
  1109. abserr1 = get_abserr(root_dir1)
  1110. abserr2 = get_abserr(root_dir2)
  1112. if paired:
  1113. test_fun = wilcoxon
  1114. # check that abserr1 and abserr2 have the same index
  1115. if not abserr1.index.equals(abserr2.index):
  1116. raise ValueError("The indices of abserr1 and abserr2 must match for paired tests.")
  1117. mmse_proportion = (abserr1['mmse_abs_error'] > abserr2['mmse_abs_error']).mean()
  1118. sob_proportion = (abserr1['sob_abs_error'] > abserr2['sob_abs_error']).mean()
  1119. else:
  1120. test_fun = mannwhitneyu
  1121. mmse_proportion = np.nan
  1122. sob_proportion = np.nan
  1125. mmse_test = test_fun(abserr1['mmse_abs_error'], abserr2['mmse_abs_error'])
  1126. sob_test = test_fun(abserr1['sob_abs_error'], abserr2['sob_abs_error'])
  1129. res = {
  1130. 'mmse_statistic': mmse_test.statistic,
  1131. 'mmse_pvalue': mmse_test.pvalue,
  1132. 'mmse_proportion': mmse_proportion,
  1133. 'sob_statistic': sob_test.statistic,
  1134. 'sob_pvalue': sob_test.pvalue,
  1135. 'sob_proportion': sob_proportion,
  1136. 'paired': paired,
  1137. f"{root_dir1}_mmse_mean": abserr1['mmse_abs_error'].mean(),
  1138. f"{root_dir2}_mmse_mean": abserr2['mmse_abs_error'].mean(),
  1139. f"{root_dir1}_sob_mean": abserr1['sob_abs_error'].mean(),
  1140. f"{root_dir2}_sob_mean": abserr2['sob_abs_error'].mean()
  1141. }
  1143. res = pd.Series(res)
  1145. return res
  1147. # ### 6.6 OASIS-3 → ADNI: Combined vs. Top 15 (Absolute errors)
  1148. get_abserr_comparison("OASIS_to_ADNI_clin_struct", "OASIS_to_ADNI_clin_struct_top15", paired=True)
  1149. # ### 6.7 ADNI → OASIS-3: Combined vs. Top 15
  1150. get_abserr_comparison("ADNI_to_OASIS_clin_struct", "ADNI_to_OASIS_clin_struct_top15", paired=True)
  1151. # ### 6.8 OASIS-3 → ADNI vs. ADNI → OASIS-3: Combined
  1152. get_abserr_comparison("OASIS_to_ADNI_clin_struct", "ADNI_to_OASIS_clin_struct", paired=False)
  1153. # ### 6.9 OASIS-3 → ADNI vs. ADNI → OASIS-3: Top 15
  1154. get_abserr_comparison("OASIS_to_ADNI_clin_struct_top15", "ADNI_to_OASIS_clin_struct_top15", paired=False)
  1156. across_results.loc[:, "clin_struct", :]
  1159. # %%
  1160. # ## 5. Subgroup Analysis
  1161. # ### 5.1 Preparation of Data for Subgroup Comparisons
  1163. subset_adni = pd.concat([adni_X_clin_filtered, get_abserr("OASIS_to_ADNI_clin_struct")], axis=1)
  1164. subset_oasis = pd.concat([oasis_X_clin_fs_filtered, get_abserr("ADNI_to_OASIS_clin_struct")], axis=1)
  1165. subset_adni_top15 = pd.concat([adni_X_clin_filtered, get_abserr("OASIS_to_ADNI_clin_struct_top15")], axis=1)
  1166. subset_oasis_top15 = pd.concat([oasis_X_clin_fs_filtered, get_abserr("ADNI_to_OASIS_clin_struct_top15")], axis=1)
  1168. # ### 5.2 Subgroup Comparison Functions and Analysis
  1170. # accessor function to make subset comparisons
  1171. def get_subset(variable, value, df, comparison_type='equal'):
  1172. if comparison_type == 'equal':
  1173. return df[df[variable] == value]
  1174. elif comparison_type == 'greater':
  1175. return df[df[variable] > value]
  1176. elif comparison_type == 'less':
  1177. return df[df[variable] < value]
  1178. elif comparison_type == 'different':
  1179. return df[df[variable] != value]
  1180. elif comparison_type == 'greater_equal':
  1181. return df[df[variable] >= value]
  1182. else:
  1183. raise ValueError("Invalid comparison type.")
  1185. def get_subset_comparison(df, variable, comparison_1, comparison_2, label=None):
  1186. (value1, comparison_type1) = comparison_1
  1187. (value2, comparison_type2) = comparison_2
  1188. subset1 = get_subset(variable, value1, df, comparison_type1)
  1189. subset2 = get_subset(variable, value2, df, comparison_type2)
  1191. sob_test = mannwhitneyu(subset1['sob_abs_error'], subset2['sob_abs_error'])
  1192. mmse_test = mannwhitneyu(subset1['mmse_abs_error'], subset2['mmse_abs_error'])
  1194. sob_mae1 = subset1['sob_abs_error'].mean()
  1195. mmse_mae1 = subset1['mmse_abs_error'].mean()
  1196. sob_mae2 = subset2['sob_abs_error'].mean()
  1197. mmse_mae2 = subset2['mmse_abs_error'].mean()
  1199. return pd.Series({
  1200. 'sob_statistic': sob_test.statistic,
  1201. 'sob_pvalue': sob_test.pvalue,
  1202. 'mmse_statistic': mmse_test.statistic,
  1203. 'mmse_pvalue': mmse_test.pvalue,
  1204. 'sob_mae1': sob_mae1,
  1205. 'mmse_mae1': mmse_mae1,
  1206. 'sob_mae2': sob_mae2,
  1207. 'mmse_mae2': mmse_mae2
  1208. }, name=label)
  1210. def subgroup_analysis(subset):
  1211. res = []
  1212. # ### 8.1 Subgroups Segregated by Diagnosis
  1213. # ##### HC (0) vs. MCI (1) and AD (2)
  1214. res.append(get_subset_comparison(subset, 'diag', (0, 'equal'), (0, 'different'), label='HC vs. MCI and AD'))
  1215. # #### MCI vs. HC and AD
  1216. res.append(get_subset_comparison(subset, 'diag', (1, 'equal'), (1, 'different'), label='MCI vs. HC and AD'))
  1217. # #### AD vs. HC and MCI
  1218. res.append(get_subset_comparison(subset, 'diag', (2, 'equal'), (2, 'different'), label='AD vs. HC and MCI'))
  1219. # #### HC vs. MCI
  1220. res.append(get_subset_comparison(subset, 'diag', (0, 'equal'), (1, 'equal'), label='HC vs. MCI'))
  1221. # #### HC vs. AD
  1222. res.append(get_subset_comparison(subset, 'diag', (0, 'equal'), (2, 'equal'), label='HC vs. AD'))
  1223. # #### MCI vs. AD
  1224. res.append(get_subset_comparison(subset, 'diag', (1, 'equal'), (2, 'equal'), label='MCI vs. AD'))
  1225. # ### 8.2 Subgroups Segregated by Sex
  1226. # #### Female (0) vs. Male (1)
  1227. res.append(get_subset_comparison(subset, 'demo_sex', (0, 'equal'), (1, 'equal'), label='Female vs. Male'))
  1228. # ### 8.3 Subgroups Segregated by APOE E4
  1229. # #### APOE E4: 0 vs. 1 and 2
  1230. res.append(get_subset_comparison(subset, 'apoe_e4count', (0, 'equal'), (0, 'greater'), label='APOE E4: 0 vs. 1 and 2'))
  1231. # #### APOE E4: 1 vs. 0 and 2
  1232. res.append(get_subset_comparison(subset, 'apoe_e4count', (1, 'equal'), (1, 'different'), label='APOE E4: 1 vs. 0 and 2'))
  1233. # #### APOE E4: 2 vs. 0 and 1
  1234. res.append(get_subset_comparison(subset, 'apoe_e4count', (2, 'equal'), (2, 'less'), label='APOE E4: 2 vs. 0 and 1'))
  1235. # #### APOE E4: 0 vs. 1
  1236. res.append(get_subset_comparison(subset, 'apoe_e4count', (0, 'equal'), (1, 'equal'), label='APOE E4: 0 vs. 1'))
  1237. # #### APOE E4: 0 vs. 2
  1238. res.append(get_subset_comparison(subset, 'apoe_e4count', (0, 'equal'), (2, 'equal'), label='APOE E4: 0 vs. 2'))
  1239. # #### APOE E4: 1 vs. 2
  1240. res.append(get_subset_comparison(subset, 'apoe_e4count', (1, 'equal'), (2, 'equal'), label='APOE E4: 1 vs. 2'))
  1241. # ### 8.4 Subgroups Segregated by Age
  1242. subset['demo_age_group'] = pd.cut(subset['demo_age'], bins=[-np.inf, 65, 70, 75, 80, np.inf], labels=['<65', '65-70', '70-75', '75-80', '>80'])
  1243. # #### AGE: <65 vs. >=65 (65-70 and 70-75 and 75-80 and and >80)
  1244. res.append(get_subset_comparison(subset, 'demo_age_group', ('65-70', 'less'), ('65-70', 'greater_equal'), label='Age: <65 vs. >=65'))
  1245. # #### AGE: 70-75 vs. <70 and >= 75
  1246. res.append(get_subset_comparison(subset, 'demo_age_group', ('70-75', 'less'), ('70-75', 'greater_equal'), label='Age: <70 vs. >= 70'))
  1247. # #### AGE: 75-80 vs. <75 and >= 80
  1248. res.append(get_subset_comparison(subset, 'demo_age_group', ('75-80', 'less'), ('75-80', 'greater_equal'), label='Age: <75 vs. >= 80'))
  1249. # #### AGE: >= 80 vs. <80
  1250. res.append(get_subset_comparison(subset, 'demo_age_group', ('>80', 'less'), ('>80', 'greater_equal'), label='Age: <80 vs. >= 80'))
  1252. res = pd.DataFrame(res)
  1254. return res
  1256. res_subset_adni = subgroup_analysis(subset_adni)[["sob_pvalue", "sob_mae1", "sob_mae2", "mmse_pvalue", "mmse_mae1", "mmse_mae2"]]
  1257. res_subset_oasis = subgroup_analysis(subset_oasis)[["sob_pvalue", "sob_mae1", "sob_mae2", "mmse_pvalue", "mmse_mae1", "mmse_mae2"]]
  1259. # ### 5.3 Results Formatting and Export
  1261. def apa_format(column):
  1262. """
  1263. Format a dataframe column in APA style.
  1264. """
  1265. # first check if the column name contains 'pvalue'
  1266. if 'pvalue' in column.name:
  1267. # format p-values
  1268. formatted = column.apply(lambda x: f"{x:.3f}" if x > 0.001 else "< .001")
  1269. else:
  1270. formatted = column.apply(lambda x: round(x, 2) if isinstance(x, float) else str(x))
  1271. return formatted
  1274. # correct for multiple comparisons using Holm-Bonferroni method
  1275. from statsmodels.stats.multitest import multipletests
  1277. res_subset_adni["sob_pvalue_corrected"] = multipletests(res_subset_adni["sob_pvalue"], method='holm')[1]
  1278. res_subset_adni["mmse_pvalue_corrected"] = multipletests(res_subset_adni["mmse_pvalue"], method='holm')[1]
  1280. res_subset_oasis["sob_pvalue_corrected"] = multipletests(res_subset_oasis["sob_pvalue"], method='holm')[1]
  1281. res_subset_oasis["mmse_pvalue_corrected"] = multipletests(res_subset_oasis["mmse_pvalue"], method='holm')[1]
  1283. # Format the results in APA style
  1284. res_subset_adni.apply(apa_format).to_csv("ADNI_subgroup_analysis.csv")
  1285. res_subset_oasis.apply(apa_format).to_csv("OASIS_subgroup_analysis.csv")
  1288. res_subset_adni.apply(apa_format)[["sob_pvalue_corrected", "mmse_pvalue_corrected"]]
  1289. res_subset_oasis.apply(apa_format)[["sob_pvalue_corrected", "mmse_pvalue_corrected"]]
  1290. # %%
  1291. # ## 6. Post-Hoc Matched Samples Analysis
  1292. # ### 6.1 Quantile-Based Sample Matching
  1294. def get_matched_samples(source_y, target_y, quantiles=np.arange(0, 1.1, 0.1)):
  1295. """
  1296. Get matched samples based on quantiles of the source dataset.
  1297. """
  1298. quant_source_sob = source_y["sob_slope"].quantile(quantiles)
  1299. quant_source_mmse = source_y["mmse_slope"].quantile(quantiles)
  1301. target_sob_cut = pd.cut(target_y["sob_slope"], quant_source_sob, include_lowest=True, duplicates="drop")
  1302. target_mmse_cut = pd.cut(target_y["mmse_slope"], quant_source_mmse, include_lowest=True, duplicates="drop")
  1304. # number of people in target with each of the categories so that it matches the proportion in source
  1305. source_sob_n = pd.cut(source_y["sob_slope"], quant_source_sob, include_lowest=True, duplicates="drop").value_counts()
  1306. source_mmse_n = pd.cut(source_y["mmse_slope"], quant_source_mmse, include_lowest=True, duplicates="drop").value_counts()
  1308. return target_sob_cut, target_mmse_cut, source_sob_n, source_mmse_n
  1311. # trained on OASIS-3, tested on ADNI
  1312. ADNI_sob_cut, ADNI_mmse_cut, OASIS_sob_n, OASIS_mmse_n = get_matched_samples(oasis_y, adni_y)
  1313. # trained on ADNI, tested on OASIS-3
  1314. OASIS_sob_cut, OASIS_mmse_cut, ADNI_sob_n, ADNI_mmse_n = get_matched_samples(adni_y, oasis_y)
  1316. subset_adni = pd.concat([subset_adni, ADNI_sob_cut.rename("sob_qcat"), ADNI_mmse_cut.rename("mmse_qcat")], axis=1)
  1317. subset_oasis = pd.concat([subset_oasis, OASIS_sob_cut.rename("sob_qcat"), OASIS_mmse_cut.rename("mmse_qcat")], axis=1)
  1318. subset_adni_top15 = pd.concat([subset_adni_top15, ADNI_sob_cut.rename("sob_qcat"), ADNI_mmse_cut.rename("mmse_qcat")], axis=1)
  1319. subset_oasis_top15 = pd.concat([subset_oasis_top15, OASIS_sob_cut.rename("sob_qcat"), OASIS_mmse_cut.rename("mmse_qcat")], axis=1)
  1322. # ### 6.2 Matched Sample Predictions and Evaluation
  1324. # ### 10.3 OASIS-3 → ADNI (Combined)
  1326. # %%
  1327. feature_names_train = oasis_X_clin_fs_filtered.columns.tolist()
  1330. # %%
  1332. # %%
  1335. def subset_results(subset, SOB_n, MMSE_n):
  1336. qcat_results = []
  1337. for _ in tqdm(range(0, 1000)):
  1338. df_sob = subset.groupby("sob_qcat", as_index=False, group_keys=False, observed=True).apply(lambda x: x.sample(n = SOB_n[x.name], replace=True), include_groups=False)
  1339. df_mmse = subset.groupby("mmse_qcat", as_index=False, group_keys=False, observed=True).apply(lambda x: x.sample(n = MMSE_n[x.name], replace=True), include_groups=False)
  1340. sob_test = df_sob['sob_slope']
  1341. mmse_test = df_mmse['mmse_slope']
  1342. sob_pred = df_sob['sob_pred']
  1343. mmse_pred = df_mmse['mmse_pred']
  1345. qcat_results.append({
  1346. 'sob_r2': sklearn.metrics.r2_score(sob_test, sob_pred),
  1347. 'sob_mse': sklearn.metrics.mean_squared_error(sob_test, sob_pred),
  1348. 'sob_mae': sklearn.metrics.mean_absolute_error(sob_test, sob_pred),
  1349. 'mmse_r2': sklearn.metrics.r2_score(mmse_test, mmse_pred),
  1350. 'mmse_mse': sklearn.metrics.mean_squared_error(mmse_test, mmse_pred),
  1351. 'mmse_mae': sklearn.metrics.mean_absolute_error(mmse_test, mmse_pred),
  1352. })
  1354. qcat_results = pd.DataFrame(qcat_results)
  1356. return qcat_results
  1358. # ### 6.3 Results for Both Combined and Top-15 Models
  1360. subset_adni_results = subset_results(subset_adni, OASIS_sob_n, OASIS_mmse_n)
  1361. subset_adni_top15_results = subset_results(subset_adni_top15, OASIS_sob_n, OASIS_mmse_n)
  1362. subset_oasis_results = subset_results(subset_oasis, ADNI_sob_n, ADNI_mmse_n)
  1363. subset_oasis_top15_results = subset_results(subset_oasis_top15, ADNI_sob_n, ADNI_mmse_n)
  1365. subset_adni_results.agg(["median"]).map(round, n_figs=2)
  1366. subset_adni_top15_results.agg(["median"]).map(round, n_figs=2)
  1368. subset_oasis_results.agg(["median"]).map(round, n_figs=2)
  1369. subset_oasis_top15_results.agg(["median"]).map(round, n_figs=2)
  1371. # %%

01_Analysis.py at commit 44acd1b, under MIT · at the source

Overview

Authors: Roya Melanie Hüppi1,2,3, Nicolas Langer1,2, Bruno Hebling Vieira1,2, for the Alzheimer’s Disease Neuroimaging Initiative
ORCID iDs: Roya Melanie Hüppi, Nicolas Langer, Bruno Hebling Vieira
  1. Methods of Plasticity Research, Department of Psychology, University of Zurich, 8050, Zurich, Zurich, Switzerland
  2. Neuroscience Center Zurich (ZNZ), University of Zurich & ETH Zurich, 8057, Zurich, Zurich, Switzerland
  3. Department of Adult Psychiatry and Psychotherapy, Psychiatric University Clinic Zurich and University of Zurich, 8032, Zurich, Zurich, Switzerland
Institutions: University of Zurich (Switzerland); ETH Zurich (Switzerland); Psychiatrische Universitätsklinik Zürich (Switzerland)
Journal: The journal of prevention of Alzheimer's disease, volume 13, issue 9, article 100646
Dates: received 14 April 2026; accepted 22 June 2026; published online 8 August 2026; in print August 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1016/j.tjpad.2026.100646 · PMID 42570468 · PMCID PMC13476611 · OpenAlex W7201989932
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: structural MRI / diffusion (modality), human (organism), Alzheimer's / dementia (population)
Methods: Statistics, Machine learning, Preprocessing
Keywords: Cognitive decline, Structural MRI, Generalizability, Machine learning, Predictive modeling
Topic: Neural Networks and Applications (Artificial Intelligence, Computer Science), according to OpenAlex
Funding: Swiss National Science Foundation (10001C_197480); UZH Postdoc (FK-23-086)
Citations: not cited yet (Europe PMC); 66 references in the paper

Abstract

Background: Predicting cognitive decline as a continuum, from healthy age-related decline to mild cognitive impairment and dementia, enables more precise individual-level predictions. However, the practical value of such models for early intervention and prevention depends on their ability to generalize to independent cohorts, a property that is often not evaluated.

Objectives: This study investigated whether adding structural magnetic resonance imaging (MRI) to non-brain data improved machine learning predictions of continuous cognitive decline and analyzed the models’ generalizability.

Design: Multi-target random forest regression models predicted annual decline in the Clinical Dementia Rating Scale Sum of Boxes (CDR-SOB) and Mini-Mental State Examination (MMSE) using non-brain data, structural MRI data, or their combination from the Alzheimer's Disease Neuroimaging Initiative (ADNI; N = 1237) and Open Access Series of Imaging Studies (OASIS-3; N = 662) datasets. Cross-site generalizability was evaluated.

Setting: Data from ADNI and OASIS-3 were used for this study.

Participants: A total of 1899 participants who had demographic, clinical, and brain imaging data from a baseline session and clinical data from at least 2 follow-up sessions were included.

Measurements: Baseline non-brain (demographics, clinical and neuropsychological scores, information on APOE genotype, cognitive diagnosis, health, and number of sessions before baseline) and/or structural MRI data were used to predict the yearly rate of change in CDR-SOB and MMSE scores.

Results: Including structural MRI data improved prediction of CDR-SOB and MMSE change, reaching respective R2 values of .41 and .33 in ADNI and .42 and .33 in OASIS-3. Model performance for across-dataset predictions was reduced (R2 between .18 and .35), unexplained by distributional shifts of target variables. Models using only top predictive features performed similarly to full models when tested externally (R2 between .18 and .34), suggesting predictor redundancy.

Conclusions: Incorporating structural MRI data enhances within-dataset prediction of continuous cognitive decline, allowing for more precise individual-level prediction and advancing towards precision medicine. Even though external validation remains limited, quantifying the generalizability gap is a crucial step towards the responsible use of ML models in clinical intervention and prevention.

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

Repository

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

rhuepp/gemacode

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 44acd1bbbf67d866b6f432c8b701b6de391fc6a8, 21 September 2026
Languages: Python (2)
Size: 4 files, 2 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (2 files), NumPy (2 files), pandas (2 files), scikit-learn (2 files), seaborn (2 files), SciPy (1 file), statsmodels (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
4 files
At the source: github.com/rhuepp/gemacode.git

Code availability

The code is available on github.com/rhuepp/gemacode (https://github.com/rhuepp/gemacode.git).

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

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:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 2 scripts, each with its path and the digest of its content;
  • 4 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

No dataset and no data link were found in the paper.

Data availability statement

ADNI data can be obtained from https://adni.loni.usc.edu/, OASIS-3 data from https://www.oasis-brains.org.

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

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 3, 28 September 2026

  • Authors: added Nicolas Langer (0000-0002-6038-9471); removed Nicolas Langer

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 5 keywords, 2 funders, 53 references.

Cite

This paper

Hüppi, R. M., Langer, N., Hebling Vieira, B., & for the Alzheimer’s Disease Neuroimaging Initiative. (2026). Quantifying generalization error in machine learning prediction of cognitive decline. The journal of prevention of Alzheimer's disease, 13(9), 100646. https://doi.org/10.1016/j.tjpad.2026.100646

BibTeX

@article{huppi2026quantifying,
author = {Hüppi, Roya Melanie and Langer, Nicolas and Hebling Vieira, Bruno and {for the Alzheimer’s Disease Neuroimaging Initiative}},
title = {{Quantifying generalization error in machine learning prediction of cognitive decline}},
journal = {The journal of prevention of Alzheimer's disease},
year = {2026},
month = aug,
volume = {13},
number = {9},
pages = {100646},
publisher = {Elsevier},
issn = {2274-5807},
doi = {10.1016/j.tjpad.2026.100646},
url = {https://doi.org/10.1016/j.tjpad.2026.100646},
pmid = {42570468},
pmcid = {PMC13476611}
}

RIS

TY - JOUR
AU - Hüppi, Roya Melanie
AU - Langer, Nicolas
AU - Hebling Vieira, Bruno
AU - for the Alzheimer’s Disease Neuroimaging Initiative
TI - Quantifying generalization error in machine learning prediction of cognitive decline
T2 - The journal of prevention of Alzheimer's disease
J2 - J Prev Alzheimers Dis
PY - 2026
DA - 2026/08/08
VL - 13
IS - 9
SP - 100646
SN - 2274-5807
PB - Elsevier
DO - 10.1016/j.tjpad.2026.100646
UR - https://doi.org/10.1016/j.tjpad.2026.100646
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.tjpad.2026.100646",
"type": "article-journal",
"title": "Quantifying generalization error in machine learning prediction of cognitive decline",
"container-title": "The journal of prevention of Alzheimer's disease",
"author": [
{
"family": "Hüppi",
"given": "Roya Melanie"
},
{
"family": "Langer",
"given": "Nicolas"
},
{
"family": "Hebling Vieira",
"given": "Bruno"
},
{
"literal": "for the Alzheimer’s Disease Neuroimaging Initiative"
}
],
"container-title-short": "J Prev Alzheimers Dis",
"volume": "13",
"issue": "9",
"page": "100646",
"DOI": "10.1016/j.tjpad.2026.100646",
"PMID": "42570468",
"PMCID": "PMC13476611",
"ISSN": "2274-5807",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.tjpad.2026.100646",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
8
]
]
}
}

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.1038/s41514-026-00390-w [code]
Mild cognitive impairment cases affect the predictive power of Alzheimer's disease diagnostic models using routine clinical variables.
Journal: npj aging
In common: seaborn, scikit-learn, pandas, 3 other tools, Alzheimer's / dementia, 2 references
[2] doi:10.1038/s42003-026-10205-z [code]
Source-space EEG alpha activity reveals brain age gaps due to neurodegeneration and disparity.
Journal: Communications biology
In common: statsmodels, seaborn, scikit-learn, 4 other tools, Alzheimer's / dementia, structural MRI / diffusion, 1 reference
[3] doi:10.21203/rs.3.rs-9914920/v1 [code]
Prediction of cognitive performance by demographics, sleep, and brain morphometry: machine learning findings from ENIGMA-Sleep Working Group
Journal: Research Square (preprint)
In common: statsmodels, seaborn, scikit-learn, 4 other tools, structural MRI / diffusion, 1 reference
[4] doi:10.1111/ejn.70480 [code]
Astrocyte Proximity Protects Synapses From Human Amyloid-Beta Induced Degeneration in a Mouse Ex Vivo Model of Early Alzheimer's Disease.
Journal: The European journal of neuroscience
In common: statsmodels, seaborn, scikit-learn, 4 other tools, Alzheimer's / dementia, 1 reference
[5] doi:10.1002/advs.202600020 [code]
Early Retinal UCHL1 Dysregulation Coupled With Synaptic Loss Reflects Alzheimer's Disease Severity.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: seaborn, scikit-learn, pandas, 3 other tools, Alzheimer's / dementia, 2 references
[6] doi:10.1162/imag.a.1208 [code]
Brain network analysis in Alzheimer's disease and mild cognitive impairment using high-density diffuse optical tomography.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: statsmodels, seaborn, scikit-learn, 4 other tools, Alzheimer's / dementia, 1 reference
[7] doi:10.1371/journal.pone.0343722 [code]
Comprehensive methodology for sample enrichment in EEG biomarker studies for Alzheimer's risk classification.
Journal: PloS one
In common: statsmodels, seaborn, scikit-learn, 4 other tools, Alzheimer's / dementia, 1 reference
[8] doi:10.3389/fnagi.2026.1847611 [code]
APOE ε4-associated hippocampal atrophy trajectories across the Alzheimer's disease continuum: a systematic review, meta-analysis, and longitudinal validation.
Journal: Frontiers in aging neuroscience
In common: statsmodels, seaborn, pandas, 3 other tools, Alzheimer's / dementia, structural MRI / diffusion, 1 reference
[9] doi:10.1038/s41467-026-75661-x [code]
A neural signature of sleep deprivation in the human brain.
Journal: Nature communications
In common: statsmodels, seaborn, scikit-learn, 4 other tools, 1 reference
[10] doi:10.1093/nc/niag029 [code]
A data-driven approach to identifying and evaluating connectivity-based neural correlates of conscious visual perception.
Journal: Neuroscience of consciousness
In common: statsmodels, seaborn, scikit-learn, 4 other tools, 1 reference

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.