demo/fmri_demo_func.R

In TCIU: Spacekime Analytics, Time Complexity and Inferential Uncertainty

# load 3d nifti data of the mask saved in the package
dim(mask)

readline("press any key to continue")
# simulate fMRI data within dimension of 64*64*40 and our saved mask
# specify the starting time points of each simulated period and the duration of it
fmri_generate = fmri_simulate_func(dim_data = c(64, 64, 40), mask = mask, 
								   ons = c(1, 21, 41, 61, 81, 101, 121, 141), 
								   dur = c(10, 10, 10, 10, 10, 10, 10, 10))
								 
readline("press any key to continue")
# for simplicity, here we use t-test to generate p-values. We also have other tests available. See more details in fmri_stimulus_detect.
# apply our p value detection function on our simulated fMRI data, choosing t-test
p_simulate_t_test_raw = fmri_stimulus_detect(fmridata= fmri_generate$fmri_data, 
											 mask = fmri_generate$mask,
											 stimulus_idx = fmri_generate$on_time,
											 method = "t-test" , 
											 ons = fmri_generate$ons, 
											 dur = fmri_generate$dur)
			
readline("press any key to continue")
# make the 2d plot of our simulated data based on the p value generated above from sagittal, coronal and axial view
for(axis in c("x", "y", "z")){
  axis_i = switch(axis, 
                  "x" = {35},
                  "y" = {30},
                  "z" = {22})
  print(fmri_2dvisual(p_simulate_t_test_raw, list(axis, axis_i), 
                      hemody_data=NULL, mask=fmri_generate$mask, 
                      p_threshold = 0.05, legend_show = TRUE, 
                      method = "scale_p",
                      color_pal = "YlOrRd", multi_pranges=TRUE))
}
			
readline("press any key to continue")
# make the 3d plot of our simulated data based on the p value generated above, using scale_p method
fmri_3d_plot = fmri_3dvisual(p_simulate_t_test_raw, fmri_generate$mask, 
							 p_threshold = 0.05, method="scale_p", multi_pranges=TRUE)

print(fmri_3d_plot$plot)

readline("press any key to continue")
# to verify the mask of the simulated data is the same as the saved p value's corresponding mask
if(all(fmri_generate$mask==mask[,,]) == TRUE){
# make comparison of the 3d plot of two set of p values. One is generated by our simulated data, 
# the other is the one saved in our package
print( fmri_pval_comparison_3d(list(p_simulate_t_test_raw, phase3_pval), mask, 
				          	   list(0.05, 0.05), list("scale_p", "scale_p"), 
				          	   multi_pranges=FALSE) )

fmri_pval_comparison_2d(list(p_simulate_t_test_raw, phase3_pval), 
						list('pval_simulated', 'pval_posthoc'),
						list(list(35, 33, 22), list(40, 26, 33)), 
						hemody_data = NULL, 
						mask = mask, p_threshold = 0.05, 
						legend_show = FALSE, method = 'scale_p',
						color_pal = "YlOrRd", multi_pranges=FALSE)
}

readline("press any key to continue")
# Here we apply real-valued fmri data to fmri_ROI_phase1(),
# and calculates p_values for the ROIs individually to test whether it is potentially activated.
ROI_phase1 = fmri_ROI_phase1(fmridata = fmri_generate$fmri_data,
                  			 label_mask = mask_label,
                  			 label_dict = mask_dict,
                  			 stimulus_idx = fmri_generate$on_time) 

readline("press any key to continue")
# Here we apply simulated real-valued fMRI data to fmri_ROI_phase2(),
# and calculate p-values for selected ROIs from phase1 by tensor-on-tensor regression.

# The detection may take a long time, and the user can directly use the sample p-values in the package to generate 3D visualization.

# ROI_phase2 = fmri_ROI_phase2(fmridata = fmri_generate$fmri_data, label_mask = label_mask, 
#                                            ROI_label_dict = ROI_detect_p_val[[2]], stimulus_idx = fmri_generate$ons,
#                                            stimulus_dur = fmri_generate$dur, rrr_rank = 3,
#                                            fmri.design_order = 2, fmri.stimulus_TR = 3, 
#                                            method = "t_test", parallel_computing = FALSE, max(detectCores()-2,1)

readline("press any key to continue")
# do the FDR correction and the spatial clustering
ROI_phase3 = fmri_post_hoc(phase2_pval , fdr_corr = "fdr",
						   spatial_cluster.thr = 0.05,
						   spatial_cluster.size = 5, 
						   show_comparison = FALSE)


readline("press any key to continue")
# Here we visualize the 3D plot of p-values for three different phases
label_index = mask_dict$index
label_name = as.character(mask_dict$name)
label_mask = mask_label
fmri_3dvisual_region(phase1_pval, label_mask, label_index, label_name, title = "phase1 p-values")

readline("press any key to continue")
# for 3D visualization, user needs to include empty region in the label
label_index = c(0, label_index)
label_name = c("empty", label_name)
fmri_3dvisual_region(phase2_pval, label_mask, label_index, label_name, title = "phase2 p-values")

readline("press any key to continue")
fmri_3dvisual_region(list(phase2_pval,phase3_pval), label_mask, label_index, label_name, title = "phase2&3 p-values")