Convergent structural brain alterations in chronic pain: a multi-metric individual participant data meta-analysis.

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] § Materials and methods › Subcortical segmentation ↔ Figure_making.R, lines 3228–3287 · score 0.94 · inferior lateral, ventral diencephalon, mid posterior, mid anterior, cerebellar white matter, cerebellar cortex
[2] § Materials and methods › MRI preprocessing—cortical parcellation ↔ Meta_analysis.R, lines 3401–3460 · score 0.80 · GausCurv, ThickAvg, MeanCurv, Gaussian curvature, SurfArea, GrayVol
[3] § Materials and methods › MRI preprocessing—cortical parcellation ↔ CP_OA.R, lines 744–774 · score 0.74 · CurvInd, FoldInd, NumVert, GausCurv
[4] § Results › IPD meta-analysis results ↔ Meta_analysis.R, lines 2084–2144 · score 0.73 · left ventral diencephalon, right cerebellar white, left pallidum, brain stem, prediction intervals, meta
[5] § Materials and methods › MRI preprocessing—cortical parcellation ↔ recon_all.sh, lines 1–54 · score 0.60 · FreeSurfer, T1 weighted, recon

Paper

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

The paper is loaded when this pane is shown.

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

The authors' code

R · 4,084 lines · 138 KB · no license · 2 matches

