OSCR

Cortical thickness changes precede high levels of amyloid by at least 7 years.

Code ↔ Paper

6 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 6 matches
  1. [1] § Methods › Participants ↔ scripts/01-prepSlopesYearsBeforeAB_simulated.r, lines 354–440 · score 0.70 · combined MRI, ageing MRI, MRI scans, PET scan, ADNI, BACS
  2. [2] § Results › Additional analyses ↔ scripts/10-rankThicknessChange.r, lines 672–712 · score 0.62 · CIs excluded zero, thickness changes, confidence, rank, derivative, interval
  3. [3] § Methods › Statistical analysis ↔ scripts/10-rankThicknessChange.r, lines 672–712 · score 0.57 · excluded zero, thickness changes, rank, gratia, derivative, CI
  4. [4] § Methods › Statistical analysis ↔ scripts/08-rankAmyloidOrder.r, lines 363–409 · score 0.53 · linear mixed models, tracer, regional, amyloid, SUVR, ADNI
  5. [5] § Methods › Statistical analysis ↔ scripts/10-rankThicknessChange.r, lines 400–463 · score 0.51 · GAMM interaction, thickness change, ICV, smooth, strength, trajectories
  6. [6] § Methods › Statistical analysis ↔ scripts/08-rankAmyloidOrder.r, lines 363–409 · score 0.51 · linear mixed models, tracer, rank, regional, SUVRs, intercept

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

