--- title: "TCIU Predictive Analytics: " subtitle: "Validation using 4D Complex (mag,phase) and Real (solely-mag reconstruction) fMRI data" author: "SOCR Team" date: "`r format(Sys.time(),'%m/%d/%Y')`" output: html_document: theme: spacelab highlight: tango includes: before_body: TCIU_header.html toc: true number_sections: true toc_depth: 2 toc_float: collapsed: false smooth_scroll: true code_folding: hide --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE,cache = FALSE, warning = FALSE) ``` # Package Loading and Data Manipulation ```{r,message=FALSE} library(TCIU) library(DT) library(R.matlab) library(AnalyzeFMRI) library(ggfortify) mat1 = readMat("subj1_run1_complex_all.mat") bigim1 = mat1$bigim[,64:1,,] bigim1_mod = Mod(bigim1) smoothmod = GaussSmoothArray(bigim1_mod,sigma = diag(3,3)) # dim(bigim1) = 64 64 40 # bigim1 contains the complex image space # dimensions are 64x*64y*40z*160t, corresponding to x,y,z,time load("mask.rda") # load the 3D nifit data of the mask load("mask_label.rda") # load the 3D nifit data of the mask with labels load("mask_dict.rda") # load the label index and label name label_index = mask_dict$index label_name = as.character(mask_dict$name) label_mask = mask_label load("hemody.rda") # load the hemodynamic contour of the brain ``` # Time-series graphs ## Interactive time-series visualization The function `fmri_time_series` is used to create four interactive time series graphs for the real, imaginary, magnitude, and phase parts for the fMRI spacetime data. We can also add a reference plotly object to the plot. This function is based on `GTSplot` function from package `TSplotly`. ### Example: X: 44, Y: 42, Z: 33 ```{r,fig.height=10,fig.width=10} reference_plot = sample[[4]] fmri_time_series(bigim1, c(44,42,33), is.4d = TRUE, ref = reference_plot) ``` # Kime-series/Kime-surfaces (spacekime analytics protocol) This Rmd notebook uses a time-series simulation to illustrate how to transform the fMRI data at a fixed voxel location into a kime-surface (kime-series). Notes: Validate all steps in this transformation protocol and finalize this 3D visualization. Try it with real fMRI data at brain voxel locations associated with the finger-tapping task. ## Pseudo-code + Randomly generate 8 $phi=\phi$ kime-phases for each of the 10 time radii. This yields an $8\times 10$ array (*phi_8_vec*) of kime phase directions. These directions can be obtained by different strategies, e.g., (1) uniform or Laplace distribution sampling over the interval $[-\pi:\pi)$, (2) randomly sampling with/without replacement from known kime-phases obtained from similar processes, etc. + Optionally, order all kime-phases (rows) from small to large for each column. + Fix the $\nu=(x,y,z)$ voxel location and extract the fMRI time-course $fMRI_{\nu}(t), \forall 1\leq t\leq 160$. + For binary stimuli (e.g., activation (ON) and rest (OFF) event-related design), we can split the 160-long fMRI series into 80 ON (Activation) and 80 OFF (rest) states, or sub-series. + Construct a data-frame with 160 rows and 4 columns; time (1:10), phases (8), states (2), and fMRI_value (Complex or Real intensity). + Convert the long DF representing fMRI_ON and fMRI_OFF from their native (old) polar coordinates to the (new) Cartesian coordinates, using polar transformations. + Check for visual (graphical) and numeric differences between the fMRI intensities during the ON vs. OFF states + Spatially smooth (blur) the matrices (in 2D) to reduce noise make them more representative of the process. May also need to temper the magnitude of the raw fMRI intensities, which can have large range. + Generate the `plot_ly` *text labels* that will be shown on mouse hover (pop-up dialogues) over each kime-surface/kime-series. These text-labels are stored in Cartesian coordinates $(-10\leq x\leq 10,-10\leq y\leq 10)$, but are computed using the polar coordinates $(1\leq t\leq 10,-\pi\leq \phi<\pi)$ and the polar-to-Cartesian transformation. This the labels are $21\times21$ arrays including the triple $(t, \phi, fMRIvalue)$. Separate text-labels are generated for each kime-surface (ONN vs. OFF stimuli). + Generate the $21\times21$ kime-domain Cartesian grid by polar-transforming the polar coordinates $(t,\phi)$ into Cartesian counterparts $(x,y)$. + Interpolate the fMRI intensities (natively anchored at $(t,\phi)$) to $fMRI(x,y), \forall -11\leq x,y \leq 11, x^2+y^2\leq 10$. + Use `plot_ly` to display in 3D the kime-series as 2D manifolds (kime-surfaces) over the Cartesian domain. ## Function main steps illustration: Kime-surfaces / kime-series construction ### Generate the kime-phases Randomly generate 8 $phi=\phi$ kime-phases for each of the 10 time radii. This yields an $8\times 10$ array (*phi_8_vec*) of kime phase directions. These directions can be obtained by different strategies, e.g., (1) uniform or Laplace distribution sampling over the interval $[-\pi:\pi)$, (2) randomly sampling with/without replacement from known kime-phases obtained from similar processes, etc. Optionally, order all kime-phases (rows) from small to large for each column. ```{r} # plot Laplacian ggfortify::ggdistribution(extraDistr::dlaplace, seq(-pi, pi, 0.01), m=0, s=0.5) # randomly generate 8 phi kime-phases for each of the 10 time radii phi_8_vec <- matrix(NA, ncol=10, nrow = 8) for (t in 1:10) { # for a given t, generate 8 new phases set.seed(t); phi_8_vec[ ,t] <- extraDistr::rlaplace(8,mu=0,sigma=0.5) # rank-order the phases for consistency # within the same foliation leaf phi_8_vec[ ,t] <- sort(phi_8_vec[ ,t]) # force phases in [-pi: pi) for (i in 1:8) { if (phi_8_vec[i,t] < -pi) phi_8_vec[i,t] <- -pi if (phi_8_vec[i,t] >= pi) phi_8_vec[i,t] <- pi } } # order all kime-phases (rows) from small to large for each column # phi_8_vec_ordered <- apply(phi_8_vec, 2, sort) ``` ### Split the fMRI series according to the event-related design For binary stimuli (e.g., activation (ON) and rest (OFF) event-related design), we can split the 160-long fMRI series into 80 ON (Activation) and 80 OFF (rest) states, or sub-series. ```{r} fMRI_ON<-bigim1_mod[40,42,33,][c(rep(TRUE,10),rep(FALSE,10))] fMRI_OFF<-bigim1_mod[40,42,33,][c(rep(FALSE,10),rep(TRUE,10))] ``` ### Preprocess the matrices Spatially smooth (blur) the matrices (in 2D) to reduce noise make them more representative of the process. May also need to temper the magnitude of the raw fMRI intensities, which can have large range. ```{r} # Convert the long DF representing fMRI_ON and fMRI_OFF from polar coordinates to Cartesian coordinates matrix_ON <- matrix(0, nrow = 21, ncol = 21) matrix_OFF <- matrix(0, nrow = 21, ncol = 21) for (t in 1:10) { for (p in 1:8) { x = 11+t*cos(phi_8_vec[p,t]) y = 11+t*sin(phi_8_vec[p,t]) matrix_ON[x,y] <- fMRI_ON[(p-1)*10 +t] matrix_OFF[x,y] <- fMRI_OFF[(p-1)*10 +t] } } # smooth/blur the matrices matrix_ON_smooth <- (1/10000)*as.matrix(spatstat::blur(spatstat::as.im(matrix_ON), sigma=0.5)) matrix_OFF_smooth <- (1/10000)*as.matrix(spatstat::blur(spatstat::as.im(matrix_OFF), sigma=0.5)) ``` ### Generate `plotly` labels Generate the `plot_ly` text labels that will be shown on mouse hover (pop-up dialogues) over each kime-surface/kime-series. These text-labels are stored in Cartesian coordinates $(-10\leq x\leq 10,-10\leq y\leq 10)$, but are computed using the polar coordinates $(1\leq t\leq 10,-\pi\leq \phi<\pi)$ and the polar-to-Cartesian transformation. This the labels are $21\times21$ arrays including the triple $(t, \phi, fMRIvalue)$. Separate text-labels are generated for each kime-surface (ONN vs. OFF stimuli). ```{r} # fix the plot_ly Text Lables x <- vector() y <- vector() i <- 1 for (t in 1:10) { for (p in 1:8) { x[i] = 11+t*cos(phi_8_vec[p,t]) y[i] = 11+t*sin(phi_8_vec[p,t]) i <- i+1 } } ``` ### Cartesian representation Generate the $21\times21$ kime-domain Cartesian grid by polar-transforming the polar coordinates $(t,\phi)$ into Cartesian counterparts $(x,y)$. ```{r} hoverText <- cbind(x=1:21, y=1:21, height=as.vector(t(matrix_ON_smooth))) # tail(mytext) custom_txt <- matrix(NA, nrow=21, ncol=21) hoverTextOFF <- cbind(x=1:21, y=1:21, height=as.vector(t(matrix_OFF_smooth))) # tail(mytext) custom_txtOFF <- matrix(NA, nrow=21, ncol=21) for (x in 1:21) { for (y in 1:21) { t = sqrt((x-11)^2 + (y-11)^2) p = atan2(y-11, x-11) custom_txt[x,y] <- paste(' fMRI: ', round(hoverText[(x-1)*21+y, 3], 3), '\n time: ', round(t, 0), '\n phi: ', round(p, 2)) custom_txtOFF[x,y] <- paste(' fMRI: ', round(hoverTextOFF[(x-1)*21+y, 3], 3), '\n time: ', round(t, 0), '\n phi: ', round(p, 2)) } } ``` ### Cartesian space interpolation Interpolate the fMRI intensities (natively anchored at $(t,\phi)$) to $fMRI(x,y), \forall -11\leq x,y \leq 11, x^2+y^2\leq 10$. ```{r} xx2 <- 11 + c(-10:10) %o% cos(seq(-pi, pi, 2*pi/20)) yy2 <- 11 + c(-10:10) %o% sin(seq(-pi, pi, 2*pi/20)) #zz2 <- as.vector(matrix_ON_smooth) %o% rep(1, 21*21) zz2 <- matrix_ON_smooth ww2 <- matrix_OFF_smooth dd2 <- matrix_ON_smooth-matrix_OFF_smooth #plot 2D into 3D and make the text of the diameter (time), height (r), and phase (phi) f <- list(family = "Courier New, monospace", size = 18, color = "black") x <- list(title = "k1", titlefont = f) y <- list(title = "k2", titlefont = f) z <- list(title = "fMRI Kime-series", titlefont = f) ``` ## Function result The following parts are the main output of our functions ### Generate a long data-frame Construct a data-frame with 160 rows and 4 columns; time (1:10), phases (8), states (2), and fMRI_value (Complex or Real intensity). Then, convert the long DF representing fMRI_ON and fMRI_OFF from their native (old) polar coordinates to the (new) Cartesian coordinates, using polar transformations. ```{r} datatable(fmri_kimesurface(bigim1_mod,c(44,42,33))[[1]]) ``` ### Display kime-surfaces Use `plot_ly` to display in 3D the kime-series as 2D manifolds (kime-surfaces) over the Cartesian domain. ```{r} fmri_kimesurface(bigim1_mod,c(44,42,33))[[2]] fmri_kimesurface(bigim1_mod,c(44,42,33))[[3]] fmri_kimesurface(bigim1_mod,c(44,42,33))[[4]] ``` # Interactive plotly Example ## Plotly method: interactive way The function `fmri_image` is used to create images for front view, side view, and top view of the fMRI image. ```{r,warning=FALSE} bigim1_mask<-bigim1_mod for (i in 1:160) { bigim1_mask[,,,i]<-bigim1_mask[,,,i]*mask } ``` Manually examples: The plot without mask restriction. ```{r,warning=FALSE,fig.height=6,fig.width=7.5} fmri_image(bigim1_mod, option = "manually", voxel_location = c(40,40,30), time = 4) ``` The plot with mask restriction. ```{r,warning=FALSE,fig.height=6,fig.width=7.5} fmri_image(bigim1_mask, option = "manually", voxel_location = c(40,40,30), time = 4) ``` ## Forcasting with time series The function `fmri_ts_forecast` will fit the ARIMA model for each location. This function is based on `TSplot_gen` function from package `TSplotly`. ```{r,fig.width=10,fig.height=4.5} fmri_ts_forecast(smoothmod,c(41,44,33)) ``` # Motor area detection ## fMRI data simulation The function `fmri_simulate_func` is used to simulate the real-valued fMRI data with the specified dimension. ```{r} 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)) ``` ## Stimulus detection An integrated function `fmri_stimulus_detect` which can apply multiple methods is used to detect motor area. "t-test" and "wilcox-test" can be applied to all real-valued fMRI data. "HotellingsT2","Wilk's Lambda" and "gLRT" method can be applied to all complex-valued fMRI data. "on_off_diff" and "HRF" method can be applied to 4D real-valued fMRI data where the first method is calculating on-off difference period polar volume to get p-values and the second method is using hemodynamic response function and linear regression. The post-hoc stat-mapping filtering can be also applied to the calculated p-values if specified in the parameter or use function `fmri_post_hoc` to do the process. ### Examples ```{r} # statistical method HRF needs parameter ons and dur pval1 = fmri_stimulus_detect(fmridata= bigim1_mod, mask = mask, stimulus_idx = c(1:160)[rep(c(TRUE,FALSE), c(10,10))], method = "HRF" , ons = c(1, 21, 41, 61, 81, 101, 121, 141), dur = c(10, 10, 10, 10, 10, 10, 10, 10) ) # statistical method t-test for real-valued fMRI data pval2 = fmri_stimulus_detect(fmridata= bigim1_mod, mask = mask, stimulus_idx = c(1:160)[rep(c(TRUE,FALSE), c(10,10))], method = "t-test") # statistical method Wilk's Lambda for complex-valued data pval3 = fmri_stimulus_detect(fmridata = bigim1, mask = mask, stimulus_idx = c(1:160)[rep(c(TRUE,FALSE), c(10,10)) ], method = "Wilks-Lambda" ) # do the fdr correction and the spatial clustering # pval4 is the pval1 after the post-hoc processing pval4 = fmri_post_hoc(pval1, fdr_corr = "fdr", spatial_cluster.thr = 0.05, spatial_cluster.size = 5, show_comparison = FALSE) ``` ```{r} summary(pval1) summary(pval2) summary(pval3) summary(pval4) ``` # Motor area visualization ## Viusalization and comparison of p-value ### 3D viusalization for p-value ```{r fig.width = 9, fig.align = "center", warning=FALSE} pval1_3d = fmri_3dvisual(pval1, mask, p_threshold = 0.05, method="scale_p", multi_pranges=TRUE, title="HRF method") pval1_3d$plot ``` ```{r fig.width = 9, fig.align = "center", warning=FALSE} pval4_3D = fmri_3dvisual(pval4, mask, p_threshold = 0.05, method="scale_p", multi_pranges=TRUE, title="HRF method after post-hoc") pval4_3D$plot ``` ### 2D viusalization for p-value Generate the 2D plot of the p-value from sagittal, coronal and axial view. Plot without hemodynamic contour. ```{r warning=FALSE} for(axis in c("x", "y", "z")){ axis_i = switch(axis, "x" = {35}, "y" = {30}, "z" = {22}) print(fmri_2dvisual(pval1, list(axis, axis_i), hemody_data=NULL, mask=mask, p_threshold = 0.05, legend_show = TRUE, method = "scale_p", color_pal = "YlOrRd", multi_pranges=TRUE)) } ``` Plot with hemodynamic contour. ```{r warning=FALSE} for(axis in c("x", "y", "z")){ axis_i = switch(axis, "x" = {35}, "y" = {30}, "z" = {22}) print(fmri_2dvisual(pval1, list(axis, axis_i), hemody_data=hemody, mask=mask, p_threshold = 0.05, legend_show = TRUE, method = "scale_p", color_pal = "YlOrRd", multi_pranges=TRUE)) } ``` ## Comparison of performance of different methods on the same fMRI data ### 3d p-value comparison `fmri_pval_comparison_3d` is to visualize two p-value data to see how they perform differently on detecting stimulated parts by 3D plot. Here we compare the difference of stimulated parts of two different fMRI data with the same mask, since our original fMRI is too big to use here for different statistical test. ```{r fig.width = 9, fig.align = "center", warning=FALSE} fmri_pval_comparison_3d(list(pval1, pval2), mask, list(0.05, 0.05), list("scale_p", "scale_p"), multi_pranges=FALSE) ``` ### 2d p-value comparison `fmri_pval_comparison_2d` is to visualize whatever number of p value (generated by different statistical tests on the same fMRI data) to see their difference by 2D plot. For simplicity here we only compare two different 3D array p value data. ```{r fig.width = 9, warning=FALSE} fmri_pval_comparison_2d(list(pval2, pval1), list('t-test', 'HRF'), 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) ``` # Three-phase ROI Analysis ## Phase1: Detect Potential Activated ROI ```{r eval=FALSE} phase1_pval = fmri_ROI_phase1(bigim1_mod, mask_label, mask_dict, stimulus_idx = c(1:160)[rep(c(TRUE,FALSE), c(10,10))])$all_ROI$pval_t ``` ## Phase2: ROI-Based Tensor-on-Tensor Regression ```{r eval=FALSE} phase2_pval = fmri_ROI_phase2(fmridata = bigim1_mod, label_mask = mask_label, label_dict = mask_dict, stimulus_idx = c(1, 21, 41, 61, 81, 101, 121, 141), stimulus_dur = c(10, 10, 10, 10, 10, 10, 10, 10), rrr_rank = 3, fmri.design_order = 2, fmri.stimulus_TR = 3, method = "t_test") ``` ## Phase3: FDR Correction and Spatial Clustering ```{r eval=FALSE} phase3_pval = fmri_post_hoc(phase2_pval , fdr_corr = "fdr", spatial_cluster.thr = 0.05, spatial_cluster.size = 5, show_comparison = FALSE) ``` ## 3D visualization based on the activated areas by regions ```{r fig.width = 9, fig.align = "center", warning=FALSE} fmri_3dvisual_region(TCIU::phase1_pval, label_mask, label_index, label_name, title = "phase1 p-values", rank = "value") ``` ```{r fig.width = 9, fig.align = "center", warning=FALSE} fmri_3dvisual_region(TCIU::phase1_pval, label_mask, label_index, label_name, 5, title = "phase1 top five p-values", rank = "value") ``` ```{r fig.width = 9, fig.align = "center", warning=FALSE} # 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(TCIU::phase2_pval, label_mask, label_index, label_name, title = "phase2 p-values") ``` ```{r fig.width = 9, fig.align = "center", warning=FALSE} fmri_3dvisual_region(list(TCIU::phase2_pval,TCIU::phase3_pval), label_mask, label_index, label_name, title = "phase2&3 p-values") ```