library(tidyverse)
library(data.table)
library(multcomp)
library(effsize)
library(writexl)
library(ggpubr)
library(plotly)
library(reshape2)
library(dplyr)
library(tidyr)
library(stringr)
library(forcats)
library(readr)
library(glmnet)
library(pROC)
library(caret)
library(pscl)
library(corrplot)
library(car)
library(readxl)
library(esvis)
library(metafor)
library(metaviz)
library(smplot2)
library(reshape2)
library(ggridges)
library(meta)
library(introdataviz)
library(smplot2)
library(gghalves)
library(ggdist)
library(meta)
library(forcats)
library(tidytext)
#Datasets used for meta analysis
#Random effects model?
migraine_summary_2_long
clbp_summary_2_long
FM_Summary_2_long
OA_DK_Summary_2_long
ptn_summary_long_2
#Need to recalculate effect size and variance
OA_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = OA_DK_Summary_2_long
)
FM_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = FM_Summary_2_long
)
CLBP_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbp_summary_2_long
)
migraine_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = migraine_summary_2_long
)
ptn_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = ptn_summary_2_long
)
fm2_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = fm2_summary_2_long
)
clbp2_s1_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbp2_s1_summary_long
)
clbp2_s2_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbp2_s2_summary_long
)
target_measures <- c("GrayVol","SurfArea","ThickAvg","MeanCurv","GausCurv")
oa_filtered <- subset(OA_escalc, Measurement %in% target_measures)
fm_filtered <- subset(FM_escalc, Measurement %in% target_measures)
clbp_filtered <- subset(CLBP_escalc, Measurement %in% target_measures)
m_filtered <- subset(migraine_escalc, Measurement %in% target_measures)
ptn_filtered <- subset(ptn_escalc, Measurement %in% target_measures)
fm2_filtered <- subset(fm2_escalc, Measurement %in% target_measures)
clbp2_s1_filtered <- subset(clbp2_s1_escalc, Measurement %in% target_measures)
clbp2_s2_filtered <- subset(clbp2_s2_escalc, Measurement %in% target_measures)
oa_filtered$study <- 'OA'
fm_filtered$study <- 'FM'
clbp_filtered$study <- 'CLBP'
m_filtered$study <- 'migraine'
ptn_filtered$study <- 'PTN'
fm2_filtered$study <- 'FM2'
clbp2_s1_filtered$study <- 'CLBP2_S1'
clbp2_s2_filtered$study <- 'CLBP2_S2'
meta_df <- rbind(oa_filtered, fm_filtered, clbp_filtered, m_filtered)
meta_results <- do.call(rbind, lapply(split(meta_df, list(meta_df$Region,
meta_df$Measurement)),
function(dfm){
if(nrow(dfm) >=2) {
res <- rma(yi = dfm$yi, vi = dfm$vi, method = 'REML')
data.frame(
Region = dfm$Region[1],
Measurement = dfm$Measurement[1],
k = res$k,
Estimate = res$b,
SE = res$se,
zval = res$zval,
pval = res$pval,
CI_lb = res$ci.lb,
CI_ub = res$ci.ub,
I2 = res$I2
)
} else {
NULL
}
}))
view(meta_results)
meta_results <- meta_results %>%
group_by(Measurement) %>%
mutate(FDR = p.adjust(pval, method = 'BH')) %>%
ungroup()
write.csv(meta_results, "~/CP_Meta_Results.csv", row.names = F)
#######################################################################
#Visualization
#######################################################################
ggplot(meta_results, aes(x = Estimate, y = reorder(paste(Region, Measurement, sep = " - "), Estimate))) +
geom_point(aes(color = FDR < 0.05), size = 3) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), height = 0.3) +
scale_color_manual(values = c("grey", "red")) +
labs(x = "Meta-analytic Hedges' g", y = "Region - Measurement", color = "Significant") +
theme_minimal()
ggplot(meta_results, aes(x = Estimate, y = Region)) +
geom_point(aes(color = FDR < 0.05), size = 3) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), height = 0.3) +
scale_color_manual(values = c("grey","red")) +
labs(x = 'Meta-analytic Hedges g', y = "Region = Measurement", color = 'Significant') +
theme_minimal() +
facet_wrap(~Measurement, ncol = 5)
ggplot(meta_results, aes(x = Estimate, y = I2)) +
geom_point(aes(color = FDR < 0.05)) +
geom_hline(yintercept = 50, linetype = "dashed", color = "gray") +
labs(x = "Effect Size (Hedges' g)", y = "I² (%)", color = "Significant") +
theme_minimal()
ggplot(meta_results, aes(x = Estimate, y = -log10(pval))) +
geom_point(aes(color = I2 > 50), size = 2) +
geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
labs(x = "Effect Size", y = "-log10(p-value)", color = "High Heterogeneity (I² > 50)") +
theme_minimal()
ggplot(meta_results, aes(x = Estimate, y = reorder(Region, Estimate))) +
geom_point(aes(color = pval < 0.05)) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), height = 0.3) +
facet_wrap(~Measurement, scales = "free_y") +
theme_bw() +
labs(x = "Meta-analytic Effect Size (Hedges' g)", y = "Region", color = "Significant")
ggplot(meta_results, aes(x = Estimate, y = reorder(Region, Estimate))) +
geom_point() +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), height = 0.2) +
facet_wrap(~Measurement, scales = "free_y", ncol = 5) +
labs(
x = "Effect Size (Hedges' g)",
y = "Region",
title = "Forest Plot by Measurement"
) +
theme_minimal()
#######################################################################
#Forest plots for meta results
#######################################################################
meta_filt <- meta_df[grepl("entorhinal", meta_df$Region), ]
meta_filt <- subset(meta_filt, Measurement %in% c('SurfArea','GrayVol'))
meta_visual <- data.frame(Hedges = meta_filt$yi, SE = sqrt(meta_filt$vi),
Measurement = meta_filt$Measurement, Study = meta_filt$study,
Region = meta_filt$Region)
meta_SA <- subset(meta_visual, Measurement == 'SurfArea')
meta_vol <- subset(meta_visual, Measurement == 'GrayVol')
viz_forest(meta_SA, group = meta_SA$Region, study_labels = meta_SA$Study,
annotate_CI = T, xlab = 'Hedges G')
viz_forest(meta_vol, group = meta_vol$Region, study_labels = meta_vol$Study,
annotate_CI = T, xlab = 'Hedges G')
viz_forest(meta_visual, group = meta_visual$Measurement, study_labels = meta_visual$Region,
annotate_CI = T, xlab = 'Hedges G')
inftemp <- meta_df[grepl("inferiortemporal", meta_df$Region), ]
inftemp <- subset(inftemp, Measurement %in% c('GausCurv'))
inftemp_vis <- data.frame(Hedges = inftemp$yi, SE = sqrt(inftemp$vi),
Measurement = inftemp$Measurement, Study = inftemp$study,
Region = inftemp$Region)
viz_forest(inftemp_vis, group = inftemp_vis$Region, study_labels = inftemp_vis$Study,
annotate_CI = T, xlab = 'Hedges G')
#######################################################################
#Add PTN and redo
#######################################################################
meta_df_2 <- rbind(oa_filtered, fm_filtered, clbp_filtered, m_filtered, ptn_filtered)
meta_results_2 <- do.call(rbind, lapply(split(meta_df_2, list(meta_df_2$Region,
meta_df_2$Measurement)),
function(dfm){
if(nrow(dfm) >=2) {
res <- rma(yi = dfm$yi, vi = dfm$vi, method = 'REML')
data.frame(
Region = dfm$Region[1],
Measurement = dfm$Measurement[1],
k = res$k,
Estimate = res$b,
SE = res$se,
zval = res$zval,
pval = res$pval,
CI_lb = res$ci.lb,
CI_ub = res$ci.ub,
I2 = res$I2
)
} else {
NULL
}
}))
view(meta_results_2)
meta_results_2 <- meta_results_2 %>%
group_by(Measurement) %>%
mutate(FDR = p.adjust(pval, method = 'BH')) %>%
ungroup()
view(meta_results_2)
meta_filt_2 <- meta_df_2[grepl("entorhinal", meta_df_2$Region), ]
meta_filt_2 <- subset(meta_filt_2, Measurement %in% c('SurfArea','GrayVol'))
meta_visual_2 <- data.frame(Hedges = meta_filt_2$yi, SE = sqrt(meta_filt_2$vi),
Measurement = meta_filt_2$Measurement, Study = meta_filt_2$study,
Region = meta_filt_2$Region)
meta_SA_2 <- subset(meta_visual_2, Measurement == 'SurfArea')
meta_vol_2 <- subset(meta_visual_2, Measurement == 'GrayVol')
viz_forest(meta_SA_2, group = meta_SA_2$Region, study_labels = meta_SA_2$Study,
annotate_CI = T, xlab = 'Hedges G')
viz_forest(meta_vol_2, group = meta_vol_2$Region, study_labels = meta_vol_2$Study,
annotate_CI = T, xlab = 'Hedges G')
write.csv(meta_results_2, "~/CP_Meta_Results_2.csv", row.names = F)
#######################################################################
#Add FM2 and redo
#######################################################################
meta_df_3 <- rbind(oa_filtered, fm_filtered, clbp_filtered, m_filtered,
ptn_filtered, fm2_filtered)
meta_results_3 <- do.call(rbind, lapply(split(meta_df_3, list(meta_df_3$Region,
meta_df_3$Measurement)),
function(dfm){
if(nrow(dfm) >=2) {
res <- rma(yi = dfm$yi, vi = dfm$vi, method = 'REML')
data.frame(
Region = dfm$Region[1],
Measurement = dfm$Measurement[1],
k = res$k,
Estimate = res$b,
SE = res$se,
zval = res$zval,
pval = res$pval,
CI_lb = res$ci.lb,
CI_ub = res$ci.ub,
I2 = res$I2
)
} else {
NULL
}
}))
view(meta_results_3)
meta_results_3 <- meta_results_3 %>%
group_by(Measurement) %>%
mutate(FDR = p.adjust(pval, method = 'BH')) %>%
ungroup()
view(meta_results_3)
meta_filt_3 <- meta_df_3[grepl("entorhinal", meta_df_3$Region), ]
meta_filt_3 <- subset(meta_filt_3, Measurement %in% c('SurfArea','GrayVol'))
meta_visual_3 <- data.frame(Hedges = meta_filt_3$yi, SE = sqrt(meta_filt_3$vi),
Measurement = meta_filt_3$Measurement, Study = meta_filt_3$study,
Region = meta_filt_3$Region)
meta_SA_3 <- subset(meta_visual_3, Measurement == 'SurfArea')
meta_vol_3 <- subset(meta_visual_3, Measurement == 'GrayVol')
viz_forest(meta_SA_3, group = meta_SA_3$Region, study_labels = meta_SA_3$Study,
annotate_CI = T, xlab = 'Hedges G')
viz_forest(meta_vol_3, group = meta_vol_3$Region, study_labels = meta_vol_3$Study,
annotate_CI = T, xlab = 'Hedges G')
ugh <- subset(meta_vol_3, Region == 'rh_entorhinal')
viz_forest(ugh, study_labels =ugh$Study, xlab = 'Hedges G',
text_size = 7)
ggplot(meta_results_3, aes(x = Estimate, y = reorder(paste(Region, Measurement, sep = " - "), Estimate))) +
geom_point(aes(color = FDR < 0.05), size = 3) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), height = 0.3) +
scale_color_manual(values = c("grey", "red")) +
labs(x = "Meta-analytic Hedges' g", y = "Region - Measurement", color = "Significant") +
theme_minimal()
ggplot(meta_results_3, aes(x = Estimate, y = Region)) +
geom_point(aes(color = FDR < 0.05), size = 3) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), height = 0.3) +
scale_color_manual(values = c("grey","red")) +
labs(x = 'Meta-analytic Hedges g', y = "Region = Measurement", color = 'Significant') +
theme_minimal() +
facet_wrap(~Measurement, ncol = 5)
ggplot(meta_results_3, aes(x = Estimate, y = I2)) +
geom_point(aes(color = FDR < 0.05)) +
geom_hline(yintercept = 50, linetype = "dashed", color = "gray") +
labs(x = "Effect Size (Hedges' g)", y = "I² (%)", color = "Significant") +
theme_minimal()
ggplot(meta_results_3, aes(x = Estimate, y = -log10(pval))) +
geom_point(aes(color = I2 > 50), size = 2) +
geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
labs(x = "Effect Size", y = "-log10(p-value)", color = "High Heterogeneity (I² > 50)") +
theme_minimal()
ggplot(meta_results_3, aes(x = Estimate, y = reorder(Region, Estimate))) +
geom_point(aes(color = pval < 0.05)) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), height = 0.3) +
facet_wrap(~Measurement, scales = "free_y", ncol = 5) +
theme_bw() +
labs(x = "Meta-analytic Effect Size (Hedges' g)", y = "Region", color = "Significant")
ggplot(meta_results_3, aes(x = Estimate, y = reorder(Region, Estimate))) +
geom_point() +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), height = 0.2) +
facet_wrap(~Measurement, scales = "free_y", ncol = 5) +
labs(
x = "Effect Size (Hedges' g)",
y = "Region",
title = "Forest Plot by Measurement"
) +
theme_minimal()
write.csv(meta_results_3, "~/CP_Meta_Results_3.csv", row.names = F)
#######################################################################
#Use M results and redo
#######################################################################
OA_m_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = OA_M_summary_long
)
oa_m_filtered <- subset(OA_m_escalc, Measurement %in% target_measures)
oa_m_filtered$study <- 'oa'
CLBP_m_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbp_m_summary_long
)
clbp_m_filtered <- subset(CLBP_m_escalc, Measurement %in% target_measures)
clbp_m_filtered$study <- 'clbp'
migraine_m_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = migraine_m_summary_long
)
migraine_m_filtered <- subset(migraine_m_escalc, Measurement %in% target_measures)
migraine_m_filtered$study <- 'migraine'
ptn_m_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = ptn_m_summary_long
)
ptn_m_filtered <- subset(ptn_m_escalc, Measurement %in% target_measures)
ptn_m_filtered$study <- 'ptn'
meta_df_m <- rbind(oa_m_filtered, clbp_m_filtered, migraine_m_filtered,
ptn_m_filtered)
meta_results_m <- do.call(rbind, lapply(split(meta_df_m, list(meta_df_m$Region,
meta_df_m$Measurement)),
function(dfm){
if(nrow(dfm) >=2) {
res <- rma(yi = dfm$yi, vi = dfm$vi, method = 'REML')
data.frame(
Region = dfm$Region[1],
Measurement = dfm$Measurement[1],
k = res$k,
Estimate = res$b,
SE = res$se,
zval = res$zval,
pval = res$pval,
CI_lb = res$ci.lb,
CI_ub = res$ci.ub,
I2 = res$I2
)
} else {
NULL
}
}))
view(meta_results_m)
meta_results_m <- meta_results_m %>%
group_by(Measurement) %>%
mutate(FDR = p.adjust(pval, method = 'BH')) %>%
ungroup()
ggplot(meta_results_m, aes(x = Estimate, y = reorder(Region, Estimate))) +
geom_point(aes(color = pval < 0.05)) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), height = 0.3) +
geom_vline(xintercept = 0, linetype = "solid", color = "black", size = 1) + # <-- bold line at 0
facet_wrap(~Measurement, scales = "free_y", ncol = 5) +
theme_bw() +
labs(x = "Meta-analytic Effect Size (Hedges' g)", y = "Region", color = "Significant")
clbp2_s1_m_filtered <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbp2_s1_dk_m_summary_long
)
clbp2_s1_m_filtered <- subset(clbp2_s1_m_filtered, Measurement %in% target_measures)
clbp2_s1_m_filtered$study <- 'cbp2_s1'
clbp2_s2_m_filtered <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbp2_s2_dk_m_summary_long
)
clbp2_s2_m_filtered <- subset(clbp2_s2_m_filtered, Measurement %in% target_measures)
clbp2_s2_m_filtered$study <- 'cbp2_s2'
meta_df_m_2 <- rbind(oa_m_filtered, clbp_m_filtered, migraine_m_filtered,
ptn_m_filtered, clbp2_s1_m_filtered, clbp2_s2_m_filtered)
meta_results_m_2 <- do.call(rbind, lapply(split(meta_df_m_2, list(meta_df_m_2$Region,
meta_df_m_2$Measurement)),
function(dfm){
if(nrow(dfm) >=2) {
res <- rma(yi = dfm$yi, vi = dfm$vi, method = 'REML')
pred <- predict(res)
data.frame(
Region = dfm$Region[1],
Measurement = dfm$Measurement[1],
k = res$k,
Estimate = res$b,
SE = res$se,
zval = res$zval,
pval = res$pval,
CI_lb = res$ci.lb,
CI_ub = res$ci.ub,
PI_lb = pred$pi.lb,
PI_ub = pred$pi.ub,
I2 = res$I2,
Q = res$QE,
pval_Q = res$QEp,
t2 = res$tau2
)
} else {
NULL
}
}))
meta_results_m_2 <- meta_results_m_2 %>%
group_by(Measurement) %>%
mutate(FDR = p.adjust(pval, method = 'BH')) %>%
ungroup()
meta_results_m_2 <- meta_results_m_2 %>%
group_by(Measurement) %>%
mutate(Q_FDR_p = p.adjust(pval_Q, method = 'BH')) %>%
ungroup()
ggplot(meta_results_m_2, aes(x = Estimate, y = reorder(Region, Estimate))) +
geom_point(aes(color = FDR < 0.05)) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), height = 0.3) +
geom_vline(xintercept = 0, linetype = "solid", color = "black", linewidth = 1) + # <-- bold line at 0
facet_wrap(~Measurement, scales = "free_y", ncol = 5) +
theme_bw() +
labs(x = "Meta-analytic Effect Size (Hedges' g)", y = "Region", color = "Significant")
#######################################################################
#Use F results and redo
#######################################################################
OA_f_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = OA_F_summary_long
)
oa_f_filtered <- subset(OA_f_escalc, Measurement %in% target_measures)
oa_f_filtered$study <- 'oa'
CLBP_f_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbp_f_summary_long
)
clbp_f_filtered <- subset(CLBP_f_escalc, Measurement %in% target_measures)
clbp_f_filtered$study <- 'clbp'
migraine_f_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = migraine_f_summary_long
)
migraine_f_filtered <- subset(migraine_f_escalc, Measurement %in% target_measures)
migraine_f_filtered$study <- 'migraine'
ptn_f_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = ptn_f_summary_long
)
ptn_f_filtered <- subset(ptn_f_escalc, Measurement %in% target_measures)
ptn_f_filtered$study <- 'ptn'
meta_df_f <- rbind(oa_m_filtered, clbp_m_filtered, migraine_m_filtered,
ptn_m_filtered, fm_filtered, fm2_filtered)
meta_results_f <- do.call(rbind, lapply(split(meta_df_f, list(meta_df_f$Region,
meta_df_f$Measurement)),
function(dfm){
if(nrow(dfm) >=2) {
res <- rma(yi = dfm$yi, vi = dfm$vi, method = 'REML')
data.frame(
Region = dfm$Region[1],
Measurement = dfm$Measurement[1],
k = res$k,
Estimate = res$b,
SE = res$se,
zval = res$zval,
pval = res$pval,
CI_lb = res$ci.lb,
CI_ub = res$ci.ub,
I2 = res$I2
)
} else {
NULL
}
}))
meta_results_f <- meta_results_f %>%
group_by(Measurement) %>%
mutate(FDR = p.adjust(pval, method = 'BH')) %>%
ungroup()
view(meta_results_f)
write.csv(meta_results_m, "~/CP_Meta_Results_M.csv", row.names = F)
write.csv(meta_results_f, "~/CP_Meta_Results_F.csv", row.names = F)
entorhinal_m <- meta_df_m[grepl("entorhinal", meta_df_m$Region), ]
entorhinal_m <- subset(entorhinal_m, Measurement %in% c('SurfArea','GrayVol'))
entorinal_m_vis <- data.frame(Hedges = entorhinal_m$yi, SE = sqrt(entorhinal_m$vi),
Measurement = entorhinal_m$Measurement, Study = entorhinal_m$study,
Region = entorhinal_m$Region)
entorhinal_sa_m <- subset(entorinal_m_vis, Measurement == 'SurfArea')
entorhinal_vol_m <- subset(entorinal_m_vis, Measurement == 'GrayVol')
viz_forest(entorhinal_sa_m, group = entorhinal_sa_m$Region, study_labels = entorhinal_sa_m$Study,
annotate_CI = T, xlab = 'Hedges G')
viz_forest(entorhinal_vol_m, group = entorhinal_vol_m$Region, study_labels = entorhinal_vol_m$Study,
annotate_CI = T, xlab = 'Hedges G')
entorhinal_f <- meta_df_f[grepl("entorhinal", meta_df_f$Region), ]
entorhinal_f <- subset(entorhinal_f, Measurement %in% c('SurfArea','GrayVol'))
entorinal_f_vis <- data.frame(Hedges = entorhinal_f$yi, SE = sqrt(entorhinal_f$vi),
Measurement = entorhinal_f$Measurement, Study = entorhinal_f$study,
Region = entorhinal_f$Region)
entorhinal_sa_f <- subset(entorinal_f_vis, Measurement == 'SurfArea')
entorhinal_vol_f <- subset(entorinal_f_vis, Measurement == 'GrayVol')
viz_forest(entorhinal_sa_f, group = entorhinal_sa_f$Region, study_labels = entorhinal_sa_f$Study,
annotate_CI = T, xlab = 'Hedges G')
viz_forest(entorhinal_vol_f, group = entorhinal_vol_f$Region, study_labels = entorhinal_vol_f$Study,
annotate_CI = T, xlab = 'Hedges G')
clbp2_s1_f_filtered <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbp2_s1_dk_f_summary_long
)
clbp2_s1_f_filtered <- subset(clbp2_s1_f_filtered, Measurement %in% target_measures)
clbp2_s1_f_filtered$study <- 'cbp2_s1'
clbp2_s2_f_filtered <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbp2_s2_dk_f_summary_long
)
clbp2_s2_f_filtered <- subset(clbp2_s2_f_filtered, Measurement %in% target_measures)
clbp2_s2_f_filtered$study <- 'cbp2_s2'
meta_df_f_2 <- rbind(oa_f_filtered, clbp_f_filtered, migraine_f_filtered,
ptn_f_filtered, fm_filtered, fm2_filtered, clbp2_s1_f_filtered,
clbp2_s2_f_filtered)
meta_results_f_2 <- do.call(rbind, lapply(split(meta_df_f_2, list(meta_df_f_2$Region,
meta_df_f_2$Measurement)),
function(dfm){
if(nrow(dfm) >=2) {
res <- rma(yi = dfm$yi, vi = dfm$vi, method = 'REML')
pred <- predict(res)
data.frame(
Region = dfm$Region[1],
Measurement = dfm$Measurement[1],
k = res$k,
Estimate = res$b,
SE = res$se,
zval = res$zval,
pval = res$pval,
CI_lb = res$ci.lb,
CI_ub = res$ci.ub,
PI_lb = pred$pi.lb,
PI_ub = pred$pi.ub,
I2 = res$I2,
Q = res$QE,
pval_Q = res$QEp,
t2 = res$tau2
)
} else {
NULL
}
}))
meta_results_f_2 <- meta_results_f_2 %>%
group_by(Measurement) %>%
mutate(FDR = p.adjust(pval, method = 'BH')) %>%
ungroup()
meta_results_f_2 <- meta_results_f_2 %>%
group_by(Measurement) %>%
mutate(Q_FDR_p = p.adjust(pval_Q, method = 'BH')) %>%
ungroup()
ggplot(meta_results_f_2, aes(x = Estimate, y = reorder(Region, Estimate))) +
geom_point(aes(color = FDR < 0.05)) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), height = 0.3) +
geom_vline(xintercept = 0, linetype = "solid", color = "black", linewidth = 1) + # <-- bold line at 0
facet_wrap(~Measurement, scales = "free_y", ncol = 5) +
theme_bw() +
labs(x = "Meta-analytic Effect Size (Hedges' g)", y = "Region", color = "Significant")
#######################################################################
#Add CLBP2 S1 and S2 and redo
#######################################################################
meta_df_4 <- rbind(oa_filtered, fm_filtered, clbp_filtered, m_filtered,
ptn_filtered, fm2_filtered, clbp2_s1_filtered, clbp2_s2_filtered)
meta_results_4 <- do.call(rbind, lapply(split(meta_df_4, list(meta_df_4$Region,
meta_df_4$Measurement)),
function(dfm){
if(nrow(dfm) >=2) {
res <- rma(yi = dfm$yi, vi = dfm$vi, method = 'REML')
data.frame(
Region = dfm$Region[1],
Measurement = dfm$Measurement[1],
k = res$k,
Estimate = res$b,
SE = res$se,
zval = res$zval,
pval = res$pval,
CI_lb = res$ci.lb,
CI_ub = res$ci.ub,
I2 = res$I2
)
} else {
NULL
}
}))
meta_results_4 <- meta_results_4 %>%
group_by(Measurement) %>%
mutate(FDR = p.adjust(pval, method = 'BH')) %>%
ungroup()
view(meta_results_4)
entorhinal_4 <- meta_df_4[grepl("entorhinal", meta_df_4$Region), ]
entorhinal_4_vol <- subset(entorhinal_4, Measurement %in% 'GrayVol')
entorhinal_4_vol_vis <- data.frame(Hedges = entorhinal_4_vol$yi, SE = sqrt(entorhinal_4_vol$vi),
Measurement = entorhinal_4_vol$Measurement, Study = entorhinal_4_vol$study,
Region = entorhinal_4_vol$Region)
entorhinal_4_vol_vis$Study[entorhinal_4_vol_vis$Study == 'OA'] <- 'Tétreault, 2016, Osteoarthritis'
entorhinal_4_vol_vis$Study[entorhinal_4_vol_vis$Study == 'FM'] <- 'Pando-Naude, 2019, Fibromyalgia'
entorhinal_4_vol_vis$Study[entorhinal_4_vol_vis$Study == 'CLBP'] <- 'Makary, 2020, Chronic Lower Back Pain'
entorhinal_4_vol_vis$Study[entorhinal_4_vol_vis$Study == 'migraine'] <- 'Seminowicz, 2020, Migraine'
entorhinal_4_vol_vis$Study[entorhinal_4_vol_vis$Study == 'PTN'] <- 'Filimonova, 2025, Primary Trigeminal Neuralgia'
entorhinal_4_vol_vis$Study[entorhinal_4_vol_vis$Study == 'FM2'] <- 'Balducci, 2022, Fibromyalgia'
entorhinal_4_vol_vis$Study[entorhinal_4_vol_vis$Study == 'CLBP2_S1'] <- 'Mano, 2018, Chronic Lower Back Pain (UK Data)'
entorhinal_4_vol_vis$Study[entorhinal_4_vol_vis$Study == 'CLBP2_S2'] <- 'Mano, 2018, Chronic Lower Back Pain (Japan Data)'
entorhinal_4_vol_vis$Region[entorhinal_4_vol_vis$Region == 'lh_entorhinal'] <- 'Left Entorhinal Cortex Volume'
entorhinal_4_vol_vis$Region[entorhinal_4_vol_vis$Region == 'rh_entorhinal'] <- 'Right Entorhinal Cortex Volume'
viz_forest(entorhinal_4_vol_vis, group = entorhinal_4_vol_vis$Region,
study_labels = entorhinal_4_vol_vis$Study,
annotate_CI = T, xlab = 'Hedges G', variant = 'rain')
ggplot(meta_results_4, aes(x = Estimate, y = Region)) +
geom_point(aes(color = pval < 0.05)) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), height = 0.3) +
geom_vline(xintercept = 0, linetype = "solid", color = "black", size = 1) + # <-- bold line at 0
facet_wrap(~Measurement, scales = "free_y", ncol = 5) +
ggplot2::theme_bw() +
labs(x = "Meta-analytic Effect Size (Hedges' g)", y = "Region", color = "Significant")
#######################################################################
#Subcortical Meta analysis
#######################################################################
OA_Sub_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = OA_subcort_summary
)
OA_Sub_escalc$study <- 'OA'
FM_Sub_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = FM_sub_summary_2
)
FM_Sub_escalc$study <- 'FM'
FM2_Sub_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = fm2_sub_summary
)
FM2_Sub_escalc$significant <- NULL
FM2_Sub_escalc$study <- 'FM2'
CLBP_Sub_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbp_sub_summary
)
CLBP_Sub_escalc$study <- 'CLBP'
CLBP_Sub_s1_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbp2_s1_sub_summary
)
CLBP_Sub_s1_escalc$study <- 'CLBP_s1'
CLBP_Sub_s2_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbp2_s2_sub_summary
)
CLBP_Sub_s2_escalc$study <- 'CLBP_s2'
migraine_sub_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = migraine_sub_summary
)
migraine_sub_escalc$study <- 'migraine'
ptn_sub_escalc <-escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = ptn_sub_summary
)
ptn_sub_escalc$study <- 'ptn'
###############################################################
#Meta subcort
###############################################################
meta_sub <- rbind(OA_Sub_escalc, FM_Sub_escalc, FM2_Sub_escalc, CLBP_Sub_escalc,
CLBP_Sub_s1_escalc, CLBP_Sub_s2_escalc, migraine_sub_escalc,
ptn_sub_escalc)
meta_sub <- meta_sub %>%
filter(!grepl("hypointensities", Region, ignore.case = T))
meta_sub <- meta_sub %>%
filter(!grepl("X5th",Region, ignore.case = T))
meta_sub_results <- do.call(rbind, lapply(split(meta_sub, list(meta_sub$Region)),
function(dfm){
if(nrow(dfm) >=2) {
res <- rma(yi = dfm$yi, vi = dfm$vi, method = 'REML')
data.frame(
Region = dfm$Region[1],
k = res$k,
Estimate = res$b,
SE = res$se,
zval = res$zval,
pval = res$pval,
CI_lb = res$ci.lb,
CI_ub = res$ci.ub,
I2 = res$I2
)
} else {
NULL
}
}))
meta_sub_results <- meta_sub_results %>%
mutate(FDR = p.adjust(pval, method = 'BH'))
ggplot(meta_sub_results, aes(x = Estimate, y = reorder(Region, Estimate))) +
geom_point(aes(color = pval < 0.05)) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), height = 0.3) +
geom_vline(xintercept = 0, linetype = "solid", color = "black", linewidth = 1) + # <-- bold line at 0
theme_bw() +
labs(x = "Meta-analytic Effect Size (Hedges' g)", y = "Region", color = "Significant")
amygdala <- meta_sub[grepl("Amygdala", meta_sub$Region), ]
amygdala_vis <- data.frame(Hedges = amygdala$yi, SE = sqrt(amygdala$vi),
Study = amygdala$study,
Region = amygdala$Region)
viz_forest(amygdala_vis, group = amygdala_vis$Region,
study_labels = amygdala_vis$Study,
annotate_CI = T, xlab = 'Hedges G', variant = 'rain')
#########################################################################
#visual - group effect size
#########################################################################
visual <- meta_results_4
visual$significant <- visual$pval < 0.05 & abs(visual$Estimate) > 0.2
ggplot(visual, aes(x = Estimate, y = -log10(pval))) +
geom_point(aes(color = significant), alpha = 0.7, size = 2) +
geom_vline(xintercept = c(-0.8, 0.8), linetype = "dashed", color = "gray") +
geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "gray") +
geom_text(aes(label = ifelse(significant, Region, "")), hjust = 1.1, vjust = 0.5, size = 3) +
scale_color_manual(values = c("grey", "red")) +
labs(x = "Effect Size (Hedge's g)",
y = "-log10(p-value)",
title = "Volcano Plot of Brain Regions",
color = "Significant") +
theme_minimal() +
theme(legend.position = "top") +
facet_wrap(~Measurement) +
xlim(-2,2)
visual_long <- visual %>%
separate(Region, into = c("hemi","region"), sep = "_") %>%
mutate(hemi = dplyr::recode(hemi, lh = 'left', rh = 'right'))
regions <- visual_long$region %>% str_to_lower()
miss <- setdiff(regions, atlas_regions)
matching <- sapply(miss, function(regions){
distances <- stringdist(regions, atlas_regions, method = 'jw')
atlas_regions[which.min(distances)]
})
recoding <- setNames(matching, miss)
recoding
visual_long <- visual_long %>%
mutate(region = str_to_lower(region),
region = ifelse(region %in% names(recoding),
recoding[region], region))
visual_long$Measurement[visual_long$Measurement == 'SurfArea'] <- 'Surface Area'
visual_long$Measurement[visual_long$Measurement == 'GrayVol'] <- 'Volume'
visual_long$Measurement[visual_long$Measurement == 'ThickAvg'] <- 'Cortical Thickness'
visual_long$Measurement[visual_long$Measurement == 'MeanCurv'] <- 'Mean Curvature'
visual_long$Measurement[visual_long$Measurement == 'GausCurv'] <- 'Gaussian Curvature'
colnames(visual_long)[which(names(visual_long) == "Estimate")] <- "Hedge's G"
visual_long$sig_outline <- ifelse(visual_long$FDR < 0.05, "sig", "ns")
ggplot(visual_long, aes(fill = -`Hedge's G`)) +
geom_brain(atlas = dk) +
scale_fill_distiller(palette = 'RdBu',limits = c(-.5,.5)) +
theme_void() +
labs(title = "Regional Effect Sizes in Individuals with Chronic Pain vs Controls") +
facet_wrap(~factor(Measurement, levels = c("Surface Area","Volume","Cortical Thickness",
"Mean Curvature",'Gaussian Curvature')), ncol = 1) +
theme(legend.position = 'bottom',
legend.text = element_text(size = 12),
legend.key.size = unit(1, 'cm'),
legend.title = element_blank())
#006f3c
ggplot(visual_long, aes(fill = -`Hedge's G`, color = sig_outline, size = sig_outline)) +
geom_brain(atlas = dk) +
scale_fill_distiller(palette = 'RdBu', limits = c(-.5, .5)) +
scale_color_manual(values = c("sig" = "black", "ns" = "grey70")) +
scale_size_manual(values = c("sig" = 0.7, "ns" = 0.2)) + # thicker red, thinner black
theme_void() +
labs(title = "Regional Effect Sizes in Individuals with Chronic Pain vs Controls") +
facet_wrap(~factor(Measurement, levels = c("Surface Area","Volume","Cortical Thickness",
"Mean Curvature",'Gaussian Curvature')), ncol = 1) +
theme(
legend.position = 'bottom',
legend.text = element_text(size = 12),
legend.key.size = unit(1, 'cm'),
legend.title = element_blank()
)
ggplot(visual_long, aes(fill = `Hedge's G`)) +
geom_brain(atlas = dk) +
scale_fill_distiller(palette = 'RdBu', limits = c(-.5, .5)) +
theme_void() +
labs(title = "Regional Effect Sizes in Individuals with Chronic Pain vs Controls") +
facet_grid(rows = vars(hemi), cols = vars(factor(Measurement, levels = c(
"Surface Area", "Volume", "Cortical Thickness", "Mean Curvature", "Gaussian Curvature"
)))) +
theme(
legend.position = c(1.05, 0.6),
strip.text = element_text(size = 9),
plot.title = element_text(hjust = 0.5)
)
gaus_meta <- subset(visual, Measurement == 'GausCurv')
gaus_meta <- subset(gaus_meta, FDR < 0.05)
gaus_meta$hemi <- sub("^(lh|rh)_.*", "\\1", gaus_meta$Region)
gaus_meta$BaseRegion <- sub("^(lh|rh)_(.*)", "\\2", gaus_meta$Region)
region_counts <- table(gaus_meta$BaseRegion, gaus_meta$hemi)
bilateral_regions <- rownames(region_counts)[rowSums(region_counts > 0) == 2]
bilateral_regions
gaus_filt <- subset(gaus_meta, BaseRegion %in% bilateral_regions)
gaus_filt_long <- gaus_filt %>%
separate(Region, into = c("hemi","region"), sep = "_") %>%
mutate(hemi = dplyr::recode(hemi, lh = 'left', rh = 'right'))
regions <- gaus_filt_long$region %>% str_to_lower()
miss <- setdiff(regions, atlas_regions)
matching <- sapply(miss, function(regions){
distances <- stringdist(regions, atlas_regions, method = 'jw')
atlas_regions[which.min(distances)]
})
recoding <- setNames(matching, miss)
recoding
gaus_filt_long <- gaus_filt_long %>%
mutate(region = str_to_lower(region),
region = ifelse(region %in% names(recoding),
recoding[region], region))
gaus_filt_long$Measurement[gaus_filt_long$Measurement == 'GausCurv'] <- 'Gaussian Curvature'
colnames(gaus_filt_long)[which(names(gaus_filt_long) == "Estimate")] <- "Hedge's G"
ggplot(gaus_filt_long, aes(fill = `Hedge's G`)) +
geom_brain(atlas = dk) +
scale_fill_distiller(palette = 'RdBu',limits = c(-1, 1)) +
theme_void() +
labs(title = "Effect sizes the Gaussian Curvature of regions significantly associated with chronic pain") +
theme(legend.position = 'bottom')
###################################################################
visual_m <- meta_results_m_2
visual_m$significant <- visual_m$pval < 0.05 & abs(visual_m$Estimate) > 0.2
ggplot(visual_m, aes(x = Estimate, y = -log10(pval))) +
geom_point(aes(color = significant), alpha = 0.7, size = 2) +
geom_vline(xintercept = c(-0.8, 0.8), linetype = "dashed", color = "gray") +
geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "gray") +
geom_text(aes(label = ifelse(significant, Region, "")), hjust = 1.1, vjust = 0.5, size = 3) +
scale_color_manual(values = c("grey", "red")) +
labs(x = "Effect Size (Hedge's g)",
y = "-log10(p-value)",
title = "Volcano Plot of Brain Regions",
color = "Significant") +
theme_minimal() +
theme(legend.position = "top") +
facet_wrap(~Measurement) +
xlim(-2,2)
meta_results_m_2
meta_results_m_2 <- meta_results_m_2 %>%
separate(Region, into = c("hemi","region"), sep = "_") %>%
mutate(hemi = dplyr::recode(hemi, lh = 'left', rh = 'right'))
regions <- meta_results_m_2$region %>% str_to_lower()
miss <- setdiff(regions, atlas_regions)
matching <- sapply(miss, function(regions){
distances <- stringdist(regions, atlas_regions, method = 'jw')
atlas_regions[which.min(distances)]
})
recoding <- setNames(matching, miss)
recoding
meta_results_m_2 <- meta_results_m_2 %>%
mutate(region = str_to_lower(region),
region = ifelse(region %in% names(recoding),
recoding[region], region))
write_xlsx(meta_results_m_2, "meta_m_results.xlsx")
visual_m <- visual_m %>%
separate(region, into = c("hemi","region"), sep = "_") %>%
mutate(hemi = dplyr::recode(hemi, lh = 'left', rh = 'right'))
regions <- visual_m$region %>% str_to_lower()
miss <- setdiff(regions, atlas_regions)
matching <- sapply(miss, function(regions){
distances <- stringdist(regions, atlas_regions, method = 'jw')
atlas_regions[which.min(distances)]
})
recoding <- setNames(matching, miss)
recoding
visual_m <- visual_m %>%
mutate(region = str_to_lower(region),
region = ifelse(region %in% names(recoding),
recoding[region], region))
visual_m$Measurement[visual_m$Measurement == 'SurfArea'] <- 'Surface Area'
visual_m$Measurement[visual_m$Measurement == 'GrayVol'] <- 'Volume'
visual_m$Measurement[visual_m$Measurement == 'ThickAvg'] <- 'Cortical Thickness'
visual_m$Measurement[visual_m$Measurement == 'MeanCurv'] <- 'Mean Curvature'
visual_m$Measurement[visual_m$Measurement == 'GausCurv'] <- 'Gaussian Curvature'
colnames(visual_m)[which(names(visual_m) == "Estimate")] <- "Hedge's G"
ggplot(visual_m_long, aes(fill = `Hedge's G`)) +
geom_brain(atlas = dk) +
scale_fill_distiller(palette = 'RdBu',limits = c(-1, 1)) +
theme_void() +
labs(title = "Regional Effect Sizes in CP vs Controls") +
facet_wrap(~factor(Measurement, levels = c("Surface Area","Volume","Cortical Thickness",
"Mean Curvature",'Gaussian Curvature')), ncol = 1) +
theme(legend.position = c(1.05,.60))
#########################################################################
visual_f <- meta_results_f_2
visual_f$significant <- visual_f$pval < 0.05 & abs(visual_f$Estimate) > 0.2
ggplot(visual_f, aes(x = Estimate, y = -log10(pval))) +
geom_point(aes(color = significant), alpha = 0.7, size = 2) +
geom_vline(xintercept = c(-0.8, 0.8), linetype = "dashed", color = "gray") +
geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "gray") +
geom_text(aes(label = ifelse(significant, region, "")), hjust = 1.1, vjust = 0.5, size = 3) +
scale_color_manual(values = c("grey", "red")) +
labs(x = "Effect Size (Hedge's g)",
y = "-log10(p-value)",
title = "Volcano Plot of Brain Regions",
color = "Significant") +
theme_minimal() +
theme(legend.position = "top") +
facet_wrap(~Measurement) +
xlim(-2,2)
meta_results_f_2
meta_results_f_2 <- meta_results_f_2 %>%
separate(Region, into = c("hemi","region"), sep = "_") %>%
mutate(hemi = dplyr::recode(hemi, lh = 'left', rh = 'right'))
regions <- meta_results_f_2$region %>% str_to_lower()
miss <- setdiff(regions, atlas_regions)
matching <- sapply(miss, function(regions){
distances <- stringdist(regions, atlas_regions, method = 'jw')
atlas_regions[which.min(distances)]
})
recoding <- setNames(matching, miss)
recoding
meta_results_f_2 <- meta_results_f_2 %>%
mutate(region = str_to_lower(region),
region = ifelse(region %in% names(recoding),
recoding[region], region))
write_xlsx(meta_results_f_2, "meta_results_f.xlsx")
visual_f_long <- visual_f %>%
separate(region, into = c("hemi","region"), sep = "_") %>%
mutate(hemi = dplyr::recode(hemi, lh = 'left', rh = 'right'))
regions <- visual_f_long$region %>% str_to_lower()
miss <- setdiff(regions, atlas_regions)
matching <- sapply(miss, function(regions){
distances <- stringdist(regions, atlas_regions, method = 'jw')
atlas_regions[which.min(distances)]
})
recoding <- setNames(matching, miss)
recoding
visual_f_long <- visual_f_long %>%
mutate(region = str_to_lower(region),
region = ifelse(region %in% names(recoding),
recoding[region], region))
visual_f$Measurement[visual_f$Measurement == 'SurfArea'] <- 'Surface Area'
visual_f$Measurement[visual_f$Measurement == 'GrayVol'] <- 'Volume'
visual_f$Measurement[visual_f$Measurement == 'ThickAvg'] <- 'Cortical Thickness'
visual_f$Measurement[visual_f$Measurement == 'MeanCurv'] <- 'Mean Curvature'
visual_f$Measurement[visual_f$Measurement == 'GausCurv'] <- 'Gaussian Curvature'
colnames(visual_f)[which(names(visual_f) == "Estimate")] <- "Hedge's G"
ggplot(visual_f_long, aes(fill = `Hedge's G`)) +
geom_brain(atlas = dk) +
scale_fill_distiller(palette = 'RdBu',limits = c(-1, 1)) +
theme_void() +
labs(title = "Regional Effect Sizes in CP vs Controls") +
facet_wrap(~factor(Measurement, levels = c("Surface Area","Volume","Cortical Thickness",
"Mean Curvature",'Gaussian Curvature')), ncol = 1) +
theme(legend.position = c(1.05,.60))
visual_f_long
identical(visual_f$region, visual_m$region)
identical(visual_f$Measurement, visual_m$Measurement)
identical(visual_f$hemi, visual_m$hemi)
x <- c(.2,-.2,-.3,.4-.68,.9)
t.test(x, mu = 0, alternative = 'two.sided')
meta_test <- data.frame(hemi = visual_m$hemi, Region = visual_m$region, Measurement = visual_m$Measurement,
estimate_m = visual_m$`Hedge's G`, estimate_f = visual_f$`Hedge's G`,
SE_m = visual_m$SE, SE_f = visual_f$SE, CI_lb_m = visual_m$CI_lb,
CI_ub_m = visual_m$CI_ub, CI_lb_f = visual_f$CI_lb, CI_ub_f = visual_f$CI_ub,
PI_lb_m = visual_m$PI_lb, PI_ub_m = visual_m$PI_ub, PI_ub_f = visual_f$PI_lb,
PI_lb_f = visual_f$PI_ub)
meta_test
meta_test$Measurement[meta_test$Measurement == 'SurfArea'] <- 'Surface Area'
meta_test$Measurement[meta_test$Measurement == 'GrayVol'] <- 'Volume'
meta_test$Measurement[meta_test$Measurement == 'ThickAvg'] <- 'Cortical Thickness'
meta_test$Measurement[meta_test$Measurement == 'MeanCurv'] <- 'Mean Curvature'
meta_test$Measurement[meta_test$Measurement == 'GausCurv'] <- 'Gaussian Curvature'
meta_test_extended <- meta_test %>%
bind_rows(tibble(
estimate_m = NA,
estimate_f = NA,
Measurement = "Legend",
Region = NA,
SE_m = NA, SE_f = NA,
CI_lb_m = NA, CI_ub_m = NA,
CI_lb_f = NA, CI_ub_f = NA,
difference = NA
))
meta_test_extended$Measurement <- factor(meta_test_extended$Measurement,
levels = c('Volume','Surface Area','Cortical Thickness',
'Mean Curvature','Gaussian Curvature', "Legend"))
ioflickhuck <- ggplot(meta_test_extended, aes(x = estimate_m, y = estimate_f)) +
geom_point(alpha = .6) +
xlim(-0.75, 0.75) +
ylim(-0.75, 0.75) +
facet_wrap(~Measurement, ncol = 3) +
geom_abline(slope = 1, intercept = 0, linetype = 'longdash',color = '#e43c40') +
geom_smooth(method = 'lm', color = '#00b2a9') +
geom_hline(yintercept = 0, linetype = 'twodash', color = 'black') +
geom_vline(xintercept = 0, linetype = 'twodash', color = 'black') +
xlab('Male estimated effect size') +
ylab('Female estiamted effect size') +
theme_bw(base_size = 20) +
theme(
strip.text = element_text( size = 20),
panel.spacing = unit(3, 'lines')
)
ioflickhuck
ggsave("sex_diff.png", plot = ioflickhuck, dpi = 500, width = 18, height = 12, units = 'in', bg = 'white')
ggplot(meta_test, aes(x = estimate_m, y = estimate_f)) +
geom_point() +
theme_bw() +
geom_smooth(method = 'lm') +
facet_wrap(~Measurement) +
xlim(-1,1) +
ylim(-1,1) +
geom_abline(slope = 1, intercept = 0, linetype = 'longdash',color = '#e43c40') +
geom_abline(slope = -1, intercept = 0, linetype = 'longdash', color = '#e43c40') +
geom_smooth(method = 'lm', color = '#00b2a9') +
geom_hline(yintercept = 0, linetype = 'twodash', color = 'black') +
geom_vline(xintercept = 0, linetype = 'twodash', color = 'black') +
theme(panel.spacing = unit(2))
ggsave("sex_diff_estimates.png", plot = ioflickhuck, dpi = 500, width = 18, height = 12, units = 'in', bg = 'white')
ggplot() +
xlim(-1,1) +
ylim(-1,1) +
theme_minimal() +
geom_hline(yintercept = 0, linetype = 'twodash', color = 'black') +
geom_vline(xintercept = 0, linetype = 'twodash', color = 'black') +
geom_abline(slope = 1, intercept = 0, linetype = 'longdash', color = '#e43c40')
meta_test$difference <- meta_test$estimate_f - meta_test$estimate_m
meta_test$z_diff <- (meta_test$estimate_m - meta_test$estimate_f) / sqrt(meta_test$SE_m^2 + meta_test$SE_f^2)
meta_test$p <- 2*pnorm(abs(meta_test$z_diff), lower.tail = FALSE)
meta_test <- meta_test %>%
group_by(Measurement) %>%
mutate(FDR = p.adjust(p, method = 'BH')) %>%
ungroup()
write_xlsx(meta_test,"meta_sex_cortical.xlsx")
meta_test$effect_direction <- case_when(
sign(meta_test$estimate_m) == sign(meta_test$estimate_f) ~ "Same Direction",
sign(meta_test$estimate_m) != sign(meta_test$estimate_f) ~ "Opposite Direction",
TRUE ~ "Undefined"
)
meta_test$z_diff_mag <- (abs(meta_test$estimate_m) - abs(meta_test$estimate_f)) /
sqrt(meta_test$SE_m^2 + meta_test$SE_f^2)
ggplot(meta_test, aes(x = Measurement, y = z_diff_mag, color = effect_direction)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(width = 0.1, alpha = 0.7, size = 2) +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_color_manual(values = c("Same Direction" = "black", "Opposite Direction" = "#e41a1c")) +
labs(
y = "Z-diff (|m| - |f|)",
title = "Effect Size Strength Difference (M > F = +)",
color = "Direction Match"
) +
theme_minimal() +
facet_wrap(~effect_direction)
ggplot(meta_test, aes(x = Measurement, y = z_diff)) +
geom_boxplot() +
geom_jitter(width = .1)+
geom_hline(yintercept = 0, linetype = "dashed") +
labs(y = "Z-difference (m - f)", title = "Z-diff by Measurement (M > F = +)") +
theme_minimal()
ggplot(meta_test, aes(x = Measurement, y = difference)) +
theme_minimal() +
geom_boxplot() +
geom_jitter(width = .1)
write.csv(meta_test, "~/sex_differences_v462.csv", row.names = F)
gc_estimate <- subset(meta_test, Measurement == 'GausCurv')
shapiro.test(gc_estimate$difference)
t.test(gc_estimate$difference, mu = 0, alternative = 'two.sided')
t.test(gc_estimate$estimate_m, gc_estimate$estimate_f, alternative = 'two.sided')
mc_estimate <- subset(meta_test, Measurement == 'MeanCurv')
shapiro.test(mc_estimate$difference)
t.test(mc_estimate$difference, mu = 0, alternative = 'two.sided')
t.test(mc_estimate$estimate_m, mc_estimate$estimate_f, alternative = 'two.sided')
gv_estimate <- subset(meta_test, Measurement == 'GrayVol')
shapiro.test(gv_estimate$difference)
t.test(gv_estimate$difference, mu = 0, alternative = 'two.sided')
t.test(gv_estimate$estimate_m, gv_estimate$estimate_f, alternative = 'two.sided')
sa_estimate <- subset(meta_test, Measurement == 'SurfArea')
shapiro.test(sa_estimate$difference)
t.test(sa_estimate$difference, mu = 0, alternative = 'two.sided')
ct_estimate <- subset(meta_test, Measurement == 'ThickAvg')
shapiro.test(ct_estimate$difference)
t.test(ct_estimate$difference, mu = 0, alternative = 'two.sided')
###########################################################################
#Sex analysis
###########################################################################
meta_df_m_2
meta_df_f_2
meta_df_m_2$sex <- 'male'
meta_df_f_2$sex <- 'female'
meta_m <- meta_df_m_2[c("Region","Measurement","yi","vi","study","sex")]
meta_f <- meta_df_f_2[c("Region","Measurement","yi","vi","study","sex")]
meta_sex <- rbind(meta_m, meta_f)
res <- rma.mv(
yi = yi, V = vi,
mods = ~ sex,
random = ~1 | study,
data = meta_sex
)
summary(res)
collapsed_combined
collapsed_by_sex
gs <- data.frame()
#use meta_sex for old answers
for (m in unique(collapsed_by_sex$Measurement)) {
df_m <- subset(collapsed_by_sex, Measurement == m)
res_mes <- rma.mv(
yi = yi,
V = vi,
mods = ~ sex,
random = ~ 1 | study,
data = df_m
)
gs <- rbind(gs, data.frame(
Measurement = m,
Effect = res_mes$b,
SE = res_mes$se,
pval = res_mes$pval,
CI_lower = res_mes$ci.lb,
CI_upper = res_mes$ci.ub,
k = res_mes$k
))
print(paste("Measurement:", m))
print(summary(res_mes))
}
measurements <- unique(meta_sex$Measurement)
group_summary <- data.frame()
group_summary_2 <- data.frame()
for (m in measurements) {
for (s in c("female", "male")) {
df_sub <- subset(meta_sex, Measurement == m & sex == s)
res_mes <- rma(
yi = yi,
vi = vi,
method = "REML",
data = df_sub
)
pred <- predict(res_mes)
group_summary_2 <- rbind(group_summary_2, data.frame(
Measurement = m,
Sex = s,
Effect = res_mes$b,
SE = res_mes$se,
CI_lower = res_mes$ci.lb,
CI_upper = res_mes$ci.ub,
PI_lb = pred$pi.lb,
PI_ub = pred$pi.ub,
tau2 = res_mes$tau2,
k = res_mes$k,
I2 = res_mes$I2,
Q = res_mes$QE,
pval_Q = res_mes$QEp
))
}
# Combined effect across both sexes (optional)
df_all <- subset(meta_sex, Measurement == m)
res_all <- rma(
yi = yi,
vi = vi,
method = "REML",
data = df_all
)
pred_all <- predict(res_all)
group_summary_2 <- rbind(group_summary_2, data.frame(
Measurement = m,
Sex = "both",
Effect = res_all$b,
SE = res_all$se,
CI_lower = res_all$ci.lb,
CI_upper = res_all$ci.ub,
PI_lb = pred_all$pi.lb,
PI_ub = pred_all$pi.ub,
tau2 = res_all$tau2,
k = res_all$k,
I2 = res_all$I2,
Q = res_all$QE,
pval_Q = res_all$QEp
))
}
group_summary$Measurement[group_summary$Measurement == 'SurfArea'] <- 'Surface Area'
group_summary$Measurement[group_summary$Measurement == 'GrayVol'] <- 'Volume'
group_summary$Measurement[group_summary$Measurement == 'ThickAvg'] <- 'Cortical Thickness'
group_summary$Measurement[group_summary$Measurement == 'MeanCurv'] <- 'Mean Curvature'
group_summary$Measurement[group_summary$Measurement == 'GausCurv'] <- 'Gaussian Curvature'
group_summary$Sex[group_summary$Sex == 'female'] <- 'Female'
group_summary$Sex[group_summary$Sex == 'male'] <- 'Male'
group_summary$Measurement <- factor(group_summary$Measurement,
levels = c("Gaussian Curvature","Mean Curvature",
"Cortical Thickness","Surface Area","Volume" ))
ggplot(group_summary, aes(x = Measurement, y = Effect, color = Sex, group = Sex)) +
geom_point(position = position_dodge(width = 0.4), size = 3) +
geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper),
position = position_dodge(width = 0.4), width = 0.3) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(title = "Group-level Effect Size by Measurement and Sex",
y = "Hedges' g",
x = "Measurement Type") +
theme_minimal() +
scale_color_manual(values = c("Female" = "#00b2a9", "Male" = "#e43c40", "both" = "black")) +
coord_flip()
group_summary_2$Measurement[group_summary_2$Measurement == 'SurfArea'] <- 'Surface Area'
group_summary_2$Measurement[group_summary_2$Measurement == 'GrayVol'] <- 'Volume'
group_summary_2$Measurement[group_summary_2$Measurement == 'ThickAvg'] <- 'Cortical Thickness'
group_summary_2$Measurement[group_summary_2$Measurement == 'MeanCurv'] <- 'Mean Curvature'
group_summary_2$Measurement[group_summary_2$Measurement == 'GausCurv'] <- 'Gaussian Curvature'
group_summary_2$Sex[group_summary_2$Sex == 'female'] <- 'Female'
group_summary_2$Sex[group_summary_2$Sex == 'male'] <- 'Male'
group_summary_2$Measurement <- factor(group_summary_2$Measurement,
levels = c("Gaussian Curvature","Mean Curvature",
"Cortical Thickness","Surface Area","Volume" ))
ggplot(group_summary_2, aes(x = Measurement, y = Effect, color = Sex, group = Sex)) +
geom_point(position = position_dodge(width = 0.4), size = 3) +
geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper),
position = position_dodge(width = 0.4), width = 0.3) +
geom_errorbar(aes(ymin = PI_lb, ymax = PI_ub),
position = position_dodge(width = 0.4), width = 0.2,
linetype = "dashed", alpha = .5, linewidth = 0.7) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(title = "Group-level Effect Size by Measurement and Sex",
y = "Hedges' g",
x = "Measurement Type") +
theme_minimal() +
scale_color_manual(values = c("Female" = "#00b2a9", "Male" = "#e43c40", "both" = "black")) +
coord_flip()
#Testing
testing_summary <- data.frame()
for (m in measurements) {
for (s in c("female", "male")) {
df_sub <- subset(meta_sex, Measurement == m & sex == s)
res_mes_stu <- rma.mv(
yi = yi,
V = vi,
method = "REML",
random = ~ 1 | study,
data = df_sub
)
testing_summary <- rbind(testing_summary, data.frame(
Measurement = m,
Sex = s,
Effect = res_mes$b,
SE = res_mes$se,
CI_lower = res_mes$ci.lb,
CI_upper = res_mes$ci.ub
))
}
# Combined effect across both sexes (optional)
df_all <- subset(meta_sex, Measurement == m)
res_all <- rma.mv(
yi = yi,
V = vi,
method = "REML",
random = ~ 1 | study,
data = df_all
)
testing_summary <- rbind(testing_summary, data.frame(
Measurement = m,
Sex = "both",
Effect = res_all$b,
SE = res_all$se,
CI_lower = res_all$ci.lb,
CI_upper = res_all$ci.ub
))
}
ggplot(testing_summary, aes(x = Measurement, y = Effect, color = Sex, group = Sex)) +
geom_point(position = position_dodge(width = 0.4), size = 3) +
geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper),
position = position_dodge(width = 0.4), width = 0.2) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(title = "Group-level Effect Size by Measurement and Sex",
y = "Hedges' g",
x = "Measurement Type") +
theme_minimal() +
scale_color_manual(values = c("female" = "#F8766D", "male" = "#00BFC4", "both" = "gray50")) +
coord_flip()
#########################################################################
#redo meta
#########################################################################
oasub_summary_2
fm1sub_summary_2
fm2sub_summary_2
clbpsub_summary_2
clbpsubs1_summary_2
clbpsubs2_summary_2
msub_summary_2
ptnsub_summary_2
oasubescalc_2 <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = oasub_summary_2
)
oasubescalc_2$study <- 'OA'
fmsubescalc_2 <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = fm1sub_summary_2
)
fmsubescalc_2$study <- 'FM1'
fm2subescalc_2 <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = fm2sub_summary_2
)
fm2subescalc_2$study <- 'FM2'
clbpsubescalc_2 <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbpsub_summary_2
)
clbpsubescalc_2$study <- 'CLBP1'
clbpsubs1escalc_2 <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbpsubs1_summary_2
)
clbpsubs1escalc_2$study <- 'CLBP_S1'
clbpsubs2escalc_2 <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbpsubs2_summary_2
)
clbpsubs2escalc_2$study <- 'CLBP_S2'
msubescalc_2 <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = msub_summary_2
)
msubescalc_2$study <- 'migraine'
ptnsubescalc_2 <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = ptnsub_summary_2
)
ptnsubescalc_2$study <- 'PTN'
meta_sub_2 <- rbind(oasubescalc_2, fmsubescalc_2, fm2subescalc_2, clbpsubescalc_2,
clbpsubs1escalc_2, clbpsubs2escalc_2, msubescalc_2,
ptnsubescalc_2)
meta_sub_3 <- rbind(oasubescalc_2, fmsubescalc_2, fm2subescalc_2, clbpsubescalc_2,
clbpsubs1escalc_2, clbpsubs2escalc_2, msubescalc_2,
ptnsubescalc_2)
meta_sub_results_2 <- do.call(rbind, lapply(split(meta_sub_2, list(meta_sub_2$Region)),
function(dfm){
if(nrow(dfm) >= 2) {
res <- rma(yi = dfm$yi, vi = dfm$vi, method = 'REML')
pred <- predict(res)
data.frame(
Region = dfm$Region[1],
k = res$k,
Estimate = res$b,
SE = res$se,
zval = res$zval,
pval = res$pval,
CI_lb = res$ci.lb,
CI_ub = res$ci.ub,
PI_lb = pred$pi.lb,
PI_ub = pred$pi.ub,
I2 = res$I2,
Q = res$QE,
pval_Q = res$QEp,
t2 = res$tau2
)
} else {
NULL
}
}))
meta_sub_results_2 <- meta_sub_results_2 %>%
mutate(FDR = p.adjust(pval, method = 'BH'))
meta_sub_results_2 <- meta_sub_results_2 %>%
mutate(Q_FDR = p.adjust(pval_Q, method = 'BH'))
write_xlsx(meta_sub_results_2, "Meta_Subcortical_Results.xlsx")
meta_sub_results_3 <- do.call(rbind, lapply(split(meta_sub_3, list(meta_sub_3$Region)),
function(dfm){
if(nrow(dfm) >= 2) {
res <- rma(yi = dfm$yi, vi = dfm$vi, method = 'REML')
pred <- predict(res)
data.frame(
Region = dfm$Region[1],
k = res$k,
Estimate = res$b,
SE = res$se,
zval = res$zval,
pval = res$pval,
CI_lb = res$ci.lb,
CI_ub = res$ci.ub,
PI_lb = pred$pi.lb,
PI_ub = pred$pi.ub,
I2 = res$I2,
Q = res$QE,
pval_Q = res$QEp,
t2 = res$tau2
)
} else {
NULL
}
}))
meta_sub_results_3 <- meta_sub_results_3 %>%
mutate(FDR = p.adjust(pval, method = 'BH'))
meta_sub_results_3 <- meta_sub_results_3 %>%
mutate(Q_FDR = p.adjust(pval_Q, method = 'BH'))
ggplot(meta_sub_results_2, aes(x = Estimate, y = reorder(Region, Estimate))) +
geom_point(aes(color = pval < 0.05)) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), height = 0.3) +
geom_vline(xintercept = 0, linetype = "solid", color = "black", linewidth = 1) + # <-- bold line at 0
theme_bw() +
labs(x = "Meta-analytic Effect Size (Hedges' g)", y = "Region", color = "Significant")
ggplot(meta_sub_results_2, aes(x = I2, y = Estimate, size = t2, color = Q_FDR < 0.05)) +
geom_point(alpha = 0.8) +
labs(title = "Effect Size vs. Heterogeneity",
x = "I² (%)", y = "Pooled Effect Size") +
scale_color_manual(values = c("gray", "red"), name = "Significant Heterogeneity") +
theme_minimal() +
ylim(-.5,.5) +
xlim(0,100) +
geom_hline(yintercept = 0, linetype = 'dashed', color = 'red') +
geom_text(aes(label = ifelse(Q_FDR < 0.05, Region, "")), hjust = 1.4, vjust = 0.5, size = 3)
meta_sub_results_2 <- meta_sub_results_2 %>%
mutate(Region = fct_reorder(Region, Estimate, .desc = TRUE))
meta_sub_results_2$Region_clean <- as.character(meta_sub_results_2$Region)
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Right.Lateral.Ventricle'] <- 'Right Lateral Ventricle'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Left.Lateral.Ventricle'] <- 'Left Lateral Ventricle'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Left.Pallidum'] <- 'Left Pallidum'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Right.Pallidum'] <- 'Right Pallidum'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Left.Cerebellum.White.Matter'] <- 'Left Cerebellum White Matter'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Right.Cerebellum.White.Matter'] <- 'Right Cerebellum White Matter'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'X3rd.Ventricle'] <- '3rd Ventricle'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'CC_Anterior'] <- 'CC Anterior'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'X4th.Ventricle'] <- '4th Ventricle'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'CC_Mid_Anterior'] <- 'CC Mid Anterior'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Left.Caudate'] <- 'Left Caudate'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Right.Caudate'] <- 'Right Caudate'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'CC_Posterior'] <- 'CC Posterior'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'CC_Mid_Posterior'] <- 'CC Mid Posterior'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'CC_Central'] <- 'CC Central'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Left.VentralDC'] <- 'Left Ventral Diencephalon'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Right.VentralDC'] <- 'Right Ventral Diencephalon'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Brain.Stem'] <- 'Brain Stem'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Left.Putamen'] <- 'Left Putamen'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Right.Putamen'] <- 'Right Putamen'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Right.Thalamus'] <- 'Right Thalamus'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Left.Thalamus'] <- 'Left Thalamus'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Left.Hippocampus'] <- 'Left Hippocampus'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Right.Hippocampus'] <- 'Right Hippocampus'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Right.Cerebellum.Cortex'] <- 'Right Cerebellum Cortex'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Left.Cerebellum.Cortex'] <- 'Left Cerebellum Cortex'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Left.Amygdala'] <- 'Left Amygdala'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Right.Amygdala'] <- 'Right Amygdala'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Right.Cerebellum.Cortex'] <- 'Right Cerebellum Cortex'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Left.Cerebellum.Cortex'] <- 'Left Cerebellum Cortex'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Left.Amygdala'] <- 'Left Amygdala'
meta_sub_results_2$Region_clean[meta_sub_results_2$Region_clean == 'Right.Amygdala'] <- 'Right Amygdala'
meta_sub_results_2 <- meta_sub_results_2 %>%
mutate(Region_clean = fct_reorder(Region_clean, Estimate, .desc = TRUE))
sub_graph <- ggplot(meta_sub_results_3, aes(y = Region_clean, x = Estimate)) +
geom_errorbarh(aes(xmin = PI_lb, xmax = PI_ub), color = "gray70", height = 0.4, size = 1) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), color = "black", height = 0.2, size = 0.8) +
geom_point(aes(color = FDR < 0.05), size = 2) +
geom_vline(xintercept = 0, linetype = "twodash") +
scale_color_manual(values = c("TRUE" = "#e41a1c", "FALSE" = "black")) +
labs(
title = "Forest Plot of Effect Sizes with Prediction Intervals",
x = "Hedges' g", y = "Region",
color = "FDR < 0.05"
) +
theme_bw(base_size = 19) +
theme(
strip.text = element_text(face = "bold", size = 19),
axis.text.y = element_text(size = 19),
legend.position = "bottom"
)
sub_graph
ggsave("sub_meta.png", plot = sub_graph, dpi = 500, width = 10, height = 12, units = 'in', bg = 'white')
ggplot(meta_sub_results_3, aes(x = Estimate, y = reorder(Region, Estimate))) +
geom_point(aes(color = FDR < 0.05)) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), height = 0.3) +
geom_vline(xintercept = 0, linetype = "solid", color = "black", linewidth = 1) + # <-- bold line at 0
theme_bw() +
labs(x = "Meta-analytic Effect Size (Hedges' g)", y = "Region", color = "Significant")
ggplot(meta_sub_results_3, aes(x = I2, y = Estimate, size = t2, color = Q_FDR < 0.05)) +
geom_point(alpha = 0.8) +
labs(title = "Effect Size vs. Heterogeneity",
x = "I² (%)", y = "Pooled Effect Size") +
scale_color_manual(values = c("gray", "red"), name = "Significant Heterogeneity") +
theme_minimal() +
ylim(-.5,.5) +
xlim(0,100) +
geom_hline(yintercept = 0, linetype = 'dashed', color = 'red') +
geom_text(aes(label = ifelse(Q_FDR < 0.05, Region, "")), hjust = 1.4, vjust = 0.5, size = 3)
meta_sub_results_3 <- meta_sub_results_3 %>%
mutate(Region = fct_reorder(Region, Estimate, .desc = TRUE))
meta_sub_results_3$Region_clean <- as.character(meta_sub_results_3$Region)
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Right.Lateral.Ventricle'] <- 'Right Lateral Ventricle'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Left.Lateral.Ventricle'] <- 'Left Lateral Ventricle'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Left.Pallidum'] <- 'Left Pallidum'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Right.Pallidum'] <- 'Right Pallidum'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Left.Cerebellum.White.Matter'] <- 'Left Cerebellum White Matter'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Right.Cerebellum.White.Matter'] <- 'Right Cerebellum White Matter'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'X3rd.Ventricle'] <- '3rd Ventricle'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'CC_Anterior'] <- 'CC Anterior'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'X4th.Ventricle'] <- '4th Ventricle'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'CC_Mid_Anterior'] <- 'CC Mid Anterior'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Left.Caudate'] <- 'Left Caudate'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Right.Caudate'] <- 'Right Caudate'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'CC_Posterior'] <- 'CC Posterior'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'CC_Mid_Posterior'] <- 'CC Mid Posterior'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'CC_Central'] <- 'CC Central'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Left.VentralDC'] <- 'Left Ventral Diencephalon'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Right.VentralDC'] <- 'Right Ventral Diencephalon'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Brain.Stem'] <- 'Brain Stem'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Left.Putamen'] <- 'Left Putamen'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Right.Putamen'] <- 'Right Putamen'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Right.Thalamus'] <- 'Right Thalamus'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Left.Thalamus'] <- 'Left Thalamus'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Left.Hippocampus'] <- 'Left Hippocampus'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Right.Hippocampus'] <- 'Right Hippocampus'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Right.Cerebellum.Cortex'] <- 'Right Cerebellum Cortex'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Left.Cerebellum.Cortex'] <- 'Left Cerebellum Cortex'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Left.Amygdala'] <- 'Left Amygdala'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Right.Amygdala'] <- 'Right Amygdala'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Left.Accumbens.area'] <- 'Left Accumbens'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Right.Accumbens.area'] <- 'Right Accumbens'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Right.Inf.Lat.Vent'] <- 'Right Inferior Lateral Ventricle'
meta_sub_results_3$Region_clean[meta_sub_results_3$Region_clean == 'Left.Inf.Lat.Vent'] <- 'Left Inferior Lateral Ventricle'
meta_sub_results_3 <- meta_sub_results_3 %>%
mutate(Region_clean = fct_reorder(Region_clean, Estimate, .desc = TRUE))
sub_graph_2 <- ggplot(meta_sub_results_3, aes(y = Region_clean, x = -Estimate)) +
geom_errorbarh(aes(xmin = -PI_ub, xmax = -PI_lb), color = "gray70", height = 0.4, size = 1) +
geom_errorbarh(aes(xmin = -CI_ub, xmax = -CI_lb), color = "black", height = 0.2, size = 0.8) +
geom_point(aes(color = FDR < 0.05), size = 2) +
geom_vline(xintercept = 0, linetype = "twodash") +
scale_color_manual(values = c("TRUE" = "#e41a1c", "FALSE" = "black")) +
labs(
title = "Forest Plot of Effect Sizes with Prediction Intervals",
x = "Hedges' g", y = "Region",
color = "FDR < 0.05"
) +
theme_bw(base_size = 19) +
theme(
strip.text = element_text(face = "bold", size = 19),
axis.text.y = element_text(size = 19),
legend.position = "bottom"
)
sub_graph_2
length(meta_sub_results_3$Region)
ggsave("sub_meta_2.png", plot = sub_graph_2, dpi = 500, width = 15, height = 12, units = 'in', bg = 'white')
meta_sub_results_3
# atlas keys we must match
atlas_keys <- aseg$data %>% distinct(hemi, region)
# canonical midline regions in aseg
cc_regions <- c("CC anterior","CC central","CC mid anterior",
"CC mid posterior","CC posterior")
other_midline <- c("3rd ventricle","4th ventricle","brain stem",
"cerebellum cortex","cerebellum white matter")
midline_rois <- c(cc_regions, other_midline)
# --- CLEAN & NORMALIZE -------------------------------------------------
clean_regions_2 <- meta_sub_results_3 %>%
mutate(
hemi = case_when(
str_detect(Region, "^Left[._ ]") ~ "left",
str_detect(Region, "^Right[._ ]") ~ "right",
TRUE ~ NA_character_
),
region_base = Region %>%
str_remove("^Left[._ ]|^Right[._ ]") %>%
str_replace_all("[._-]", " ") %>%
str_squish() %>%
str_to_lower()
) %>%
# map to aseg spellings
mutate(
region_base = str_replace(region_base, "^third ventricle$", "3rd ventricle"),
region_base = str_replace(region_base, "^fourth ventricle$", "4th ventricle"),
region_base = str_replace(region_base, "^x?3rd ventricle$", "3rd ventricle"),
region_base = str_replace(region_base, "^x?4th ventricle$", "4th ventricle"),
region_base = str_replace(region_base, "^brainstem$", "brain stem"),
region_base = str_replace(region_base, "^thalamus$", "thalamus proper"),
region_base = str_replace(region_base, "^ventral ?dc$", "ventral DC"),
region_base = str_replace(region_base, "^(nucleus )?accumbens( area)?$", "accumbens area"),
region_base = str_replace(region_base, "^cerebel+um cortex$|^cerebel+ar cortex$", "cerebellum cortex"),
region_base = str_replace(region_base, "^cerebel+um white matter$|^cerebel+ar white matter$", "cerebellum white matter"),
region_base = str_replace(region_base, "^cc ?ant(erior)?$", "CC anterior"),
region_base = str_replace(region_base, "^cc ?cent(ral)?$", "CC central"),
region_base = str_replace(region_base, "^cc ?mid ?ant(erior)?$", "CC mid anterior"),
region_base = str_replace(region_base, "^cc ?mid ?post(erior)?$", "CC mid posterior"),
region_base = str_replace(region_base, "^cc ?post(erior)?$", "CC posterior"),
# final region string
region = region_base %>%
# ensure "CC " stays caps and the rest is exactly as needed
(\(x) ifelse(str_starts(x, "cc "), str_replace(x, "^cc ", "CC "), x))() %>%
str_squish()
) %>%
# set hemi="midline" for all midline ROIs (CC + others)
mutate(
hemi = if_else(region %in% midline_rois, "midline", hemi),
hemi = str_squish(hemi)
) %>%
select(region, hemi, Estimate)
# --- DIAGNOSE: what still fails to match? ------------------------------
still_unmatched <- clean_regions %>%
anti_join(atlas_keys, by = c("hemi","region"))
if (nrow(still_unmatched)) {
message("Rows not matching atlas (hemi, region):")
print(still_unmatched)
}
# --- PLOT --------------------------------------------------------------
L <- max(abs(clean_regions_2$Estimate), na.rm = TRUE)
ggplot() +
geom_brain(
atlas = aseg,
data = clean_regions_2,
aes(fill = -Estimate),
colour = "grey60", # outline color
linewidth = 0.2 # or size = 0.2 if on older ggplot2
) +
scale_fill_distiller(
palette = "RdBu",
limits = c(-L, L),
na.value = "grey90",
name = "Hedges g"
) +
theme_void() +
theme(
legend.position = "bottom",
legend.text = element_text(size = 8)
) +
guides(fill = guide_colorbar(barwidth = 10, barheight = 2))
sub_colors <- data.frame(Region = meta_sub_results_2$Region,
estimate = meta_sub_results_2$Estimate)
min_val <- -0.4
max_val <- 0.4
# Step 2: Define custom symmetric palette
my_palette <- colorRampPalette(c("#2569AE", "white", "#B31D2C"))
palette_n <- 100
custom_colors <- my_palette(palette_n)
# Step 3: Clip your effect sizes to stay within bounds
sub_colors$clipped_estimate <- pmax(pmin(sub_colors$estimate, max_val), min_val)
# Step 4: Convert effect sizes to indices (1 to 100)
sub_colors$scaled_index <- round(scales::rescale(sub_colors$clipped_estimate, to = c(1, palette_n), from = c(min_val, max_val)))
# Step 5: Get hex color
sub_colors$hex_code <- custom_colors[sub_colors$scaled_index]
# ✅ Done!
head(sub_colors)
RPall <- meta_sub_2[grepl("Right.Pallidum", meta_sub_2$Region), ]
LPall <- meta_sub_2[grepl("Left.Pallidum", meta_sub_2$Region), ]
Lamy <- meta_sub_2[grepl("Left.Amygdala", meta_sub_2$Region),]
Lamy_viz <- data.frame(Hedges = Lamy$yi, SE = sqrt(Lamy$vi), Study = Lamy$study,
Region = Lamy$Region)
rpall_viz <- data.frame(Hedges = RPall$yi, SE = sqrt(RPall$vi), Study = RPall$study,
Region = RPall$Region)
lpall_viz <- data.frame(Hedges = LPall$yi, SE = sqrt(LPall$vi), Study = LPall$study,
Region = LPall$Region)
rpall_viz$Study[rpall_viz$Study == 'OA'] <- 'Osteoarthritis: Tétreault, 2016'
rpall_viz$Study[rpall_viz$Study == 'FM1'] <- 'Fibromyalgia: Pando-Naude, 2019'
rpall_viz$Study[rpall_viz$Study == 'CLBP1'] <- 'Chronic Lower Back Pain: Makary, 2020'
rpall_viz$Study[rpall_viz$Study == 'migraine'] <- 'Migraine: Seminowicz, 2020'
rpall_viz$Study[rpall_viz$Study == 'PTN'] <- 'Primary Trigeminal Neuralgia: Filimonova, 2025'
rpall_viz$Study[rpall_viz$Study == 'FM2'] <- 'Balducci, 2022, Fibromyalgia'
rpall_viz$Study[rpall_viz$Study == 'CLBP_S1'] <- 'Chronic Lower Back Pain: Mano, 2018 (UK Data)'
rpall_viz$Study[rpall_viz$Study == 'CLBP_S2'] <- 'Chronic Lower Back Pain: Mano, 2018 (Japan Data)'
rpall_viz$Study[rpall_viz$Study == 'Balducci, 2022, Fibromyalgia'] <- 'Fibromyalgia: Balducci, 2022'
lpall_viz$Study[lpall_viz$Study == 'OA'] <- 'Osteoarthritis: Tétreault, 2016'
lpall_viz$Study[lpall_viz$Study == 'FM1'] <- 'Fibromyalgia: Pando-Naude, 2019'
lpall_viz$Study[lpall_viz$Study == 'CLBP1'] <- 'Chronic Lower Back Pain: Makary, 2020'
lpall_viz$Study[lpall_viz$Study == 'migraine'] <- 'Migraine: Seminowicz, 2020'
lpall_viz$Study[lpall_viz$Study == 'PTN'] <- 'Primary Trigeminal Neuralgia: Filimonova, 2025'
lpall_viz$Study[lpall_viz$Study == 'FM2'] <- 'Balducci, 2022, Fibromyalgia'
lpall_viz$Study[lpall_viz$Study == 'CLBP_S1'] <- 'Chronic Lower Back Pain: Mano, 2018 (UK Data)'
lpall_viz$Study[lpall_viz$Study == 'CLBP_S2'] <- 'Chronic Lower Back Pain: Mano, 2018 (Japan Data)'
lpall_viz$Study[lpall_viz$Study == 'Balducci, 2022, Fibromyalgia'] <- 'Fibromyalgia: Balducci, 2022'
res_lamy <- metagen(TE = Hedges,
seTE = SE,
studlab = Study,
data = Lamy_viz,
method.tau = "REML",
prediction = TRUE,
random = TRUE)
forest(res_lamy,
main = "Left Amygdala Volume",
prediction = TRUE,
col.predict = "gray50",
col.random = "navy",
xlab = "Hedges' g",
common = F,
print.common = F,
leftlabs = c("Study","Estimated effect size","Standard Error"))
res_rpall <- metagen(TE = Hedges,
seTE = SE,
studlab = Study,
data = rpall_viz,
method.tau = "REML",
prediction = TRUE,
random = TRUE)
res_lpall <- metagen(TE = Hedges,
seTE = SE,
studlab = Study,
data = lpall_viz,
method.tau = "REML",
prediction = TRUE,
random = TRUE)
forest(res_rpall,
main = "Right Pallidum Volume",
prediction = TRUE,
col.predict = "gray50",
col.random = "navy",
xlab = "Hedges' g",
common = F,
print.common = F,
leftlabs = c("Study","Estimated effect size","Standard Error"))
forest(res_lpall,
main = "Left Pallidum Volume",
prediction = TRUE,
col.predict = "gray50",
col.random = "navy",
xlab = "Hedges' g",
common = F,
print.common = F,
leftlabs = c("Study","Estimated effect size","Standard Error"))
#########################################################################
#meta with sex
#########################################################################
oasub_m_summary_2
clbpsub_m_summary_2
clbpsubs1_m_summary_2
clbpsubs2_m_summary_2
msub_m_summary_2
ptnsub_m_summary_2
oasub_m_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = oasub_m_summary_2
)
oasub_m_escalc$study <- 'OA'
clbpsub_m_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbpsub_m_summary_2
)
clbpsub_m_escalc$study <- 'CLBP1'
clbpsubs1_m_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbpsubs1_m_summary_2
)
clbpsubs1_m_escalc$study <- 'CLBP_S1'
clbpsubs2_m_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbpsubs2_m_summary_2
)
clbpsubs2_m_escalc$study <- 'CLBP_S2'
msub_m_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = msub_m_summary_2
)
msub_m_escalc$study <- 'migraine'
ptnsub_m_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = ptnsub_m_summary_2
)
ptnsub_m_escalc$study <- 'PTN'
meta_sub_m <- rbind(oasub_m_escalc, clbpsub_m_escalc,
clbpsubs1_m_escalc, clbpsubs2_m_escalc, msub_m_escalc,
ptnsub_m_escalc)
meta_sub_m_results <- do.call(rbind, lapply(split(meta_sub_m, list(meta_sub_m$Region)),
function(dfm){
if(nrow(dfm) >=2) {
res <- rma(yi = dfm$yi, vi = dfm$vi, method = 'REML')
pred <- predict(res)
data.frame(
Region = dfm$Region[1],
k = res$k,
Estimate = res$b,
SE = res$se,
zval = res$zval,
pval = res$pval,
CI_lb = res$ci.lb,
CI_ub = res$ci.ub,
PI_lb = pred$pi.lb,
PI_ub = pred$pi.ub,
I2 = res$I2,
Q = res$QE,
pval_Q = res$QEp,
t2 = res$tau2
)
} else {
NULL
}
}))
meta_sub_m_results <- meta_sub_m_results %>%
mutate(FDR = p.adjust(pval, method = 'BH'))
meta_sub_m_results <- meta_sub_m_results %>%
mutate(Q_FDR_p = p.adjust(pval_Q, method = 'BH'))
ggplot(meta_sub_m_results, aes(x = -Estimate, y = reorder(Region, Estimate))) +
geom_point(aes(color = FDR < 0.05)) +
geom_errorbarh(aes(xmax = -CI_lb, xmin = -CI_ub), height = 0.3) +
geom_vline(xintercept = 0, linetype = "solid", color = "black", linewidth = 1) + # <-- bold line at 0
theme_bw() +
labs(x = "Meta-analytic Effect Size (Hedges' g)", y = "Region", color = "Significant")
################################################################in F now
oasub_f_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = oasub_f_summary_2
)
oasub_f_escalc$study <- 'OA'
clbpsub_f_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbpsub_f_summary_2
)
clbpsub_f_escalc$study <- 'CLBP1'
clbpsubs1_f_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbpsubs1_f_summary_2
)
clbpsubs1_f_escalc$study <- 'CLBP_S1'
clbpsubs2_f_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = clbpsubs2_f_summary_2
)
clbpsubs2_f_escalc$study <- 'CLBP_S2'
msub_f_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = msub_f_summary_2
)
msub_f_escalc$study <- 'migraine'
ptnsub_f_escalc <- escalc(
measure = 'SMDH',
m1i = Mean_Group1,
sd1i = SD_Group1,
n1i = N_Group1,
m2i = Mean_Group2,
sd2i = SD_Group2,
n2i = N_Group2,
data = ptnsub_f_summary_2
)
ptnsub_f_escalc$study <- 'PTN'
fmsubescalc_2$sex <- 'female'
fm2subescalc_2$sex <- 'female'
msub_f_escalc$sex <- 'female'
meta_sub_f <- rbind(oasub_f_escalc, clbpsub_f_escalc,
clbpsubs1_f_escalc, clbpsubs2_f_escalc, msub_f_escalc,
ptnsub_f_escalc, fmsubescalc_2, fm2subescalc_2)
meta_sub_f_results <- do.call(rbind, lapply(split(meta_sub_f, list(meta_sub_f$Region)),
function(dfm){
if(nrow(dfm) >=2) {
res <- rma(yi = dfm$yi, vi = dfm$vi, method = 'REML')
pred <- predict(res)
data.frame(
Region = dfm$Region[1],
k = res$k,
Estimate = res$b,
SE = res$se,
zval = res$zval,
pval = res$pval,
CI_lb = res$ci.lb,
CI_ub = res$ci.ub,
PI_lb = pred$pi.lb,
PI_ub = pred$pi.ub,
I2 = res$I2,
Q = res$QE,
pval_Q = res$QEp,
t2 = res$tau2
)
} else {
NULL
}
}))
meta_sub_f_results <- meta_sub_f_results %>%
mutate(FDR = p.adjust(pval, method = 'BH'))
meta_sub_f_results <- meta_sub_f_results %>%
mutate(Q_FDR_p = p.adjust(pval_Q, method = 'BH'))
ggplot(meta_sub_f_results, aes(x = Estimate, y = reorder(Region, Estimate))) +
geom_point(aes(color = FDR < 0.05)) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), height = 0.3) +
geom_vline(xintercept = 0, linetype = "solid", color = "black", linewidth = 1) + # <-- bold line at 0
theme_bw() +
labs(x = "Meta-analytic Effect Size (Hedges' g)", y = "Region", color = "Significant")
########################################################################
meta_sub_m_results
meta_sub_f_results
sub_sex_result <- data.frame(Region = meta_sub_m_results$Region,
estimate_m = meta_sub_m_results$Estimate, estimate_f = meta_sub_f_results$Estimate,
SE_m = meta_sub_m_results$SE, SE_f = meta_sub_f_results$SE, CI_lb_m = meta_sub_m_results$CI_lb,
CI_ub_m = meta_sub_m_results$CI_ub, CI_lb_f = meta_sub_f_results$CI_lb, CI_ub_f = meta_sub_f_results$CI_ub)
ggplot(sub_sex_result, aes(x = estimate_m, y = estimate_f)) +
geom_point() +
theme_minimal() +
xlim(-1,1) +
ylim(-1,1) +
geom_hline(yintercept = 0, linetype = 'twodash', color = 'black') +
geom_vline(xintercept = 0, linetype = 'twodash', color = 'black') +geom_abline(slope = 1, intercept = 0, linetype = 'longdash',color = '#e43c40') +
geom_abline(slope = -1, intercept = 0, linetype = 'longdash', color = '#e43c40') +
sm_statCorr(method = 'lm', color = '#00b2a9') +
ylab("Female estimated effect size") +
xlab("Male estimated effect size")
###############################################################################
meta_sub_m
meta_sub_f
meta_sub_f_results
meta_sub_m_results
write_xlsx(meta_sub_f_results, "meta_sub_f.xlsx")
write_xlsx(meta_sub_m_results, "meta_sub_m.xlsx")
identical(meta_sub_f_results$Region, meta_sub_m_results$Region)
meta_sub_df <- data.frame(Region = meta_sub_m_results$Region, estimate_m = meta_sub_m_results$Estimate,
estimate_f = meta_sub_f_results$Estimate, SE_m = meta_sub_m_results$SE,
SE_f = meta_sub_f_results$SE, CI_l_m = meta_sub_m_results$CI_lb,
CI_u_m = meta_sub_m_results$CI_ub, CI_l_f = meta_sub_f_results$CI_lb,
CI_u_f = meta_sub_f_results$CI_ub, PI_lb_m = meta_sub_m_results$PI_lb,
PI_ub_m = meta_sub_m_results$PI_ub, PI_lb_f = meta_sub_f_results$PI_lb,
PI_ub_f = meta_sub_f_results$PI_ub)
meta_sub_df$z_diff <- (meta_sub_df$estimate_m - meta_sub_df$estimate_f) / sqrt(meta_sub_df$SE_m^2 + meta_sub_df$SE_f^2)
write_xlsx(meta_sub_df, "meta_sub_sex_analysis.xlsx")
meta_sub_df$effect_direction <- case_when(
sign(meta_sub_df$estimate_m) == sign(meta_sub_df$estimate_f) ~ "Same Direction",
sign(meta_sub_df$estimate_m) != sign(meta_sub_df$estimate_f) ~ "Opposite Direction",
TRUE ~ "Undefined"
)
meta_sub_df
meta_sub_m$sex <- 'male'
meta_sub_f$sex <- 'female'
sub_m <- meta_sub_m[c("Region","yi","vi","study","sex")]
sub_f <- meta_sub_f[c("Region","yi","vi","study","sex")]
meta_sub_sex <- rbind(sub_m, sub_f)
meta_sub_sex$sex <- factor(meta_sub_sex$sex)
meta_sub_sex$study <- factor(meta_sub_sex$study)
res_sex <- rma(yi = yi, vi = vi, mods = ~ sex, data = meta_sub_sex, method = "REML")
summary(res_sex)
res_sex_region <- rma(yi, vi, mods = ~ sex + Region, data = meta_sub_sex, method = "REML")
summary(res_sex_region)
res_3lvl <- rma.mv(yi = yi,
V = vi,
mods = ~ sex,
random = ~ 1 | study,
data = meta_sub_sex,
method = "REML")
summary(res_3lvl)
res_3lvl
group_summary_sub <- data.frame()
for (s in c("female", "male")) {
df_sub <- subset(meta_sub_sex, sex == s)
res_sub <- rma(
yi = yi,
vi = vi,
method = "REML",
data = df_sub
)
group_summary_sub <- rbind(group_summary_sub, data.frame(
Sex = s,
Effect = res_sub$b,
SE = res_sub$se,
CI_lower = res_sub$ci.lb,
CI_upper = res_sub$ci.ub,
))
}
# Combined effect across both sexes (optional)
sub_collapsed <- meta_sub_sex %>%
group_by(study, sex) %>%
summarise(
yi = sum(yi / vi) / sum(1 / vi),
vi = 1 / sum(1 / vi),
.groups = "drop"
)
res_all_sub <- rma(
yi = yi,
vi = vi,
method = "REML",
data = meta_sub_sex
)
res_sub2 <- rma(
yi = yi,
vi = vi,
method = "REML",
data = sub_collapsed
)
res_sub2_df <- rbind(res_sub2, data.frame(
Effect = res_sub2$b,
SE = res_sub2$se,
CI_lower = res_sub2$ci.lb,
CI_upper = res_sub2$ci.ub
))
group_summary_sub <- rbind(group_summary_sub, data.frame(
Sex = "both",
Effect = res_all_sub$b,
SE = res_all_sub$se,
CI_lower = res_all_sub$ci.lb,
CI_upper = res_all_sub$ci.ub
))
ggplot(group_summary_sub, aes(x = Sex, y = Effect, color = Sex, group = Sex)) +
geom_point(position = position_dodge(width = 0.4), size = 3) +
geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper),
position = position_dodge(width = 0.4), width = 0.3) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(title = "Group-level Effect Size by Sex for Subcortical Volumes",
y = "Hedges' g",
x = "Sex") +
theme_minimal() +
scale_color_manual(values = c("female" = "#00b2a9", "male" = "#e43c40", "both" = "black")) +
coord_flip() +
ylim(-.1,.25)
##################################################################
#Visualizing Meta results
##################################################################
heatmap_data2 <- meta_results_4[, c("Region", "Measurement", "Estimate")]
heatmap_data2 <- dcast(heatmap_data2, Region ~ Measurement, value.var = "Estimate")
# Convert to long format for ggplot
long_heatmap <- melt(heatmap_data2, id.vars = "Region")
# Plot
ggplot(long_heatmap, aes(x = variable, y = Region, fill = value)) +
geom_tile() +
scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0,
name = "Effect Size") +
theme_minimal() +
labs(x = "Measurement", y = "Region", title = "Effect Sizes by Region and Measurement") +
theme(axis.text.y = element_text(size = 6)) +
coord_flip() +
theme(axis.text.x = element_text(angle = 60, vjust = .5))
ggplot(meta_results_4, aes(x = Estimate, y = Measurement, fill = Measurement)) +
geom_density_ridges(alpha = 0.8, scale = 1.2) +
theme_ridges() +
theme(legend.position = "none") +
labs(title = "Distribution of Effect Sizes by Measurement",
x = "Effect Size Estimate",
y = "Measurement")
ggplot(meta_results_4, aes(x = Estimate, y = I2, color = Measurement)) +
geom_point(alpha = 0.8, size = 3) +
theme_minimal() +
labs(title = "Effect Size vs. Heterogeneity (I²)",
x = "Effect Size Estimate",
y = "I² (%)") +
geom_vline(xintercept = 0, linetype = "dashed") +
scale_color_brewer(palette = "Set1") +
facet_wrap(~Measurement, ncol = 5) +
ylim(0,100) +
geom_hline(yintercept = c(25,50,75), linetype = 'dashed', color = 'red')
ggplot(meta_results_4, aes(x = Estimate, y = I2, color = Measurement)) +
ggdist::stat_halfeye(
adjust = .5,
width = .6,
.width = 0,
justification = -.3
) +
geom_boxplot(width = .25, outlier.shape = NA) +
geom_point(alpha = 0.8, size = 3,
position = position_jitter(seed = 1, width = .1)) +
theme_minimal() +
labs(title = "Effect Size vs. Heterogeneity (I²)",
x = "Effect Size Estimate",
y = "I² (%)") +
geom_vline(xintercept = 0, linetype = "dashed") +
facet_wrap(~Measurement, ncol = 5) +
geom_hline(yintercept = 50, linetype = "dashed", color = 'red')
funnel(metagen(TE = meta_results_4$Estimate, seTE = meta_results_4$SE))
meta_gaus <- subset(meta_results_4, Measurement %in% "GausCurv")
meta_gaus$significant <- meta_gaus$FDR < 0.05
ggplot(meta_gaus, aes(x = Estimate, y = I2)) +
geom_point(aes(color = significant), alpha = .7, size = 2) +
geom_hline(yintercept = c(25,50,75), linetype = 'dashed', color = 'red') +
geom_vline(xintercept = 0, linetype = 'solid', color = 'black') +
geom_text(aes(label = ifelse(significant, Region, "")) , hjust = 1.1, vjust = 0.5, size = 3) +
labs(x = "Effect Size (Hedge's g)",
y = "Heterogeneity I2 (%)",
title = "Effect size vs heterogeneity of Gaussian curve measurements",
color = "Significant") +
theme_minimal() +
xlim(-.75,.75) +
ylim(0,100)
meta_df_4
meta_results_5 <- do.call(rbind, lapply(split(meta_df_4, list(meta_df_4$Region,
meta_df_4$Measurement)),
function(dfm){
if(nrow(dfm) >= 2) {
res <- rma(yi = dfm$yi, vi = dfm$vi, method = 'REML')
pred <- predict(res)
data.frame(
Region = dfm$Region[1],
Measurement = dfm$Measurement[1],
k = res$k,
Estimate = res$b,
SE = res$se,
zval = res$zval,
pval = res$pval,
CI_lb = res$ci.lb,
CI_ub = res$ci.ub,
PI_lb = pred$pi.lb,
PI_ub = pred$pi.ub,
I2 = res$I2,
Q = res$QE,
pval_Q = res$QEp,
t2 = res$tau2
)
} else {
NULL
}
}))
meta_results_5 <- meta_results_5 %>%
group_by(Measurement) %>%
mutate(FDR = p.adjust(pval, method = 'BH')) %>%
ungroup()
meta_results_5 <- meta_results_5 %>%
group_by(Measurement) %>%
mutate(FDR_Q = p.adjust(pval_Q, method = 'BH')) %>%
ungroup()
meta_results_5 <- meta_results_5 %>%
mutate(Region_facet = paste(Measurement, Region, sep = "__")) %>% # unique within facet
group_by(Measurement) %>%
mutate(Region_facet = fct_reorder(Region_facet, Estimate, .desc = TRUE)) %>%
ungroup()
m5v <- meta_results_5
flipped <- data.frame(Region = m5v$Region, Measurement = m5v$Measurement, k = m5v$k,
Estimate = -m5v$Estimate, SE = m5v$SE, zval = m5v$zval,
pval = m5v$pval, CI_lb = -m5v$CI_ub, CI_ub = -m5v$CI_lb,
PI_lb = -m5v$PI_ub, PI_ub = -m5v$PI_lb, I2 = m5v$I2, Q = m5v$Q,
pval_Q = m5v$pval_Q, t2 = m5v$t2, FDR = m5v$FDR, Region_facet = m5v$Region_facet,
FDR_Q = m5v$FDR_Q)
m5v$Measurement[m5v$Measurement == 'SurfArea'] <- 'Surface Area'
m5v$Measurement[m5v$Measurement == 'GrayVol'] <- 'Gray Matter Volume'
m5v$Measurement[m5v$Measurement == 'ThickAvg'] <- 'Average Thickness'
m5v$Measurement[m5v$Measurement == 'GausCurv'] <- 'Gaussian Curvature'
m5v$Measurement[m5v$Measurement == 'MeanCurv'] <- 'Mean Curvature'
flipped$Measurement[flipped$Measurement == 'SurfArea'] <- 'Surface Area'
flipped$Measurement[flipped$Measurement == 'GrayVol'] <- 'Gray Matter Volume'
flipped$Measurement[flipped$Measurement == 'ThickAvg'] <- 'Average Thickness'
flipped$Measurement[flipped$Measurement == 'GausCurv'] <- 'Gaussian Curvature'
flipped$Measurement[flipped$Measurement == 'MeanCurv'] <- 'Mean Curvature'
flipped <- flipped %>%
separate(Region, into = c("hemi","region"), sep = "_") %>%
mutate(hemi = dplyr::recode(hemi, lh = 'Left', rh = 'Right'))
regions <- flipped$region %>% str_to_lower()
miss <- setdiff(regions, atlas_regions)
matching <- sapply(miss, function(regions){
distances <- stringdist(regions, atlas_regions, method = 'jw')
atlas_regions[which.min(distances)]
})
recoding <- setNames(matching, miss)
recoding
flipped <- flipped %>%
mutate(region = str_to_lower(region),
region = ifelse(region %in% names(recoding),
recoding[region], region))
write_xlsx(flipped, "Cortical_meta_results.xlsx")
metaplot <- ggplot(m5v, aes(y = Region_facet, x = Estimate)) +
geom_errorbarh(aes(xmin = PI_lb, xmax = PI_ub), color = "gray70", height = 0.4, size = 1) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), color = "black", height = 0.2, size = 0.8) +
geom_point(aes(color = FDR < 0.05), size = 2) +
geom_vline(xintercept = 0, linetype = "twodash") +
facet_wrap(~factor(Measurement, levels = c("Surface Area","Gray Matter Volume","Average Thickness",
"Gaussian Curvature",'Mean Curvature')),
scales = "free_y", ncol = 5) +
scale_color_manual(values = c("TRUE" = "#e41a1c", "FALSE" = "black")) +
scale_y_discrete(labels = function(x) gsub(".*__", "", x)) +
labs(
title = "Forest Plot of Effect Sizes with Prediction Intervals",
x = "Hedges' g", y = "Region",
color = "FDR < 0.05"
) +
ggplot2::theme_bw(base_size = 12) +
theme(
strip.text = element_text(face = "bold", size = 12),
axis.text.y = element_text(size = 12),
legend.position = "bottom"
)
metaplot
metaplot2 <- ggplot(flipped, aes(y = Region_facet, x = Estimate)) +
geom_errorbarh(aes(xmin = PI_lb, xmax = PI_ub), color = "gray70", height = 0.4, size = 1) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), color = "black", height = 0.2, size = 0.8) +
geom_point(aes(color = FDR < 0.05), size = 2) +
geom_vline(xintercept = 0, linetype = "twodash") +
facet_wrap(~factor(Measurement, levels = c("Surface Area","Gray Matter Volume","Average Thickness",
"Gaussian Curvature",'Mean Curvature')),
scales = "free_y", ncol = 5) +
scale_color_manual(values = c("TRUE" = "#e41a1c", "FALSE" = "black")) +
scale_y_discrete(labels = function(x) gsub(".*__", "", x)) +
labs(
title = "Forest Plot of Effect Sizes with Prediction Intervals",
x = "Hedges' g", y = "Region",
color = "FDR < 0.05"
) +
ggplot2::theme_bw(base_size = 12) +
theme(
strip.text = element_text(face = "bold", size = 12),
axis.text.y = element_text(size = 12),
legend.position = "bottom"
)
metaplot2
ggsave("meta_forest.png", plot = metaplot, dpi = 500, width = 26, height = 15, units = 'in', bg = 'white')
flipped_gc <- subset(flipped, Measurement %in% 'Gaussian Curvature')
gc_long <- flipped_gc %>%
separate(Region, into = c("hemi","region"), sep = "_") %>%
mutate(hemi = dplyr::recode(hemi, lh = 'Left', rh = 'Right'))
regions <- gc_long$region %>% str_to_lower()
miss <- setdiff(regions, atlas_regions)
matching <- sapply(miss, function(regions){
distances <- stringdist(regions, atlas_regions, method = 'jw')
atlas_regions[which.min(distances)]
})
recoding <- setNames(matching, miss)
recoding
gc_long <- gc_long %>%
mutate(region = str_to_lower(region),
region = ifelse(region %in% names(recoding),
recoding[region], region))
gc_long <- gc_long %>%
mutate(region = reorder_within(region, Estimate, hemi))
gc_meta <- ggplot(gc_long, aes(y = region, x = Estimate)) +
geom_errorbarh(aes(xmin = PI_lb, xmax = PI_ub), color = "gray70", height = 0.4, size = 3) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), color = "black", height = 0.2, size = 2) +
geom_point(aes(color = FDR < 0.05), size = 5) +
geom_vline(xintercept = 0, linetype = "twodash") +
scale_color_manual(values = c("TRUE" = "#e41a1c", "FALSE" = "black")) +
scale_y_discrete(labels = function(x) gsub("__.+", "", x)) + # clean label if needed
labs(
title = "Forest Plot of Effect Sizes with Prediction Intervals",
x = "Hedges' g", y = "Region",
color = "FDR < 0.05"
) +
ggplot2::theme_bw(base_size = 28) +
theme(
strip.text = element_text(face = "bold", size = 28),
axis.text.y = element_text(size = 28),
legend.position = "bottom"
) +
facet_wrap(~hemi, scales = "free_y", ncol = 1) +
tidytext::scale_y_reordered(position = 'left')
gc_meta
ggplot(gc_long, aes(y = region, x = Estimate)) +
geom_errorbarh(aes(xmin = PI_lb, xmax = PI_ub), color = "gray70", height = 0.4, size = 3) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), color = "black", height = 0.2, size = 2) +
geom_point(aes(color = FDR < 0.05), size = 5) +
geom_vline(xintercept = 0, linetype = "twodash") +
scale_color_manual(values = c("TRUE" = "#e41a1c", "FALSE" = "black")) +
scale_y_discrete(labels = function(x) gsub("__.+", "", x)) + # clean label if needed
labs(
title = "Forest Plot of Effect Sizes with Prediction Intervals",
x = "Hedges' g", y = "Region",
color = "FDR < 0.05"
) +
ggplot2::theme_bw(base_size = 12) +
theme(
strip.text = element_text(face = "bold", size = 12),
axis.text.y = element_text(size = 12),
legend.position = "bottom"
) +
facet_wrap(~hemi, scales = "free_y", ncol = 2) +
tidytext::scale_y_reordered(position = 'right')
ggsave("gc_meta.png", plot = gc_meta, dpi = 500, width = 8, height = 26, units = 'in', bg = 'white')
flipped_mc <- subset(flipped, Measurement %in% 'Mean Curvature')
mc_long <- flipped_mc %>%
separate(Region, into = c("hemi","region"), sep = "_") %>%
mutate(hemi = dplyr::recode(hemi, lh = 'Left', rh = 'Right'))
regions <- mc_long$region %>% str_to_lower()
miss <- setdiff(regions, atlas_regions)
matching <- sapply(miss, function(regions){
distances <- stringdist(regions, atlas_regions, method = 'jw')
atlas_regions[which.min(distances)]
})
recoding <- setNames(matching, miss)
recoding
mc_long <- mc_long %>%
mutate(region = str_to_lower(region),
region = ifelse(region %in% names(recoding),
recoding[region], region))
mc_long <- mc_long %>%
mutate(region = reorder_within(region, Estimate, hemi))
mc_meta <- ggplot(mc_long, aes(y = region, x = Estimate)) +
geom_errorbarh(aes(xmin = PI_lb, xmax = PI_ub), color = "gray70", height = 0.4, size = 3) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), color = "black", height = 0.2, size = 2) +
geom_point(aes(color = FDR < 0.05), size = 5) +
geom_vline(xintercept = 0, linetype = "twodash") +
scale_color_manual(values = c("TRUE" = "#e41a1c", "FALSE" = "black")) +
scale_y_discrete(labels = function(x) gsub("__.+", "", x)) + # clean label if needed
labs(
title = "Forest Plot of Effect Sizes with Prediction Intervals",
x = "Hedges' g", y = "Region",
color = "FDR < 0.05"
) +
ggplot2::theme_bw(base_size = 28) +
theme(
strip.text = element_text(face = "bold", size = 28),
axis.text.y = element_text(size = 28),
legend.position = "bottom"
) +
facet_wrap(~hemi, scales = "free_y", ncol = 1) +
tidytext::scale_y_reordered(position = 'left')
mc_meta
ggsave("mc_meta.png", plot = mc_meta, dpi = 500, width = 8, height = 26, units = 'in', bg = 'white')
flipped_gmv <- subset(flipped, Measurement %in% 'Gray Matter Volume')
gm_long <- flipped_gmv %>%
separate(Region, into = c("hemi","region"), sep = "_") %>%
mutate(hemi = dplyr::recode(hemi, lh = 'Left', rh = 'Right'))
regions <- gm_long$region %>% str_to_lower()
miss <- setdiff(regions, atlas_regions)
matching <- sapply(miss, function(regions){
distances <- stringdist(regions, atlas_regions, method = 'jw')
atlas_regions[which.min(distances)]
})
recoding <- setNames(matching, miss)
recoding
gm_long <- gm_long %>%
mutate(region = str_to_lower(region),
region = ifelse(region %in% names(recoding),
recoding[region], region))
gm_long <- gm_long %>%
mutate(region = reorder_within(region, Estimate, hemi))
gmv_meta <- ggplot(gm_long, aes(y = region, x = Estimate)) +
geom_errorbarh(aes(xmin = PI_lb, xmax = PI_ub), color = "gray70", height = 0.4, size = 3) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), color = "black", height = 0.2, size = 2) +
geom_point(aes(color = FDR < 0.05), size = 5) +
geom_vline(xintercept = 0, linetype = "twodash") +
scale_color_manual(values = c("TRUE" = "#e41a1c", "FALSE" = "black")) +
scale_y_discrete(labels = function(x) gsub("__.+", "", x)) + # clean label if needed
labs(
title = "Forest Plot of Effect Sizes with Prediction Intervals",
x = "Hedges' g", y = "Region",
color = "FDR < 0.05"
) +
ggplot2::theme_bw(base_size = 28) +
theme(
strip.text = element_text(face = "bold", size = 28),
axis.text.y = element_text(size = 28),
legend.position = "bottom"
) +
facet_wrap(~hemi, scales = "free_y", ncol = 1) +
tidytext::scale_y_reordered(position = 'left')
gmv_meta
ggsave("gmv_meta.png", plot = gmv_meta, dpi = 500, width = 8, height = 26, units = 'in', bg = 'white')
flipped_sa <- subset(flipped, Measurement %in% 'Surface Area')
sa_long <- flipped_sa %>%
separate(Region, into = c("hemi","region"), sep = "_") %>%
mutate(hemi = dplyr::recode(hemi, lh = 'Left', rh = 'Right'))
regions <- sa_long$region %>% str_to_lower()
miss <- setdiff(regions, atlas_regions)
matching <- sapply(miss, function(regions){
distances <- stringdist(regions, atlas_regions, method = 'jw')
atlas_regions[which.min(distances)]
})
recoding <- setNames(matching, miss)
recoding
sa_long <- sa_long %>%
mutate(region = str_to_lower(region),
region = ifelse(region %in% names(recoding),
recoding[region], region))
sa_long <- sa_long %>%
mutate(region = reorder_within(region, Estimate, hemi))
sa_meta <- ggplot(sa_long, aes(y = region, x = Estimate)) +
geom_errorbarh(aes(xmin = PI_lb, xmax = PI_ub), color = "gray70", height = 0.4, size = 3) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), color = "black", height = 0.2, size = 2) +
geom_point(aes(color = FDR < 0.05), size = 5) +
geom_vline(xintercept = 0, linetype = "twodash") +
scale_color_manual(values = c("TRUE" = "#e41a1c", "FALSE" = "black")) +
scale_y_discrete(labels = function(x) gsub("__.+", "", x)) + # clean label if needed
labs(
title = "Forest Plot of Effect Sizes with Prediction Intervals",
x = "Hedges' g", y = "Region",
color = "FDR < 0.05"
) +
ggplot2::theme_bw(base_size = 28) +
theme(
strip.text = element_text(face = "bold", size = 28),
axis.text.y = element_text(size = 28),
legend.position = "bottom"
) +
facet_wrap(~hemi, scales = "free_y", ncol = 1) +
tidytext::scale_y_reordered(position = 'right')
sa_meta
ggsave("gmv_meta.png", plot = gmv_meta, dpi = 500, width = 8, height = 26, units = 'in', bg = 'white')
flipped_ct <- subset(flipped, Measurement %in% 'Average Thickness')
ct_long <- flipped_ct %>%
separate(Region, into = c("hemi","region"), sep = "_") %>%
mutate(hemi = dplyr::recode(hemi, lh = 'Left', rh = 'Right'))
regions <- ct_long$region %>% str_to_lower()
miss <- setdiff(regions, atlas_regions)
matching <- sapply(miss, function(regions){
distances <- stringdist(regions, atlas_regions, method = 'jw')
atlas_regions[which.min(distances)]
})
recoding <- setNames(matching, miss)
recoding
ct_long <- ct_long %>%
mutate(region = str_to_lower(region),
region = ifelse(region %in% names(recoding),
recoding[region], region))
ct_long <- ct_long %>%
mutate(region = reorder_within(region, Estimate, hemi))
ct_meta <- ggplot(ct_long, aes(y = region, x = Estimate)) +
geom_errorbarh(aes(xmin = PI_lb, xmax = PI_ub), color = "gray70", height = 0.4, size = 3) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), color = "black", height = 0.2, size = 2) +
geom_point(aes(color = FDR < 0.05), size = 5) +
geom_vline(xintercept = 0, linetype = "twodash") +
scale_color_manual(values = c("TRUE" = "#e41a1c", "FALSE" = "black")) +
scale_y_discrete(labels = function(x) gsub("__.+", "", x)) + # clean label if needed
labs(
title = "Forest Plot of Effect Sizes with Prediction Intervals",
x = "Hedges' g", y = "Region",
color = "FDR < 0.05"
) +
ggplot2::theme_bw(base_size = 28) +
theme(
strip.text = element_text(face = "bold", size = 28),
axis.text.y = element_text(size = 28),
legend.position = "bottom"
) +
facet_wrap(~hemi, scales = "free_y", ncol = 1) +
tidytext::scale_y_reordered(position = 'left')
ct_meta
ggsave("ct_meta.png", plot = ct_meta, dpi = 500, width = 8, height = 26, units = 'in', bg = 'white')
left_ent <- ent_test %>%
filter(grepl("Left", Region))
gmv_ent <- m5v_gmv %>%
filter(grepl("entorhinal", Region))
gmv_ent$Region[gmv_ent$Region == 'lh_entorhinal'] <- 'Left Entorhinal Cortex'
gmv_ent$Region[gmv_ent$Region == 'rh_entorhinal'] <- 'Right Entorhinal Cortex'
ent_plot <- ggplot(gmv_ent, aes(y = Region, x = Estimate)) +
geom_errorbarh(aes(xmin = PI_lb, xmax = PI_ub), color = "gray70", height = 0.4, size = 3) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), color = "black", height = 0.2, size = 2) +
geom_point(aes(color = FDR < 0.05), size = 5) +
geom_vline(xintercept = 0, linetype = "twodash", size = 1) +
scale_color_manual(values = c("TRUE" = "#e41a1c", "FALSE" = "black")) +
scale_y_discrete(labels = function(x) gsub(".*__", "", x)) +
labs(
title = "Forest Plot of Effect Sizes with Prediction Intervals",
x = "Hedges' g", y = "Region",
color = "FDR < 0.05"
) +
ggplot2::theme_bw(base_size = 20) +
theme(
strip.text = element_text(face = "bold", size = 20),
axis.text.y = element_text(size = 20),
legend.position = "bottom",
legend.background = element_rect(fill = "transparent", color = NA),
panel.background = element_rect(fill = "transparent", color = NA),
plot.background = element_rect(fill = "transparent", color = NA),
panel.grid = element_blank()
) +
xlim(-.3,.6)
ent_plot
ggsave(
"forest_plot_transparent.png",
plot = ent_plot,
width = 10,
height = 4,
dpi = 300,
bg = "transparent"
)
ent_test <- entorhinal_4_vol_vis
ent_test$Study[ent_test$Study == 'Tétreault, 2016, Osteoarthritis'] <- 'Osteoarthritis: Tétreault, 2016'
ent_test$Study[ent_test$Study == 'Pando-Naude, 2019, Fibromyalgia'] <- 'Fibromyalgia: Pando-Naude, 2019'
ent_test$Study[ent_test$Study == 'Makary, 2020, Chronic Lower Back Pain'] <- 'Chronic Lower Back Pain: Makary, 2020'
ent_test$Study[ent_test$Study == 'Seminowicz, 2020, Migraine'] <- 'Migraine: Seminowicz, 2020'
ent_test$Study[ent_test$Study == 'Filimonova, 2025, Primary Trigeminal Neuralgia'] <- 'Primary Trigeminal Neuralgia: Filimonova, 2025'
ent_test$Study[ent_test$Study == 'Balducci, 2022, Fibromyalgia'] <- 'Balducci, 2022, Fibromyalgia'
ent_test$Study[ent_test$Study == 'Mano, 2018, Chronic Lower Back Pain (UK Data)'] <- 'Chronic Lower Back Pain: Mano, 2018 (UK Data)'
ent_test$Study[ent_test$Study == 'Mano, 2018, Chronic Lower Back Pain (Japan Data)'] <- 'Chronic Lower Back Pain: Mano, 2018 (Japan Data)'
ent_test$Study[ent_test$Study == 'Balducci, 2022, Fibromyalgia'] <- 'Fibromyalgia: Balducci, 2022'
left_ent <- ent_test %>%
filter(grepl("Left", Region))
right_ent <- ent_test %>%
filter(grepl("Right", Region))
left_ent <- left_ent %>%
separate(Study, into = c("Condition","Author"), sep = ":")
right_ent <- right_ent %>%
separate(Study, into = c("Condition","Author"), sep = ":")
res_left <- metagen(TE = Hedges,
seTE = SE,
studlab = Author,
data = left_ent,
method.tau = "REML",
prediction = TRUE,
random = TRUE)
res_right <- metagen(TE = Hedges,
seTE = SE,
studlab = Author,
data = right_ent,
method.tau = "REML",
prediction = TRUE,
random = TRUE)
forest(res_left,
prediction = TRUE,
col.predict = "gray50",
col.random = "navy",
xlab = "Hedges' g",
common = FALSE,
print.common = FALSE,
leftcols = c("Condition","studlab", "TE", "seTE"),
leftlabs = c("Condition","Author", "Effect Size", "SE"),
sortvar = Condition)
forest(res_right,
prediction = TRUE,
col.predict = "gray50",
col.random = "navy",
xlab = "Hedges' g",
common = FALSE,
print.common = FALSE,
leftcols = c("Condition","studlab", "TE", "seTE"),
leftlabs = c("Condition","Author", "Effect Size", "SE"),
sortvar = Condition)
#Visualize lh precentral gyrus gaus curv
precentral <- meta_df_4[grepl("precentral", meta_df_4$Region), ]
precentral <- subset(precentral, Measurement %in% 'GausCurv')
precentral_viz <- data.frame(Hedges = precentral$yi, SE = sqrt(precentral$vi),
Measurement = precentral$Measurement, Study = precentral$study,
Region = precentral$Region)
left_precentral <- precentral_viz %>%
filter(grepl("lh", Region))
right_precentral <- precentral_viz %>%
filter(grepl("rh", Region))
res_lhpre <- metagen(TE = Hedges,
seTE = SE,
studlab = Study,
data = left_precentral,
method.tau = "REML",
prediction = TRUE,
random = TRUE)
res_rhpre <- metagen(TE = Hedges,
seTE = SE,
studlab = Study,
data = right_precentral,
method.tau = "REML",
prediction = TRUE,
random = TRUE)
forest(res_lhpre,
main = "Left Precentral Gyrus Gaussian Curve",
prediction = TRUE,
col.predict = "gray50",
col.random = "navy",
xlab = "Hedges' g",
common = F,
print.common = F,
leftlabs = c("Study","Estimated effect size","Standard Error"))
forest(res_rhpre,
main = "Right Precentral Gyrus Gaussian Curve",
prediction = TRUE,
col.predict = "gray50",
col.random = "navy",
xlab = "Hedges' g",
common = F,
print.common = F,
leftlabs = c("Study","Estimated effect size","Standard Error"))
subsetting_meta5 <- meta_results_5
subsetting_meta5$Measurement[subsetting_meta5$Measurement == 'SurfArea'] <- 'Surface Area'
subsetting_meta5$Measurement[subsetting_meta5$Measurement == 'GrayVol'] <- 'Gray Matter Volume'
subsetting_meta5$Measurement[subsetting_meta5$Measurement == 'ThickAvg'] <- 'Cortical Thickness'
subsetting_meta5$Measurement[subsetting_meta5$Measurement == 'GausCurv'] <- 'Gaussian Curvature'
subsetting_meta5$Measurement[subsetting_meta5$Measurement == 'MeanCurv'] <- 'Mean Curvature'
ggplot(subsetting_meta5, aes(x = I2, y = Estimate, size = t2, color = FDR_Q < 0.05)) +
geom_point(alpha = 0.8) +
labs(title = "Effect Size vs. Heterogeneity",
x = "I² (%)", y = "Pooled Effect Size") +
scale_color_manual(values = c("gray", "red"), name = "Significant Heterogeneity") +
facet_wrap(~factor(Measurement, levels = c("Surface Area","Gray Matter Volume","Cortical Thickness",
"Mean Curvature",'Gaussian Curvature')), ncol = 1) +
theme_minimal() +
ylim(-.5,.5) +
xlim(0,100) +
geom_hline(yintercept = 0, linetype = 'solid', color = 'black') +
geom_vline(xintercept = c(25,50,75), linetype = 'twodash',color = 'red') +
geom_text(aes(label = ifelse(FDR_Q < 0.05, Region, "")), hjust = 1.4, vjust = 0.5, size = 3)
ggplot(meta_results_5, aes(x = Estimate, y = -log10(FDR_Q), color = FDR_Q < 0.05)) +
geom_point() +
geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
labs(title = "Effect Size vs. Heterogeneity Significance",
x = "Effect Estimate", y = "-log10(FDR_Q)")+
facet_wrap(~Measurement, ncol = 5) +
theme_minimal()
ggplot(meta_results_5, aes(x = t2)) +
geom_histogram(bins = 20, fill = "steelblue", color = "white") +
labs(title = "Distribution of Between-Study Variance (τ²)") +
facet_wrap(~Measurement) +
theme_minimal()
testing <- subset(meta_results_5, Measurement %in% "ThickAvg")
densityPlot(testing$t2)
#########################################################################
#Group level summary again
#########################################################################
meta_sex
rma(yi, vi, mods = ~ sex + Region, data = meta_sex, method = "REML")
collapsed_by_sex <- meta_sex %>%
group_by(study, sex, Measurement) %>%
summarise(
yi = sum(yi / vi) / sum(1 / vi),
vi = 1 / sum(1 / vi),
.groups = "drop"
)
collapsed_combined <- collapsed_by_sex %>%
group_by(study, Measurement) %>%
summarise(
yi = if(n() == 2) sum(yi / vi) / sum(1 / vi) else yi,
vi = if(n() == 2) 1 / sum(1 / vi) else vi,
.groups = "drop"
)
group_summary_3 <- data.frame()
for (m in measurements) {
for (s in c("female", "male")) {
df_sub <- subset(collapsed_meta, Measurement == m & sex == s)
res_mes <- rma(yi = yi, vi = vi, method = "REML", data = df_sub)
pred <- predict(res_mes)
group_summary_3 <- rbind(group_summary_3, data.frame(
Measurement = m,
Sex = s,
Effect = res_mes$b,
SE = res_mes$se,
CI_lower = res_mes$ci.lb,
CI_upper = res_mes$ci.ub,
PI_lb = pred$pi.lb,
PI_ub = pred$pi.ub,
tau2 = res_mes$tau2,
k = res_mes$k,
I2 = res_mes$I2,
Q = res_mes$QE,
pval_Q = res_mes$QEp
))
}
# Combined (both sexes)
df_all <- subset(collapsed_combined, Measurement == m)
res_all <- rma(yi = yi, vi = vi, method = "REML", data = df_all)
pred_all <- predict(res_all)
group_summary_3 <- rbind(group_summary_3, data.frame(
Measurement = m,
Sex = "both",
Effect = res_all$b,
SE = res_all$se,
CI_lower = res_all$ci.lb,
CI_upper = res_all$ci.ub,
PI_lb = pred_all$pi.lb,
PI_ub = pred_all$pi.ub,
tau2 = res_all$tau2,
k = res_all$k,
I2 = res_all$I2,
Q = res_all$QE,
pval_Q = res_all$QEp
))
}
group_summary_3$Measurement[group_summary_3$Measurement == 'SurfArea'] <- 'Surface Area'
group_summary_3$Measurement[group_summary_3$Measurement == 'GrayVol'] <- 'Gray Matter Volume'
group_summary_3$Measurement[group_summary_3$Measurement == 'ThickAvg'] <- 'Cortical Thickness'
group_summary_3$Measurement[group_summary_3$Measurement == 'GausCurv'] <- 'Gaussian Curvature'
group_summary_3$Measurement[group_summary_3$Measurement == 'MeanCurv'] <- 'Mean Curvature'
group_summary_3$Measurement <- factor(group_summary_3$Measurement,
levels = c('Gaussian Curvature','Mean Curvature','Surface Area',
'Cortical Thickness','Gray Matter Volume'))
pfaff <- ggplot(group_summary_3, aes(x = Measurement, y = Effect, color = Sex, group = Sex)) +
geom_point(position = position_dodge(width = 0.7), size = 3) +
geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper),
position = position_dodge(width = 0.7), width = 0.6) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(title = "Group-level Effect Size by Measurement and Sex",
y = "Hedges' g",
x = "Measurement Type") +
scale_color_manual(values = c("female" = "#00b2a9", "male" = "#e43c40", "both" = "black")) +
coord_flip() +
ggplot2::theme_bw(base_size = 16) +
theme(
strip.text = element_text(face = "bold", size = 16),
axis.text.y = element_text(size = 16),
legend.position = "bottom"
)
pfaff
ggsave("meta_sex_groupestimate.png", plot = pfaff, dpi = 500, width = 15, height = 6, units = 'in', bg = 'white')
meta_sex$sex <- factor(meta_sex$sex)
meta_sex$Region <- factor(meta_sex$Region)
meta_function <- function(data) {
rma.mv(yi, vi, random = ~1 | study/Region, method = 'REML', data = data)
}
meta_test_results <- meta_sex %>%
group_by(Measurement, sex) %>%
filter(n() >= 2) %>% # only analyze if you have enough data
group_modify(~ {
tryCatch({
model <- rma.mv(yi, vi, random = ~1 | study/Region, method = "REML", data = .x)
tibble(g = model$b[1], se = model$se, ci.lb = model$ci.lb, ci.ub = model$ci.ub)
}, error = function(e) {
tibble(g = NA, se = NA, ci.lb = NA, ci.ub = NA)
})
}) %>%
ungroup()
meta_test_results$Measurement[meta_test_results$Measurement == 'SurfArea'] <- 'Surface Area'
meta_test_results$Measurement[meta_test_results$Measurement == 'GrayVol'] <- 'Gray Matter Volume'
meta_test_results$Measurement[meta_test_results$Measurement == 'ThickAvg'] <- 'Cortical Thickness'
meta_test_results$Measurement[meta_test_results$Measurement == 'GausCurv'] <- 'Gaussian Curvature'
meta_test_results$Measurement[meta_test_results$Measurement == 'MeanCurv'] <- 'Mean Curvature'
meta_test_results$Measurement <- factor(meta_test_results$Measurement,
levels = c('Gaussian Curvature','Mean Curvature','Cortical Thickness',
'Surface Area','Gray Matter Volume'))
pfaff <- ggplot(meta_test_results, aes(x = g, y = Measurement, color = sex)) +
geom_point(position = position_dodge(width = 0.5), size = 3) +
geom_errorbar(aes(xmin = ci.lb, xmax = ci.ub),
position = position_dodge(width = 0.5), width = 0.2) +
geom_vline(xintercept = 0, linetype = "dashed") +
scale_color_manual(values = c("red", "turquoise3")) +
labs(
x = "Hedges' g",
y = "Measurement Type",
title = "Group-level Effect Size by Measurement and Sex"
) +
theme_minimal(base_size = 14)
pfaff
############################################################################
#Redo more shit
############################################################################
measurement_types <- unique(meta_sex$Measurement)
# Create empty list to store results
results_list <- list()
# Loop
for (m in measurement_types) {
cat("\n\n-------------------\nAnalyzing:", m, "\n-------------------\n")
# Subset data
sub_df <- subset(meta_sex, Measurement == m)
# Run 3-level meta-regression with sex as moderator, random intercept by study
model <- rma.mv(yi = yi,
V = vi,
mods = ~ sex,
random = ~ 1 | study,
data = sub_df,
method = "REML")
print(summary(model))
# Optionally store result
results_list[[as.character(m)]] <- model
}
######################################################################
oa_vol_summary
fm_vol_summary
fm2_vol_summary
clbp_vol_summary
clbps1_vol_summary
clbps2_vol_summary
mvol_summary
pvol_summary
oa_vol_summary$study <- 'OA'
fm_vol_summary$study <- 'FM'
fm2_vol_summary$study <- 'FM2'
clbp_vol_summary$study <- 'CLBP'
clbps1_vol_summary$study <- 'CLBP_S1'
clbps2_vol_summary$study <- 'CLBP_S2'
mvol_summary$study <- 'migraine'
pvol_summary$study <- 'PTN'
meta_vol_df <- rbind(oa_vol_summary, fm_vol_summary, fm2_vol_summary,
clbp_vol_summary, clbps1_vol_summary, clbps2_vol_summary,
mvol_summary, pvol_summary)
meta_vol_df$yi <- meta_vol_df$Hedges_g
meta_vol_df$vi <- meta_vol_df$SE_Hedges_g^2
meta_vol_results <- do.call(rbind, lapply(split(meta_vol_df, meta_vol_df$Region),
function(dfm){
if(nrow(dfm) >= 2) {
res <- rma(yi = dfm$yi, vi = dfm$vi, method = 'REML')
pred <- predict(res)
data.frame(
Region = dfm$Region[1],
k = res$k,
Estimate = res$b,
SE = res$se,
zval = res$zval,
pval = res$pval,
CI_lb = res$ci.lb,
CI_ub = res$ci.ub,
PI_lb = pred$pi.lb,
PI_ub = pred$pi.ub,
I2 = res$I2,
Q = res$QE,
pval_Q = res$QEp,
t2 = res$tau2
)
} else {
NULL
}
}))
meta_vol_results <- meta_vol_results %>%
mutate(FDR = p.adjust(pval, method = 'BH'))
oasub_summary
fm1sub_summary
fm2sub_summary
clbpsub_summary
clbpsubs1_summary
clbpsubs2_summary
msub_summary
ptnsub_summary
oasub_summary$study <- 'OA'
fm1sub_summary$study <- 'FM'
fm2sub_summary$study <- 'FM2'
clbpsub_summary$study <- 'CLBP'
clbpsubs1_summary$study <- 'CLBP_S1'
clbpsubs2_summary$study <- 'CLBP_S2'
msub_summary$study <- 'migraine'
ptnsub_summary$study <- 'PTN'
meta_sub_df <- rbind(oasub_summary, fm1sub_summary, fm2sub_summary,
clbpsub_summary, clbpsubs1_summary, clbpsubs2_summary,
msub_summary, ptnsub_summary)
meta_sub_df$yi <- meta_sub_df$Hedges_g
meta_sub_df$vi <- meta_sub_df$SE_Hedges_g^2
meta_sub_results <- do.call(rbind, lapply(split(meta_sub_df, meta_sub_df$Region),
function(dfm){
if(nrow(dfm) >= 2) {
res <- rma(yi = dfm$yi, vi = dfm$vi, method = 'REML')
pred <- predict(res)
data.frame(
Region = dfm$Region[1],
k = res$k,
Estimate = res$b,
SE = res$se,
zval = res$zval,
pval = res$pval,
CI_lb = res$ci.lb,
CI_ub = res$ci.ub,
PI_lb = pred$pi.lb,
PI_ub = pred$pi.ub,
I2 = res$I2,
Q = res$QE,
pval_Q = res$QEp,
t2 = res$tau2
)
} else {
NULL
}
}))
meta_sub_results <- meta_sub_results %>%
mutate(FDR = p.adjust(pval, method = 'BH'))
meta_sub_results <- meta_sub_results %>%
mutate(FDR_Q = p.adjust(pval_Q, method = 'BH'))
meta_sub_results$Region_clean <- as.character(meta_sub_results$Region)
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Right.Lateral.Ventricle'] <- 'Right Lateral Ventricle'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Left.Lateral.Ventricle'] <- 'Left Lateral Ventricle'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Left.Pallidum'] <- 'Left Pallidum'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Right.Pallidum'] <- 'Right Pallidum'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Left.Cerebellum.White.Matter'] <- 'Left Cerebellum White Matter'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Right.Cerebellum.White.Matter'] <- 'Right Cerebellum White Matter'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'X3rd.Ventricle'] <- '3rd Ventricle'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'CC_Anterior'] <- 'CC Anterior'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'X4th.Ventricle'] <- '4th Ventricle'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'CC_Mid_Anterior'] <- 'CC Mid Anterior'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Left.Caudate'] <- 'Left Caudate'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Right.Caudate'] <- 'Right Caudate'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'CC_Posterior'] <- 'CC Posterior'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'CC_Mid_Posterior'] <- 'CC Mid Posterior'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'CC_Central'] <- 'CC Central'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Left.VentralDC'] <- 'Left Ventral Diencephalon'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Right.VentralDC'] <- 'Right Ventral Diencephalon'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Brain.Stem'] <- 'Brain Stem'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Left.Putamen'] <- 'Left Putamen'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Right.Putamen'] <- 'Right Putamen'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Right.Thalamus'] <- 'Right Thalamus'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Left.Thalamus'] <- 'Left Thalamus'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Left.Hippocampus'] <- 'Left Hippocampus'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Right.Hippocampus'] <- 'Right Hippocampus'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Right.Cerebellum.Cortex'] <- 'Right Cerebellum Cortex'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Left.Cerebellum.Cortex'] <- 'Left Cerebellum Cortex'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Left.Amygdala'] <- 'Left Amygdala'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Right.Amygdala'] <- 'Right Amygdala'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Right.Cerebellum.Cortex'] <- 'Right Cerebellum Cortex'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Left.Cerebellum.Cortex'] <- 'Left Cerebellum Cortex'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Left.Amygdala'] <- 'Left Amygdala'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Right.Amygdala'] <- 'Right Amygdala'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Left.Inf.Lat.Vent'] <- 'Left Inferior Lateral Ventricle'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Right.Inf.Lat.Vent'] <- 'Right Inferior Lateral Ventricle'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Left.Accumbens.area'] <- 'Left Accumbens'
meta_sub_results$Region_clean[meta_sub_results$Region_clean == 'Right.Accumbens.area'] <- 'Right Accumbens'
meta_sub_results <- meta_sub_results %>%
mutate(Region_clean = fct_reorder(Region_clean, Estimate, .desc = TRUE))
sub_graph_adjusted <- ggplot(meta_sub_results, aes(y = Region_clean, x = -Estimate)) +
geom_errorbarh(aes(xmin = -PI_ub, xmax = -PI_lb), color = "gray70", height = 0.4, size = 1) +
geom_errorbarh(aes(xmin = -CI_ub, xmax = -CI_lb), color = "black", height = 0.2, size = 0.8) +
geom_point(aes(color = FDR < 0.05), size = 2) +
geom_vline(xintercept = 0, linetype = "twodash") +
scale_color_manual(values = c("TRUE" = "#e41a1c", "FALSE" = "black")) +
labs(
title = "Forest Plot of Effect Sizes with Prediction Intervals",
x = "Hedges' g", y = "Region",
color = "FDR < 0.05"
) +
theme_bw(base_size = 24) +
theme(
strip.text = element_text(face = "bold", size = 24),
axis.text.y = element_text(size = 24),
legend.position = "bottom"
) +
tidytext::scale_y_reordered(position = 'right')
sub_graph_adjusted
ggsave("sub_meta_adjusted.png", plot = sub_graph_adjusted, dpi = 500, width = 9, height = 12, units = 'in', bg = 'white')
identical(meta_sub_results$Region, meta_sub_results_462$Region)
meta_diff <- data.frame(Region = meta_sub_results$Region_clean, unadjusted_estimate = meta_sub_results_462$Estimate,
adjusted = meta_sub_results$Estimate)
meta_diff$difference <- meta_diff$unadjusted_estimate - meta_diff$adjusted
meta_diff$Region_name <- meta_sub_results$Region
ilovegraphs <- ggplot(meta_diff, aes(x = difference, y = Region)) +
geom_point(size = 5) +
theme_bw() +
geom_vline(xintercept = 0, linetype = 'dotdash') +
xlim(-0.02, 0.02) +
ggplot2::theme_bw(base_size = 28) +
theme(
strip.text = element_text(face = "bold", size = 28),
axis.text.y = element_text(size = 28),
legend.position = "bottom"
)
ilovegraphs
ggsave("submeta_diff.png", plot = ilovegraphs, dpi = 500, width = 9, height = 12, units = 'in', bg = 'white')
clean_regions <- meta_diff %>%
mutate(
hemi = case_when(
str_detect(Region_name, "^Left[._ ]") ~ "left",
str_detect(Region_name, "^Right[._ ]") ~ "right",
TRUE ~ NA_character_
),
region_base = Region_name %>%
str_remove("^Left[._ ]|^Right[._ ]") %>%
str_replace_all("[._-]", " ") %>%
str_squish() %>%
str_to_lower()
) %>%
# map to aseg spellings
mutate(
region_base = str_replace(region_base, "^third ventricle$", "3rd ventricle"),
region_base = str_replace(region_base, "^fourth ventricle$", "4th ventricle"),
region_base = str_replace(region_base, "^x?3rd ventricle$", "3rd ventricle"),
region_base = str_replace(region_base, "^x?4th ventricle$", "4th ventricle"),
region_base = str_replace(region_base, "^brainstem$", "brain stem"),
region_base = str_replace(region_base, "^thalamus$", "thalamus proper"),
region_base = str_replace(region_base, "^ventral ?dc$", "ventral DC"),
region_base = str_replace(region_base, "^(nucleus )?accumbens( area)?$", "accumbens area"),
region_base = str_replace(region_base, "^cerebel+um cortex$|^cerebel+ar cortex$", "cerebellum cortex"),
region_base = str_replace(region_base, "^cerebel+um white matter$|^cerebel+ar white matter$", "cerebellum white matter"),
region_base = str_replace(region_base, "^cc ?ant(erior)?$", "CC anterior"),
region_base = str_replace(region_base, "^cc ?cent(ral)?$", "CC central"),
region_base = str_replace(region_base, "^cc ?mid ?ant(erior)?$", "CC mid anterior"),
region_base = str_replace(region_base, "^cc ?mid ?post(erior)?$", "CC mid posterior"),
region_base = str_replace(region_base, "^cc ?post(erior)?$", "CC posterior"),
# final region string
region = region_base %>%
# ensure "CC " stays caps and the rest is exactly as needed
(\(x) ifelse(str_starts(x, "cc "), str_replace(x, "^cc ", "CC "), x))() %>%
str_squish()
) %>%
# set hemi="midline" for all midline ROIs (CC + others)
mutate(
hemi = if_else(region %in% midline_rois, "midline", hemi),
hemi = str_squish(hemi)
) %>%
select(region, hemi, difference)
# --- DIAGNOSE: what still fails to match? ------------------------------
still_unmatched <- clean_regions %>%
anti_join(atlas_keys, by = c("hemi","region"))
if (nrow(still_unmatched)) {
message("Rows not matching atlas (hemi, region):")
print(still_unmatched)
}
# --- PLOT --------------------------------------------------------------
L <- max(abs(clean_regions$difference), na.rm = TRUE)
ggplot() +
geom_brain(
atlas = aseg,
data = clean_regions,
aes(fill = difference),
colour = "grey60" # or size = 0.2 if on older ggplot2
) +
scale_fill_distiller(
palette = "PRGn",
limits = c(-L,L),
na.value = "grey90"
) +
theme(
legend.position = "bottom",
legend.text = element_text(size = 8)
) +
guides(fill = guide_colorbar(barwidth = 2, barheight = 10)) +
theme_void()
head(meta_vol_results)
metavol_flip <- data.frame(Region = meta_vol_results$Region, k = meta_vol_results$k,
Estimate = -meta_vol_results$Estimate, SE = meta_vol_results$SE,
zval = meta_vol_results$zval, pval = meta_vol_results$pval,
CI_lb = -meta_vol_results$CI_ub, CI_ub = -meta_vol_results$CI_lb,
PI_lb = -meta_vol_results$PI_ub, PI_ub = -meta_vol_results$PI_lb,
I2 = meta_vol_results$I2, Q = meta_vol_results$Q,
pval_Q = meta_vol_results$pval_Q, t2 = meta_vol_results$t2,
FDR = meta_vol_results$FDR)
metavol_flip <- metavol_flip %>%
mutate(FDR_Q = p.adjust(pval_Q, method = 'BH'))
gm_long2 <- metavol_flip %>%
separate(Region, into = c("hemi","region","measurement"), sep = "_") %>%
mutate(hemi = dplyr::recode(hemi, lh = 'Left', rh = 'Right'))
regions <- gm_long2$region %>% str_to_lower()
miss <- setdiff(regions, atlas_regions)
matching <- sapply(miss, function(regions){
distances <- stringdist(regions, atlas_regions, method = 'jw')
atlas_regions[which.min(distances)]
})
recoding <- setNames(matching, miss)
recoding
gm_long2 <- gm_long2 %>%
mutate(region = str_to_lower(region),
region = ifelse(region %in% names(recoding),
recoding[region], region))
gm_long2 <- gm_long2 %>%
mutate(region = reorder_within(region, Estimate, hemi))
gmv_meta2 <- ggplot(gm_long2, aes(y = region, x = Estimate)) +
geom_errorbarh(aes(xmin = PI_lb, xmax = PI_ub), color = "gray70", height = 0.4, size = 3) +
geom_errorbarh(aes(xmin = CI_lb, xmax = CI_ub), color = "black", height = 0.2, size = 2) +
geom_point(aes(color = FDR < 0.05), size = 5) +
geom_vline(xintercept = 0, linetype = "twodash") +
scale_color_manual(values = c("TRUE" = "#e41a1c", "FALSE" = "black")) +
scale_y_discrete(labels = function(x) gsub("__.+", "", x)) + # clean label if needed
labs(
title = "Forest Plot of Effect Sizes with Prediction Intervals",
x = "Hedges' g", y = "Region",
color = "FDR < 0.05"
) +
ggplot2::theme_bw(base_size = 28) +
theme(
strip.text = element_text(face = "bold", size = 28),
axis.text.y = element_text(size = 28),
legend.position = "bottom"
) +
facet_wrap(~hemi, scales = "free_y", ncol = 1) +
tidytext::scale_y_reordered(position = 'right')
gmv_meta2
ggsave("gmv_meta2.png", plot = gmv_meta2, dpi = 500, width = 8, height = 26, units = 'in', bg = 'white')
gm_long2$unadjusted_estimate <- flipped_gmv$Estimate
gm_long2$difference <- gm_long2$unadjusted_estimate - gm_long2$Estimate
gm_long2$percent_reduction <- abs((gm_long2$unadjusted_estimate - gm_long2$Estimate))/abs(gm_long2$unadjusted_estimate)
head(gm_long2)
gm_figure <- data.frame(hemi = gm_long2$hemi, region = flipped_gmv$region, estimate = gm_long2$Estimate,
unadjusted_estimate = gm_long2$unadjusted_estimate)
head(gm_figure)
rawdata <- rawdata %>%
pivot_longer(cols = !c(subject, cohort, sex, study),
names_to = c('hemi','region','measurement'),
names_sep = "_",
values_to = 'value'
)
gm_figure <- gm_figure %>%
pivot_longer(cols = c(estimate, unadjusted_estimate),
names_to = 'type',
values_to = 'estimate')
gm_figure <- gm_figure %>%
pivot_wider(names_from = type,
values_from = estimate)
gm_figure$difference <- gm_figure$unadjusted_estimate - gm_figure$estimate
gm_figure <- gm_figure %>%
pivot_longer(cols = c(estimate, unadjusted_estimate),
names_to = 'type',
values_to = 'estimate')
gmv_difference <- ggplot(gm_figure, aes(x = difference, y = region)) +
geom_point(size = 5) +
theme_bw() +
facet_wrap(~hemi) +
geom_vline(xintercept = 0, linetype = 'dotdash') +
xlim(-0.1, 0.1) +
ggplot2::theme_bw(base_size = 28) +
theme(
strip.text = element_text(face = "bold", size = 28),
axis.text.y = element_text(size = 28),
legend.position = "bottom"
)
gmv_difference
ggsave("gmv_difference.png", plot = gmv_difference, dpi = 500, width = 10, height = 12, units = 'in', bg = 'white')
gm_figure$hemi[gm_figure$hemi == 'Left'] <- 'left'
gm_figure$hemi[gm_figure$hemi == 'Right'] <- 'right'
ggplot(gm_figure, aes(fill = difference)) +
geom_brain(atlas = dk, position = position_brain(hemi~side)) +
scale_fill_distiller(palette = 'PRGn',
limits = c(-.1,.1)) +
theme(legend.text = element_text(size = 12), plot.title = element_text(size = 20)) +
theme_void()