R · 1,114 lines · 39 KB · MIT · 3 matches

  1. #========================================================================================#
  2. # Author: James M Roe, Ph.D.
  3. # Center for Lifespan Changes in Brain and Cognition, University of Oslo
  4. #
  5. # Purpose: Compute the rank order of thickness changes across all converter MRIs and correlate with the rank order of regional Aβ deposition.
  6. # Uses linear time-to-Aβ trajectories or SILA input, and computes time-to-Aβ thickness trajectories in regions of high-v-low Aβ.
  7. # Script requires individual-level data as input and is not executable
  8. #========================================================================================#
  11. #---load packages
  12. loadPackages = function() {
  13. packages = c("here", "tidyverse","magrittr","gamm4","itsadug","numDeriv","gratia","mgcv","viridis","wesanderson","asbio","broom","cowplot","data.table","stringi","tictoc","lmerTest","effects","ggpubr")
  14. new.packages = packages[!(packages %in% installed.packages()[,"Package"])]
  15. if(length(new.packages)) {
  16. install.packages(new.packages)
  17. }
  18. print(sapply(packages, require, character.only = T))
  19. print(sapply(packages, function(p) as.character(packageVersion(p))))
  20. }
  21. loadPackages()
  22. here()
  26. #---set dir
  27. b = "/cluster/projects/p274/projects/p040-ad_change/Berkeley"
  28. # b = here()
  29. setwd(b)
  30. #---set dir
  33. #---make dirstruct
  34. plotdir = "plots"; if (! dir.exists(plotdir)) { dir.create(plotdir)}
  35. resdir = "results"; if (! dir.exists(resdir)) { dir.create(resdir)}
  38. #---load data
  39. savefigs=F
  40. nTime=2
  41. agecut=30
  42. saveres=T
  45. load(file.path(b, "reproduce/data/DF_LONG_4570.Rda"))
  46. dim(DF); length(unique(DF$subject_id))
  50. # load converter/nonconverter data
  51. load(file.path(b, "reproduce/data/converters_all_negfirst_ADNINC_UPDATE.Rda"))
  52. load(file.path(b, "reproduce/data/converters_all_negfirst_BACS_UPDATE.Rda"))
  53. load(file.path(b, "reproduce/data/converters_all_negfirst_LCBC_UPDATE_REPRO_corthresh.Rda"))
  54. load(file.path(b, "reproduce/data/converters_data_for_plot_BACS_UPDATE_REPRO.Rda"))
  55. load(file.path(b, "reproduce/data/converters_data_for_plot_ADNINC_UPDATE_REPRO.Rda"))
  56. load(file.path(b, "reproduce/data/converters_data_for_plot_LCBC_UPDATE_REPRO_corthresh.Rda"))
  59. converters_allMRI_negfirst_LCBC = left_join(converters_allMRI_negfirst_LCBC,
  60. converters_allMRI_negfirst_LCBC_plotdat %>% select(imageLink, slope, intercept, age_at_threshold, contains("CL_at_thresh")) %>%
  61. rename(slope_centiloid = slope,
  62. intercept_centiloid = intercept)
  63. )
  65. converters_allMRI_negfirst_ADNINC = left_join(converters_allMRI_negfirst_ADNINC,
  66. converters_allMRI_negfirst_ADNINC_plotdat %>% select(imageLink, slope, intercept, age_at_threshold, contains("CL_at_thresh")) %>%
  67. rename(slope_centiloid = slope,
  68. intercept_centiloid = intercept)
  69. )
  71. converters_allMRI_negfirst_BACS = left_join(converters_allMRI_negfirst_BACS,
  72. converters_allMRI_negfirst_BACS_plotdat %>% select(imageLink, slope, intercept, age_at_threshold, contains("CL_at_thresh")) %>%
  73. rename(slope_centiloid = slope,
  74. intercept_centiloid = intercept)
  75. )
  79. allFeat = readLines(file.path(b, "reproduce/data/allFeatures364.txt"))
  80. adnioutlier = "029_S_0845"
  81. rois=allFeat[grepl("thickness", allFeat)]
  82. nrois=length(rois)
  83. subset.size=nrois; jj = 1
  84. N = ceiling(nrois/subset.size)
  85. print(N)
  86. start = (jj*subset.size)-subset.size+1
  88. if (jj == N) {
  89. end = nrois
  90. loopend = length(rois[start:end])
  91. } else {
  92. end = jj*subset.size
  93. loopend = subset.size
  94. }
  95. print(paste("subsetting cols", start, "-", end))
  96. rois = rois[start:end]
  97. print(rois)
  98. ROIs = rois
  101. pb = txtProgressBar(min=2, max=end, style=3)
  102. Usubs = length(unique(DF$subject_id))
  105. # load sila outputs
  106. osila_bacs = fread("/cluster/projects/p274/projects/p040-ad_change/Berkeley/scripts/SILA-AD-Biomarker/demo/output/testBACS.csv")
  107. osila_adni = fread("/cluster/projects/p274/projects/p040-ad_change/Berkeley/scripts/SILA-AD-Biomarker/demo/output/testADNINC.csv")
  108. osila_lcbc = fread("/cluster/projects/p274/projects/p040-ad_change/Berkeley/scripts/SILA-AD-Biomarker/demo/output2/testLCBC.csv")
  111. # load sila inputs (ids get changed in sila modelling)
  112. isila_all = fread("/cluster/projects/p274/projects/p040-ad_change/Berkeley/scripts/SILA-AD-Biomarker/demo/df_silo_amyloidtimeCorrect.csv")
  113. isila_adni = isila_all[isila_all$cohort == "ADNINC",] %>% rename(subject_id = subid) %>% rename(subid = subjid)
  114. isila_bacs = isila_all[isila_all$cohort == "BACS",] %>% rename(subject_id = subid) %>% rename(subid = subjid)
  115. isila_lcbc = isila_all[isila_all$cohort == "LCBC",] %>% rename(subject_id = subid) %>% rename(subid = subjid)
  118. range(osila_adni$subid)
  119. range(isila_adni$subid)
  120. range(osila_bacs$subid)
  121. range(isila_bacs$subid)
  122. range(isila_lcbc$subid)
  123. range(osila_lcbc$subid)
  127. osila_adni = left_join(osila_adni, isila_adni)
  128. osila_bacs = left_join(osila_bacs, isila_bacs)
  129. osila_lcbc = left_join(osila_lcbc, isila_lcbc)
  133. mytheme = theme(
  134. plot.background = element_rect(fill = "white"),
  135. panel.background = element_rect(fill = "white"),
  136. panel.grid.major = element_blank(),
  137. panel.grid.minor = element_blank(),
  138. title = element_text(size=17),
  139. text = element_text(color = "black", size = 18, family="Nimbus Sans Narrow"),
  140. plot.title = element_text(hjust = 0.5),
  141. # axis.ticks = element_blank(),
  142. axis.title.y = element_text(color = "black", size = 22, vjust =-1, margin = margin(0,20,0,0)),
  143. axis.title.x = element_text(color = "black", size = 22, vjust = -2, margin = margin(0,20,20,0)),
  144. axis.text = element_text(color = "black", size = 18),
  145. legend.key.size = unit(1,"cm"))
  148. pal = wes_palette("FantasticFox1", n = 5)
  149. plotSila = function(dat, cohort) {
  151. # dat = osila_adni
  152. (p_sila1 =
  153. dat %>%
  154. ggplot(.) +
  155. geom_line(data=dat,aes(x=age,val,group=subid, col = factor(conv)),alpha=0.6, size=0.5) +
  156. geom_point(data=dat,aes(x=age,val,group=subid, col = factor(conv)),stat="identity",alpha=1, size=0.5) +
  157. scale_color_manual(values = c(pal[2], pal[5])) +
  158. # geom_smooth(method = "gam", col = "black", se = F) +
  159. ggtitle(cohort) +
  160. labs(x = "Age") +
  161. theme_classic() + mytheme)
  163. (p_sila2 =
  164. dat %>%
  165. ggplot(.) +
  166. geom_line(data=dat,aes(x=estdtt0,val,group=subid, col = factor(conv)),alpha=0.6, size=0.5) +
  167. geom_point(data=dat,aes(x=estdtt0,val,group=subid, col = factor(conv)),stat="identity",alpha=1, size=0.5) +
  168. scale_color_manual(values = c(pal[2], pal[5])) +
  169. geom_hline(yintercept = dat$valt0, linetype = 2, col = "black") +
  170. ggtitle(cohort) +
  171. labs(x = "Years to Aβ+ (SILA)") +
  172. theme_classic() + mytheme)
  174. return(list(p_sila1 = p_sila1, p_sila2 = p_sila2, threshold = dat$valt0[1]))
  176. }
  179. p_sila_adni = plotSila(osila_adni, "ADNI")
  180. p_sila_bacs = plotSila(osila_bacs, "BACS")
  181. p_sila_lcbc = plotSila(osila_lcbc, "LCBC")
  182. p_sila_adni$p_sila2
  183. p_sila_bacs$p_sila2
  184. p_sila_lcbc$p_sila2
  187. if (savefigs) {
  188. ggsave(filename = "/cluster/projects/p274/projects/p040-ad_change/Berkeley/paper2/figs_yearsBeforeAB/p_sila_adni.pdf",
  189. plot = p_sila_adni$p_sila2 + theme(legend.position = "none"),
  190. width = 13,
  191. height = 13,
  192. dpi = 600,
  193. units = "cm",
  194. device = cairo_pdf
  195. )
  196. ggsave(filename = "/cluster/projects/p274/projects/p040-ad_change/Berkeley/paper2/figs_yearsBeforeAB/p_sila_bacs.pdf",
  197. plot = p_sila_bacs$p_sila2 + theme(legend.position = "none"),
  198. width = 13,
  199. height = 13,
  200. dpi = 600,
  201. units = "cm",
  202. device = cairo_pdf
  203. )
  204. ggsave(filename = "/cluster/projects/p274/projects/p040-ad_change/Berkeley/paper2/figs_yearsBeforeAB/p_sila_lcbc_threshcorrect.pdf",
  205. plot = p_sila_lcbc$p_sila2 + theme(legend.position = "none"),
  206. width = 13,
  207. height = 13,
  208. dpi = 600,
  209. units = "cm",
  210. device = cairo_pdf
  211. )
  212. }
  215. DF_SILA = rbind(osila_adni,
  216. osila_bacs,
  217. osila_lcbc)
  218. DF_SILA$visit_age = DF_SILA$age
  219. DF_SILA$subject_id %in% DF$subject_id
  220. DF_SILA$subject_id[!DF_SILA$subject_id %in% DF$subject_id]
  221. DF_SILA$age_at_sila_threshold = DF_SILA$age - DF_SILA$estdtt0
  224. DF_SILA %>%
  225. filter(subject_id == "002_S_4213") %>%
  226. pull(age_at_sila_threshold) %>%
  227. dput()
  229. # fix four subjects that have very slightly different age at threshold across their observations
  230. (checkSubs = DF_SILA %>%
  231. group_by(subject_id) %>%
  232. summarise(n_unique = n_distinct(round(age_at_sila_threshold, 4))) %>%
  233. filter(n_unique != 1)
  234. )
  235. DF_SILA[DF_SILA$subject_id == "021_S_4276",]
  236. DF_SILA[DF_SILA$subject_id == "031_S_4021",]
  237. DF_SILA[DF_SILA$subject_id == "036_S_4491",]
  238. DF_SILA[DF_SILA$subject_id == "1100591",]
  240. # by taking their last estimate
  241. DF_SILA$age_at_sila_threshold[DF_SILA$subject_id == "021_S_4276"] = DF_SILA$age_at_sila_threshold[DF_SILA$subject_id == "021_S_4276"][4]
  242. DF_SILA$age_at_sila_threshold[DF_SILA$subject_id == "031_S_4021"] = DF_SILA$age_at_sila_threshold[DF_SILA$subject_id == "031_S_4021"][3]
  243. DF_SILA$age_at_sila_threshold[DF_SILA$subject_id == "036_S_4491"] = DF_SILA$age_at_sila_threshold[DF_SILA$subject_id == "036_S_4491"][4]
  244. DF_SILA$age_at_sila_threshold[DF_SILA$subject_id == "1100591"] = DF_SILA$age_at_sila_threshold[DF_SILA$subject_id == "1100591"][2]
  246. DF_SILA$age_at_sila_threshold = round(DF_SILA$age_at_sila_threshold, 3)
  247. DF_SILA %>%
  248. group_by(subject_id) %>%
  249. summarise(n_unique = n_distinct(age_at_sila_threshold)) %>%
  250. filter(n_unique != 1)
  253. DF = left_join(DF,
  254. DF_SILA %>% select(subject_id, age_at_sila_threshold, conv) %>% distinct()
  255. )
  257. DF$visit_age
  258. DF$time_from_sila_threshold = DF$age_at_sila_threshold - DF$visit_age
  259. DF$time_from_sila_threshold_flip = DF$time_from_sila_threshold*-1
  260. head(DF %>% select(matches("age|time", ignore.case = TRUE)))
  264. # NB! this is correctly 76 (not 77) due to the negative centiloid slope BACS person -------
  265. # and this person had 10 MRI scans
  266. # hence the difference between 477 scans (in this analysis) and max 487 MRI scans in converters (in paper) is correct
  267. load("/cluster/projects/p274/projects/p040-ad_change/Berkeley/reproduce/data/DF.convallMRI.rda")
  268. length(unique(DF.convallMRI$subject_id)); dim(DF.convallMRI)
  269. DF.convallMRI$diff_mriAge_predABpos
  272. # choose to estimate via original method (linear estimates)
  273. # or SILA
  274. estType = "SILA"
  275. convonly = 1
  276. estType = "ORIG"
  278. if (estType != "SILA") {
  280. # if not SILA analysis (original with converters only)
  281. DF = DF.convallMRI
  282. dim(DF)
  283. length(unique(DF$subject_id))
  284. converters_allMRI_negfirst_BACS[converters_allMRI_negfirst_BACS$SID == "B16-220",] %>% dim()
  285. convonly = 0
  287. } else if (estType == "SILA") {
  289. # if SILA analysis (review)
  290. DF = DF %>% filter(!is.na(time_from_sila_threshold))
  291. dim(DF)
  292. length(unique(DF$subject_id))
  293. length(unique(DF$subject_id[DF$conv == 1]))
  295. # if testing SILA across only converter group
  296. if (convonly) {
  297. dim(DF %>% filter(conv == 1))
  298. DF %<>% filter(conv == 1)
  299. }
  300. }
  305. # high ab region map
  306. frontal_regions_lh <- c(
  307. "lh_superiorfrontal_thickness.aparcnative71",
  308. "lh_rostralmiddlefrontal_thickness.aparcnative71",
  309. "lh_caudalmiddlefrontal_thickness.aparcnative71",
  310. "lh_parsopercularis_thickness.aparcnative71",
  311. "lh_parstriangularis_thickness.aparcnative71",
  312. "lh_parsorbitalis_thickness.aparcnative71",
  313. "lh_lateralorbitofrontal_thickness.aparcnative71",
  314. "lh_medialorbitofrontal_thickness.aparcnative71",
  315. "lh_frontalpole_thickness.aparcnative71"
  316. )
  317. frontal_regions_rh <- c(
  318. "rh_superiorfrontal_thickness.aparcnative71",
  319. "rh_rostralmiddlefrontal_thickness.aparcnative71",
  320. "rh_caudalmiddlefrontal_thickness.aparcnative71",
  321. "rh_parsopercularis_thickness.aparcnative71",
  322. "rh_parstriangularis_thickness.aparcnative71",
  323. "rh_parsorbitalis_thickness.aparcnative71",
  324. "rh_lateralorbitofrontal_thickness.aparcnative71",
  325. "rh_medialorbitofrontal_thickness.aparcnative71",
  326. "rh_frontalpole_thickness.aparcnative71"
  327. )
  328. parietalregions = c(
  329. "lh_precuneus_thickness.aparcnative71",
  330. "rh_precuneus_thickness.aparcnative71",
  331. "lh_inferiorparietal_thickness.aparcnative71",
  332. "rh_inferiorparietal_thickness.aparcnative71",
  333. "lh_supramarginal_thickness.aparcnative71",
  334. "rh_supramarginal_thickness.aparcnative71"
  335. )
  336. cingulate_regions = c(
  337. "lh_posteriorcingulate_thickness.aparcnative71",
  338. "rh_posteriorcingulate_thickness.aparcnative71",
  339. "lh_isthmuscingulate_thickness.aparcnative71",
  340. "rh_isthmuscingulate_thickness.aparcnative71",
  341. "lh_rostralanteriorcingulate_thickness.aparcnative71",
  342. "rh_rostralanteriorcingulate_thickness.aparcnative71",
  343. "lh_caudalanteriorcingulate_thickness.aparcnative71",
  344. "rh_caudalanteriorcingulate_thickness.aparcnative71"
  345. )
  346. temporal_regions = c(
  347. "lh_middletemporal_thickness.aparcnative71",
  348. "rh_middletemporal_thickness.aparcnative71"
  349. )
  352. ab_regions_high = c(frontal_regions_lh, frontal_regions_rh, parietalregions, temporal_regions, cingulate_regions)
  355. # low ab regions - everything else
  356. ROIs[!ROIs %in% ab_regions_high]
  357. ab_regions_low = ROIs[!ROIs %in% ab_regions_high]
  358. ab_regions_low = ab_regions_low[!grepl("Mean", ab_regions_low)]
  361. # FDR regions in thickness analysis + homologues
  362. FDR_mask = c(frontal_regions_lh, frontal_regions_rh,
  363. "lh_precentral_thickness.aparcnative71",
  364. "rh_precentral_thickness.aparcnative71",
  365. "lh_paracentral_thickness.aparcnative71",
  366. "rh_paracentral_thickness.aparcnative71",
  367. "lh_insula_thickness.aparcnative71",
  368. "rh_insula_thickness.aparcnative71",
  369. "lh_supramarginal_thickness.aparcnative71",
  370. "rh_supramarginal_thickness.aparcnative71"
  371. )
  374. # minus frontal poles which were not FDR sig
  375. FDR_mask = FDR_mask[!grepl("pole", FDR_mask)]
  376. ROIs = FDR_mask
  379. # reverse years to AB to be correct (-years to AB)
  380. if (estType != "SILA") {
  381. DF$diff_mriAge_predABpos_flip = DF$diff_mriAge_predABpos*-1
  382. } else if (estType == "SILA") {
  383. DF$diff_mriAge_predABpos_flip = DF$time_from_sila_threshold_flip
  384. }
  385. DF$brainvarLow = rowMeans(DF[,ab_regions_low])
  386. DF$brainvarHigh = rowMeans(DF[,ab_regions_high])
  387. DF$subject_id = as.factor(DF$subject_id)
  388. fullDF = DF
  391. if (estType != "SILA") {
  392. facSmooth = T
  393. } else {
  394. facSmooth = F
  395. }
  398. if (facSmooth) {
  400. # ordered factor approach to get test statistics for GAMM interaction
  401. dat = DF
  402. dat$brainvarfac = dat$brainvarHigh
  403. stackdat = rbind(
  404. dat %>% mutate(highlow = "high"),
  405. dat %>% mutate(highlow = "low")
  406. )
  407. stackdat$brainvarfac[stackdat$highlow == "low"] = dat$brainvarLow
  408. stackdat = mutate(stackdat,
  409. ohighlow = factor(highlow, levels = c("low","high"),ordered = T))
  412. gamm.trajectories = gamm4(brainvarfac ~ s(diff_mriAge_predABpos_flip, by = as.factor(highlow)) + as.factor(highlow) + visit_age + subject_sex + cohort + scanStrength + ICV,
  413. data = stackdat, random = ~ (1 |subject_id))
  414. gamm.sum = summary(gamm.trajectories$gam)
  415. g = gamm.trajectories$gam
  416. plot.gam(g)
  419. # estimate smooth for set reflevel and a smoothed difference between ref and other levels
  420. ogamm.trajectories = gamm4(brainvarfac ~ as.factor(highlow) + s(diff_mriAge_predABpos_flip) + s(diff_mriAge_predABpos_flip, by = ohighlow) + visit_age + subject_sex + cohort + scanStrength + ICV,
  421. data = stackdat, random = ~ (1 |subject_id))
  422. ogamm.sum = summary(ogamm.trajectories$gam)
  423. plot.gam(ogamm.trajectories$gam)
  426. # thickness trajectory in high Aβ regions was significantly different than low Aβ regions
  427. ogamm.sum
  428. }
  432. # three analyses to run through - select which here
  433. regionTest = "high"
  434. regionTest = "low"
  435. regionTest = "cortex"
  439. if (regionTest == "high") {
  440. loopend = 1
  441. } else if (regionTest == "low") {
  442. loopend = 1
  443. } else if (regionTest == "cortex") {
  444. loopend = length(ROIs)
  445. }
  449. for (i in 1:loopend) {
  450. # tic()
  451. print(paste(i,"/",length(ROIs)))
  452. if (i == 1) {
  453. derivMat = gratiaderivMat = gratiaderivCIlwrMat = gratiaderivCIuprMat = matrix(NA, nrow = 1000, ncol = length(ROIs))
  454. fitMat = matrix(NA, nrow = 100, ncol = length(ROIs))
  455. RR = list()
  456. p_derivs = p_derivs_se = list()
  457. yearsBeforepredAB_on_ci_exclusion = c()
  458. yearsBeforepredAB_on_se_exclusion = c()
  459. maxaccels = c()
  460. }
  461. DF = fullDF
  462. setTxtProgressBar(pb,i)
  463. set.seed(123)
  466. # ab regions low / high
  467. if (regionTest == "high") {
  468. ROI = "high AB composite"
  469. DF$brainvar = rowMeans(DF[,ab_regions_high])
  470. palcol = pal[5]
  471. anaTitle = "Aβ high"
  472. }
  474. if (regionTest == "low") {
  475. ROI = "low AB composite"
  476. DF$brainvar = rowMeans(DF[,ab_regions_low])
  477. palcol = pal[2]
  478. anaTitle = "Aβ low"
  479. }
  481. if (regionTest == "cortex") {
  482. ROI = ROIs[i]
  483. DF$brainvar = DF[[ROI]]
  484. palcol = "darkgrey"
  485. anaTitle = ROI
  486. }
  489. g_convall = gamm4(brainvar ~ s(diff_mriAge_predABpos_flip) + visit_age + subject_sex + cohort + scanStrength + ICV, data = DF, random = ~ (1 | subject_id))
  490. g_convall_mgcv = gam(
  491. brainvar ~
  492. s(diff_mriAge_predABpos_flip) +
  493. s(subject_id, bs = "re") +
  494. visit_age + subject_sex + cohort + scanStrength + ICV,
  495. data = DF,
  496. method = "REML"
  497. )
  498. summary(g_convall$gam)
  499. g_sum = summary(g_convall_mgcv)
  500. RR[[i]] = g_sum
  503. lmm_convall = lmer(
  504. brainvar ~ diff_mriAge_predABpos_flip + visit_age + subject_sex + cohort + scanStrength + ICV + (1 | subject_id),
  505. data = DF
  506. )
  507. # plot.gam(g_convall$gam, residuals = T)
  508. summary(lmm_convall)
  511. # clearly the covariates are important
  512. ggplot(DF, aes(y=brainvar, x = diff_mriAge_predABpos_flip)) + geom_point(aes(group = subject_id)) + geom_smooth(method = "gam")
  513. plot(effect("diff_mriAge_predABpos_flip", lmm_convall, residuals=TRUE)) #warning is due to scaling of Y and X
  516. predictions <- DF %>%
  517. mutate(visit_age = mean(DF$visit_age), subject_sex = "Female", cohort = "ADNINC", scanStrength = "1-5T", ICV=0) %>%
  518. select(diff_mriAge_predABpos_flip, visit_age,subject_sex,ICV,mri_info_site_name,scanStrength, cohort) %>%
  519. predict(g_convall$gam, newdata = ., se.fit = T)
  522. residualsg <- residuals(g_convall$mer)
  523. DF$partial_residuals = predictions$fit + residualsg
  524. DF$fit = predictions$fit
  525. DF$sefit = predictions$se.fit
  526. DF$cifit = predictions$se.fit*1.96
  530. #colour palette ---
  531. pal = wesanderson::wes_palettes$FantasticFox1
  532. pointcol = "#6faca8"
  533. #colour palette ---
  535. if (estType != "SILA") {
  536. tmpDF = DF %>% filter(diff_mriAge_predABpos >= 1)
  537. } else {
  538. tmpDF = DF %>% filter(time_from_sila_threshold >= 1)
  539. }
  540. (fig1 =
  541. DF %>% #filter(subject_id != adnioutlier) %>%
  542. ggplot(.) +
  543. geom_line(data=tmpDF,aes(x=diff_mriAge_predABpos_flip,partial_residuals,group=subject_id),color=pointcol,alpha=0.6, size=0.5) +
  544. geom_point(data=tmpDF,aes(x=diff_mriAge_predABpos_flip,partial_residuals,group=subject_id),color=pointcol,stat="identity",alpha=1, size=0.5) +
  545. # geom_ribbon(data=dug,aes(x=diff_mriAge_predABpos_flip,ymin=fit-CI,ymax=fit+CI),alpha=.7,show.legend=F,fill="dark grey") +
  546. geom_line(data=tmpDF,aes(x=diff_mriAge_predABpos_flip,y=fit),col="black") +
  547. ggtitle(ROI) +
  548. labs(x = "Age") +
  549. theme_classic() + mytheme)
  553. # Extract plotting data without rendering the plot
  554. plot_data = plot.gam(g_convall$gam, select = 1, se = TRUE, rug = FALSE, shade = T, pages = 0)
  557. # get gam intercept
  558. g_sum = summary(g_convall$gam)
  559. g_sum$p.coeff
  562. # The output is a list, extract x, fit, and se
  563. df = data.frame(
  564. x = plot_data[[1]]$x,
  565. fit = plot_data[[1]]$fit + g_sum$p.coeff[1],
  566. se = plot_data[[1]]$se
  567. )
  570. # Calculate upper and lower confidence intervals
  571. df$upper = df$fit + 1 * df$se #NB! help(plot.gam) shows it is already using CI
  572. df$lower = df$fit - 1 * df$se
  575. ggplot(df, aes(x = x, y = fit)) +
  576. geom_line() +
  577. geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2) +
  578. labs(
  579. ) +
  580. theme_minimal()
  583. # Approximate first derivative
  584. df$deriv = c(NA, diff(df$fit) / diff(df$x))
  586. # Find index of minimum (i.e., where uptick starts)
  587. min_idx = which.min(df$fit)
  589. # Optionally, find first point where derivative turns positive
  590. first_uptick_idx = which(df$deriv > 0 & seq_along(df$deriv) > min_idx)[1]
  592. # Get corresponding x-value
  593. uptick_point = df$x[first_uptick_idx]
  594. uptick_pointy = df$fit[first_uptick_idx]
  597. fitMat[,i] = df$fit
  600. (pderiv0 = ggplot(df, aes(x = x, y = fit)) +
  601. geom_line() +
  602. labs(
  603. x = "Years to Aβ+",
  604. y = "Thickness"
  605. ) +
  606. # geom_point(aes(x = uptick_point, y = uptick_pointy), col = "black", size = 5) +
  607. theme_minimal())
  610. (pderiv0_ci = ggplot(df, aes(x = x, y = fit)) +
  611. geom_line(col = palcol) +
  612. ggtitle(ROI) +
  613. geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.1, col = palcol, fill = palcol) +
  614. labs(
  615. x = "Years to Aβ+",
  616. y = "Thickness"
  617. ) +
  618. # geom_point(aes(x = uptick_point, y = uptick_pointy), col = "black", size = 5) +
  619. theme_minimal())
  622. # first derivative
  623. df$deriv = c(NA, diff(df$fit) / diff(df$x))
  624. df$deriv_lower = c(NA, diff(df$lower) / diff(df$x))
  625. df$deriv_upper = c(NA, diff(df$upper) / diff(df$x))
  628. (pderiv1 = ggplot(df, aes(x = x, y = deriv)) +
  629. geom_line() +
  630. geom_hline(yintercept = 0, linetype = "dashed") +
  631. labs(x = "Years to Aβ+", y = "Rate of change") +
  632. theme_minimal())
  635. (pderiv1_se = ggplot(df, aes(x = x, y = deriv)) +
  636. geom_line() +
  637. geom_ribbon(aes(ymin = deriv_lower, ymax = deriv_upper), alpha = 0.2) +
  638. geom_hline(yintercept = 0, linetype = "dashed") +
  639. labs(x = "Years to Aβ+", y = "Rate of change") +
  640. theme_minimal())
  643. derivMat[,i] = df$deriv
  646. # Compute second derivative (numerical)
  647. second_derivative = diff(diff(df$fit)) / diff(df$x[-1])
  649. # Align x-axis (midpoints between x values)
  650. second_x = df$x[-c(1, 2)] + diff(df$x)[-1] / 2
  652. #point of maximum accerelated change
  653. maxaccel_y = second_derivative[which(second_derivative == max(second_derivative))]
  654. maxaccel_x = second_x[which(second_derivative == max(second_derivative))]
  656. if (length(maxaccel_x) > 1) {
  657. maxaccel_y = NA
  658. maxaccel_x = NA
  659. }
  660. if (is.na(maxaccel_x)) {
  661. maxaccels[i] = NA
  662. } else {
  663. maxaccels[i] = maxaccel_x
  665. pderiv2 = ggplot(data.frame(x = second_x, second_derivative = second_derivative), aes(x = x, y = second_derivative)) +
  666. geom_line() +
  667. geom_hline(yintercept = 0, linetype = "dashed") +
  668. labs(y = "Acceleration", x = "Years to Aβ+") +
  669. geom_point(aes(x = maxaccel_x[i], y = maxaccel_y[i]), col = "black", size = 2) +
  670. theme_minimal()
  672. # check all derivatives
  673. (pcheckderiv = cowplot::plot_grid(pderiv0, pderiv1, pderiv2, nrow = 3))
  675. }
  679. # NB! se of calculated derivative not ideal - using gratia instead
  680. (psumse1 = cowplot::plot_grid(pderiv0_ci, pderiv1_se, nrow = 3))
  683. # Compute first derivative of the smooth term
  684. d1 = gratia::derivatives(g_convall_mgcv, term = "s(diff_mriAge_predABpos_flip)", interval = "confidence", n = 1000)
  685. d1 = as.data.frame(d1)
  686. d1$.lower_ci[which(d1$.lower_ci > 0)]
  689. # point at which the derivative CI excluded zero used to estimate the rank order of thickness changes
  690. cross_point = which(d1$.lower_ci > 0)[1]
  693. data.frame(d1$diff_mriAge_predABpos_flip, d1$.lower_ci, logical = d1$.lower_ci > 0)
  696. if (is.na(cross_point)) {
  697. print("CIs do not exclude zero")
  698. yearsBeforepredAB_on_ci_exclusion[i] = "never"
  699. exclude0 = 0
  700. } else if (cross_point == 1) {
  701. print("CIs always exclude zero")
  702. yearsBeforepredAB_on_ci_exclusion[i] = "always"
  703. cross_y = 0
  704. cross_x = d1$diff_mriAge_predABpos_flip[cross_point]
  705. exclude0 = 1
  706. } else {
  707. print("CIs exclude zero")
  708. cross_y = d1$.lower_ci[cross_point]
  709. cross_x = d1$diff_mriAge_predABpos_flip[cross_point]
  710. yearsBeforepredAB_on_ci_exclusion[i] = cross_x
  711. exclude0 = 1
  712. }
  715. #repeat for SE (since not all cross)
  716. d1$upper_se = d1$.derivative + d1$.se
  717. d1$lower_se = d1$.derivative - d1$.se
  721. # point at which the derivative SE excluded zero
  722. cross_point_se = which(d1$lower_se > 0)[1]
  724. # two instances where the deriv excluded zero on the downward trajectory - fixed to first point on upward trajectory
  725. data.frame(d1$diff_mriAge_predABpos_flip, d1$lower_se, logical = d1$lower_se > 0)
  726. # if (estType == "SILA" & convonly == 0) { if (i == 12) { cross_point_se = 420 } } # fix the one across full group
  727. # if (estType == "SILA" & convonly == 0) { if (i == 2) { cross_point_se = 320 } } # fix the one across full group
  728. # if (estType == "SILA" & convonly == 0) { if (i == 9) { cross_point_se = 537 } } # fix the one across full group
  729. # if (estType == "SILA" & convonly == 0) { if (i == 10) { cross_point_se = 322 } } # fix the one across full group
  730. # if (estType == "SILA" & convonly == 1) { if (i == 4) { cross_point_se = 262 } } # fix if only converters
  731. # if (estType == "SILA" & convonly == 1) { if (i == 2) { cross_point_se = 281 } } # fix the one across full group
  732. # if (estType == "SILA" & convonly == 1) { if (i == 15) { cross_point_se = 296 } } # fix the one across full group
  733. # if (estType == "SILA" & convonly == 1) { if (i == 21) { cross_point_se = NA } } # fix the one across full group
  735. if (is.na(cross_point_se)) {
  736. print("SEs do not exclude zero")
  737. yearsBeforepredAB_on_se_exclusion[i] = "never"
  738. } else if (cross_point_se == 1) {
  739. print("SEs always exclude zero")
  740. cross_yse = 0
  741. cross_xse = d1$diff_mriAge_predABpos_flip[cross_point_se]
  742. yearsBeforepredAB_on_se_exclusion[i] = "always"
  743. } else {
  744. print("SEs exclude zero")
  745. cross_yse = d1$lower_se[cross_point_se]
  746. cross_xse = d1$diff_mriAge_predABpos_flip[cross_point_se]
  747. yearsBeforepredAB_on_se_exclusion[i] = cross_xse
  748. }
  751. # check crossing points
  752. if (!is.na(cross_point)) {
  753. ggplot(d1, aes(x = diff_mriAge_predABpos_flip, y = .derivative)) +
  754. geom_line() +
  755. geom_ribbon(aes(ymin = .lower_ci, ymax = .upper_ci), alpha = 0.1) +
  756. geom_hline(yintercept = 0, linetype = "dashed") +
  757. geom_point(aes(x = cross_x, y = cross_y), col = "black", size = 5)
  758. }
  759. if (!is.na(cross_point_se) & cross_point_se != 1) {
  760. ggplot(d1, aes(x = diff_mriAge_predABpos_flip, y = .derivative)) +
  761. geom_line() +
  762. geom_ribbon(aes(ymin = lower_se, ymax = upper_se), alpha = 0.1) +
  763. geom_hline(yintercept = 0, linetype = "dashed") +
  764. geom_point(aes(x = cross_xse, y = cross_yse), col = "black", size = 5)
  765. }
  769. (pderiv1_gratia_ci = ggplot(d1, aes(x = diff_mriAge_predABpos_flip, y = .derivative)) +
  770. geom_line(col = palcol) +
  771. geom_ribbon(aes(ymin = .lower_ci, ymax = .upper_ci), alpha = 0.1, colour = palcol, fill = palcol) +
  772. geom_hline(yintercept = 0, linetype = "dashed") +
  773. labs(
  774. x = "Years to Aβ+",
  775. y = "First derivative"
  776. ) +
  777. # geom_point(aes(x = cross_x, y = cross_y), col = "black", size = 1) +
  778. theme_classic() + mytheme
  779. )
  782. (pderiv1_gratia_se = ggplot(d1, aes(x = diff_mriAge_predABpos_flip, y = .derivative)) +
  783. geom_line() +
  784. geom_ribbon(aes(ymin = lower_se, ymax = upper_se), alpha = 0.2) +
  785. geom_hline(yintercept = 0, linetype = "dashed") +
  786. labs(
  787. x = "Years to Aβ+",
  788. y = "First derivative"
  789. ) +
  790. # geom_point(aes(x = cross_xse, y = cross_yse), col = "black", size = 1) +
  791. theme_classic() + mytheme
  792. )
  796. if (regionTest == "cortex") {
  797. if (!is.na(cross_point)) {
  798. # add crossing point to plot
  799. pderiv1_gratia_ci = pderiv1_gratia_ci + geom_point(aes(x = cross_x, y = cross_y), size = 6, fill = "#cbac09", col = "#cbac09")
  800. }
  801. if (!is.na(cross_point_se)) {
  802. pderiv1_gratia_se = pderiv1_gratia_se + geom_point(aes(x = cross_xse, y = cross_yse), size = 6, fill = "#cbac09", col = "#cbac09")
  803. }
  804. }
  807. # save gratia outputs
  808. gratiaderivMat[,i] = d1$.derivative
  809. gratiaderivCIlwrMat[,i] = d1$.lower_ci
  810. gratiaderivCIuprMat[,i] = d1$.upper_ci
  813. # main plot
  814. (
  815. p_combine_fitderiv_ci = ggpubr::ggarrange(
  816. pderiv0_ci + theme_classic() + mytheme + ggtitle(anaTitle) +
  817. pderiv1_gratia_ci,
  818. nrow = 1,
  819. align = "hv"
  820. )
  821. )
  823. (
  824. p_combine_fitderiv_se = ggpubr::ggarrange(
  825. pderiv0_ci + theme_classic() + mytheme + ggtitle(anaTitle) +
  826. pderiv1_gratia_se,
  827. nrow = 1,
  828. align = "hv"
  829. )
  830. )
  833. if (regionTest == "cortex") {
  834. (
  835. p_combine_fitderiv_ci = ggpubr::ggarrange(
  836. pderiv0_ci + theme_classic() + mytheme + ggtitle(anaTitle) + theme(axis.ticks = element_blank(), axis.line = element_blank()),
  837. pderiv1_gratia_ci + theme(axis.ticks = element_blank(), axis.line = element_blank()),
  838. nrow = 1,
  839. align = "hv"
  840. )
  841. )
  843. }
  844. p_derivs[[i]] = p_combine_fitderiv_ci
  845. p_derivs_se[[i]] = p_combine_fitderiv_se
  848. if (savefigs == 1) {
  849. if (regionTest == "high") {
  850. print("saving high plot")
  851. # ggsave(plot = p_combine_fitderiv_ci,
  852. # filename = paste0("/cluster/projects/p274/projects/p040-ad_change/Berkeley/paper2/figs_yearsBeforeAB/p_thickTraj_ABhigh.pdf"),
  853. # width=20, height=9, units="cm", dpi=600, device = cairo_pdf
  854. # )
  855. } else if (regionTest == "low") {
  856. print("saving low plot")
  857. # ggsave(plot = p_combine_fitderiv_ci,
  858. # filename = paste0("/cluster/projects/p274/projects/p040-ad_change/Berkeley/paper2/figs_yearsBeforeAB/p_thickTraj_ABlow.pdf"),
  859. # width=20, height=9, units="cm", dpi=600, device = cairo_pdf
  860. # )
  861. } else {
  862. # ggsave(plot = p_combine_fitderiv_ci,
  863. # filename = paste0("/cluster/projects/p274/projects/p040-ad_change/Berkeley/paper2/figs_yearsBeforeAB/p_thickTraj_", ROI, ".pdf"),
  864. # width=20, height=9, units="cm", dpi=600, device = cairo_pdf
  865. # )
  866. }
  867. }
  869. }
  871. # check cross point on all 24 roi derivs
  872. p_derivs[[1]]
  873. p_derivs[[2]]
  874. p_derivs[[3]]
  875. p_derivs[[4]]
  876. p_derivs[[5]]
  877. p_derivs[[6]]
  878. p_derivs[[7]]
  879. p_derivs[[8]]
  880. p_derivs[[9]]
  881. p_derivs[[10]]
  882. p_derivs[[11]]
  883. p_derivs[[12]]
  884. p_derivs[[13]]
  885. p_derivs[[14]]
  886. p_derivs[[15]]
  887. p_derivs[[16]]
  888. p_derivs[[17]]
  889. p_derivs[[18]]
  890. p_derivs[[19]]
  891. p_derivs[[20]]
  892. p_derivs[[21]]
  893. p_derivs[[22]]
  894. p_derivs[[23]]
  895. p_derivs[[24]]
  898. p_derivs_se[[1]]
  899. p_derivs_se[[2]]
  900. p_derivs_se[[3]]
  901. p_derivs_se[[4]]
  902. p_derivs_se[[5]]
  903. p_derivs_se[[6]]
  904. p_derivs_se[[7]]
  905. p_derivs_se[[8]]
  906. p_derivs_se[[9]]
  907. p_derivs_se[[10]]
  908. p_derivs_se[[11]]
  909. p_derivs_se[[12]]
  910. p_derivs_se[[13]]
  911. p_derivs_se[[14]]
  912. p_derivs_se[[15]]
  913. p_derivs_se[[16]]
  914. p_derivs_se[[17]]
  915. p_derivs_se[[18]]
  916. p_derivs_se[[19]]
  917. p_derivs_se[[20]]
  918. p_derivs_se[[21]]
  919. p_derivs_se[[22]]
  920. p_derivs_se[[23]]
  921. p_derivs_se[[24]]
  922. # all confirmed correct
  925. # rank order on thickness change
  926. data.frame(yearsBeforepredAB_on_ci_exclusion,
  927. yearsBeforepredAB_on_se_exclusion
  928. )
  930. rankThickTraj = data.frame(rois = ROIs,
  931. rankedThickTraj = yearsBeforepredAB_on_ci_exclusion,
  932. rankedThickTrajSE = yearsBeforepredAB_on_se_exclusion,
  933. rankedThickTrajMaxAccel = maxaccels)
  935. # set as 0 if CI always excludes 0
  936. rankThickTraj$rankedThickTraj[rankThickTraj$rankedThickTraj == "always"] = 0
  937. rankThickTraj$rankedThickTraj[rankThickTraj$rankedThickTraj == "never"] = NA
  938. rankThickTraj$rankedThickTrajSE[rankThickTraj$rankedThickTrajSE == "always"] = 0
  939. rankThickTraj$rankedThickTrajSE[rankThickTraj$rankedThickTrajSE == "never"] = NA
  940. rankThickTraj$rankedThickTraj = as.numeric(rankThickTraj$rankedThickTraj)
  941. rankThickTraj$rankedThickTrajSE = as.numeric(rankThickTraj$rankedThickTrajSE)
  945. # rank and reorder ROIs based on the point at which the CI / SE of the derivative crosses 0
  946. # NB! no need to reverse as yearsBeforepredAB was flipped in model (diff_mriAge_predABpos_flip)
  947. # CI
  948. rankThickTraj$rankedThickTrajRev <- ifelse(
  949. is.na(rankThickTraj$rankedThickTraj),
  950. NA,
  951. ifelse(rankThickTraj$rankedThickTraj == 0,
  952. 0,
  953. rank(rankThickTraj$rankedThickTraj, ties.method = "first"))
  954. )
  955. rankThickTraj %<>% arrange(rankedThickTrajRev)
  958. # SE
  959. rankThickTraj$rankedThickTrajSERev <- ifelse(
  960. is.na(rankThickTraj$rankedThickTrajSE),
  961. NA,
  962. ifelse(rankThickTraj$rankedThickTrajSE == 0,
  963. 0,
  964. rank(rankThickTraj$rankedThickTrajSE, ties.method = "first"))
  965. )
  966. rankThickTraj %<>% arrange(rankedThickTrajSERev)
  970. # load in amyloid order
  971. rankAB = fread(file.path(b, "reproduce/data/rankABpredTraj.csv"))
  972. rankAB %<>% rename(rankedABTraj = rankedTraj,
  973. rankedABTrajRev = revRankedTraj)
  974. rankThickTraj$rois = paste0("CTX_", toupper(gsub("_thickness.aparcnative71", "", rankThickTraj$rois)), "_SUVR")
  975. rankThickTraj$rois %in% rankAB$rois
  977. rankedBoth = merge(
  978. rankThickTraj, rankAB
  979. )
  981. # filter data where there is no rank order for thickness (i.e. ROI derivative did not exclude 0 and thus could not be ranked)
  982. rankedBoth_cut = rankedBoth %>% filter(!is.na(rankedThickTrajRev))
  985. # make ranking plots
  986. rankedBoth$roi = gsub("CTX_", "", rankedBoth$rois)
  987. rankedBoth$roi = gsub("_SUVR", "", rankedBoth$roi)
  988. rankedBoth$roi = tolower(rankedBoth$roi)
  989. orderplot = arrange(rankedBoth, rankedABTrajRev) %>% select(roi)
  990. rankedBoth$roi = factor(rankedBoth$roi, levels = rev(orderplot$roi))
  991. rankRange = range(rankedBoth$rankedABTrajRev, na.rm= T)
  993. table(rankedBoth$rankedABTrajRev)
  994. rankedBoth$rankedABTrajRev = rankedBoth$rankedABTrajRev-1 # make 0 indexed (so colours match up with thickness rank plot)
  995. (p_rankAB_FDRregions = ggplot(rankedBoth %>%
  996. arrange(rankedABTrajRev),
  997. aes(x = rankedABTrajRev, y = factor(roi, levels = rev(unique(roi))))) +
  998. geom_bar(aes(x=rankedABTrajRev, y =roi), stat = "identity", colour = "grey", alpha = 0.5, width = 0.01, size=0.5) +
  999. geom_point(aes(colour = rankedABTrajRev), fill = "white", shape = 21, stroke = 2, size = 5) +
  1000. scale_colour_viridis(option = "E", direction = -1, name = "Rank", limits = c(0, rankRange[2]), oob = scales::squish) +
  1001. theme_classic() +
  1002. mytheme +
  1003. labs(x = "Rank (Aβ)",y = NULL) + theme(legend.position = "none")
  1004. )
  1007. # ggsave(plot = p_rankAB_FDRregions,
  1008. # filename = paste0("/cluster/projects/p274/projects/p040-ad_change/Berkeley/paper2/figs_yearsBeforeAB/p_rankAB_FDRregions.png"),
  1009. # width=7, height=26, units="cm", dpi=600)
  1012. # ensure rank plot for thickness is ordered the same
  1013. (p_rankThick_frontal = ggplot(rankedBoth,
  1014. aes(x = rankedThickTrajRev, y = factor(roi, levels = rev(unique(roi))))) +
  1015. geom_bar(aes(x=rankedThickTrajRev, y =roi), stat = "identity", colour = "grey", alpha = 0.5, width = 0.01, size=0.5) +
  1016. geom_point(aes(colour = rankedThickTrajRev), fill = "white", shape = 21, stroke = 2, size = 5) +
  1017. scale_colour_viridis(option = "E", direction = -1, name = "Rank", limits = c(0, rankRange[2]), oob = scales::squish) +
  1018. theme_classic() +
  1019. mytheme +
  1020. labs(x = "Rank (Thickness)",y = NULL) + theme(legend.position = "none")
  1021. )
  1023. table(rankedBoth$rankedThickTrajSERev)
  1024. (p_rankThick_frontal_bySE = ggplot(rankedBoth,
  1025. aes(x = rankedThickTrajSERev, y = factor(roi, levels = rev(unique(roi))))) +
  1026. geom_bar(aes(x=rankedThickTrajSERev, y =roi), stat = "identity", colour = "grey", alpha = 0.5, width = 0.01, size=0.5) +
  1027. geom_point(aes(colour = rankedThickTrajSERev), fill = "white", shape = 21, stroke = 2, size = 5) +
  1028. scale_colour_viridis(option = "E", direction = -1, name = "Rank", limits = c(0, rankRange[2]), oob = scales::squish) +
  1029. theme_classic() +
  1030. mytheme +
  1031. labs(x = "Rank (Thickness)",y = NULL) + theme(legend.position = "none")
  1032. )
  1034. cp1 = cowplot::plot_grid(p_rankAB_FDRregions, p_rankThick_frontal + theme(axis.text.y = element_text(color = "transparent"), axis.line.y = element_blank(), axis.ticks.y = element_blank()))
  1035. cp2 = cowplot::plot_grid(p_rankAB_FDRregions, p_rankThick_frontal_bySE + theme(axis.text.y = element_text(color = "transparent"), axis.line.y = element_blank(), axis.ticks.y = element_blank()))
  1037. # ggsave(plot = p_rankAB_FDRregions,
  1038. # filename = paste0("/cluster/projects/p274/projects/p040-ad_change/Berkeley/paper2/figs_yearsBeforeAB/p_rankAB_FDRregions.pdf"),
  1039. # width=12, height=26, units="cm", dpi=600, device = cairo_pdf)
  1041. # ggsave(plot = p_rankThick_frontal,
  1042. # filename = paste0("/cluster/projects/p274/projects/p040-ad_change/Berkeley/paper2/figs_yearsBeforeAB/p_rankRThick_FDRregions.pdf"),
  1043. # width=12, height=26, units="cm", dpi=600, device = cairo_pdf)
  1047. table(rankedBoth_cut$rankedThickTrajRev)
  1048. pear1 = cor.test(rankedBoth_cut$rankedThickTrajRev, rankedBoth_cut$rankedABTrajRev)
  1049. pear2 = cor.test(rankedBoth_cut$rankedThickTrajSERev, rankedBoth_cut$rankedABTrajRev)
  1051. pval1 = format(
  1052. tidy(pear1)$p.value[1],
  1053. scientific = TRUE, digits = 2)
  1054. bval1 = round(
  1055. tidy(pear1)$estimate[1],
  1056. digits = 2)
  1057. pval2 = format(
  1058. tidy(pear2)$p.value[1],
  1059. scientific = TRUE, digits = 2)
  1060. bval2 = round(
  1061. tidy(pear2)$estimate[1],
  1062. digits = 2)
  1064. table(rankedBoth_cut$rankedThickTrajRev)
  1067. # correlation plot based on CI (fig 5g)
  1068. (p_rankCI = ggplot(rankedBoth_cut, aes(x = rankedThickTrajRev, y = rankedABTrajRev)) +
  1069. geom_point(aes(col = rankedThickTrajRev), size = 3) +
  1070. scale_colour_viridis(option = "E", direction = -1, name = "Rank", limits = c(0, rankRange[2]), oob = scales::squish) +
  1071. theme_classic() +
  1072. mytheme +
  1073. # coord_fixed() +
  1074. geom_smooth(method = "lm", se = T, col = "black", alpha = .1) +
  1075. labs(x = "Rank (Thickness)",y = "Rank (Aβ)") + theme(legend.position = "none") +
  1076. ggtitle("order of CIs excluding 0") +
  1077. annotate("text", x = 1, y = max(rankedBoth_cut$rankedABTrajRev)+2, label = paste("p =", pval1), color = "black", hjust = 0) +
  1078. annotate("text", x = 1, y = max(rankedBoth_cut$rankedABTrajRev)+3, label = paste("B =", bval1), color = "black", hjust = 0)
  1079. )
  1081. # ggsave(plot = p_rankCI,
  1082. # filename = paste0("/cluster/projects/p274/projects/p040-ad_change/Berkeley/paper2/figs_yearsBeforeAB/p_rankCI.pdf"),
  1083. # width=9, height=12, units="cm", dpi=600, device = cairo_pdf)
  1085. # ggsave(plot = p_rankCI,
  1086. # filename = paste0("/cluster/projects/p274/projects/p040-ad_change/Berkeley/paper2/figs_yearsBeforeAB/p_rankCI_", estType, "convonly", convonly, ".pdf"),
  1087. # width=9, height=12, units="cm", dpi=600, device = cairo_pdf)
  1090. table(rankedBoth_cut$rankedThickTrajSERev)
  1093. # correlation plot based on SE (fig 5h)
  1094. (p_rankSE = ggplot(rankedBoth_cut, aes(x = rankedThickTrajSERev, y = rankedABTrajRev)) +
  1095. geom_point(aes(col = rankedThickTrajSERev), size = 3) +
  1096. scale_colour_viridis(option = "E", direction = -1, name = "Rank", limits = c(0, rankRange[2]), oob = scales::squish) +
  1097. theme_classic() +
  1098. mytheme +
  1099. # coord_fixed() +
  1100. geom_smooth(method = "lm", se = T, col = "black", alpha = .2) +
  1101. labs(x = "Rank (Thickness)",y = "Rank (Aβ)") + theme(legend.position = "none") +
  1102. ggtitle("order of SEs excluding 0") +
  1103. annotate("text", x = 1, y = max(rankedBoth_cut$rankedABTrajRev)+2, label = paste("p =", pval2), color = "black", hjust = 0) +
  1104. annotate("text", x = 1, y = max(rankedBoth_cut$rankedABTrajRev)+3, label = paste("B =", bval2), color = "black", hjust = 0)
  1105. )
  1107. # ggsave(plot = p_rankSE,
  1108. # filename = paste0("/cluster/projects/p274/projects/p040-ad_change/Berkeley/paper2/figs_yearsBeforeAB/p_rankSE.pdf"),
  1109. # width=9, height=12, units="cm", dpi=600, device = cairo_pdf)
  1111. # ggsave(plot = p_rankSE,
  1112. # filename = paste0("/cluster/projects/p274/projects/p040-ad_change/Berkeley/paper2/figs_yearsBeforeAB/p_rankSE_", estType, "convonly", convonly, ".pdf"),
  1113. # width=9, height=12, units="cm", dpi=600, device = cairo_pdf)

