Plotting peak voxels from DBM analysis - CoBrALab/documentation GitHub Wiki
You will want to identify peak-voxels and plot them in either a boxplot or scatter-plot with lines, depending on whether your data are cross-sectional or longitudinal. This is especially helpful for understanding the direction of interaction terms or time terms (is volume decreasing, or is it just reduced in growth rate compared to a control group?)
You can identify peak voxels one of two ways. The first is to open the stats files on the average, as described in detail here. You can hover a specific voxel with your cursor and see the coordinates in the bottom left of the Display window. The world voxels you will need are Xw, Yw, and Zw. Then, you would use an atlas to identify where the voxel is located. This continues to be the gold standard to ensure that you know what you are looking at and you should always put eyes on a paper atlas.
The second is to use a label-file appropriate to your age group to automatically identify and label peak voxels, described here. This method is good to confirm what you identified from the first method, but should always be supplementary.
Once you have identified your peak voxels, you want to plot them properly, following the recommendations here.
To properly plot the voxel, you will have to rerun an LMER at precisely that voxel by identifying it with world coordinates so that you can extract the parameters of the model. For example, I used the following function to supply world coordinates and a region name and return the plot:
plot_peak_voxel <- function(Xw, Yw, Zw, region){
data$voxel <-mincGetWorldVoxel(filenames = data$file , v1 = Xw, v2 = Yw, v3 = Zw)
model_voxel = lmer(voxel ~ treatment * poly(I(age-24),2)*sex + (1|ID) + (1|mom_id), data = data)
fitdata_voxel = as.data.frame(Effect(c("treatment", "age","sex"),model_voxel, xlevels=list(age=seq(20,90,5))))
p <- ggplot(data = data, aes(y=voxel,x=age, color=treatment, group=treatment)) +
geom_point() +
geom_line(data = fitdata_voxel, aes(y=fit),size = 2) +
geom_ribbon(data = fitdata_voxel, aes(y=fit, ymin=lower, ymax=upper), alpha=0.2) +
xlab("Age") +
ylab(bquote('Absolute Jacobians')) +
theme_classic()+
labs(fill = "Treatment", group="Treatment", color="Treatment", title = paste0("Peak Voxel in ", region))+
scale_color_manual('Treatment', values = c('SAL' = '#420085', 'THC' = '#018700'))+
scale_size_manual(values=c(2,2))+
theme(plot.title = element_text(hjust = 0.5,size = 16, face = "bold"))+
theme(axis.text=element_text(size=18), axis.title=element_text(size=20,face="bold"),
legend.text=element_text(size=16), legend.title=element_text(size=18))+
scale_x_continuous(breaks=c(20,40,60,90))+
facet_wrap(~sex)
}In this case, data is the data frame, where file is the Jacobians. Having it in a function, you can then call it later in a loop if you used an automated approach.
Remember at this point the result of that model does not matter, you have already identified the peak voxels.
# plotting peak voxel in Perirhinal area
data$voxel <- mincGetWorldVoxel(filenames = data$relative_jacobian , v1 = 4.3, v2 = 0.3, v3 = -0.7)
baseline_mean_rPERI <- data %>%
filter(timepoint == 1) %>%
dplyr::group_by(injection) %>%
dplyr::summarise(
baseline_mean_rPERI = mean(voxel, na.rm = TRUE),
.groups = "drop"
)
# Merge baseline means back into the original dataset
data <- data %>%
left_join(baseline_mean_rPERI, by = "injection") %>%
mutate(voxel_centered_rPERI = voxel - baseline_mean_rPERI)
data$dpi <- as.numeric(data$dpi)
peak_voxel_model_mlm <- lmer(voxel_centered_rPERI ~ dpi * injection + sex + (1|ear_ID), data = data)
fitdata_voxel_mlm = as.data.frame(Effect(c("dpi", "injection"),peak_voxel_model_mlm, xlevels=list(dpi=seq(-12,122,1))))
ggplot(data=data, aes(y=voxel_centered_rPERI, x=dpi, group=injection, color=injection, fill=injection)) +
geom_point(aes(shape=sex),size=2) +
geom_line(data=fitdata_voxel_mlm, aes(y=fit)) +
geom_ribbon(data=fitdata_voxel_mlm, aes(y=fit, ymin=lower, ymax=upper, color=NULL), alpha=0.4) + #NULL color = no outline color for ribbon
scale_color_manual(values = c("#619CFF","#B2182B"))+
scale_fill_manual(values = c("#619CFF","#B2182B"))+
#scale_x_continuous(breaks=c(-7,30,90,120,130)) +
scale_y_continuous(
name = "Volume Change (%)",
labels = function(x) {
#exp() takes the exponent base e of the logged Jacobian to give the Jacobian which is 0-1 for decreases and 1 to infinity for increases
# a log of 1.02 indicates a 2% change, so we subtract 1 for .02 and multiply 100
# a log of .98 is a -2% change so we subtract 1 for -.02 and multiply by 100.
percent_change <- (exp(x) - 1) * 100
# Add a label formatted to 1 decimal place.
sprintf("%.1f%%", percent_change) # Format as "10.5%"
}) +
ylab("% Volume Change from Baseline") + xlab ("Days Post-Injection") +
ggtitle("Peak Voxel in Right Perirhinal Area") +
#theme_dark(base_size = 14)
theme(text = element_text(size=12), plot.title = element_text(hjust = 0.5))
ggsave(filename = "plot-percent-injection_dpi_rightPERI.png", plot=last_plot(), width = 5.5, height = 4.5, units = "in", dpi = 1000)