Meta_analysis.R at commit d81138d, no license · at the source

Overview

Authors: Ryan W J Loke^1,2, Oscar Ortiz^2,3, Sylvia M Gustin^4,5, Michèle Hubli^6,7, Clas Linnman⁸, Abigail Livny^9,10, Yann Quidé^4,5, Paulina S Scheuren^1,2, John L K Kramer^1,2

ORCID iDs: Abigail Livny, Yann Quidé, Paulina S Scheuren

Institutions: University of British Columbia (Canada); International Collaboration On Repair Discoveries (Canada); UNSW Sydney (Australia); Neuroscience Research Australia (Australia); University of Zurich (Switzerland); ETH Zurich (Switzerland); Universitätsklinik Balgrist (Switzerland); Harvard University (United States); Massachusetts General Hospital (United States); Mass General Brigham (United States); Tel Aviv University (Israel); Sheba Medical Center (Israel)

Journal: Brain communications, volume 8, issue 3, article fcag146

Dates: received 22 October 2025; accepted 21 April 2026; published online 24 April 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1093/braincomms/fcag146 · PMID 42099305 · PMCID PMC13148768 · OpenAlex W4414077135

Open access: gold, a free copy (OpenAlex)

Preprint: doi:10.1101/2025.09.04.25335117

Status: code verified

