OSCR

ComBat-Predict Enhances Generalizability of Neuroimaging Models to New Sites.

A correction to this paper has been published: the notice, 42670026, from Europe PMC.

Code ↔ Paper

5 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 5 matches
  1. [1] § Methods › Generalizability of Normative Models to a New Site ↔ ComBat-Predict_section1-5_share.R, lines 812–868 · score 0.67 · Rank biserial correlation, Wilcoxon rank sum, centile score, LBCC, ADNI, harmonization
  2. [2] § Results › Harmonization Aligns Normative Scores Between Reference and Test Cohorts ↔ ComBat-Predict_section1-5_share.R, lines 616–657 · score 0.64 · harmonized ADNI CN, unharmonized ADNI, LBCC control, centile scores, LMCI, healthy
  3. [3] § Results › CB‐Predict Harmonization Reduces Heterogeneity Between Healthy Control Groups ↔ ComBat-Predict_section1-5_share.R, lines 812–868 · score 0.62 · rank biserial correlation, Wilcoxon rank sum, centile score, LBCC, ADNI, Predict
  4. [4] § Methods › Generalizability of Normative Models to a New Site ↔ ComBat-Predict_section1-5_share.R, lines 487–549 · score 0.60 · residual standard deviations, centile score, spline, sex, ComBat, formulation
  5. [5] § Methods › Generalizability of Normative Models to a New Site ↔ ComBat-Predict_section1-5_share.R, lines 400–485 · score 0.57 · Normative centile scores, healthy control, CB Predict, ComBat, trained, LBCC

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 · 996 lines · 37 KB · no license · 5 matches

  1. ##################################################################
  2. ####project: ComBat-predict for out-of-sample harmonization and ##
  3. #### improved normative scoring application ##
  4. ####author: Yao Xin ##
  5. ####Created: 08/01/2025 ##
  6. ##################################################################
  9. # load libraries
  10. #analysis
  11. library(ComBatFamily) #ComBatFamily package
  12. library(gam)
  13. library(gamlss)
  14. library(mgcv)
  15. #data cleaning
  16. library(dplyr)
  17. library(tidyr)
  18. library(tidyverse)
  19. #visualization
  20. library(ggplot2)
  21. library(patchwork) #plotting multiple plots together
  22. library(RColorBrewer) #color palettes
  24. ## section 1 data preparation and visualization ####################
  26. # import cleaned data here - our dataset is called adni.lite
  28. ### 1.1 visualization of site effects ########
  29. #### 1.1.1 arrange the sites by mean thickness ####
  30. ordered_sites <- adni.lite %>%
  31. group_by(site) %>%
  32. dplyr::summarize(mean_thickness = mean(thickness.left.caudal.anterior.cingulate)) %>%
  33. arrange(mean_thickness) %>%
  34. pull(site)
  36. adni.lite$site <- factor(adni.lite$site, levels = ordered_sites)
  38. # Create the boxplot in the order of sites, arranged by mean
  39. ggplot(adni.lite,
  40. aes(x = factor(site, levels = ordered_sites), y = thickness.left.caudal.anterior.cingulate,
  41. fill = as.factor(site), color = as.factor(site))) +
  42. geom_boxplot() +
  43. labs(title = "Boxplot of Cortical Thickness Outcome by Site - left caudal anterior cingulate",
  44. x = "Site",
  45. y = "Cortical Thickness",
  46. fill = "Site",
  47. color = "Site") +
  48. theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotate x-axis labels for better readability
  50. #### 1.1.2 arrange the sites by variance of thickness #####
  51. ordered_sites <- adni.lite %>%
  52. dplyr::group_by(site) %>%
  53. summarize(mean_thickness = sd(thickness.left.caudal.anterior.cingulate)) %>%
  54. arrange(mean_thickness) %>%
  55. pull(site)
  57. adni.lite$site <- factor(adni.lite$site, levels = ordered_sites)
  59. # Create the boxplot in the order of sites, arranged by variance
  60. ggplot(adni.lite,
  61. aes(x = factor(site, levels = ordered_sites), y = thickness.left.caudal.anterior.cingulate,
  62. fill = as.factor(site), color = as.factor(site))) +
  63. geom_boxplot() +
  64. labs(title = "Boxplot of Cortical Thickness Outcome by Site - left caudal anterior cingulate",
  65. x = "Site",
  66. y = "Cortical Thickness",
  67. fill = "Site",
  68. color = "Site") +
  69. theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotate x-axis labels for better readability
  72. #### 1.1.3 arrange the sites by **sample size** ####
  73. ordered_sites <- adni.lite %>%
  74. group_by(site) %>%
  75. summarise(freq = n()) %>%
  76. arrange(freq) %>%
  77. pull(site)
  79. adni.lite$site <- factor(adni.lite$site, levels = ordered_sites)
  81. # Create the boxplot in the order of sites
  82. ggplot(adni.lite,
  83. aes(x = factor(site, levels = ordered_sites), y = thickness.left.caudal.anterior.cingulate,
  84. fill = as.factor(site), color = as.factor(site))) +
  85. geom_boxplot() +
  86. labs(title = "Boxplot of Cortical Thickness Outcome by Site - left caudal anterior cingulate",
  87. x = "Site",
  88. y = "Cortical Thickness",
  89. fill = "Site",
  90. color = "Site") +
  91. theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotate x-axis labels for better readability
  93. #### print the sample sizes
  94. adni.lite %>%
  95. group_by(site) %>%
  96. summarise(freq = n()) %>%
  97. arrange(freq) %>% print(n=64)
  101. ## section 2 - Out-of-sample harmonization performance####################
  103. ### 2.1 full sample harmonization - for in-sample prediction ##########
  104. # Fit comfam on the full sample - all batches
  105. fit2 <- comfam(data = adni.lite[,9:70], #data
  106. bat = as.factor(adni.lite$batch), #batch
  107. covar=adni.lite[,(3:5)], # covariates
  108. lm, #model
  109. y ~ AGE + as.factor(SEX) + as.factor(DIAGNOSIS), #formula
  110. eb = TRUE, #empirical bayes,
  111. robust.LS = FALSE,
  112. ref.batch = 64
  113. )
  115. #apply the comfam fit on all batches, get outcomes of in-sample harmonization
  116. fit2_pred_in <- predict(fit2,
  117. newdata = adni.lite[,9:70],
  118. newbat = adni.lite$batch,
  119. newcovar = adni.lite[,(3:5)],
  120. robust.LS = FALSE,
  121. eb=TRUE
  122. )
  124. # in-sample-harmonized version of the data (for all 1~64 batches)
  125. full_sample_fitted <- cbind(adni.lite$ID,adni.lite$batch,
  126. as.data.frame(fit2_pred_in$dat.combat ) )
  127. colnames(full_sample_fitted)[1:2] <- c("ID", "batch")
  129. # the resulting harmonzied data frame is full_sample_fitted, with
  130. # 62 columns for outcomes. 505 rows for subjects
  133. ### 2.2 Evaluating out-of-sample prediction by leaving one site out #######
  135. #create empty objects to store results
  136. Nbatch <- 64
  137. n.train <- n.test <- rep(NA, Nbatch)
  138. MSE_batch <- error_list <- pred_list <- vector(mode='list', length=Nbatch)
  140. for (i in 1:Nbatch) {
  141. # split data to training and testing :
  142. # train with all batches except batch i , test with batch i
  143. train_data <- adni.lite[adni.lite$batch %in% seq(1:Nbatch)[-i],]
  144. test_data <- adni.lite[adni.lite$batch==i,] #batch i is being tested
  145. row_index <- which(adni.lite$batch==i)
  146. n.train[i] <- nrow(train_data)
  147. n.test[i] <- nrow(test_data)
  149. # fit ComBat model on training data
  150. fit_train <- comfam(data = train_data[,9:70], #cortical thickness columns
  151. bat = as.factor(train_data$batch), #batch
  152. covar=train_data[,(3:5)], # covariates: AGE, SEX, DIAGNOSIS
  153. lm, # type of model
  154. y ~ AGE + as.factor(SEX) + as.factor(DIAGNOSIS), #formula
  155. eb = TRUE,
  156. robust.LS = FALSE,
  157. ref.batch = max(as.numeric(train_data$batch)) #use the largest batch as ref
  158. )
  160. # apply the comfam fit on the testing data (batch i)
  161. out_pred <- predict(fit_train,
  162. newdata=test_data[,9:70], #new data
  163. newbat=as.factor(test_data$batch), #new batch
  164. newcovar = test_data[,(3:5)],
  165. robust.LS = FALSE, eb = TRUE
  166. )
  168. # store the out-of-sample harmonized outcomes for batch i in the i_ith element of pred_list
  169. pred_list[[i]] <- out_pred$dat.combat
  172. #error for each subject in batch i: store in the i_th element of error_list
  173. error_list[[i]] <- as.matrix(pred_list[[i]]) - # out-of-sample harmonized outcomes
  174. as.matrix(full_sample_fitted[row_index,-c(1:2)] ) # subtracting: in-sample harmonized outcomes
  176. MSE_batch[[i]] <- colMeans(error_list[[i]]^2) #MSE for all subjects in batch i
  178. }
  180. result <- list(MSE_batch = MSE_batch,
  181. error_list = error_list,
  182. pred_list = pred_list)
  183. #save(result, file = "path/result_name.RData")
  186. ### 2.3 Examining&visualizing the difference of out-of-sample and in-sample -##########
  187. ##load in saved result
  188. # load("path/xxx.RData")
  189. Nbatch <- 64
  190. MSE_batch <- result$MSE_batch
  191. error_list <- result$error_list
  192. pred_list <- result$pred_list
  194. ### re-organize MSE into new data frame, and derive RMSE
  195. batch_index <- c()
  196. MSE_vec <- c()
  198. for (j in 1:Nbatch) {
  199. MSE_vec <- rbind(MSE_vec, as.data.frame(as.numeric(MSE_batch[[j]])) )
  200. batch_index <- c(batch_index, rep(j, length(MSE_batch[[j]]) ) )
  201. }
  203. MSE_df <- as.data.frame(cbind(batch_index, MSE_vec))
  204. colnames(MSE_df) <- c("batch", "MSE")
  206. MSE_df <- MSE_df %>% mutate(RMSE = sqrt(MSE))
  209. #### how low are the RMSE values?
  210. summary(MSE_df$RMSE)
  212. #### [Figure 1a] Plot RMSE in each batch, arranged by batch sample size #####
  213. batch_table_sub1 <- batch_table_sub %>%
  214. mutate(batch = batch_order ) %>%
  215. dplyr::select(batch, freq)
  217. MSE_df <- MSE_df %>%
  218. left_join(batch_table_sub1, by = "batch") %>%
  219. mutate(freq=as.factor(freq) ) %>%
  220. arrange(freq) #merge in sample size of each batch
  222. #different sample sizes - categorized for plot
  223. category_breaks <- c(1, 9, 16, 23, 29,
  224. 37, 42, 48, 52, 54,
  225. 57, 59, 61, 62, 63)-0.5
  226. category_labels <- c("3", "4", "5", "6", "7",
  227. "8", "9", "10", "11","12",
  228. "13", "14", "17", "20", "23")
  230. # set color palette
  231. mycolors <- colorRampPalette(brewer.pal(8, "Set2"))(15)
  233. # set up a uniform theme for publication
  234. publication_theme <- function(base_size = 10, base_family = "Liberation Sans") {
  235. theme_minimal(base_size = base_size, base_family = base_family) +
  236. theme(
  237. axis.text = element_text(size = base_size, family = base_family),
  238. axis.title = element_text(size = base_size + 2, family = base_family),
  239. plot.title = element_text(size = base_size + 4, family = base_family, hjust = 0.5),
  240. plot.subtitle = element_text(size = base_size + 2, family = base_family),
  241. legend.text = element_text(size = base_size),
  242. legend.title = element_text(size = base_size + 1),
  243. axis.text.x = element_text(angle = 45, hjust = 1), # For your specific case
  244. panel.grid.minor = element_blank(), # Clean look for publication
  245. plot.margin = margin(10, 10, 10, 10) # Consistent margins
  246. )
  247. }
  250. ggplot(MSE_df[MSE_df$batch != 64,],
  251. aes(x = factor(batch), y = (RMSE),
  252. fill = as.factor(freq), color = as.factor(freq) )) +
  253. geom_boxplot() +
  254. stat_summary(fun = median, geom = "crossbar",
  255. width = 1, fatten = 0.4,
  256. color = "white", size = 1) +
  257. labs(
  258. subtitle = "a)",
  259. x = "Sites - ordered by sample size",
  260. y = "RMSE (mm)",
  261. fill = "Sample Size",
  262. color = "Sample Size") +
  263. publication_theme() +
  264. theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7),
  265. plot.subtitle = element_text(hjust = 0)) + # Add this line # Add dashed lines between the categories
  266. geom_vline(xintercept = category_breaks,
  267. linetype = "dashed", color = "grey60") +
  268. scale_fill_manual(values = mycolors,
  269. labels = category_labels) +
  270. scale_color_manual(values = mycolors)
  272. #### [Figure 1b] RMSE by AD% - plots #########
  274. # arrange the sites by **diagnosis%**
  276. ordered_batch <- adni.lite %>%
  277. group_by(batch) %>%
  278. summarise(ADprop = mean(DIAGNOSIS == "AD")) %>%
  279. arrange(ADprop) %>%
  280. pull(batch)
  282. MSE_df$batch <- factor(MSE_df$batch, levels = ordered_batch)
  284. # Create the boxplot in the order of batches by AD%
  285. ggplot(MSE_df[MSE_df$batch != 64,],
  286. aes(x = factor(batch, levels = ordered_batch),
  287. y = RMSE,
  288. fill = as.factor(batch), color = as.factor(batch))) +
  289. geom_boxplot() +
  290. stat_summary(fun = median, geom = "crossbar",
  291. width = 1, fatten = 0.4,
  292. color = "white", size = 1) +
  293. labs(#title = "Boxplot of RMSE of CB-Predict vs ComBat/d-Combat",
  294. #subtitle = "Averaged across subjects in each site",
  295. subtitle="b)",
  296. x = "Sites - ordered by AD%",
  297. y = "RMSE (mm)",
  298. fill = "Batch (ordered by AD%, 0.00% - 66.7%)",
  299. color = "Batch (ordered by AD%, 0.00% - 66.7%)") +
  300. publication_theme() +
  301. theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7) , # smaller font on xaxis labels
  302. legend.position = "none",
  303. plot.subtitle = element_text(hjust = 0))
  307. ## section 3 - Out-of-sample prediction performance####################
  308. #comparing the errors of linear model prediction, from non-harmonized vs. harmonized
  310. ### 3.1 Residuals from linear model prediction, with and without ComBat harmonization ##########
  311. Nbatch <- 64
  312. n.train <- n.test <- rep(NA, Nbatch)
  313. pred_list <- error_old_list <- error_new_list <- vector(mode='list', length=Nbatch)
  314. MSE_new_batch<- MSE_old_batch <- vector(mode='list', length=Nbatch)
  316. for (i in 1:Nbatch) {
  317. # similar to 2.2, split data to training and testing by leaving site i out
  318. train_data <- adni.lite[adni.lite$batch!=i,]
  319. test_data <- adni.lite[adni.lite$batch==i,] #batch i is being tested
  320. row_index <- which(adni.lite$batch==i)
  321. n.train[i] <- nrow(train_data)
  322. n.test[i] <- nrow(test_data)
  324. # fit ComBat model on training data
  325. fit_train <- comfam(data = train_data[,9:70], # CT outcomes from 62 regions
  326. bat = as.factor(train_data$batch), #sites
  327. covar=train_data[,(3:5)], # covariates
  328. lm, #model
  329. y ~ AGE + as.factor(SEX) + as.factor(DIAGNOSIS), #formula
  330. eb = TRUE,
  331. robust.LS = FALSE,
  332. ref.batch = max(as.numeric(train_data$batch)) #use the largest batch as ref
  333. )
  335. # harmonized version of the training data
  336. train_data_post <- as.data.frame(cbind(train_data[,1:8],fit_train$dat.combat) )
  338. # Multivariate linear model
  339. fit_LM_train <- lm( fit_train$dat.combat ~
  340. AGE + as.factor(SEX) + as.factor(DIAGNOSIS),
  341. data=train_data_post )
  343. # Applying the linear model - predicted CTs of testing data
  344. fit_LM_pred <- predict(fit_LM_train,
  345. newdata=test_data)
  347. # out-of-sample harmonization of the testing data
  348. fit_test <- predict(
  349. fit_train, # using the fitted ComBat model from training data
  350. newdata = test_data[,9:70], # harmonizing based on observed CTs in testing data
  351. newbat = as.factor(test_data$batch), #sites in the testing
  352. newcovar = test_data[,(3:5)],
  353. robust.LS = FALSE,
  354. eb=TRUE
  355. )
  357. test_data_harmonized <- fit_test$dat.combat # harmonized CTs in testing
  359. ## prediction error _ harmonized vs predicted
  360. error_new_list[[i]] <- as.matrix(fit_LM_pred) - as.matrix(test_data_harmonized) # for all subjects
  361. MSE_new_batch[[i]] <- colMeans(error_new_list[[i]]^2) #MSE for all subjects in batch i
  363. ### prediction error _ non-harmonized vs predicted
  364. error_old_list[[i]] <- as.matrix(fit_LM_pred) - as.matrix(test_data[,c(9:70)] )# for all subjects
  365. MSE_old_batch[[i]] <- colMeans(error_old_list[[i]]^2) #MSE for all subjects in batch i
  367. }
  369. result2 <- list(MSE_new_batch=MSE_new_batch,
  370. MSE_old_batch = MSE_old_batch,
  371. error_new_list = error_new_list,
  372. error_old_list = error_old_list,
  373. pred_list = pred_list)
  374. save(result2, file = "path/result2.Rdata")
  377. ### 3.2 Compare prediction errors of linear models for ComBat-Pred vs. non-harmonized #####
  378. ##load in elements in result2
  379. load("path/result2.Rdata")
  380. MSE_new_batch <- result2$MSE_new_batch
  381. MSE_old_batch <- result2$MSE_old_batch
  382. error_new_list <- result2$error_new_list
  383. error_old_list <- result2$error_old_list
  384. pred_list <- result2$pred_list
  386. ### re-organize the two sets of MSEs into new data frame
  387. MSE_new_vec<- MSE_old_vec <- c()
  388. batch_index <- c()
  389. Nbatch <-64
  391. for (j in 1:Nbatch) {
  392. MSE_new_vec <- rbind(MSE_new_vec, as.data.frame(as.numeric(MSE_new_batch[[j]])) )
  393. MSE_old_vec <- rbind(MSE_old_vec, as.data.frame(as.numeric(MSE_old_batch[[j]])) )
  394. batch_index <- c(batch_index, rep(j, length(MSE_new_batch[[j]]) ) )
  395. }
  397. MSE_oldnew_df <- as.data.frame(cbind(batch_index, MSE_new_vec, MSE_old_vec))
  398. colnames(MSE_oldnew_df) <- c("batch", "MSE_new","MSE_old")
  400. # add RMSE column for each set of MSE
  401. MSE_oldnew_df <- MSE_oldnew_df %>%
  402. mutate(difference = MSE_old - MSE_new,
  403. RMSE_new = sqrt(MSE_new),
  404. RMSE_old = sqrt(MSE_old)
  405. )
  408. #### [Figure 2] boxplot of new and old MSE by batch ########
  410. ggplot(MSE_oldnew_df[MSE_oldnew_df$batch != 64,],
  411. aes(x = factor(batch))) +
  412. geom_boxplot(aes(y = RMSE_new, fill = "Harmonized"),color = "green4", alpha = 0.8) + # New MSEs in light blue
  413. geom_boxplot(aes(y = RMSE_old, fill = "Non-harmonized"),color = "grey", alpha = 0.4) + # Old MSEs in grey
  414. scale_fill_manual(values = c("Harmonized" = "green4", "Non-harmonized" = "grey")) + # Define fill colors
  415. #scale_color_manual(values = c(ComBat = "purple", LM = "grey")) + # Define border colors
  416. labs(#title = "Boxplot of Out-of-Sample Prediction RMSE with vs. without CB-Predict Harmonization",
  417. #subtitle = "Averaged across subjects in each batch",
  418. x = "Sites - ordered by sample size",
  419. y = "RMSE (mm)",
  420. fill = "Model",
  421. color = "Model") +
  422. publication_theme() + # use the same publication theme as defined in fig1.1
  423. theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 6) )
  431. ## Section 4 - normative centile scores ###########
  432. #with healthy subjects in training, generate normative scores for healthy subjects in testing
  434. #this section is substituted by section 5
  436. ## Section 5 - centile scores, with LBCC healthy dataset as reference for normative scoring ###########
  438. ### 5.1 prepare the reference data for analysis --------
  440. # read in cleaned larger CN reference data - our imported dataset is called CN_lite
  442. # obtain colnames of regions
  443. CT_names1 <- colnames(CN_lite)[6:67]
  445. #### 5.1.1 reorder ADNI1 lite data in the same region order with CN_lite #########
  447. ## obtain left and right region names in ADNI
  448. ref_names <- colnames(adni.lite)
  449. ref_names_left <- ref_names[grepl("thickness.left", ref_names)]
  450. ref_names_right <- ref_names[grepl("thickness.right", ref_names)]
  452. #order alphabetically
  453. ref_names_left <- ref_names_left[order(ref_names_left)]
  454. ref_names_right <- ref_names_right[order(ref_names_right)]
  455. ref_names <- c(ref_names_left, ref_names_right)
  456. # compare with CN data ordering
  457. print(cbind(CT_names1,ref_names )) # manually check
  459. ## reorder columns in adni.lite by alphabetical order of regions
  460. adni.lite.reorder <- adni.lite %>%
  461. select(1:8, all_of(ref_names))
  463. #### 5.1.2 subset eligible batches in reordered ADNI lite #####
  465. # select out eligible batches: with >= 3 healthy controls
  466. batches_CN <- adni.lite.reorder %>%
  467. filter(DIAGNOSIS=="CN") %>% # filter out
  468. group_by(batch) %>%
  469. summarise(freq=n()) %>%
  470. filter(freq >= 3) %>%
  471. ungroup()
  473. # testing_batch_IDs contains the index of eligible batches - with at least 3 healthy controls
  474. testing_batch_IDs<-batches_CN$batch
  476. # subset of all subjects in the eligible batches
  477. adni.lite.reorder.sub <- adni.lite.reorder %>%
  478. filter( batch %in% testing_batch_IDs)
  480. # subset the healthy controls in the eligible batches
  481. adni.lite.reorder.CN_sub <- adni.lite.reorder %>%
  482. filter(DIAGNOSIS=="CN"& batch %in% testing_batch_IDs)
  485. ### 5.2 harmonizing the CN_lite healthy control data ######
  487. # function for harmonizing data and converting it into centile scores
  488. harmonize_qscore_4CN_lite <- function(dataset){
  490. #harmonizing the healthy subjects
  491. fit_harmonize <- comfam(data = dataset[,6:67], # data columns of outcomes
  492. bat = as.factor(dataset$batch), #batch
  493. covar=dataset[,(2:3)], # covariates: age and sex
  494. model = lm, #model
  495. formula = y ~ AGE + as.factor(SEX), #formula
  496. eb = TRUE,
  497. robust.LS = FALSE,
  498. ref.batch = max(dataset$batch)
  499. )
  501. #harmonized version of dataset
  502. df_post <- as.data.frame(cbind(dataset[,1:5],
  503. fit_harmonize$dat.combat) )
  505. #create empty matrix for storing centile score ver. of the same data
  506. df_pscores <- df_post
  507. df_pscores[,] <- NA
  508. #create empty list for storing each lm fit and sigma_hat
  509. fit_LM_post_list <- list()
  510. sigma2_hat_list <- list()
  512. message("Fitting individual models and calculating p-scores...")
  514. for (k in 6:ncol(df_post) ) { # for column 6 to the last col
  515. col_name <- names(df_post)[k]
  517. message(paste("Processing column:", col_name, "(", k, "/", ncol(df_post), ")"))
  519. # fit individual lm <--- can include a spline term, in here
  520. fit_LM_post_k <- lm(as.matrix(df_post[,k]) ~
  521. AGE + as.factor(SEX), data=df_post)
  523. #----produce centile scores based on est. distribution of yhat -----
  525. # variance of Yhat|X
  526. sigma_hat_k <- summary(fit_LM_post_k)$sigma
  528. # Calculate the centile (p-score) for each observation in column k
  529. # using the fitted mean and the residual standard deviation for model k
  530. df_pscores[, k] <- pnorm(df_post[, k],
  531. mean = fitted(fit_LM_post_k),
  532. sd = sigma_hat_k)
  534. # Store the model fit and sigma_hat in their respective lists
  535. # Use the column name as the list element name for easy access
  536. fit_LM_post_list[[col_name]] <- fit_LM_post_k
  537. sigma2_hat_list[[col_name]] <- sigma_hat_k
  539. }
  541. message("Function finished.")
  542. return(list(
  543. df_post = df_post, # The harmonized data
  544. df_pscores = df_pscores, # df of percentile scores
  545. fit_harmonize = fit_harmonize, # The original object returned by comfam
  546. fit_LM_post_list = fit_LM_post_list, # List containing LM for each region
  547. sigma2_hat_list = sigma2_hat_list # List containing each sigma_hat
  548. ))
  549. }
  552. # Harmonize the CN_lite data, AND store the results
  553. pscore_CN <- harmonize_qscore_4CN_lite(CN_lite)
  555. # save result, load in result for use later
  556. # save(pscore_CN,file="yourpath/pscore_CN.Rdata")
  558. # load("yourpath/pscore_CN.Rdata")
  560. #load in elements in pscore_CN
  561. CN_fit <- pscore_CN$fit_harmonize #the comfam model fit
  562. CN_scores <- pscore_CN$df_pscores #centile scores
  563. CN_LM <- pscore_CN$fit_LM_post_list # List containing LM for each region
  564. CN_sigma2 <- pscore_CN$sigma2_hat_list # List containing each sigma_hat
  567. ### 5.3 Applying the model to harmonize+normalize adni.lite data ######
  569. #### 5.3.1 The function for applying the CN model to harmonized/unharmonized test data #####
  571. test_data_qscore_fn <- function(test_data_CN, test_data){
  573. # Harmonizing out-of-sample testing data, with the comfam fit based on CN_lite
  574. test_harmonized_applied <- predict(CN_fit,
  575. newdata = test_data[,9:70],
  576. newbat = as.factor(test_data$batch),
  577. newcovar = test_data[,3:4],
  578. robust.LS = FALSE,
  579. eb=TRUE )
  581. # harmonized ver. of testing data
  582. test_raw_values <- test_harmonized_applied$dat.combat # test data harmonized
  583. test_raw_values2 <- as.matrix(test_data[,9:70]) # test data un-harmonized(Raw)
  585. #empty dfs for storing centile score ver. of the testing data
  586. test_pscores_harmonized <-
  587. test_pscores_unharmonized <-
  588. matrix(NA, nrow = nrow(test_data), ncol = 62)
  590. # centile scores from harmonized data
  591. for (j in 1:62) {
  592. test_means <- predict(CN_LM[[j]], # the LM for region j
  593. newdata = test_data[,3:4]) # fitted mean for testing data
  594. test_sigma <- CN_sigma2[[j]] # the fitted sigma_hat for region j
  596. #obtain the p-score of the
  597. test_pscores_harmonized[,j] <- pnorm(test_raw_values[,j],
  598. mean = test_means,
  599. sd = test_sigma )
  600. }
  603. # centile scores from NON-harmonized data
  604. for (j in 1:62) {
  605. test_pscores_unharmonized[,j] <- pnorm(test_raw_values2[,j],
  606. mean = test_means,
  607. sd = test_sigma)
  608. }
  610. # output both datasets - harmonized centile scores, non-harmonized centile scores
  611. return(list(test_pscores_harmonized = test_pscores_harmonized,
  612. test_pscores_unharmonized = test_pscores_unharmonized)
  613. )
  614. }
  616. #### 5.3.2 Applying the function on the ADNI data (as testing) ######
  618. # dataset of harmonized centile scores of our testing data
  619. test_pscores_harmonized <-
  620. test_data_qscore_fn(test_data_CN = adni.lite.reorder.CN_sub,
  621. # Using the healthy subset in testing to estimate site effects
  622. test_data = adni.lite.reorder.sub
  623. # To harmonize all subjects in testing data
  624. )$test_pscores_harmonized # the harmonized centile scores
  626. # dataset of NON-harmonized centile scores of our testing data
  627. test_pscores_unharmonized <-
  628. test_data_qscore_fn(test_data_CN = adni.lite.reorder.CN_sub,
  629. test_data = adni.lite.reorder.sub
  630. )$test_pscores_unharmonized # UNharmonized centile scores
  634. #### 5.3.3 Comparing the centile scores in CN_lite, and in ADNI.lite.sub #######
  636. ##### 5.3.3.1 save the centile scores of the four groups - LBCC_control, CN, LMCI, AD ########
  638. # dataset for plotting LBCC vs. three groups in testing data - harmonized
  639. levels_diagnosis <- c( "LBCC_control","CN","LMCI", "AD" ) # add LBCC control as one of the levels
  641. # harmonized ADNI centile scores
  642. centile_scores_3groups <-
  643. cbind(adni.lite.reorder.sub[,1:5] ,
  644. test_pscores_harmonized) %>%
  645. as.data.frame() %>%
  646. mutate(DIAGNOSIS = factor(DIAGNOSIS,
  647. levels = c("LBCC_control", "CN", "LMCI", "AD"),
  648. labels = c("LBCC(control)", "ADNI - CN", "ADNI - LMCI", "ADNI - AD")))
  650. # unharmonized ADNI centile scores
  651. centile_scores_3groups_raw <-
  652. cbind(adni.lite.reorder.sub[,1:5] ,
  653. test_pscores_unharmonized) %>%
  654. as.data.frame() %>%
  655. mutate(DIAGNOSIS = factor(DIAGNOSIS,
  656. levels = c("LBCC_control", "CN", "LMCI", "AD"),
  657. labels = c("LBCC(control)", "ADNI - CN", "ADNI - LMCI", "ADNI - AD")))
  659. # LBCC centile scores
  660. CN_scores_1 <- CN_scores %>%
  661. as.data.frame() %>%
  662. mutate( DIAGNOSIS = "LBCC_control",
  663. DIAGNOSIS = factor(DIAGNOSIS, #add the levels
  664. levels = c("LBCC_control", "CN", "LMCI", "AD"),
  665. labels = c("LBCC(control)", "ADNI - CN", "ADNI - LMCI", "ADNI - AD")),
  666. batch=factor(batch)) %>%
  667. rename(ID = participant) %>%
  668. select(ID, batch, AGE, SEX, DIAGNOSIS, 6:67)
  671. # append the ADNI + LBCC centile scores
  673. colnames(centile_scores_3groups) <- colnames(CN_scores_1)
  674. colnames(centile_scores_3groups_raw) <- colnames(CN_scores_1)
  675. centile_scores_4groups <- dplyr::bind_rows(centile_scores_3groups, CN_scores_1) # harmonized
  676. centile_scores_4groups_raw <- dplyr::bind_rows(centile_scores_3groups_raw, CN_scores_1) # unharmonized
  679. ##### 5.3.3.2. plot centile scores in 1 group of LBCC + 3 groups of testing data======== #####
  680. #### [Figure 4a-b] ######
  681. # Define color palette for consistency
  682. mycolors <- c("grey65", "#4059ad", "#97d8c4", "#f4b942")
  684. ## figure 4A - boxplot of centile scores in 4 groups, without harmonization
  685. plot_a <- centile_scores_4groups_raw %>%
  686. # pivot by region for plotting each region in one box
  687. tidyr::pivot_longer(cols = 6:67, names_to = "region", values_to = "centile") %>%
  688. arrange(region) %>%
  689. ggplot(aes(x = region, y = centile, color = DIAGNOSIS, fill = DIAGNOSIS)) +
  690. geom_boxplot(outlier.size = 0.5, size = 0.3) +
  691. stat_summary(fun = median, geom = "crossbar",
  692. width = 1, fatten = 0.4,
  693. color = "white", size = 1) +
  694. labs(#title=str_wrap("Boxplot of Fitted Centile Scores by Region and Group with and without CB-Predict Harmonization", width = 60),
  695. subtitle = "a) Testing Data Not Harmonized",
  696. x = "Regions 1–62",
  697. y = "Fitted Centile Scores",
  698. color = "Group", # Legend title
  699. fill = "Group") + # Legend title
  700. facet_wrap(~ DIAGNOSIS, ncol = 4) +
  701. theme_minimal(base_size = 10, base_family = "Liberation Sans") +
  702. theme(
  703. axis.text = element_text(size = 10),
  704. axis.title = element_text(size = 12),
  705. plot.title = element_text(size = 16, hjust = 0, vjust=1),
  706. plot.subtitle = element_text(size = 12,hjust=0),
  707. legend.text = element_text(size = 10),
  708. legend.title = element_text(size = 11),
  709. strip.text = element_text(size = 10),
  710. axis.text.x = element_blank(),
  711. panel.grid.minor = element_blank(), # Clean look for publication
  712. plot.margin = margin(10, 10, 10, 10) # Consistent margins
  713. ) +
  714. ylim(0, 1) +
  715. scale_color_manual(values = mycolors) +
  716. scale_fill_manual(values = mycolors) +
  717. geom_hline(yintercept = 0.5, linetype = "longdash", linewidth = 0.8, color = "grey")
  719. ## figure 4B - boxplot of centile scores in 4 groups, with harmonization
  720. plot_b <- centile_scores_4groups %>%
  721. tidyr::pivot_longer(cols = 6:67, names_to = "region", values_to = "centile") %>%
  722. arrange(region) %>%
  723. ggplot(aes(x = region, y = centile, color = DIAGNOSIS, fill = DIAGNOSIS)) +
  724. geom_boxplot(outlier.size = 0.5, size = 0.3) +
  725. stat_summary(fun = median, geom = "crossbar",
  726. width = 1, fatten = 0.4,
  727. color = "white", size = 1) +
  728. labs(#title = "Boxplot of Fitted Centile Scores by Region and Group with and without CB-Predict Harmonization",
  729. subtitle = "b) Testing Data Harmonized",
  730. x = "Regions 1–62",
  731. y = "Fitted Centile Scores",
  732. color = "Group", # Legend title
  733. fill = "Group") + # Legend title
  734. facet_wrap(~ DIAGNOSIS, ncol = 4) +
  735. theme_minimal(base_size = 10, base_family = "Liberation Sans") +
  736. theme(
  737. axis.text = element_text(size = 10),
  738. axis.title = element_text(size = 12),
  739. plot.title = element_text(size = 16, hjust = 0, vjust=1),
  740. plot.subtitle = element_text(size = 12,hjust=0),
  741. legend.text = element_text(size = 10),
  742. legend.title = element_text(size = 11),
  743. strip.text = element_text(size = 10),
  744. axis.text.x = element_blank(),
  745. panel.grid.minor = element_blank(), # Clean look for publication
  746. plot.margin = margin(10, 10, 10, 10) # Consistent margins
  747. ) +
  748. ylim(0, 1) +
  749. scale_color_manual(values = mycolors) +
  750. scale_fill_manual(values = mycolors) +
  751. geom_hline(yintercept = 0.5, linetype = "longdash", linewidth = 0.8, color = "grey")
  753. # Combine plots using patchwork
  754. combined_plot <- plot_a / plot_b +
  755. plot_layout(heights = c(1, 1), guides = "collect") &
  756. theme(legend.position = "bottom")
  758. # Print the combined plot
  759. combined_plot
  762. ### 5.4 Harmonization improves the comparability of normative scores #######
  764. #### 5.4.1 Non-parametric test of the differences between ADNI and LBCC CN groups: ####
  765. # (Wilcoxon rank sum test )
  767. # we are going to compare the CN in ADNI vs. the CN in LBCC
  768. centile_scores_2groups_CN <-
  769. subset(centile_scores_4groups,
  770. DIAGNOSIS %in% c("LBCC(control)", "ADNI - CN"))
  772. ##### 5.4.1.1 function for conducting wilcoxon test and calculating effect size ######
  773. fn_wilcox <- function(datacolumn, datasetname) {
  774. # build formula
  775. form <- as.formula(paste(datacolumn, "~ DIAGNOSIS"))
  777. # run Wilcoxon rank‑sum (two‐sample) test
  778. wtest <- wilcox.test(form,
  779. data = datasetname,
  780. exact = FALSE,
  781. correct = TRUE)
  783. # extract W and p
  784. W <- unname(wtest$statistic)
  785. p_value <- wtest$p.value
  787. # group sizes
  788. grp_tbl <- table(datasetname$DIAGNOSIS)
  789. n1 <- as.numeric(grp_tbl[1])
  790. n2 <- as.numeric(grp_tbl[2])
  791. N <- n1 + n2
  793. # compute Z and r
  794. mu_W <- n1 * n2 / 2
  795. sigma_W <- sqrt(n1 * n2 * (N + 1) / 12)
  796. Z <- (W - mu_W) / sigma_W
  797. r <- Z / sqrt(N)
  799. # Bonferroni‐adjusted significance
  800. sig <- ifelse(p_value < 0.05/62, "*", "")
  802. return(list(
  803. W = round(W, 1),
  804. Z = round(Z, 3),
  805. r = round(r, 3),
  806. p_value = round(p_value, 8),
  807. sig = sig
  808. ))
  809. }
  812. ##### 5.4.1.2 Conduct the wilcoxon test for each region - harmonized, comparing ADNI.CN vs. LBCC.CN ######
  813. # prepare result table
  814. wilcox_table <- data.frame(
  815. region = colnames(centile_scores_2groups_CN)[6:67],
  816. W = NA_real_, #Wilcoxon rank sum test stat
  817. Z = NA_real_, #standardized score
  818. r = NA_real_, #rank-biserial correlation
  819. p_value = NA_real_, #p-value
  820. significance = "" #compared with bonferroni corrected
  821. )
  823. # loop through each column
  824. for (i in seq_along(wilcox_table$region)) {
  825. colname <- wilcox_table$region[i]
  826. res <- fn_wilcox(datacolumn = colname,
  827. datasetname = centile_scores_2groups_CN)
  829. wilcox_table[i, c("W","Z","r","p_value","significance")] <-
  830. unlist(res)
  831. }
  833. # wilcox_table stores the wilcoxon test statistics
  834. #print the test statistics:
  835. # for better visuals, put left and right side to side
  836. # Combine the two halves column-wise
  837. cbind(
  838. wilcox_table[1:31, , drop = FALSE],
  839. wilcox_table[32:62, , drop = FALSE]
  840. ) %>%
  841. knitr::kable("html",
  842. col.names = c("Region-left", "W", "Z", "r", "P-value", "Significance",
  843. "Region-right", "W", "Z", "r", "P-value", "Significance")
  844. )
  847. ##### 5.4.1.3 Conduct the wilcoxon test for each region - NON-harmonized, comparing ADNI.CN vs. LBCC.CN ######
  848. centile_scores_2groups_raw <-
  849. subset(centile_scores_4groups_raw,
  850. DIAGNOSIS %in% c("LBCC(control)", "ADNI - CN"))
  852. wilcox_table_raw <- data.frame(
  853. region = colnames(centile_scores_2groups_raw)[6:67],
  854. W = NA_real_, #Wilcoxon rank sum test stat
  855. Z = NA_real_, #standardized score
  856. r = NA_real_, #rank-biserial correlation
  857. p_value = NA_real_, #p-value
  858. significance = "" #compared with bonferroni corrected
  859. )
  861. for (i in seq_along(wilcox_table_raw$region)) {
  862. colname <- wilcox_table_raw$region[i]
  863. res <- fn_wilcox(datacolumn = colname,
  864. datasetname = centile_scores_2groups_raw)# the un-harmonized CN groups
  866. wilcox_table_raw[i, c("W","Z","r","p_value","significance")] <-
  867. unlist(res)
  868. }
  870. # wilcox_table_raw stores the wilcoxon test statistics
  871. # Print the wilcoxon test statistics:
  872. # for better visuals, put left and right side to side
  873. # Combine the two halves column-wise
  874. cbind(
  875. wilcox_table_raw[1:31, , drop = FALSE],
  876. wilcox_table_raw[32:62, , drop = FALSE]
  877. ) %>%
  878. knitr::kable("html",
  879. col.names = c("Region-left", "W", "Z", "r", "P-value", "Significance",
  880. "Region-right", "W", "Z", "r", "P-value", "Significance")
  881. )
  884. #### 5.4.2 visualization of the 2 sets of wilcoxon tests ######
  886. # Rename columns so harmonized vs unharmonized line up
  887. colnames(wilcox_table) <- c(
  888. "outcome",
  889. "W_harmonized", "Z_harmonized", "r_harmonized",
  890. "p_value_harmonized", "sig_harmonized"
  891. )
  892. colnames(wilcox_table_raw) <- c(
  893. "outcome",
  894. "W_unharmonized", "Z_unharmonized", "r_unharmonized",
  895. "p_value_unharmonized", "sig_unharmonized"
  896. )
  898. # Wide‐format join for printing side‑by‑side HTML table
  899. wilcox_wide <- wilcox_table %>%
  900. left_join(wilcox_table_raw, by = "outcome") %>%
  901. select(
  902. outcome,
  903. Z_unharmonized, r_unharmonized,
  904. Z_harmonized, r_harmonized
  905. )
  907. # Long format of wilcoxon tables for plotting
  908. wt_h <- wilcox_table %>%
  909. select(outcome, Z_harmonized, r_harmonized) %>%
  910. mutate(source = "harmonized") %>%
  911. # rename to match
  912. rename(
  913. Z = Z_harmonized,
  914. r = r_harmonized
  915. )
  917. wt_u <- wilcox_table_raw %>%
  918. select(outcome, Z_unharmonized, r_unharmonized) %>%
  919. mutate(source = "unharmonized") %>%
  920. # rename to match
  921. rename(
  922. Z = Z_unharmonized,
  923. r = r_unharmonized
  924. )
  926. results_long <- bind_rows(wt_h, wt_u)
  929. #### [Figure 5a] lollipop plots of the rank biseral --------
  931. # Define color palette
  932. mycolors <- c("#4059ad", "#f4b942") # For harmonized and unharmonized
  934. # Lollipop chart
  935. results_dodged <- results_long %>%
  936. mutate(
  937. r_abs = abs(as.numeric(r)),
  938. dodge_pos = as.numeric(factor(outcome)) +
  939. ifelse(source == unique(source)[2], 0.2, -0.2) # Manual dodge offset
  940. )
  942. # Lollipop chart with vertical stems
  943. ggplot(results_dodged, aes(x = dodge_pos, y = r_abs, color = source)) +
  944. geom_segment(aes(xend = dodge_pos, yend = 0),
  945. linewidth = 0.5) + # Vertical stems from dodged position
  946. geom_point(size = 2) + # Points at dodged position
  947. labs(
  948. #title = "Rank-Biserial Correlation of Two-Sample Wilcoxon Rank Sum Test",
  949. #subtitle = "Harmonized vs. Unharmonized",
  950. x = "Regions 1–62",
  951. y = "Absolute Rank-Biserial Correlation",
  952. color = ""
  953. ) +
  954. publication_theme() +
  955. theme(
  956. axis.text.x = element_text(size = 8),
  957. #axis.text.x = element_blank(),
  958. #axis.ticks.x = element_blank(),
  959. legend.position = "bottom"
  960. ) +
  961. scale_y_continuous(
  962. breaks = seq(0, 1, by = 0.1),
  963. minor_breaks = NULL
  964. ) +
  965. scale_color_manual(values = mycolors) +
  966. scale_x_continuous(breaks = 1:62, labels = 1:62) # Map dodged positions back to regions
  970. #### [Figure 5b] Plot standardized Z by region -------
  971. results_long %>%
  972. ggplot(aes(x = outcome, y = as.numeric(Z),
  973. color = source, group = source)) +
  974. geom_point(size = 1) +
  975. labs(
  976. title = "Standardized Z by Region",
  977. subtitle = "Harmonized vs. Unharmonized",
  978. x = "Region (1–62)",
  979. y = "Z statistic"
  980. ) +
  981. theme(
  982. axis.text.x = element_blank(),
  983. axis.ticks.x = element_blank()
  984. )+
  985. scale_color_manual(values = c("#4059ad", "#f4b942"))+
  986. scale_fill_manual(values = c("#4059ad", "#f4b942"))