10-rankThicknessChange.r at commit 8606e75, under MIT · at the source

Overview

Authors: James M. Roe1,2, William J. Jagust3, Susan M. Landau3, Theresa M. Harrison3, Håkon Grydeland1, Maksim Slivka1, José-Luis Alatorre-Warren1, Pablo F. Garrido1, Øystein Sørensen1, Edvard O. S. Grødem1,2, Tyler J. Ward3, Esten H. Leonardsen1,4, Alice Murphy3, JiaQie Lee3, Tormod Fladby5,6, Atle Bjørnerud1,2, Kristine B. Walhovd1,2, Anders M. Fjell1,2, Didac Vidal-Piñeiro1, Yunpeng Wang1
ORCID iDs: James M. Roe, William J. Jagust, Maksim Slivka, Pablo F. Garrido, Øystein Sørensen, Tyler J. Ward, Tormod Fladby, Kristine B. Walhovd, Anders M. Fjell, Didac Vidal-Piñeiro, Yunpeng Wang
  1. Center for Lifespan Changes in Brain and Cognition (LCBC), Department of Psychology, University of Oslo,Oslo, Norway
  2. Computational Radiology and Artificial Intelligence, Department of Radiology and Nuclear Medicine, Oslo University Hospital,Oslo, Norway
  3. Department of Neuroscience, University of California, Berkeley,Berkeley, CA USA
  4. Centre for Precision Psychiatry, Oslo University Hospital & Institute of Clinical Medicine, University of Oslo,Oslo, Norway
  5. Department of Neurology, Akershus University Hospital,Lørenskog, Norway
  6. Institute for Clinical Medicine, University of Oslo,Oslo, Norway