Categories: structural MRI / diffusion (modality), human (organism), pain (population)

Methods: Statistics, Preprocessing, fMRI & imaging

Keywords: chronic pain, brain morphology, individual participant data meta-analysis, curvature

Topic: Musculoskeletal pain and rehabilitation (Pharmacology, Medicine), according to OpenAlex

Funding: Natural Sciences and Engineering Research Council of Canada

Citations: cited by 3 papers (Europe PMC); 61 references in the paper

Abstract

Chronic pain is a leading contributor to all-cause morbidity and disability, encompassing numerous biopsychosocial dimensions that persistently engage complex networks of brain regions. Meta-analyses have advanced our understanding of structural brain differences in chronic pain but rely exclusively on summary statistics which may introduce heterogeneity related to completeness of reporting and differences in methodological approaches. To address these limitations, we conducted the first individual participant data (IPD) meta-analysis of brain structure alterations in chronic pain. Using traditional morphometric measures (i.e. volume, cortical thickness, and surface area) and differential-geometric shape metrics (i.e. intrinsic and extrinsic curvature), we aimed to reveal alterations in brain structure convergent across chronic pain conditions. We hypothesized that chronic pain would be associated with region-specific grey matter reductions in regions previously implicated in chronic pain (e.g. parahippocampal gyrus and insula) and explored whether curvature metrics would reveal additional structural changes. Anatomical MRI images from eight publicly available datasets spanning five conditions and 401 individuals with chronic pain (and 245 age- and sex- matched healthy controls) were analysed: (i) knee osteoarthritis, (ii) chronic low back pain, (iii) fibromyalgia, (iv) migraine, and (v) primary trigeminal neuralgia. FreeSurfer was used to parcellate T1-weighted anatomical images, and metrics for cortical and subcortical regions were extracted. Meta-analysis revealed a range of structural changes in the brain associated with chronic pain. Cortical thinning and volume loss were small and localized to the temporo-occipital regions, including bilateral volumetric reductions in the entorhinal cortex in individuals with chronic pain. Increases in intrinsic curvature were widespread, involving 49 out of 68 cortical regions. No significant alterations were detected in subcortical volumes. Intrinsic curvature and subcortical volumetric estimates had higher levels of inter-study heterogeneity compared to other metrics, reflecting potential condition and sample-specific variability. Leveraging harmonized processing across a large sample size, our novel IPD meta-analysis highlights both widespread and region-specific structural remodelling of chronic pain-related neuroanatomy.

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

