2 Analysis Setup
2.1 Defining Variables and Folders
First, we set up some global variables that will be used throughout the analyses. This includes the identifiers of the subjects to include in the analysis, specifics of the design such as the number of days and events per day and the regions of interest for representational similarity analysis.
#-------- DEFINE GLOBAL VARIABLES ---------
# list of subject IDs(excluded: 58 and 68 because of bad memory performance and MRI acquisition problems)
subjects <- c("031", "032", "033", "034", "035", "036", "037", "038", "039", "054",
"055", "056", "057", "059", "061", "062", "063", "064", "065", "066",
"069", "070", "071", "072", "073", "074", "075", "076")
n_subs <- length(subjects)
# PVT scanning runs and blocks within each run
runs <- c("pre", "post")
n_runs <- length(runs)
n_blocks <- 10
# design parameters
n_days <- 4
n_events_day <- 5
n_pics <- n_days*n_events_day
# Main regions of interest for analyses
rois <- c("aHPC_lr", "alEC_lr")
# Regions of interest in MNI space
rois_ec = c("alEC_lr")
rois_hpc = c("aHPC_lr")
rois_mni <- c(rois_hpc, rois_ec)
roi_colors <- c("#dd4731", "#079cd6")
# Regions of interest to get from Freesurfer
rois_fs = c("hpc_lr", "ec_lr")
# numeric Freesurfer labels
#(https://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/AnatomicalROI/FreeSurferColorLUT)
labels_fs = list(c(17,53), c(1006,2006))
roi_colors_fs <- wes_palette(n = length(rois_fs), name = "FantasticFox1")
# Define colors to use in plots
event_colors <- scico::scico(n=5, begin = 0, end = 0.6, palette = "bamako")
time_colors <- c("#26588E", "#33602D", "#C1AA6E")
aHPC_colors <- c("#dd4731", # main
"#26588E", # within main, from scico::scico(n=5, begin = 0.1, end = 0.7, palette = "devon")[3],
"#e31a1c", # within low
"#800026", # within high
"#A54440", # across main, from scico::scico(n=5, begin = 0.3, end = 0.9, palette = "lajolla")[3]
"#feb24c", # across low
"#fc4e2a") # across high
names(aHPC_colors) <- c("main", "within_main", "within_low", "within_high", "across_main", "across_low", "across_high")
alEC_colors <- c("#855C85",
"#225ea8",
"#1d91c0",
"#0c2c84",
"#7fcdbb",
"#c7e9b4",
"#41b6c4")
names(alEC_colors) <- c("main", "within_main", "within_low", "within_high", "across_main", "across_low", "across_high")
day_time_int_color <- "#F5DF4D"
time_within_across_color <- unname(alEC_colors["main"])
ultimate_gray <- "#939597"Folder Structure
In a second step, we create some folders that will be used during the analysis. These are folders that contain raw data as well as folders into which processed data is written. CAVE: Some Markdown files still create their own folders.
#-------- SET UP FOLDERS ---------
dirs <- c()
# directory with logs from picture viewing task
dirs$pvt_log_dir <- here("data", "behavior", "logs", "picture_viewing_task")
# directory with logs from day learning task
dirs$dlt_log_dir <- here("data", "behavior", "logs", "day_learning")
# directory with logs from the timeline task
dirs$timeline_dat_dir <- here("data", "behavior", "timeline")
# directories for freesurfer ROIs
dirs$rois_fs_dirs = here("data", "mri", "rois", rois_fs)
# directories for final ROIs in analysis space
dirs$rois_ss_dirs <- here("data", "mri", "rois", rois, "samespace")
# directories for the MNI ROIs in analysis space
dirs$rois_mni_ss_dirs <- here("data", "mri", "rois", rois_mni, "samespace")
# directory where preprocessed data lies for each run and block
dirs$feat_dir <- here("data", "mri", "processed", "functionalDataPerBlock")
# directory with MRI data in the analysis space (samespace)
dirs$samespace_dir <- here("data", "mri", "processed", "samespace")
# base directory for RSA
dirs$rsa_dir <- here("data", "mri", "rsa")
# directories for RSA correlation matrices
dirs$rsa_roi_corr_mat_dirs <- here("data", "mri", "rsa", "correlation_matrices", rois)
# directories for RSA pattern similarity change
dirs$rsa_roi_ps_change_dirs <- here("data", "mri", "rsa", "pattern_similarity_change", rois)
# directories for the cleaned timeseries data
dirs$rsa_roi_clean_timeseries_dirs <- here("data", "mri", "rsa", "clean_roi_timeseries", rois)
# directories for the relevant volumes of each ROI
dirs$rsa_roi_rel_vol_dirs <- here("data", "mri", "rsa", "relevant_roi_volumes", rois)
# directory for data on which to run RSA
dirs$rsa_dat_dir <- here("data", "mri", "rsa", "data_for_rsa")
dirs$mask_dir <- here("data", "mri", "processed", "group_masks")
# base directory for RSA searchlights
dirs$searchlight <- here("data", "mri", "rsa", "searchlight")
# data for final analysis (this will be shared)
dirs$data4analysis <- here("data_for_analysis")
# directory for source data for nature communications
dirs$source_dat_dir <-here("source_data_natcomms")
if(dir.exists(dirs$source_dat_dir)){unlink(dirs$source_dat_dir, recursive = TRUE)} # delete if already existing to make sure to start with empty directory
# directory to save figures to
fig_dir = here("figures")
if(!dir.exists(fig_dir)){dir.create(fig_dir)}
# create all the directories
dirs_created <- lapply(unlist(dirs), function(x) if(!dir.exists(x)) dir.create(x, recursive=TRUE))2.2 Define Analysis Functions
Analysis Logic
We run both the behavioral analysis and the RSA using two approaches: A summary statistics approach and with linear mixed effect models.
In the analysis of the timeline task, we test whether virtual time explains remembered times when competing for variance with order and real time in seconds.
In order to assess the change in pattern similarity that is due to the learning task, we will later calculate the difference of the Fisher transformed mean correlation coefficients for every pair of scenes between the pre and the post picture viewing tasks (see this section). We will analyze how these difference values relate to various predictor variables derived from the learning task, such as the temporal distance between pairs of scenes within a day.
The two analysis approaches are briefly outlined below.
Summary Statistics
The summary statistics approach is based on permutation testing.
# control number of permutations used throughout the analyses
n_perm <- 10000We will use 10000 random permutations throughout the analyses.
In the summary statistics approach, we use the different time metrics as predictors for the remembered times in the timeline task. We will thus run one GLM per participant. In RSA, we set up a GLM with the given variable from the learning task as a predictor and the pairwise RSA difference values as criterion for every participant.
The resulting model coefficients are then compared to a null distribution obtained from shuffling the dependent variable of the linear model (i.e. pattern similarity change) for a large number of times. This results in a p-value for each coefficient, which is transformed to a Z-score. The Z-scores are then taken to the second level for group-level statistics.
This is also described in the methods section of the manuscript:
For the summary statistics approach, we ran a multiple regression analysis for each participant with virtual time, sequence position (order), and real time since the first event of a day as predictors of responses in the timeline task. To test whether virtual time indeed explained participants’ responses even when competing for variance with order and real time, included in the model as control predictors of no interest, we compared the participant-specific t-values of the resulting regression coefficients against null distributions obtained from shuffling the remembered times against the predictors 10000 times. We converted the resulting p-values to Z-values and tested these against zero using a permutation-based t-test (10000 random sign-flips, Figure 2E).
In the summary statistics approach, we used the different time metrics as predictors for pattern similarity change. We set up a GLM with the given variable from the day learning task as a predictor and the pairwise representational change values as the criterion for every participant. The t-values of the resulting model coefficients were then compared to a null distribution obtained from shuffling the dependent variable of the linear model (i.e. pattern similarity change) 10000 times. This approach to permutation-testing of regression coefficients controls Type I errors even under situations of collinear regressors (Anderson and Legendre, 1999). Resulting p-values for each coefficient were transformed to a Z-score. The Z-scores were then used for group-level inferential statistics.
We start by defining the function that permutes the dependent variable of the linear model. This approach was described e.g. by Manly (1997) and is referred to as permutation of raw data in Anderson & Legendre (1999), who compare different ways to implement permutation tests for (partial) regression coefficients. Their simulations show that the chosen approach does well in terms of controlling type I errors and power, even under situations of collinear regressors.
# define function that calculates z-value for permutation
lm_perm_jb <- function(in_dat = df, lm_formula = lm_formula, nsim = 1000){
# run the model for original data and store observed t-values
lm_fit <- lm(formula = lm_formula, data=in_dat)
obs <- coef(summary(lm_fit))[,"t value"]
# extract the dependent variable from the formula
dv <- str_extract(lm_formula, "[^~]+")
dv <- str_replace_all(dv, fixed(" "), "")
if(!(dv %in% colnames(in_dat))){stop("Cannot find dependent variable in input data");}
# initialize df for permutation
data_perm <- in_dat
# set aside space for results
res <- matrix(nrow = nsim, ncol = length(obs))
for (i in 1:nsim) {
# scramble response value
perm <- sample(nrow(in_dat))
dv_dat <- in_dat[dv]
data_perm[dv] <- dv_dat[perm,]
# compute linear model and store the t-value of predictor
lm_fit_perm <- lm(formula = lm_formula, data=data_perm)
res[i,] <- coef(summary(lm_fit_perm))[,"t value"]
}
# append the observed value to the list of results
res <- rbind(res,obs)
# calculate p-value for each coefficient and transform to z
p <- rep(0,length(obs))
z <- rep(0,length(obs))
for (i_coef in 1:length(obs)){
p[i_coef] <- sum(res[,i_coef] >= obs[i_coef])/nrow(res)
z[i_coef] <- qnorm(1-p[i_coef])
}
return(z)
}The Z-values resulting from this first-level permutation are then analyzed on the group level. For t-tests, we use random sign-flips (c.f. the one-sample t-test in FSL Randomise) to non-parametrically test against 0 or assess within-participant differences between conditions. For this, we use the function defined below that is a reduced version of this Matlab implementation.
Group-level statistics were carried out using permutation-based procedures. For t-tests, we compared the observed t-values against a surrogate distribution obtained from 10000 random sign-flips to non-parametrically test against 0 or to assess within-participant differences between conditions. Permutation-based repeated measures ANOVAs were carried out using the permuco package (Frossard and Renaud, 2019).
# function for permutation-based one-sample t-test against 0
# requires tidyverse and broom
# diff should be a vector of values (e.g. differences between paired samples)
# returns the data frame from broom::tidy(t.test(diff)) with an additional column
# for the p-value from n_perm random sign flips
# this function is based on the mult_comp_perm_t1 function for Matlab by David Groppe)
paired_t_perm_jb <- function(in_dat, n_perm=10000){
# make sure input is a vector
in_dat <- unname(unlist(in_dat))
# number of observations
n <- length(in_dat)
# get t of unpermuted data
t_stats <- t.test(in_dat) %>% broom::tidy()
# run n_perm iterations
t_perm <- vector("numeric", n_perm)
for(i in 1:n_perm){
# randomly shuffle sign of each observation and get t-value
dat_perm <- in_dat*sample(c(-1,1),n,replace=TRUE)
t_perm[i] <- abs(t.test(dat_perm)$statistic)
}
#add the negative of all values since we assume null hypothesis distribution is symmetric
t_perm <-c(t_perm, -t_perm)
# add observed t-value so p cannot be 0
t_perm <-c(t_perm, t_stats$statistic)
# calculate two-tailed p-value
p_perm <- mean(t_perm >= abs(t_stats$statistic))*2
t_stats <- t_stats %>%
tibble::add_column(p_perm = p_perm, .after = "p.value")
return(t_stats)
}Linear Mixed Effects
Linear Mixed Effects models consist of fixed and random effects. In our case the fixed case are of interest and consist of the temporal relationships that could explain remembered times or pattern similarity change (the dependent variables in the behavioral analysis and RSA, respectively). Random effects account for the fact that data points from different participants enter the estimation of the fixed effects.
Following the recommendation for maximal random effect structures by Barr et al. (JML, 2013), we first attempt to fit a model with a random effects structure including random intercepts and random slopes for participants. If the model does not converge or results in singular fits, we reduce the random effects structure, attempting to always at least keep random slopes for the fixed effect of interest in the model as these are crucial to avoid anti-conservativity.
Statistical inference about a model is done using a likelihood ratio tests against a nested, reduced model. The reduced model is identical to the full model, only the fixed effect of interest is removed.
LME assumptions
The code below defines a function to generate three diagnostic plots to visually assess the assumptions of a mixed model.
- Linearity & Homoscedasticity: Residual plot
- Normality of residuals: QQ-Plot and histogram of the residuals
It returns a ggplot object based on the input LME model. Probably these diagnostic plots will not end up in the manuscript because space is limited.
lmm_diagplots_jb <- function(lmm = lmm_full){
# residual plot to inspect homoscedasticity
resids_gg <- ggplot() +
geom_point(aes(x = fitted(lmm), y = residuals(lmm)),
size = 1, shape = 16, alpha = 0.5) +
geom_smooth(aes(x = fitted(lmm), y = residuals(lmm)),
formula = "y ~ x", method='glm', se = TRUE, color = "darkred") +
ylab('residuals') +
xlab('fitted values') +
ggtitle("Residual Plot")+
theme_cowplot()
# data frame for QQ plot and histogram
df <- data.frame(r = residuals(lmm))
# QQ-Plot look at residuals of the model
qqplot_gg <- ggplot(df, aes(sample = r)) +
stat_qq(size = 1, shape = 16, alpha = 0.5) +
stat_qq_line(color = "darkred") +
ggtitle("QQ-Plot")+
theme_cowplot()
# Histogram of residuals
hist_gg <-ggplot(df, aes(x=r)) +
geom_histogram(aes(y=..density..), colour="black", fill="darkgrey", bins = 50) +
geom_density(color="darkred") +
xlab("magnitude of residual") +
ggtitle("Histogram of Residuals")+
theme_cowplot()
# collect for final figure
diag_fig <- resids_gg + qqplot_gg + hist_gg &
theme(text = element_text(size=10),
axis.text = element_text(size=8),
legend.position = 'auto',
aspect.ratio = 1,
plot.title = element_text(hjust = 0.5)) &
plot_annotation(tag_levels = 'A',
caption = paste0(deparse(formula(lmm)), collapse=""))
return (diag_fig)
}LME Summary Tables
To summarize final models we create tables inspired by the best practice guidelines by Meteyard & Davies (JML, 2020). Examples can be found on their OSF page, in particular the example reporting table on page 4 of this document.
To create these tables we rely on the broom.mixed package to get the LME model summary in a tidy format. The tidy dataframes are then converted to huxtables, which can be nicely formatted.
The function below takes as an input the tidy data frames for fixed effects, random effects as well as the ANOVA results from the comparison of the full model against a reduced model without the fixed effect of interest. All are merged into one table to limit the number of tables.
make_lme_huxtable <- function(fix_df, ran_df, aov_mdl, fe_terms=NULL, re_terms=NULL, re_groups=NULL, lme_form=NULL, caption="Summary of Linear Mixed Effects Model"){
########### FIXED EFFECTS
# create huxtable for fixed effects and format it
fix_hux <- fix_df %>% huxtable::huxtable(., add_colnames = TRUE) %>%
# set standard error and t-value column names
set_contents(row = 1, col = 4, value = "SE") %>%
set_contents(row = 1, col = 5, value = "t-value") %>%
# merge the confidence column title and reset column name
merge_cells(row = 1, col = 6:7) %>%
set_contents(row = 1, col = 6, value = "95% CI") %>%
set_align(row = 1, col = 6, value = "center") %>%
# align the contents
set_align(row = 1, col = 2:5, "center") %>%
set_align(row = 1, col = 1, "center") %>%
set_align(row = 2:nrow(.), col = 3:7, "center") %>%
# how many digits to print?
set_number_format(row = 2:nrow(.), col = 3:7, 6) %>%
set_number_format(row = 2:nrow(.), col = 5, 2) %>%
# add header row for fixed effects
dplyr::select(-effect) %>% # remove column to the left because not needed
huxtable::insert_row(.,c("fixed effects", rep("",ncol(.)-1)), after=0) %>%
set_align(row=1,col=1, value="center") %>%
set_header_cols(col=1, value=TRUE) %>%
# bottom border
set_bottom_border(row=nrow(.), value=0.5)
# if names for fixed effect terms are supplied use them
if(!is.null(fe_terms)){fix_hux$term[3:nrow(fix_hux)] <- fe_terms}
########### RANDOM EFFECTS
# create huxtable for random effects and format it
ran_hux <- ran_df %>% huxtable::huxtable(., add_colnames = TRUE) %>%
# how many digits to print?
set_number_format(row = 2:nrow(.), col = 4, 6) %>%
# align header row
set_align(row = 1, col = c(1,2), "left") %>%
# add header row for random effects
select(-effect) %>% # remove column to the left because not needed
huxtable::insert_row(.,c("random effects", rep("",ncol(.)-1)), after=0) %>%
set_align(row=1,col=1, value="center") %>%
set_header_cols(col=1, value=TRUE) %>%
# center estimate data
set_align(row=2:nrow(.),col=ncol(.), value="center") %>%
# bottom border
set_bottom_border(row=nrow(.), value=0.5)
# if names for random effect grouping factor are supplied use them
if(!is.null(re_groups)){ran_hux$group[3:nrow(ran_hux)] <- re_groups}
# if names for random effect terms are supplied use them
if(!is.null(re_terms)){ran_hux$term[3:nrow(ran_hux)] <- re_terms}
########### MODEL COMPARISON
# determine how many digits of p-value to print
if(aov_mdl$`Pr(>Chisq)`[2]>=0.001)
{num_fmt<-3} # 3 digits if p>0.001
else{
num_fmt<-NA # use default -> scientific notation
aov_mdl$`Pr(>Chisq)`[2]<-formatC(aov_mdl$`Pr(>Chisq)`[2], format = "e", digits = 2)}
aov_hux <- aov_mdl %>% as.data.frame() %>%
tibble::rownames_to_column(var ="term") %>%
huxtable(add_colnames = TRUE) %>%
# align contents
set_align(row = 1:nrow(.), col = c(2:ncol(.)), "center") %>%
set_align(row = 1, col = 1, "left") %>%
# add header row for model comparison
huxtable::insert_row(.,c("model comparison", rep("",ncol(.)-1)), after=0) %>%
#merge_cells(1:nrow(.), col=1) %>%
set_align(row=1,col=1, value="center") %>%
set_header_cols(col=1, value=TRUE) %>%
# add model names
set_contents(row = 2:nrow(.), col = "term", value = c("model", "reduced model", "full model")) %>%
set_align(row = 2, col = 2, "left") %>%
# change some column names
set_contents(row = 2, col="Pr(>Chisq)", value = "p") %>%
set_contents(row = 2, col="logLik", value = "LL") %>%
set_contents(row = 2, col="Df", value = "df") %>%
set_contents(row = 2, col="Chisq", value = "X2") %>%
#set_contents(row = 2, col="Chisq", value = expression("$\\chi"^"2")) %>%
# how many digits to print? (do at the end because affected by adding columns)
set_number_format(row = 3:4, col = 3:7, 2) %>%
set_number_format(row = 4, col = 9, num_fmt) %>%
# remove the deviance and BIC column
remove_cols(deviance) %>%
remove_cols(BIC) %>%
# bottom border
set_bottom_border(row=nrow(.), value=0.5)
####### MERGE THE THREE HUXTABLES
# to be able to merge the tables, they need to have the same number of columns.
# the AOV table has 7 columns, so we add 1 empty columns to the fixed effects
# table and 4 empty columns to the random effects table. We merge these with
# existing columns immediately
ran_hux <-ran_hux %>% mutate(a=NA, b=NA, c=NA, d=NA, .after = "term") %>%
merge_across(row = 1:nrow(.), col = 2:6)
fix_hux <-fix_hux %>% mutate(a=NA, .after = "term") %>%
merge_across(row = 1:nrow(.), col = 1:2)
# merge
lmm_hux <- huxtable::add_rows(fix_hux, ran_hux) %>%
huxtable::add_rows(., aov_hux) %>%
# set the header and bottom border
set_caption(caption) %>%
set_caption_pos("topleft") %>%
set_bottom_border(row=nrow(.), value=0.5) %>%
add_footnote(sprintf("model: %s; \nSE: standard error, CI: confidence interval, SD: standard deviation, npar: number of parameters, LL: log likelihood, df: degrees of freedom, corr.: correlation",
lme_form), border = 0.5)
return(lmm_hux)
}The resulting huxtables nicely summarize the mixed models in the HTML documentation. To collect the tables in a word file that accompanies the manuscript, we convert them to flextables. These can be written to Word documents using the officer package. The function below does the conversion plus some touching up to end up with nicely formatted tables in Word.
convert_huxtable_to_flextable <-function(ht = lmm_hux){
# define border style to apply to selected cells
def_cell_l=officer::fp_cell(border.right = fp_border(), border.top = fp_border())
def_cell_t=officer::fp_cell(border.bottom = fp_border(), border.top = fp_border())
# define text style to apply to selected cells
def_par=officer::fp_par(text.align = "center", padding=3)
def_tex=officer::fp_text(bold=TRUE)
# find the header rows (where new sections of the table begin)
head_rows <- match(c("term", "group", "model"), ht$term)
# find the cell where we want to add Chi2
x2_cell <- which(ht=="X2", arr.ind = TRUE)
# create the flextable
ft <- ht %>% huxtable::as_flextable() %>%
# add border at the bottom to the rows that have the names
flextable::style(i=head_rows, pr_c = def_cell_t, pr_p = def_par, pr_t = def_tex) %>%
flextable::style(i=head_rows-1, pr_c = def_cell_t, pr_p = def_par, pr_t = def_tex) %>%
flextable::bg(i=head_rows-1, bg="lightgrey", part="all") %>%
# left-align the first and second column (apart from no. of params in model comparison) and add consistent padding
flextable::align(., j=c(1,2), align="left") %>%
flextable::align(., i = c(head_rows[3], head_rows[3]+1, head_rows[3]+2), j=2, align="center") %>%
flextable::padding(., padding = 3, part="all") %>%
# set font style
flextable::fontsize(size=10, part="all") %>%
flextable::font(fontname = font2use) %>%
# replace X2 with greek letter chi and superscript 2
flextable::compose(i=x2_cell[1], j=x2_cell[2], value = as_paragraph("\U03C7", as_sup("2"))) %>%
# autofit to page
flextable::set_table_properties(layout="autofit", width=1)
# set caption style
ft <- flextable::set_caption(ft,caption = ft$caption$value, style="Normal")
return(ft)
}Here we open the word file that we want to write the tables to. It is based on a .docx-file where the themes for headings and text was manually modified to match the style of the manuscript.
stables_docx <- officer::read_docx(here("virtem_code", "modified_headings.docx"))
stables_docx <- stables_docx %>%
officer::body_add_par("Supplemental Tables", style = "heading 1", pos = "on")Brain Plots in ggplot
These functions are based on the ggBrain package. Plus, there are two simple custom functions to transform between 1mm MNI and matrix coordinates.
ggBrain functions to get template background
The first function returns a ggplot object of the template brain (MNI 1mm in our case).
# function from ggBrain that returns a ggplot object of the template brain
getggTemplate<-function(col_template,row_ind,col_ind, ...){
templateFrame<-getBrainFrame(row_ind=row_ind, col_ind=col_ind, ...)
n<-length(col_template)
if(n>1) col_cut<-as.numeric(cut(templateFrame$value,n))
if(n==1) col_cut=1
p<-ggplot()+facet_just_unique(row_ind,col_ind)
for(i in 1:n){
if(all(col_cut!=i)) next
drop_ind<-which(names(templateFrame)=='value')
#so it doesn't conflict with mappings to "value" later on
templateFrame_col<-templateFrame[col_cut==i,-drop_ind]
p<- p + geom_tile(data=templateFrame_col,aes(x=row,y=col),fill=col_template[i])
}
p
}
# another internal function from ggBrain used to plot the template
facet_just_unique<-function(row_ind,col_ind){
if( all(row_ind==row_ind[1]) & all(col_ind==col_ind[1]))
out<-NULL
if( all(row_ind==row_ind[1]) & !all(col_ind==col_ind[1]))
out<-facet_grid(.~col_ind)
if(!all(row_ind==row_ind[1]) & all(col_ind==col_ind[1]))
out<-facet_grid(row_ind~.)
if(!all(row_ind==row_ind[1]) & !all(col_ind==col_ind[1]))
out<-facet_grid(row_ind~col_ind)
return(out)
}Coordinate transforms
The code section below defines two helper functions that transform between the MNI coordinate system and the matrix coordinates of R. CAVE: Only use with 1mm MNI space.
# to get accurate labels of panels in MNI space, define functions to convert coords
mni2vox <- function(mni_coords, output_zero_based = FALSE){
x<- mni_coords[1] * -1 + 90
y<- mni_coords[2] * 1 + 126
z<- mni_coords[3] * 1 + 72
vox_coords <- c(x,y,z)
# add +1 if output coordinates should not be zero-based
if(!output_zero_based){vox_coords <- vox_coords+1}
return(vox_coords)
}
# convert from voxel coordinates to MNI with option for whether voxel coordinates are from
# R (i.e. not zero-based) or FSL (i.e. zero-based)
vox2mni <- function(vox_coords, input_zero_based = FALSE){
# if input coordinates are from R, i.e. not zero-based, subtract -1
if(!input_zero_based){vox_coords <- vox_coords-1}
# subtract voxel space coordinates of origin (MNI=0,0,0 at 90,126,72 in FSL voxel coordinates)
mni_x <- -1*(vox_coords[1]-90)
mni_y <- vox_coords[2]-126
mni_z <- vox_coords[3]-72
mni_coords <- c(mni_x,mni_y,mni_z)
return(mni_coords)
}Waiting for Cluster Jobs
The function below is a convenient function to pause the execution of a script until a batch of condor jobs has finished. The function takes the ID of a condor batch as an input. It uses condor_q via the system-command. If 0 is returned, it terminates and moves on. Else the function sleeps for (a default of) 60s and then tries again. No checks are performed whether the jobs finished without errors. This function merely pauses the execution of the code to wait for long jobs to finish.
pause_until_batch_done<- function(batch_id, wait_interval = 60){
print(sprintf("Monitoring batch job %s. Checking status every %d seconds.",
batch_id, wait_interval))
batch_running <- TRUE
tictoc::tic()
while (batch_running){
# sleep for the specified number of seconds
Sys.sleep(wait_interval)
# check if batch finished (condor_q returns character(0))
status <- system(sprintf("condor_q -l %s", batch_id), intern=TRUE)
if(identical(status, character(0))){
print("Batch jobs finished. Moving on.")
tictoc::toc()
batch_running <- FALSE
} else{
print(sprintf("Batch jobs still running. Waiting another %d seconds.",
wait_interval))
}
}
}