Institutions: University of Oslo (Norway); Oslo University Hospital (Norway); University of California, Berkeley (United States); Akershus University Hospital (Norway)
Journal: Nature neuroscience, volume 29, issue 9, pages 2164-2175
Dates: received 24 August 2025; accepted 9 June 2026; published online 19 August 2026; in print 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41593-026-02363-4 · PMID 42618753 · PMCID PMC13533846 · OpenAlex W7203752607
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: structural MRI / diffusion (modality), PET / SPECT (modality), human (organism), Alzheimer's / dementia (population)
Methods: Statistics, Connectivity, Smoothing, state filtering, decompositions, Preprocessing
Keywords: Neural ageing, Predictive markers
MeSH: Alzheimer Disease*, Amyloid beta-Peptides*, Cerebral Cortex*, Aged, Aged, 80 and over, Disease Progression, Female, Humans, Longitudinal Studies, Magnetic Resonance Imaging, Male, Positron-Emission Tomography (* major topic)
Topic: Dementia and Cognitive Impairment Research (Psychiatry and Mental health, Medicine), according to OpenAlex
Funding: University of California, Berkeley | Peder Sather Center for Advanced Study (Peder Sather Center, UC Berkley) (Markers of brain activity in early Alzheimer's Disease: The Berkeley-Oslo alliance); EC | EU Framework Programme for Research and Innovation H2020 | H2020 Priority Excellent Science | H2020 European Research Council (H2020 Excellent Science - European Research Council) (313440, 725025)
Citations: cited by 1 paper (Europe PMC); 68 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

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

jamesmroe/yearsBeforeAB

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 8606e753805b25245b52466dbf4d2c26cb94a43a, 7 May 2026
Languages: R (14), Shell (1)
Size: 264 files, 15 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (requirements.txt), documentation
Not found: CITATION.cff, tests, continuous integration
Tools: tidyverse (14 files), data.table (12 files), broom (8 files), cowplot (7 files), ggpubr (6 files), mgcv (5 files), patchwork (4 files), easystats (3 files), ggseg (3 files), ggplot2 (2 files), lmerTest (2 files), lme4 (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
17 files
At the source: github.com/jamesmroe/yearsBeforeAB

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s41593-026-02363-4.

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

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

  • it says that the data are available on request

Read it in the paper: doi.org/10.1038/s41593-026-02363-4.

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

  • Publisher: n/a → Nature Portfolio

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 20 authors, 2 keywords, 12 MeSH terms, 2 funders, 65 references.

Cite

This paper

Roe, J. M., Jagust, W. J., Landau, S. M., Harrison, T. M., Grydeland, H., Slivka, M., Alatorre-Warren, J.-L., Garrido, P. F., Sørensen, Ø., Grødem, E. O. S., Ward, T. J., Leonardsen, E. H., Murphy, A., Lee, J., Fladby, T., Bjørnerud, A., Walhovd, K. B., Fjell, A. M., Vidal-Piñeiro, D., & Wang, Y. (2026). Cortical thickness changes precede high levels of amyloid by at least 7 years. Nature neuroscience, 29(9), 2164-2175. https://doi.org/10.1038/s41593-026-02363-4

BibTeX

@article{roe2026cortical,
author = {Roe, James M. and Jagust, William J. and Landau, Susan M. and Harrison, Theresa M. and Grydeland, Håkon and Slivka, Maksim and Alatorre-Warren, José-Luis and Garrido, Pablo F. and Sørensen, Øystein and Grødem, Edvard O. S. and Ward, Tyler J. and Leonardsen, Esten H. and Murphy, Alice and Lee, JiaQie and Fladby, Tormod and Bjørnerud, Atle and Walhovd, Kristine B. and Fjell, Anders M. and Vidal-Piñeiro, Didac and Wang, Yunpeng},
title = {{Cortical thickness changes precede high levels of amyloid by at least 7 years}},
journal = {Nature neuroscience},
year = {2026},
month = aug,
volume = {29},
number = {9},
pages = {2164--2175},
publisher = {Nature Portfolio},
issn = {1097-6256},
doi = {10.1038/s41593-026-02363-4},
url = {https://doi.org/10.1038/s41593-026-02363-4},
pmid = {42618753},
pmcid = {PMC13533846}
}

RIS

TY - JOUR
AU - Roe, James M.
AU - Jagust, William J.
AU - Landau, Susan M.
AU - Harrison, Theresa M.
AU - Grydeland, Håkon
AU - Slivka, Maksim
AU - Alatorre-Warren, José-Luis
AU - Garrido, Pablo F.
AU - Sørensen, Øystein
AU - Grødem, Edvard O. S.
AU - Ward, Tyler J.
AU - Leonardsen, Esten H.
AU - Murphy, Alice
AU - Lee, JiaQie
AU - Fladby, Tormod
AU - Bjørnerud, Atle
AU - Walhovd, Kristine B.
AU - Fjell, Anders M.
AU - Vidal-Piñeiro, Didac
AU - Wang, Yunpeng
TI - Cortical thickness changes precede high levels of amyloid by at least 7 years
T2 - Nature neuroscience
J2 - Nat Neurosci
PY - 2026
DA - 2026/08/19
VL - 29
IS - 9
SP - 2164
EP - 2175
SN - 1097-6256
PB - Nature Portfolio
DO - 10.1038/s41593-026-02363-4
UR - https://doi.org/10.1038/s41593-026-02363-4
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41593-026-02363-4",
"type": "article-journal",
"title": "Cortical thickness changes precede high levels of amyloid by at least 7 years",
"container-title": "Nature neuroscience",
"author": [
{
"family": "Roe",
"given": "James M."
},
{
"family": "Jagust",
"given": "William J."
},
{
"family": "Landau",
"given": "Susan M."
},
{
"family": "Harrison",
"given": "Theresa M."
},
{
"family": "Grydeland",
"given": "Håkon"
},
{
"family": "Slivka",
"given": "Maksim"
},
{
"family": "Alatorre-Warren",
"given": "José-Luis"
},
{
"family": "Garrido",
"given": "Pablo F."
},
{
"family": "Sørensen",
"given": "Øystein"
},
{
"family": "Grødem",
"given": "Edvard O. S."
},
{
"family": "Ward",
"given": "Tyler J."
},
{
"family": "Leonardsen",
"given": "Esten H."
},
{
"family": "Murphy",
"given": "Alice"
},
{
"family": "Lee",
"given": "JiaQie"
},
{
"family": "Fladby",
"given": "Tormod"
},
{
"family": "Bjørnerud",
"given": "Atle"
},
{
"family": "Walhovd",
"given": "Kristine B."
},
{
"family": "Fjell",
"given": "Anders M."
},
{
"family": "Vidal-Piñeiro",
"given": "Didac"
},
{
"family": "Wang",
"given": "Yunpeng"
}
],
"container-title-short": "Nat Neurosci",
"volume": "29",
"issue": "9",
"page": "2164-2175",
"DOI": "10.1038/s41593-026-02363-4",
"PMID": "42618753",
"PMCID": "PMC13533846",
"ISSN": "1097-6256",
"publisher": "Nature Portfolio",
"URL": "https://doi.org/10.1038/s41593-026-02363-4",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
19
]
]
}
}

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.1162/imag.a.1242 [code]
Stable individual differences dominate adult brain volume variation until later life.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: structural MRI / diffusion, 2 references, 4 authors
[2] doi:10.1162/imag.a.1196 [code]
Failure to detect entorhinal grid-like signals in a passive navigation human fMRI study.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: lmerTest, tidyverse, 3 authors
[3] doi:10.1016/j.tjpad.2026.100622 [code]
The time interval from amyloid to tau PET positivity varies by age, sex and APOE-ε4 status.
Journal: The journal of prevention of Alzheimer's disease
In common: broom, lmerTest, lme4, 3 other tools, PET / SPECT, Alzheimer's / dementia, 5 references
[4] doi:10.1002/alz.71609
Cortical gray-white matter contrast alterations precede amyloid-β positivity and macrostructural changes in older adults without dementia.
Journal: Alzheimer's & dementia : the journal of the Alzheimer's Association
In common: PET / SPECT, Alzheimer's / dementia, structural MRI / diffusion, 9 references
[5] doi:10.1002/alz.71567 [code]
Associations of dementia polyexposure scores to Alzheimer's disease endophenotypes in a diverse population.
Journal: Alzheimer's & dementia : the journal of the Alzheimer's Association
In common: easystats, broom, cowplot, 4 other tools, PET / SPECT, Alzheimer's / dementia, 4 references
[6] doi:10.1038/s41398-026-04010-9 [code]
Bullying victimization and brain development: a longitudinal structural magnetic resonance imaging study from adolescence to early adulthood.
Journal: Translational psychiatry
In common: ggseg, easystats, broom, 7 other tools, structural MRI / diffusion
[7] doi:10.1093/brain/awaf413 [code]
Estimating the time course of biomarker changes in Alzheimer's disease.
Journal: Brain : a journal of neurology
In common: PET / SPECT, Alzheimer's / dementia, structural MRI / diffusion, 9 references
[8] doi:10.1093/braincomms/fcag176 [code]
Tau topography subtypes account for clinical heterogeneity and longitudinal trajectories in early-onset Alzheimer's disease.
Journal: Brain communications
In common: easystats, lmerTest, lme4, 4 other tools, PET / SPECT, Alzheimer's / dementia, 3 references
[9] doi:10.1162/imag.a.1235 [code]
Intracranial volume: To adjust or not to adjust? It is not a matter of if, but how.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: ggseg, broom, lmerTest, 6 other tools, structural MRI / diffusion, 1 reference
[10] doi:10.1002/hbm.70605 [code]
BrainEnrich: Revealing Biological Insights for Imaging-Derived Phenotypes Through Transcriptomic Enrichment.
Journal: Human brain mapping
In common: ggseg, easystats, broom, 6 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.