Repository

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

lokeryan/ChronicPainIPD

License: none: the authors keep all their rights

State: the link answers, verified on 29 September 2026

Evidence: files inventoried

Commit: d81138d3bc2bd0455c701187ca4ff41e7ad54b97, 22 January 2026

Languages: R (8), Shell (2)

Size: 12 files, 10 scripts

Software Heritage: not archived

Found in: “Data availability”

Holds: README

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

Tools: tidyverse (8 files), data.table (7 files), ggpubr (7 files), multcomp (7 files), Plotly (7 files), reshape2 (7 files), car (6 files), caret (6 files), glmnet (6 files), pROC (6 files), ggplot2 (2 files), broom (1 file), FreeSurfer (1 file), ggseg (1 file), metafor (1 file), psych (1 file)

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

29 September 2026: the link answers

11 files

CP_OA.R, R, 1,623 lines, 1 match
FM_2.R, R, 405 lines
FM_Exploring.R, R, 887 lines
Figure_making.R, R, 3,715 lines, 1 match
Meta_analysis.R, R, 4,084 lines, 2 matches
PTN.R, R, 911 lines
clbp.R, R, 2,213 lines
extract_ICV.sh, Shell, 32 lines
migraine.R, R, 828 lines
recon_all.sh, Shell, 150 lines, 1 match
README.md, Text, 4 lines

