OSCR

Biological brain aging, cognitive-motor decline and vascular risk: a multivariate imaging analysis of 40,579 individuals.

Code ↔ Paper

3 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 3 matches
  1. [1] § Materials and methods › Statistics › Partial least square correlation analysis ↔ pyls/base.py, lines 438–527 · score 0.77 · singular vector, standard error, bootstrap resampling, bootstrap ratio, replacement, weight
  2. [2] § Materials and methods › Statistics › Partial least square correlation analysis ↔ pyls/structures.py, lines 299–326 · score 0.68 · confidence intervals, singular vector, latent variable, permuting, resampling, correlation
  3. [3] § Results › Imaging markers of biological brain aging are associated with cognitive and motor function ↔ pyls/types/behavioral.py, lines 172–227 · score 0.59 · confidence interval, cross validation, bootstrap ratio, squares, scores, PLS

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 · 759 lines · 31 KB · GPL-2.0 · 1 match

  1. # -*- coding: utf-8 -*-
  3. import gc
  4. import warnings
  5. import numpy as np
  6. from sklearn.utils.validation import check_random_state
  7. from . import compute, structures, utils
  10. def gen_permsamp(groups, n_cond, n_perm, seed=None, verbose=True):
  11. """
  12. Generates permutation arrays for PLS permutation testing
  14. Parameters
  15. ----------
  16. groups : (G,) list
  17. List with number of subjects in each of `G` groups
  18. n_cond : int
  19. Number of conditions, for each subject. Default: 1
  20. n_perm : int
  21. Number of permutations for which to generate resampling arrays
  22. seed : {int, :obj:`numpy.random.RandomState`, None}, optional
  23. Seed for random number generation. Default: None
  24. verbose : bool, optional
  25. Whether to print status updates as permutations are generated.
  26. Default: True
  28. Returns
  29. -------
  30. permsamp : (S, P) `numpy.ndarray`
  31. Subject permutation arrays, where `S` is the number of subjects and `P`
  32. is the requested number of permutations (i.e., `P = n_perm`)
  33. """
  34. Y = utils.dummy_code(groups, n_cond)
  35. permsamp = np.zeros(shape=(len(Y), n_perm), dtype=int)
  36. subj_inds = np.arange(np.sum(groups), dtype=int)
  37. rs = check_random_state(seed)
  38. warned = False
  40. # calculate some variables for permuting conditions within subject
  41. # do this here to save on calculation time
  42. indices, grps = np.where(Y)
  43. grp_conds = np.split(indices, np.where(np.diff(grps))[0] + 1)
  44. to_permute = [np.vstack(grp_conds[i:i + n_cond]) for i in
  45. range(0, Y.shape[-1], n_cond)]
  46. splitinds = np.cumsum(groups)[:-1]
  47. check_grps = utils.dummy_code(groups).T.astype(bool)
  49. for i in utils.trange(n_perm, verbose=verbose, desc='Making permutations'):
  50. count, duplicated = 0, True
  51. while duplicated and count < 500:
  52. count, duplicated = count + 1, False
  53. # generate conditions permuted w/i subject
  54. inds = np.hstack([utils.permute_cols(i, seed=rs) for i
  55. in to_permute])
  56. # generate permutation of subjects across groups
  57. perm = rs.permutation(subj_inds)
  58. # confirm subjects *are* mixed across groups
  59. if len(groups) > 1:
  60. for grp in check_grps:
  61. if np.all(np.sort(perm[grp]) == subj_inds[grp]):
  62. duplicated = True
  63. # permute conditions w/i subjects across groups and stack
  64. perminds = np.hstack([f.flatten('F') for f in
  65. np.split(inds[:, perm].T, splitinds)])
  66. # make sure permuted indices are not a duplicate sequence
  67. dupe_seq = perminds[:, None] == permsamp[:, :i]
  68. if dupe_seq.all(axis=0).any():
  69. duplicated = True
  70. # if we broke out because we tried 500 permutations and couldn't
  71. # generate a new one, just warn that we're using duplicate
  72. # permutations and give up
  73. if count == 500 and not warned:
  74. warnings.warn('WARNING: Duplicate permutations used.')
  75. warned = True
  76. # store the permuted indices
  77. permsamp[:, i] = perminds
  79. return permsamp
  82. def gen_bootsamp(groups, n_cond, n_boot, seed=None, verbose=True):
  83. """
  84. Generates bootstrap arrays for PLS bootstrap resampling
  86. Parameters
  87. ----------
  88. groups : (G,) list
  89. List with number of subjects in each of `G` groups
  90. n_cond : int
  91. Number of conditions, for each subject. Default: 1
  92. n_boot : int
  93. Number of boostraps for which to generate resampling arrays
  94. seed : {int, :obj:`numpy.random.RandomState`, None}, optional
  95. Seed for random number generation. Default: None
  96. verbose : bool, optional
  97. Whether to print status updates as bootstrap samples are genereated.
  98. Default: True
  100. Returns
  101. -------
  102. bootsamp : (S, B) `numpy.ndarray`
  103. Subject bootstrap arrays, where `S` is the number of subjects and `B`
  104. is the requested number of bootstraps (i.e., `B = n_boot`)
  105. """
  106. Y = utils.dummy_code(groups, n_cond)
  107. bootsamp = np.zeros(shape=(len(Y), n_boot), dtype=int)
  108. subj_inds = np.arange(np.sum(groups), dtype=int)
  109. rs = check_random_state(seed)
  110. warned = False
  111. min_subj = int(np.ceil(Y.sum(axis=0).min() * 0.5))
  113. # calculate some variables for ensuring we resample with replacement
  114. # subjects across all their conditions. do this here to save on
  115. # calculation time
  116. indices, grps = np.where(Y)
  117. grp_conds = np.split(indices, np.where(np.diff(grps))[0] + 1)
  118. inds = np.hstack([np.vstack(grp_conds[i:i + n_cond]) for i
  119. in range(0, len(grp_conds), n_cond)])
  120. splitinds = np.cumsum(groups)[:-1]
  121. check_grps = utils.dummy_code(groups).T.astype(bool)
  123. for i in utils.trange(n_boot, verbose=verbose, desc='Making bootstraps'):
  124. count, duplicated = 0, True
  125. while duplicated and count < 500:
  126. count, duplicated = count + 1, False
  127. # empty container to store current bootstrap attempt
  128. boot = np.zeros(shape=(subj_inds.size), dtype=int)
  129. # iterate through and resample from w/i groups
  130. for grp in check_grps:
  131. curr_grp, all_same = subj_inds[grp], True
  132. while all_same:
  133. num_subj = curr_grp.size
  134. boot[curr_grp] = np.sort(rs.choice(curr_grp,
  135. size=num_subj,
  136. replace=True),
  137. axis=0)
  138. # make sure bootstrap has enough unique subjs
  139. if np.unique(boot[curr_grp]).size >= min_subj:
  140. all_same = False
  141. # resample subjects (with conditions) and stack groups
  142. bootinds = np.hstack([f.flatten('F') for f in
  143. np.split(inds[:, boot].T, splitinds)])
  144. # make sure bootstrap is not a duplicated sequence
  145. for grp in check_grps:
  146. curr_grp = subj_inds[grp]
  147. check = bootinds[curr_grp, None] == bootsamp[curr_grp, :i]
  148. if check.all(axis=0).any():
  149. duplicated = True
  150. # if we broke out because we tried 500 bootstraps and couldn't
  151. # generate a new one, just warn that we're using duplicate
  152. # bootstraps and give up
  153. if count == 500 and not warned:
  154. warnings.warn('WARNING: Duplicate bootstraps used.')
  155. warned = True
  156. # store the bootstrapped indices
  157. bootsamp[:, i] = bootinds
  159. return bootsamp
  162. def gen_splits(groups, n_cond, n_split, seed=None, test_size=0.5):
  163. """
  164. Generates splitting arrays for PLS split-half resampling and CV
  166. Parameters
  167. ----------
  168. groups : (G,) list
  169. List with number of subjects in each of `G` groups
  170. n_cond : int
  171. Number of conditions, for each subject. Default: 1
  172. n_split : int
  173. Number of splits for which to generate resampling arrays
  174. seed : {int, :obj:`numpy.random.RandomState`, None}, optional
  175. Seed for random number generation. Default: None
  176. test_size : (0, 1) float, optional
  177. Percent of subjects to include in the split halves. Default: 0.5
  179. Returns
  180. -------
  181. splitsamp : (S, I) `numpy.ndarray`
  182. Subject split arrays, where `S` is the number of subjects and `I`
  183. is the requested number of splits (i.e., `I = n_split`)
  184. """
  185. Y = utils.dummy_code(groups, n_cond)
  186. splitsamp = np.zeros(shape=(len(Y), n_split), dtype=bool)
  187. subj_inds = np.arange(np.sum(groups), dtype=int)
  188. rs = check_random_state(seed)
  189. warned = False
  191. # calculate some variables for permuting conditions within subject
  192. # do this here to save on calculation time
  193. indices, grps = np.where(Y)
  194. grp_conds = np.split(indices, np.where(np.diff(grps))[0] + 1)
  195. inds = np.hstack([np.vstack(grp_conds[i:i + n_cond]) for i
  196. in range(0, len(grp_conds), n_cond)])
  197. splitinds = np.cumsum(groups)[:-1]
  198. check_grps = utils.dummy_code(groups).T.astype(bool)
  200. for i in range(n_split):
  201. count, duplicated = 0, True
  202. while duplicated and count < 500:
  203. count, duplicated = count + 1, False
  204. # empty containter to store current split half attempt
  205. split = np.zeros(shape=(subj_inds.size), dtype=bool)
  206. # iterate through and split each group separately
  207. for grp in check_grps:
  208. curr_grp = subj_inds[grp]
  209. take = rs.choice([np.ceil, np.floor])
  210. num_subj = int(take(curr_grp.size * (1 - test_size)))
  211. splinds = rs.choice(curr_grp,
  212. size=num_subj,
  213. replace=False)
  214. split[splinds] = True
  215. # split subjects (with conditions) and stack groups
  216. half = np.hstack([f.flatten('F') for f in
  217. np.split(((inds + 1).astype(bool)
  218. * [split[None]]).T,
  219. splitinds)])
  220. # make sure split half is not a duplicated sequence
  221. dupe_seq = half[:, None] == splitsamp[:, :i]
  222. if dupe_seq.all(axis=0).any():
  223. duplicated = True
  224. if count == 500 and not warned:
  225. warnings.warn('WARNING: Duplicate split halves used.')
  226. warned = True
  227. splitsamp[:, i] = half
  229. return splitsamp
  232. class BasePLS():
  233. """
  234. Base PLS class to be subclassed
  236. Contains most of the math required for PLS, leaving a few functions for PLS
  237. subclasses to implement. This will not run without those implementations.
  239. Parameters
  240. ----------
  241. {input_matrix}
  242. {groups}
  243. {conditions}
  244. **kwargs : optional
  245. Additional key-value pairs; see :obj:`pyls.structures.PLSInputs` for
  246. more info
  248. References
  249. ----------
  251. {references}
  252. """.format(**structures._pls_input_docs)
  254. def __init__(self, X, Y=None, groups=None, n_cond=1, **kwargs):
  255. # if groups aren't provided or are provided wrong, fix them
  256. if groups is None:
  257. groups = [len(X) // n_cond]
  258. elif not isinstance(groups, (list, np.ndarray)):
  259. groups = [groups]
  261. # coerce groups to integers
  262. groups = [int(g) for g in groups]
  264. # check that data matrices and groups + n_cond inputs jibe
  265. n_samples = sum([g * n_cond for g in groups])
  266. if len(X) != n_samples:
  267. raise ValueError('Number of samples specified by `groups` and '
  268. '`n_cond` does not match number of samples in '
  269. 'input array(s).\n'
  270. ' EXPECTED: {}\n'
  271. ' ACTUAL: {} (groups: {} * n_cond: {})'
  272. .format(len(X), n_samples, groups, n_cond))
  274. if Y is not None and len(X) != len(Y):
  275. raise ValueError('Provided `X` and `Y` matrices must have the '
  276. 'same number of samples. Provided matrices '
  277. 'differed: X: {}, Y: {}'.format(len(X), len(Y)))
  279. self.inputs = structures.PLSInputs(X=X, Y=Y, groups=groups,
  280. n_cond=n_cond, **kwargs)
  281. # store dummy-coded array of groups / conditions (save on computation)
  282. self.dummy = utils.dummy_code(groups, n_cond)
  283. self.rs = check_random_state(self.inputs.get('seed'))
  285. # check for parallel processing desire
  286. n_proc = self.inputs.get('n_proc')
  287. if n_proc is not None and n_proc != 1 and not utils.joblib_avail:
  288. self.inputs.n_proc = None
  289. warnings.warn('Setting n_proc > 1 requires the joblib module. '
  290. 'Considering installing joblib and re-running this '
  291. 'if you would like parallelization. Resetting '
  292. 'n_proc to 1 for now.')
  294. def gen_covcorr(self, X, Y, groups=None):
  295. """
  296. Should generate cross-covariance array to be used in `self._svd()`
  298. Must accept the listed parameters and return one array
  300. Parameters
  301. ----------
  302. X : (S, B) array_like
  303. Input data matrix, where `S` is observations and `B` is features
  304. Y : (S, T) array_like
  305. Input data matrix, where `S` is observations and `T` is features
  306. groups : (G,) array_like
  307. Array with number of subjects in each of `G` groups
  309. Returns
  310. -------
  311. crosscov : np.ndarray
  312. Covariance array for decomposition
  313. """
  315. raise NotImplementedError
  317. def gen_distrib(self, X, Y, groups=None, original=None):
  318. """
  319. Should generate behavioral correlations or contrast for bootstrap
  321. Parameters
  322. ----------
  323. X : (S, B) array_like
  324. Input data matrix, where `S` is observations and `B` is features
  325. Y : (S, T) array_like
  326. Input data matrix, where `S` is observations and `T` is features
  327. groups : (S, J) array_like
  328. Dummy coded array, where `S` is observations and `J` corresponds to
  329. the number of different groups x conditions represented in `X` and
  330. `Y`. A value of 1 indicates that an observation belongs to a
  331. specific group or condition
  333. Returns
  334. -------
  335. distrib : (T, L)
  336. Behavioral correlations or contrast for single bootstrap resample
  337. """
  339. raise NotImplementedError
  341. def run_pls(self, X, Y):
  342. """
  343. Runs PLS analysis
  345. Parameters
  346. ----------
  347. X : (S, B) array_like
  348. Input data matrix, where `S` is observations and `B` is features
  349. Y : (S, T) array_like
  350. Input data matrix, where `S` is observations and `T` is features
  352. Returns
  353. -------
  354. results : :obj:`pyls.structures.PLSResults`
  355. Results of PLS (not including PLS type-specific outputs)
  356. """
  358. # initate results structure
  359. self.res = res = structures.PLSResults(inputs=self.inputs)
  361. # get original singular vectors / values
  362. res['x_weights'], res['singvals'], res['y_weights'] = \
  363. self.svd(X, Y, seed=self.rs)
  364. res['x_scores'] = X @ res['x_weights']
  366. if self.inputs.n_perm > 0:
  367. # compute permutations and get statistical significance of LVs
  368. d_perm, ucorrs, vcorrs = self.permutation(X, Y, seed=self.rs)
  369. res['permres']['pvals'] = compute.perm_sig(res['singvals'], d_perm)
  370. res['permres']['permsamples'] = self.permsamp
  372. if self.inputs.n_split is not None:
  373. # get ucorr / vcorr (via split half resampling) for original,
  374. # unpermuted `X` and `Y` arrays
  375. di = np.linalg.inv(res['singvals'])
  376. orig_ucorr, orig_vcorr = self.split_half(X, Y,
  377. res['x_weights'] @ di,
  378. res['y_weights'] @ di,
  379. seed=self.rs)
  380. # get p-values for ucorr/vcorr
  381. ucorr_prob = compute.perm_sig(np.diag(orig_ucorr), ucorrs)
  382. vcorr_prob = compute.perm_sig(np.diag(orig_vcorr), vcorrs)
  384. # get confidence intervals for ucorr/vcorr
  385. ucorr_ll, ucorr_ul = compute.boot_ci(ucorrs, ci=self.inputs.ci)
  386. vcorr_ll, vcorr_ul = compute.boot_ci(vcorrs, ci=self.inputs.ci)
  388. # update results object with split-half resampling results
  389. res['splitres'].update(dict(ucorr=orig_ucorr,
  390. vcorr=orig_vcorr,
  391. ucorr_pvals=ucorr_prob,
  392. vcorr_pvals=vcorr_prob,
  393. ucorr_lolim=ucorr_ll,
  394. vcorr_lolim=vcorr_ll,
  395. ucorr_uplim=ucorr_ul,
  396. vcorr_uplim=vcorr_ul))
  398. return res
  400. def svd(self, X, Y, groups=None, seed=None):
  401. """
  402. Runs SVD on cross-covariance matrix computed from `X` and `Y`
  404. Parameters
  405. ----------
  406. X : (S, B) array_like
  407. Input data matrix, where `S` is observations and `B` is features
  408. Y : (S, T) array_like
  409. Input data matrix, where `S` is observations and `T` is features
  410. groups : (S, J) array_like
  411. Dummy coded array, where `S` is observations and `J` corresponds to
  412. the number of different groups x conditions represented in `X` and
  413. `Y`. A value of 1 indicates that an observation belongs to a
  414. specific group or condition
  415. seed : {int, :obj:`numpy.random.RandomState`, None}, optional
  416. Seed for random number generation. Default: None
  418. Returns
  419. -------
  420. U : (B, L) `numpy.ndarray`
  421. Left singular vectors from singular value decomposition
  422. d : (L, L) `numpy.ndarray`
  423. Diagonal array of singular values from singular value decomposition
  424. V : (J, L) `numpy.ndarray`
  425. Right singular vectors from singular value decomposition
  426. """
  428. # make dummy-coded grouping array if not provided
  429. if groups is None:
  430. groups = utils.dummy_code(self.inputs.groups, self.inputs.n_cond)
  432. # generate cross-covariance matrix and determine # of components
  433. crosscov = self.gen_covcorr(X, Y, groups=groups)
  434. U, d, V = compute.svd(crosscov, seed=seed)
  436. return U, d, V
  438. def bootstrap(self, X, Y, seed=None):
  439. """
  440. Bootstraps `X` and `Y` (w/replacement) and recomputes SVD
  442. Parameters
  443. ----------
  444. X : (S, B) array_like
  445. Input data matrix, where `S` is observations and `B` is features
  446. Y : (S, T) array_like
  447. Input data matrix, where `S` is observations and `T` is features
  448. seed : {int, :obj:`numpy.random.RandomState`, None}, optional
  449. Seed for random number generation. Default: None
  451. Returns
  452. -------
  453. distrib : (T, L) numpy.ndarray
  454. Either behavioral correlations or group x condition contrast;
  455. depends on PLS type
  456. u_sum : (B, L) numpy.ndarray
  457. Sum of the left singular vectors across all bootstraps
  458. u_square : (B, L) numpy.ndarray
  459. Sum of the squared left singular vectors across all bootstraps
  460. """
  462. # generate bootstrap resampled indices (unless already provided)
  463. self.bootsamp = self.inputs.get('bootsamples', None)
  464. if self.bootsamp is None:
  465. self.bootsamp = gen_bootsamp(self.inputs.groups,
  466. self.inputs.n_cond,
  467. self.inputs.n_boot,
  468. seed=seed,
  469. verbose=self.inputs.verbose)
  471. # make empty arrays to store bootstrapped singular vectors
  472. # these will be used to calculate the standard error later on for
  473. # creation of bootstrap ratios
  474. u_sum = np.zeros_like(self.res['x_weights'])
  475. u_square = np.zeros_like(self.res['x_weights'])
  477. # `distrib` corresponds either to the behavioral correlations (if
  478. # running a behavioral PLS) or to the group/condition contrast (if
  479. # running a mean-centered PLS); we'll just extend it and then stack
  480. # all the individual matrices together later (they're quite small so we
  481. # don't need to be too worried about memory usage, here)
  482. distrib = []
  484. # determine the number of bootstraps we'll run each iteration
  485. iters = 1 if self.inputs.n_proc is None else self.inputs.n_proc
  486. gen = utils.trange(self.inputs.n_boot, verbose=self.inputs.verbose,
  487. desc='Running bootstraps')
  489. with utils.get_par_func(self.inputs.n_proc,
  490. self.__class__._single_boot) as (par, func):
  491. boots = 0
  492. while boots < self.inputs.n_boot:
  493. # determine number of bootstraps to run this round
  494. # we don't want to overshoot the requested number, so make
  495. # sure to cut it off if that's what wold happen
  496. top = boots + iters
  497. if top >= self.inputs.n_boot:
  498. top = self.inputs.n_boot
  500. # run the bootstraps
  501. d, usu = zip(*par(func(self, X=X, Y=Y,
  502. inds=self.bootsamp[..., i],
  503. groups=self.dummy,
  504. original=self.res['x_weights'],
  505. seed=i)
  506. for i in range(boots, top)))
  508. # sum bootstrapped singular vectors and store
  509. u_sum += np.sum(usu, axis=0)
  510. u_square += np.sum(np.square(usu), axis=0)
  511. distrib.extend(d)
  513. # force garbage collection
  514. # this is only really needed when parallelizing bootstraps
  515. # the `usu` variable can get REALLY GIANT if either `X` or `Y`
  516. # is large and `n_proc` is > 1, so we really don't want to keep
  517. # it around for any longer than absolutely necessary
  518. if self.inputs.n_proc is not None:
  519. del usu
  520. gc.collect()
  522. # update progress bar and # of bootstraps already run
  523. gen.update(top - boots)
  524. boots = top
  525. gen.close()
  527. return np.stack(distrib, axis=-1), u_sum, u_square
  529. def _single_boot(self, X, Y, inds, groups=None, original=None, seed=None):
  530. """
  531. Bootstraps `X` and `Y` (w/replacement) and recomputes SVD
  533. Parameters
  534. ----------
  535. X : (S, B) array_like
  536. Input data matrix, where `S` is observations and `B` is features
  537. Y : (S, T) array_like
  538. Input data matrix, where `S` is observations and `T` is features
  539. groups : (S, J) array_like
  540. Dummy coded input array, where `S` is observations and `J`
  541. corresponds to the number of different groups x conditions. A value
  542. of 1 indicates that an observation belongs to a specific group or
  543. condition.
  544. original : (B, L) array_like
  545. Left singular vector from original decomposition of `X` and `Y`.
  546. Used to perform Procrustes rotation on permuted singular vectors
  547. seed : {int, :obj:`numpy.random.RandomState`, None}, optional
  548. Seed for random number generation. Default: None
  550. Returns
  551. -------
  552. distrib : np.ndarray
  553. Either behavioral correlations or contrast, depending on PLS type;
  554. generated with self.gen_distrib() which should be specified by the
  555. PLS subclass
  556. U_sum : (B, L) array_like
  557. Left singular vectors from decomposition of bootstrap resampled `X`
  558. and `Y`
  559. """
  561. # make sure we have original (non-bootstrapped) singular vectors
  562. # these are required for the procrustes rotation to ensure our
  563. # singular vectors are all in the same orientation
  564. if original is None:
  565. original = self.svd(X, Y, groups=groups, seed=seed)[0]
  567. # perform SVD of bootstrapped arrays and rotate left singular vectors
  568. U, d = self.svd(X[inds], Y[inds], groups=groups, seed=seed)[:-1]
  569. U_boot = compute.procrustes(original, U, d)
  571. # get contrast / behavcorrs (this function should be specified by the
  572. # subclass)
  573. distrib = self.gen_distrib(X[inds], Y[inds], original, groups)
  575. return distrib, U_boot
  577. def make_permutation(self, X, Y, perminds):
  578. """
  579. Permutes `Y` according to `perminds`, leaving `X` un-permuted
  581. Parameters
  582. ----------
  583. X : (S, B) array_like
  584. Input data matrix, where `S` is observations and `B` is features
  585. Y : (S, T) array_like
  586. Input data matrix, where `S` is observations and `T` is features
  587. perminds : (S,) array_like
  588. Array by which to permute `Y`
  590. Returns
  591. -------
  592. Xp : (S, B) array_like
  593. Identical to `X`
  594. Yp : (S, T) array_like
  595. `Y`, permuted according to `perminds`
  596. """
  598. return X, Y[perminds]
  600. def permutation(self, X, Y, seed=None):
  601. """
  602. Permutes `X` (w/o replacement) and recomputes SVD
  604. Parameters
  605. ----------
  606. X : (S, B) array_like
  607. Input data matrix, where `S` is observations and `B` is features
  608. Y : (S, T) array_like
  609. Input data matrix, where `S` is observations and `T` is features
  610. seed : {int, :obj:`numpy.random.RandomState`, None}, optional
  611. Seed for random number generation. Default: None
  613. Returns
  614. -------
  615. d_perm : (L, P) `numpy.ndarray`
  616. Permuted singular values, where `L` is the number of singular
  617. values and `P` is the number of permutations
  618. ucorrs : (L, P) `numpy.ndarray`
  619. Split-half correlations of left singular values. Only set if
  620. `self.inputs.n_split != 0`
  621. vcorrs : (L, P) `numpy.ndarray`
  622. Split-half correlations of right singular values. Only set if
  623. `self.inputs.n_split != 0`
  624. """
  626. # generate permuted indices (unless already provided)
  627. self.permsamp = self.inputs.get('permsamples')
  628. if self.permsamp is None:
  629. self.permsamp = gen_permsamp(self.inputs.groups,
  630. self.inputs.n_cond,
  631. self.inputs.n_perm,
  632. seed=seed,
  633. verbose=self.inputs.verbose)
  635. # get permuted values (parallelizing as requested)
  636. gen = utils.trange(self.inputs.n_perm, verbose=self.inputs.verbose,
  637. desc='Running permutations')
  638. with utils.get_par_func(self.inputs.n_proc,
  639. self.__class__._single_perm) as (par, func):
  640. out = par(func(self, X=X, Y=Y, inds=self.permsamp[:, i],
  641. groups=self.dummy, original=self.res['y_weights'],
  642. seed=i)
  643. for i in gen)
  644. d_perm, ucorrs, vcorrs = [np.stack(o, axis=-1) for o in zip(*out)]
  646. return d_perm, ucorrs, vcorrs
  648. def _single_perm(self, X, Y, inds, groups=None, original=None, seed=None):
  649. """
  650. Permutes `X` (w/o replacement) and recomputes SVD
  652. Parameters
  653. ----------
  654. X : (S, B) array_like
  655. Input data matrix, where `S` is observations and `B` is features
  656. Y : (S, T) array_like
  657. Input data matrix, where `S` is observations and `T` is features
  658. inds : (S,) array_like
  659. Permutation resampling array
  660. original : (J, L) array_like
  661. Right singular vector from original decomposition of `X` and `Y`.
  662. Used to perform Procrustes rotation on permuted singular values,
  663. if desired
  664. seed : {int, :obj:`numpy.random.RandomState`, None}, optional
  665. Seed for random number generation. Default: None
  667. Returns
  668. -------
  669. d_perm : (L,) `numpy.ndarray`
  670. Permuted singular values, where `L` is the number of singular
  671. values
  672. ucorrs : (L,) `numpy.ndarray`
  673. Split-half correlations of left singular values. Only set if
  674. `self.inputs.n_split != 0`
  675. vcorrs : (L,) `numpy.ndarray`
  676. Split-half correlations of right singular values. Only set if
  677. `self.inputs.n_split != 0`
  678. """
  680. # calculate SVD of permuted matrices
  681. Xp, Yp = self.make_permutation(X, Y, inds)
  682. U, d, V = self.svd(Xp, Yp, groups=groups, seed=seed)
  684. # optionally get rotated/rescaled singular values
  685. if self.inputs.rotate:
  686. if original is None:
  687. original = self.svd(X, Y, groups=groups, seed=seed)[-1]
  688. ssd = np.sqrt(np.sum(compute.procrustes(original, V, d)**2,
  689. axis=0))
  690. else:
  691. ssd = np.diag(d)
  693. # get ucorr/vcorr if split-half resampling requested
  694. if self.inputs.n_split is not None:
  695. di = np.linalg.inv(d)
  696. ucorr, vcorr = self.split_half(Xp, Yp, U @ di, V @ di,
  697. groups=groups, seed=seed)
  698. else:
  699. ucorr, vcorr = None, None
  701. return ssd, ucorr, vcorr
  703. def split_half(self, X, Y, ud=None, vd=None, groups=None, seed=None):
  704. """
  705. Parameters
  706. ----------
  707. X : (S, B) array_like
  708. Input data matrix, where `S` is observations and `B` is features
  709. Y : (S, T) array_like
  710. Input data matrix, where `S` is observations and `T` is features
  711. ud : (B, L) array_like
  712. Left singular vectors, scaled by singular values
  713. vd : (J, L) array_like
  714. Right singular vectors, scaled by singular values
  715. seed : {int, :obj:`numpy.random.RandomState`, None}, optional
  716. Seed for random number generation. Default: None
  718. Returns
  719. -------
  720. ucorr : (L,) `numpy.ndarray`
  721. Average correlation of left singular vectors across split-halves
  722. vcorr : (L,) `numpy.ndarray`
  723. Average correlation of right singular vectors across split-halves
  724. """
  726. # generate splits
  727. splitsamp = gen_splits(self.inputs.groups,
  728. self.inputs.n_cond,
  729. self.inputs.n_split,
  730. seed=seed,
  731. test_size=0.5).astype(bool)
  733. # make dummy-coded grouping array if not provided
  734. if groups is None:
  735. groups = utils.dummy_code(self.inputs.groups, self.inputs.n_cond)
  737. # generate original singular vectors if not provided
  738. if ud is None or vd is None:
  739. U, d, V = self.svd(X, Y, groups=groups, seed=seed)
  740. di = np.linalg.inv(d)
  741. ud, vd = U @ di, V @ di
  743. # empty arrays to hold split-half correlations
  744. ucorr = np.zeros(shape=(ud.shape[-1], self.inputs.n_split))
  745. vcorr = np.zeros(shape=(vd.shape[-1], self.inputs.n_split))
  747. for i in range(self.inputs.n_split):
  748. # calculate cross-covariance matrix for both splits
  749. spl = splitsamp[:, i]
  750. D1 = self.gen_covcorr(X[spl], Y[spl], groups=groups[spl])
  751. D2 = self.gen_covcorr(X[~spl], Y[~spl], groups=groups[~spl])
  753. # project cross-covariance matrices onto original SVD to obtain
  754. # left & right singular vector and correlate between split halves
  755. ucorr[:, i] = compute.efficient_corr(D1.T @ vd, D2.T @ vd)
  756. vcorr[:, i] = compute.efficient_corr(D1 @ ud, D2 @ ud)
  758. # return average correlations for singular vectors across `n_split`
  759. return np.mean(ucorr, axis=-1), np.mean(vcorr, axis=-1)

base.py at commit d8a19d5, under GPL-2.0 · at the source

Overview

Authors: Marvin Petersen1, Moritz A. Link1, Carola Mayer1, Felix L. Nägele1, Maximilian Schell1, Märit Jensen1, Eckhard Schlemm1, Jens Fiehler2, Jürgen Gallinat3, Simone Kühn3, Raphael Twerenbold4,5,6,7, Amir Omidvarnia8,9, Felix Hoffstaedter8,9, Kaustubh R. Patil8,9, Simon B. Eickhoff8,9, Götz Thomalla1, Bastian Cheng1
  1. Department of Neurology, University Medical Center Hamburg-Eppendorf, Hamburg, Germany
  2. Department of Diagnostic and Interventional Neuroradiology, University Medical Center Hamburg-Eppendorf, Hamburg, Germany
  3. Department of Psychiatry and Psychotherapy, University Medical Center Hamburg-Eppendorf, Hamburg, Germany
  4. Department of General and Interventional Cardiology, University Heart and Vascular Center, Hamburg, Germany
  5. Epidemiological Study Center, University Medical Center Hamburg-Eppendorf, Hamburg, Germany
  6. German Center for Cardiovascular Research (DZHK), Partner Site Hamburg/Kiel/Luebeck, Hamburg, Germany
  7. University Center of Cardiovascular Science, University Heart and Vascular Center, Hamburg, Germany
  8. Faculty of Medicine, Institute of Systems Neuroscience, Heinrich Heine University Düsseldorf, Düsseldorf, Germany
  9. Institute of Neuroscience and Medicine, Brain and Behaviour (INM-7), Research Center Jülich, Jülich, Germany
Institutions: University Medical Center Hamburg-Eppendorf (Germany); German Centre for Cardiovascular Research (Germany); Forschungszentrum Jülich (Germany); Heinrich Heine University Düsseldorf (Germany)
Journal: Frontiers in aging neuroscience, volume 18, article 1789408
Dates: received 16 January 2026; accepted 10 April 2026; published online 1 May 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.3389/fnagi.2026.1789408 · PMID 42147461 · PMCID PMC13176138 · OpenAlex W7159953708
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: structural MRI / diffusion (modality), human (organism), stroke (population), cognitive (subfield)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, fMRI & imaging
Keywords: brain age, cardiovascular risk, cerebral small vessel disease, cognitive function, motor function, neuroimaging, white matter hyperintensities
Topic: Dementia and Cognitive Impairment Research (Psychiatry and Mental health, Medicine), according to OpenAlex
Citations: not cited yet (Europe PMC); 54 references in the paper

Abstract

Introduction: Age-related declines in cognitive and motor functions show highly variable trajectories. To better understand the underlying mechanisms, we investigated multivariate associative effects between modifiable vascular risk factors, biological brain aging, cognitive, and motor performance in 40,579 individuals from the population-based UK Biobank and Hamburg City Health Study.

Methods: We employed partial least squares correlation analysis (PLS) to model associations between multi-domain cognitive and motor test scores and three distinct MRI-derived markers of biological brain aging: relative brain age (from morphometric brain imaging), white matter hyperintensity load, and peak width of skeletonized mean diffusivity. Furthermore, we conducted mediation analyses to assess if these markers mediate the impact of vascular risk on functional decline.

Results: PLS identified a single dominant latent dimension explaining 94.7% of the shared variance between neuroimaging and behavior. This dimension linked higher biological brain aging markers – with relative brain age showing the strongest contribution – to poorer cognitive and motor performance, particularly in executive function and processing speed. Mediation analysis revealed that biological brain aging acts as a partial mediator for the negative effects of blood pressure, glucose, waist-hip ratio, and smoking load on cognitive and motor function. Notably, this mediating effect was not observed for cholesterol levels. These results were consistent across both cohorts.

Discussion: Our study illustrates the associative interplay between vascular health, biological brain aging, and cognitive and motor performance, emphasizing the need for preventive strategies to maintain late-life independence in aging populations.

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 3 matches between paragraphs and lines of code.

rmarkello/pyls

License: GPL-2.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: d8a19d564cc5804249527b68c937df3a5fd8c7cc, 4 November 2019
Languages: Python (32)
Size: 65 files, 32 scripts
Software Heritage: archived
Found in: the end of the paper
Holds: README, license file, environment (requirements.txt, setup.cfg, setup.py, docs/requirements.txt), tests, continuous integration, documentation
Not found: CITATION.cff
Tools: NumPy (18 files), scikit-learn (4 files), h5py (2 files), pandas (2 files), SciPy (2 files), seaborn (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
34 files
At the source: github.com/rmarkello/pyls

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;
  • 32 scripts, each with its path and the digest of its content;
  • 3 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

Data availability statement

UK Biobank data can be obtained via its standardized data access procedure (https://www.ukbiobank.ac.uk/). HCHS participant data used in this analysis is not publicly available for privacy reasons, but access can be established via request to the HCHS steering committee.

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 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 17 authors, 7 keywords, 54 references.

Cite

This paper

Petersen, M., Link, M. A., Mayer, C., Nägele, F. L., Schell, M., Jensen, M., Schlemm, E., Fiehler, J., Gallinat, J., Kühn, S., Twerenbold, R., Omidvarnia, A., Hoffstaedter, F., Patil, K. R., Eickhoff, S. B., Thomalla, G., & Cheng, B. (2026). Biological brain aging, cognitive-motor decline and vascular risk: a multivariate imaging analysis of 40,579 individuals. Frontiers in aging neuroscience, 18, 1789408. https://doi.org/10.3389/fnagi.2026.1789408

BibTeX

@article{petersen2026biological,
author = {Petersen, Marvin and Link, Moritz A. and Mayer, Carola and Nägele, Felix L. and Schell, Maximilian and Jensen, Märit and Schlemm, Eckhard and Fiehler, Jens and Gallinat, Jürgen and Kühn, Simone and Twerenbold, Raphael and Omidvarnia, Amir and Hoffstaedter, Felix and Patil, Kaustubh R. and Eickhoff, Simon B. and Thomalla, Götz and Cheng, Bastian},
title = {{Biological brain aging, cognitive-motor decline and vascular risk: a multivariate imaging analysis of 40,579 individuals}},
journal = {Frontiers in aging neuroscience},
year = {2026},
month = may,
volume = {18},
pages = {1789408},
publisher = {Frontiers Media SA},
issn = {1663-4365},
doi = {10.3389/fnagi.2026.1789408},
url = {https://doi.org/10.3389/fnagi.2026.1789408},
pmid = {42147461},
pmcid = {PMC13176138}
}

RIS

TY - JOUR
AU - Petersen, Marvin
AU - Link, Moritz A.
AU - Mayer, Carola
AU - Nägele, Felix L.
AU - Schell, Maximilian
AU - Jensen, Märit
AU - Schlemm, Eckhard
AU - Fiehler, Jens
AU - Gallinat, Jürgen
AU - Kühn, Simone
AU - Twerenbold, Raphael
AU - Omidvarnia, Amir
AU - Hoffstaedter, Felix
AU - Patil, Kaustubh R.
AU - Eickhoff, Simon B.
AU - Thomalla, Götz
AU - Cheng, Bastian
TI - Biological brain aging, cognitive-motor decline and vascular risk: a multivariate imaging analysis of 40,579 individuals
T2 - Frontiers in aging neuroscience
J2 - Front Aging Neurosci
PY - 2026
DA - 2026/05/01
VL - 18
SP - 1789408
SN - 1663-4365
PB - Frontiers Media SA
DO - 10.3389/fnagi.2026.1789408
UR - https://doi.org/10.3389/fnagi.2026.1789408
LA - en
ER -

CSL-JSON

{
"id": "10.3389/fnagi.2026.1789408",
"type": "article-journal",
"title": "Biological brain aging, cognitive-motor decline and vascular risk: a multivariate imaging analysis of 40,579 individuals",
"container-title": "Frontiers in aging neuroscience",
"author": [
{
"family": "Petersen",
"given": "Marvin"
},
{
"family": "Link",
"given": "Moritz A."
},
{
"family": "Mayer",
"given": "Carola"
},
{
"family": "Nägele",
"given": "Felix L."
},
{
"family": "Schell",
"given": "Maximilian"
},
{
"family": "Jensen",
"given": "Märit"
},
{
"family": "Schlemm",
"given": "Eckhard"
},
{
"family": "Fiehler",
"given": "Jens"
},
{
"family": "Gallinat",
"given": "Jürgen"
},
{
"family": "Kühn",
"given": "Simone"
},
{
"family": "Twerenbold",
"given": "Raphael"
},
{
"family": "Omidvarnia",
"given": "Amir"
},
{
"family": "Hoffstaedter",
"given": "Felix"
},
{
"family": "Patil",
"given": "Kaustubh R."
},
{
"family": "Eickhoff",
"given": "Simon B."
},
{
"family": "Thomalla",
"given": "Götz"
},
{
"family": "Cheng",
"given": "Bastian"
}
],
"container-title-short": "Front Aging Neurosci",
"volume": "18",
"page": "1789408",
"DOI": "10.3389/fnagi.2026.1789408",
"PMID": "42147461",
"PMCID": "PMC13176138",
"ISSN": "1663-4365",
"publisher": "Frontiers Media SA",
"URL": "https://doi.org/10.3389/fnagi.2026.1789408",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
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.1371/journal.pbio.3003856 [code]
Aging and metabolism contribute separately to brain-body health.
Journal: PLoS biology
In common: seaborn, scikit-learn, pandas, 2 other tools, ukbiobank.ac.uk/media/0xsbmfmw, structural MRI / diffusion, 6 references
[2] doi:10.1038/s41467-026-71271-9 [code]
Exposome-wide patterns predict brain health in aging.
Journal: Nature communications
In common: seaborn, scikit-learn, pandas, 2 other tools, 5 references
[3] doi:10.7554/elife.108109 [code]
Multimodal MRI marker of cognition explains the association between cognition and mental health in the UK Biobank.
Journal: eLife
In common: seaborn, scikit-learn, pandas, 2 other tools, cognitive, structural MRI / diffusion, 5 references
[4] doi:10.1016/j.ebiom.2026.106312 [code]
Subgingival microbiota composition is associated with brain health in the general population-the PAROMIND study.
Journal: EBioMedicine
In common: scikit-learn, pandas, SciPy, 1 other tool, cognitive, 4 references
[5] 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: h5py, seaborn, scikit-learn, 3 other tools, structural MRI / diffusion, 4 references
[6] doi:10.1038/s41514-026-00456-9 [code]
Exploring the link between body physiology and cognition: the role of the brain and aging.
Journal: npj aging
In common: seaborn, scikit-learn, pandas, 2 other tools, cognitive, 4 references
[7] doi:10.1186/s40708-026-00316-y [code]
Generalizable and explainable deep learning for brain MRI: a multi-cohort evaluation of 3D architectures for age and sex prediction.
Journal: Brain informatics
In common: seaborn, scikit-learn, pandas, 2 other tools, structural MRI / diffusion, 3 references
[8] doi:10.1162/imag.a.1242 [code]
Stable individual differences dominate adult brain volume variation until later life.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: seaborn, pandas, SciPy, 1 other tool, structural MRI / diffusion, 4 references
[9] doi:10.1038/s41593-026-02359-0 [code]
The cross-site reproducibility of MRI morphometric phenotypes in psychiatric disorders.
Journal: Nature neuroscience
In common: scikit-learn, pandas, SciPy, 1 other tool, structural MRI / diffusion, 5 references
[10] doi:10.1016/j.isci.2026.117180 [code]
Developmental changes in similarity between neural representations of mental arithmetic and artificial neural networks.
Journal: iScience
In common: h5py, seaborn, scikit-learn, 3 other tools, 3 references

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.