ComBat-Predict_section1-5_share.R at commit f7341bc, no license · at the source

Overview

Authors: Yao Xin1, Margaret Gardner2,3, Nicholas J. Tustison4, Philip Cook5, James Gee5, Andreana Benitez6, Jens H. Jensen7, Alzheimer's Disease Neuroimaging Initiative, Lifespan Brain Chart Consortium, Richard Bethlehem8, Jakob Seidlitz2,9,10,11, Aaron F. Alexander‐Bloch2,8,9,10,11, Andrew A. Chen1
ORCID iDs: Yao Xin, Margaret Gardner, Andreana Benitez, Jens H. Jensen
  1. Department of Public Health Sciences Medical University of South Carolina Charleston South Carolina USA
  2. Brain‐Gene‐Development Lab The Children's Hospital of Philadelphia and Penn Medicine Philadelphia Pennsylvania USA
  3. Neuroscience Graduate Group, Perelman School of Medicine University of Pennsylvania Philadelphia Pennsylvania USA
  4. Department of Radiology and Medical Imaging University of Virginia Charlottesville Virginia USA
  5. Department of Radiology University of Pennsylvania Philadelphia Pennsylvania USA
  6. Department of Neurology Medical University of South Carolina Charleston South Carolina USA
  7. Department of Neuroscience Medical University of South Carolina Charleston South Carolina USA
  8. Department of Psychology University of Cambridge Cambridge UK
  9. Department of Psychiatry University of Pennsylvania Philadelphia Pennsylvania USA
  10. Department of Child and Adolescent Psychiatry and Behavioral Sciences The Children's Hospital of Philadelphia Philadelphia Pennsylvania USA
  11. Lifespan Brain Institute of The Children's Hospital of Philadelphia and Penn Medicine Philadelphia Pennsylvania USA