At the source: github.com/lokeryan/ChronicPainIPD

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:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
10 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

Beyond current analyses, we developed an RShiny application (https://lokeryan.shinyapps.io/chronic_pain/) that enables interactive visualization of the raw data, study-level comparisons, and meta-analyses. This tool can allow other researchers to explore analyses in detail and facilitate new lines of inquiry beyond the scope of the present study. All MRI images analysed for this study were obtained from OpenNeuro (https://openneuro.org/), an open-science neuroinformatics database, or OpenPain (https://openpain.org/), an open access data sharing platform for brain imaging studies of human pain. All study-level and sex analysis results can be found in the Supplementary material and can be explored in the RShiny application. In-house bash scripts used to process the MRI images have been uploaded to our GitHub repository, along with the generated R scripts used to analyse all datasets. Additionally, code generated for the purpose of data processing and statistical analysis is included in the following GitHub: https://github.com/lokeryan/ChronicPainIPD. Upon reasonable request, the corresponding author can provide additional information and data to interested researchers for the purpose of academic research and further scientific investigations. If interested, other investigators may contribute their data to our study by reaching out to the corresponding author, which we can implement into our RShiny application. All study-level and meta-analysis results can be found in the supplementary. Please refer to the supplementary table dictionary for descriptions of all sheet names and column variables.

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

Recorded: type, language, journal, volume, issue, pages, dates, 9 authors, 4 keywords, 1 funder, 56 references.

Cite

This paper

Loke, R. W. J., Ortiz, O., Gustin, S. M., Hubli, M., Linnman, C., Livny, A., Quidé, Y., Scheuren, P. S., & Kramer, J. L. K. (2026). Convergent structural brain alterations in chronic pain: a multi-metric individual participant data meta-analysis. Brain communications, 8(3), fcag146. https://doi.org/10.1093/braincomms/fcag146

BibTeX

@article{loke2026convergent, author = {Loke, Ryan W J and Ortiz, Oscar and Gustin, Sylvia M and Hubli, Michèle and Linnman, Clas and Livny, Abigail and Quidé, Yann and Scheuren, Paulina S and Kramer, John L K}, title = {{Convergent structural brain alterations in chronic pain: a multi-metric individual participant data meta-analysis}}, journal = {Brain communications}, year = {2026}, month = apr, volume = {8}, number = {3}, pages = {fcag146}, publisher = {Oxford University Press}, issn = {2632-1297}, doi = {10.1093/braincomms/fcag146}, url = {https://doi.org/10.1093/braincomms/fcag146}, pmid = {42099305}, pmcid = {PMC13148768} }

RIS

TY - JOUR AU - Loke, Ryan W J AU - Ortiz, Oscar AU - Gustin, Sylvia M AU - Hubli, Michèle AU - Linnman, Clas AU - Livny, Abigail AU - Quidé, Yann AU - Scheuren, Paulina S AU - Kramer, John L K TI - Convergent structural brain alterations in chronic pain: a multi-metric individual participant data meta-analysis T2 - Brain communications J2 - Brain Commun PY - 2026 DA - 2026/04/24 VL - 8 IS - 3 SP - fcag146 SN - 2632-1297 PB - Oxford University Press DO - 10.1093/braincomms/fcag146 UR - https://doi.org/10.1093/braincomms/fcag146 LA - en ER -

CSL-JSON

{ "id": "10.1093/braincomms/fcag146", "type": "article-journal", "title": "Convergent structural brain alterations in chronic pain: a multi-metric individual participant data meta-analysis", "container-title": "Brain communications", "author": [ { "family": "Loke", "given": "Ryan W J" }, { "family": "Ortiz", "given": "Oscar" }, { "family": "Gustin", "given": "Sylvia M" }, { "family": "Hubli", "given": "Michèle" }, { "family": "Linnman", "given": "Clas" }, { "family": "Livny", "given": "Abigail" }, { "family": "Quidé", "given": "Yann" }, { "family": "Scheuren", "given": "Paulina S" }, { "family": "Kramer", "given": "John L K" } ], "container-title-short": "Brain Commun", "volume": "8", "issue": "3", "page": "fcag146", "DOI": "10.1093/braincomms/fcag146", "PMID": "42099305", "PMCID": "PMC13148768", "ISSN": "2632-1297", "publisher": "Oxford University Press", "URL": "https://doi.org/10.1093/braincomms/fcag146", "language": "en", "issued": { "date-parts": [ [ 2026, 4, 24 ] ] } }

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.1126/sciadv.aec9291 [code]: Computational mechanisms of perception in autism revealed using games inspired by rodent operant tasks.
Journal: Science advances
In common: pROC, glmnet, caret, 8 other tools
[2] doi:10.1093/braincomms/fcag121 [code]: Anterior insular co-activation patterns associated with stress markers in chronic primary pain.
Journal: Brain communications
In common: multcomp, psych, car, 5 other tools, pain, 2 references
[3] doi:10.1093/bioinformatics/btag592 [code]: Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: multcomp, pROC, car, 7 other tools
[4] doi:10.1016/j.isci.2026.115657 [code]: Integration of machine learning to develop a disulfidptosis model for predicting glioma prognosis, immunotherapy response, and drug.
Journal: iScience
In common: pROC, glmnet, psych, 7 other tools
[5] doi:10.1093/braincomms/fcag236 [code]: Dynamic, state-dependent characteristics of cognitive fluctuations in Lewy body dementia: a magnetoencephalography study.
Journal: Brain communications
In common: ggseg, pROC, caret, 6 other tools
[6] doi:10.1080/20002297.2026.2705667 [code]: Oral microbiota dysbiosis related to the cortical thinning and cognitive impairment in cerebral small vessel disease.
Journal: Journal of oral microbiology
In common: ggseg, pROC, glmnet, 5 other tools, structural MRI / diffusion
[7] doi:10.3390/ijms27104466 [code]: Uncovering the Key Circuit FOSL2/FOS/EGR3/EGR1, Contributing to the Hyperexcitability of Excitatory Neurons in the Epileptic Temporal Cortex and Hippocampus.
Journal: International journal of molecular sciences
In common: pROC, glmnet, caret, 6 other tools
[8] doi:10.1038/s42003-026-10957-8 [code]: Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: glmnet, caret, broom, 6 other tools
[9] doi:10.1038/s41398-026-04131-1 [code]: Multimodal phenotypic classification of generalized anxiety and panic using structural MRI data and psychosocial factors: machine learning results from the German National Cohort (NAKO) study.
Journal: Translational psychiatry
In common: pROC, caret, psych, 5 other tools, structural MRI / diffusion
[10] doi:10.1038/s41467-026-73262-2 [code]: Robust but independent sex differences in human brain function, structure, and behavior.
Journal: Nature communications
In common: caret, car, broom, 5 other tools, structural MRI / diffusion, 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.