Institutions: Medical University of South Carolina (United States); Children's Hospital of Philadelphia (United States); University of Pennsylvania Health System (United States); University of Pennsylvania (United States); University of Virginia (United States); University of Cambridge (United Kingdom)
Journal: Human brain mapping, volume 47, issue 8, article e70546
Dates: received 4 October 2025; accepted 5 May 2026; published online 19 May 2026; in print June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1002/hbm.70546 · PMID 42157534 · PMCID PMC13581068 · OpenAlex W7161767811
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: structural MRI / diffusion (modality), human (organism), Alzheimer's / dementia (population), computational (subfield)
Methods: Statistics, Machine learning, Preprocessing
Keywords: ADNI, ComBat, cortical thickness, neuroimaging, normative modeling, site effects
MeSH: Alzheimer Disease*, Brain*, Brain Cortical Thickness*, Cerebral Cortex*, Magnetic Resonance Imaging*, Neuroimaging*, Aging, Female, Humans, Male (* major topic)
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: NIH (R01MH123550, R01MH132934, R01MH133843, R01MH134896)
Citations: cited by 4 papers (Europe PMC); 40 references in the paper
Notices: A correction to this paper has been published (42670026, from Europe PMC)

Abstract

Neuroimaging is vital in quantifying brain atrophy due to typical aging and due to neurodegenerative diseases. To collect large samples necessary to model lifespan brain development, research consortiums aggregate images acquired across multiple study sites. Previous studies have demonstrated that this multi‐site study design can lead to site‐related bias, necessitating harmonization of these “site effects.” However, current methodologies are unable to generalize to new sites outside the original harmonized sample, limiting translation to new sites or clinical practice. Here, we propose a method called ComBat‐Predict (CB‐Predict) building upon the ComBat method for site effect adjustment, which extends to data from a new site with smaller sample sizes and unknown site effects. In data from the Alzheimer's Disease Neuroimaging Initiative, our proposed method mitigates bias and yields high accuracy in predicting cortical thickness measures when generalizing the model to new data. Furthermore, we demonstrate that our proposed harmonization method can reduce site‐related variance in centile scores estimated using data from the Lifespan Brain Chart Consortium. Altogether, our results demonstrate that CB‐Predict effectively harmonizes new sites and thereby enables effective translation of neuroimaging models to additional samples.

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

Repositories

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

andy1764/ComBatFamily

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: e0a8af655de110997af722800b5a1d4090f35480, 1 July 2026
Languages: R (7)
Size: 24 files, 7 scripts
Software Heritage: not checked
Found in: the text, “CB‐Predict Extension”
Holds: README, environment (DESCRIPTION), documentation, 2 notebooks
Not found: license file, CITATION.cff, tests, continuous integration
Tools: mgcv (4 files), lme4 (2 files), neuroCombat (2 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
8 files
At the source: github.com/andy1764/ComBatFamily

ntustison/CrossLong

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 6abd90eb1a7c02c7a9e42afcba6daaff553bef10, 11 August 2026
Languages: R (52), Stan (15)
Size: 280 files, 67 scripts
Software Heritage: not checked
Found in: the text, “ADNI Data Processing”
Holds: README, 38 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Stan (28 files), tidyverse (14 files), ggplot2 (4 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
68 files
At the source: github.com/ntustison/CrossLong/

brainchart/Lifespan

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 4a6faa19b8fb7e11e2173fe11998ab3b87c8c5b4, 21 February 2025
Languages: R (21)
Size: 361 files, 21 scripts
Software Heritage: not checked
Found in: “Data Availability Statement”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: ggplot2 (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
22 files
At the source: github.com/brainchart/Lifespan

Munchkin-233/ComBat-Predict_eval

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: f7341bc63f00b5d0380ec087aeaad6ffa22206c9, 5 September 2025
Languages: R (1)
Size: 2 files, 1 script
Software Heritage: not checked
Found in: “Data Availability Statement”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: ggplot2 (1 file), mgcv (1 file), patchwork (1 file), tidyverse (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
2 files
At the source: github.com/Munchkin-233/ComBat-Predict_eval

The paper's code and data availability statement is in the Data section.

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:

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

Data were obtained from the publicly available ADNI database (adni.loni.usc.edu) and the LBCC consortium dataset (brainchart.shinyapps.io/brainchart), both accessed under their respective data use agreements. Links to open datasets of LBCC are listed on https://github.com/brainchart/Lifespan. The investigators within the ADNI contributed to the design and implementation of the ADNI and/or provided data but did not participate in the analysis or writing of this report. A complete listing of ADNI investigators can be found at https://adni.loni.usc.edu/wp‐content/uploads/how_to_apply/ADNI_Acknowledgement_List.pdf (https://adni.loni.usc.edu/wp-content/uploads/how_to_apply/ADNI_Acknowledgement_List.pdf). ComBat‐Predict is implemented as part of the ComBatFamily R package as predict.comfam and publicly accessible via Github: https://github.com/andy1764/ComBatFamily. Code and data analysis scripts are available at: https://github.com/Munchkin‐233/ComBat‐Predict_eval (https://github.com/Munchkin-233/ComBat-Predict_eval).

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, issue, pages, dates, 13 authors, 6 keywords, 10 MeSH terms, 1 funder, 40 references, 1 integrity notice.

Cite

This paper

Xin, Y., Gardner, M., Tustison, N. J., Cook, P., Gee, J., Benitez, A., Jensen, J. H., Alzheimer's Disease Neuroimaging Initiative, Lifespan Brain Chart Consortium, Bethlehem, R., Seidlitz, J., Alexander‐Bloch, A. F., & Chen, A. A. (2026). ComBat-Predict Enhances Generalizability of Neuroimaging Models to New Sites. Human brain mapping, 47(8), e70546. https://doi.org/10.1002/hbm.70546

BibTeX

@article{xin2026combat,
author = {Xin, Yao and Gardner, Margaret and Tustison, Nicholas J. and Cook, Philip and Gee, James and Benitez, Andreana and Jensen, Jens H. and {Alzheimer's Disease Neuroimaging Initiative} and {Lifespan Brain Chart Consortium} and Bethlehem, Richard and Seidlitz, Jakob and Alexander‐Bloch, Aaron F. and Chen, Andrew A.},
title = {{ComBat-Predict Enhances Generalizability of Neuroimaging Models to New Sites}},
journal = {Human brain mapping},
year = {2026},
month = jun,
volume = {47},
number = {8},
pages = {e70546},
publisher = {Wiley},
issn = {1065-9471},
doi = {10.1002/hbm.70546},
url = {https://doi.org/10.1002/hbm.70546},
pmid = {42157534},
pmcid = {PMC13581068}
}

RIS

TY - JOUR
AU - Xin, Yao
AU - Gardner, Margaret
AU - Tustison, Nicholas J.
AU - Cook, Philip
AU - Gee, James
AU - Benitez, Andreana
AU - Jensen, Jens H.
AU - Alzheimer's Disease Neuroimaging Initiative
AU - Lifespan Brain Chart Consortium
AU - Bethlehem, Richard
AU - Seidlitz, Jakob
AU - Alexander‐Bloch, Aaron F.
AU - Chen, Andrew A.
TI - ComBat-Predict Enhances Generalizability of Neuroimaging Models to New Sites
T2 - Human brain mapping
J2 - Hum Brain Mapp
PY - 2026
DA - 2026/06/01
VL - 47
IS - 8
SP - e70546
SN - 1065-9471
PB - Wiley
DO - 10.1002/hbm.70546
UR - https://doi.org/10.1002/hbm.70546
LA - en
ER -

CSL-JSON

{
"id": "10.1002/hbm.70546",
"type": "article-journal",
"title": "ComBat-Predict Enhances Generalizability of Neuroimaging Models to New Sites",
"container-title": "Human brain mapping",
"author": [
{
"family": "Xin",
"given": "Yao"
},
{
"family": "Gardner",
"given": "Margaret"
},
{
"family": "Tustison",
"given": "Nicholas J."
},
{
"family": "Cook",
"given": "Philip"
},
{
"family": "Gee",
"given": "James"
},
{
"family": "Benitez",
"given": "Andreana"
},
{
"family": "Jensen",
"given": "Jens H."
},
{
"literal": "Alzheimer's Disease Neuroimaging Initiative"
},
{
"literal": "Lifespan Brain Chart Consortium"
},
{
"family": "Bethlehem",
"given": "Richard"
},
{
"family": "Seidlitz",
"given": "Jakob"
},
{
"family": "Alexander‐Bloch",
"given": "Aaron F."
},
{
"family": "Chen",
"given": "Andrew A."
}
],
"container-title-short": "Hum Brain Mapp",
"volume": "47",
"issue": "8",
"page": "e70546",
"DOI": "10.1002/hbm.70546",
"PMID": "42157534",
"PMCID": "PMC13581068",
"ISSN": "1065-9471",
"publisher": "Wiley",
"URL": "https://doi.org/10.1002/hbm.70546",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
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.1162/imag.a.1290 [code]
Calibration of MRI-based reference intervals to new samples.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Alzheimer's / dementia, structural MRI / diffusion, 11 references, 3 authors
[2] doi:10.1038/s41467-026-73072-6 [code]
Mapping the spatiotemporal continuum of structural connectivity development across the human connectome in youth.
Journal: Nature communications
In common: neuroCombat, mgcv, lme4, 3 other tools, structural MRI / diffusion, 2 references
[3] doi:10.1002/hbm.70559 [code]
Replicability of Functional Brain Networks: A Study Through the Lens of Seven Resting-State Networks.
Journal: Human brain mapping
In common: neuroCombat, lme4, ggplot2, 1 other tool, 4 references
[4] doi:10.1038/s41593-026-02363-4 [code]
Cortical thickness changes precede high levels of amyloid by at least 7 years.
Journal: Nature neuroscience
In common: mgcv, lme4, patchwork, 2 other tools, Alzheimer's / dementia, structural MRI / diffusion, 2 references
[5] doi:10.1162/imag.a.1269 [code]
From early to contemporary normative modeling: Mapping individual differences in neurophysiological signals.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: computational, 7 references
[6] doi:10.1016/j.celrep.2026.117505 [code]
Impaired spatial coding and neuronal hyperactivity in the medial entorhinal cortex of aged APP knock-in mice.
Journal: Cell reports
In common: Stan, mgcv, lme4, 3 other tools, Alzheimer's / dementia
[7] doi:10.1162/imag.a.1337 [code]
Data quality biases normative models derived from fetal brain MRI.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: mgcv, patchwork, ggplot2, 1 other tool, structural MRI / diffusion, 3 references
[8] doi:10.1371/journal.pone.0355165 [code]
Pupillary dynamics during hands-off L2 driving and transitions of control under high cognitive load.
Journal: PloS one
In common: Stan, mgcv, lme4, 3 other tools
[9] doi:10.1093/braincomms/fcag343 [code]
Long-term brain volume trajectories and lifestyle associations in cognitively normal adults: the BRAIN-STRIDE study.
Journal: Brain communications
In common: lme4, patchwork, ggplot2, 1 other tool, structural MRI / diffusion, 3 references
[10] doi:10.7554/elife.103097 [code]
Canonical neurodevelopmental trajectories of structural and functional manifolds.
Journal: eLife
In common: mgcv, ggplot2, tidyverse, structural MRI / diffusion, 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.