SOCR ≫ | TCIU Website ≫ | TCIU GitHub ≫ |
This Spacekime TCIU Learning Module presents the core elements of spacekime analytics including:
TCIU and other R
package dependencies …
In this case, we are just loading some exemplary fMRI data, which is available here.
# pathname <- file.path("./test_data", "subj1_run1_complex_all.mat")
mat1 = readMat("./test_data/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("./test_data/mask.rda") # load the 3D nifti data of the mask
load("./test_data/mask_label.rda") # load the 3D nifti data of the mask with labels
load("./test_data/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("./test_data/hemody.rda") # load the hemodynamic contour of the brain
The TCIU
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 the GTSplot
function from package
TSplotly
.
This example uses a time-series simulation to illustrate how to transform the fMRI time-series at a fixed voxel location into a kime-surface (kime-series).
Notes:
Validate all steps in this time-series to kime-surfaces transformation protocol of simulated data, and finalize this 3D visualization.
Try it with real fMRI data at brain voxel locations associated with the finger-tapping task or musical genre study.
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. 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).plot_ly
to display in 3D the kime-series as 2D
manifolds (kime-surfaces) over the Cartesian domain.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.
# plot Laplacian
# ggfortify::ggdistribution(extraDistr::dlaplace, seq(-pi, pi, 0.01), m=0, s=0.5)
t <- seq(-pi, pi, 0.01)
f_t <- extraDistr::dlaplace(t, mu=0, sigma=0.5)
plot_ly(x=t, y=f_t, type="scatter", mode="lines", name="Laplace Disdtribution") %>%
layout(title="Laplace Disdtribution: Kime-Phase Prior")
# 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)
Here should be included any study-specific data preprocessing. In the case of fMRI series, we may need to split off the individual repeated measurement fMRI time-series from the master (single) fMRI time-series according to the specific event-related fMRI design.
For simplicity, consider a simulated binary stimulus paradigm (e.g., activation (ON) and rest (OFF) event-related design). We can split the \(160\)-timepoint fMRI series into \(80\) ON (Activation) and \(80\) OFF (rest) states, or sub-series, consisting of \(8\) repeats, each of length \(10\) time points, where each time point corresponds to about \(2.5\ s\) of clock time.
In practice, some spatial smoothing (blurring) the 1D time-series or their 2D array (tensor representations as 2D matrices (\(row=repear\), \(column=time\)) to reduce noise and make the data more natural (low-pass filtering, avoiding high-pitch noise). Sometimes, we may also need to temper the magnitude of the raw time-series intensities, which can have a large range.
# Convert the long DF representing fMRI_ON and fMRI_OFF from polar coordinates to Cartesian coordinates
library(spatstat)
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(blur(as.im(matrix_ON), sigma=0.5))
matrix_OFF_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_OFF), sigma=0.5))
plotly
labelsGenerate the plot_ly
text labels that will be shown over
the graph, upon mouse hovering (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. The labels are \(21\times 21\) arrays including the triple
\((t, \phi, fMRIvalue)\). Separate
text-labels are generated for each kime-surface (ON vs. OFF
stimuli).
Generate the \(21\times 21\) kime-domain Cartesian grid by polar-transforming the polar coordinates \((t,\phi)\) into Cartesian counterparts \((x,y)\).
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))
}
}
Interpolate the fMRI intensities, natively anchored at polar (kime) coordinates) \(\left (\underbrace{t}_{time},\underbrace{\phi}_{repeat} \right)\), into Cartesian coordinates \(fMRI(x,y), \forall -11\leq x,y \leq 11, x^2+y^2\leq 10\). Note that for specificity, we hard-coded polar-grid parameterization to time \(t\in\{1,2,3, \cdots,10\},\ |T|=10\) and phase \(\phi=seq\left (from=-\pi, to=\pi, step=\frac{2\pi}{20}\right )\in\{-3.1415927, -2.8274334, \cdots, 0.0,\cdots, 2.8274334, 3.1415927 \},\ |\Phi|=21\).
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)
Convert the entire dataset into a long data-frame, i.e., construct a data-frame (DF) with \(\underbrace{160}_{total}=\underbrace{2}_{binary\\ On/Off}\times \underbrace{8}_{repeats}\times \underbrace{10}_{timepoints}\) rows and \(4\) columns representing time (\(1:10\)), phases (\(8\)), states (\(2\)), and fMRI_value (Complex or Real time-series intensity). Then, transform the long DF representing the fMRI_ON and fMRI_OFF time-sources from their native (old) polar coordinates to the (new) Cartesian coordinates, using polar transformations.
Use plot_ly
to display in 3D the kime-series as 2D
manifolds (kime-surfaces) over the Cartesian domain. In this specific
binary case-study, we demonstrate the following \(3\) kime-surface reconstructions:
First plot the simulated On (stimulus) and Off (rest) fMRI time-series at a fixed voxel location \((44,42,33)\in\mathbb{R}^3\), along with the averaged (pooled) On and Off signal over all \(8\) repeats in the single run (epoch) of \(160\) time-points.
# Extract just the On kime-series at voxel (44,42,33), each time-series has 8 repeats!
onFMRISeries <-
bigim1_mod[44,42,33, c(1:10, 21:30, 41:50, 61:70, 81:90, 101:110, 121:130, 141:150)]
# the corresponding Off kime-series at voxel (44,42,33) will be the temporal complement
offFMRISeries <-
bigim1_mod[44,42,33, -c(1:10, 21:30, 41:50, 61:70, 81:90, 101:110, 121:130, 141:150)]
t_indx <- seq(1, 80, 1) # index of time (not actual time value)
f_On <- onFMRISeries
f_Off <- offFMRISeries
vline <- function(x=0, color = "lightgray") {
list(type="line", y0=0, y1=1, yref="paper", name="time break points",
opacity = 0.5, x0=x, x1=x, line=list(color=color))
}
hline <- function(y = 0, color = "blue") {
list(type="line", x0=0, x1=1, xref="paper", name="intensity break points",
opacity = 0.3, y0=y, y1=y, line=list(color=color))
}
plot_ly() %>%
# On fMRI time-series
add_trace(x=t_indx, y=f_On, type="scatter", mode="lines", name="On time-series") %>%
# On fMRI time-series
add_trace(x=t_indx, y=f_Off, type="scatter", mode="lines", name="Off time-series") %>%
# Repeated measurement break points
# add_trace(x=c(10,20,30,40,50,60,70,80), y=f_Off, type="scatter", mode="lines", name="Off time-series") %>%
layout(title="3D fMRI Simulation: On & Off Time-series at Voxel(44,42,33)",
shapes = list(
vline(10),vline(20),vline(30),vline(40),vline(50),vline(60),vline(70),vline(80)),
legend = list(orientation='h', y=-0.2))
# Compute and plot against each other the Average On and Average Off time-series
seriesAvg <- function(f=0, period=0) {
tempAvg <- rep(0,period) # initialize avg signal
for (i in c(1:period)) {
counter =0 # emperically count the number of repeated samples (here it's 8)
for (j in c(1:length(f))) {
if (j %% period == i-1) { # e.g., 159 %% 10 # [1] 9
tempAvg[i] = tempAvg[i] + f[j]
counter = counter + 1
}
}
tempAvg[i] = tempAvg[i]/counter
}
return(tempAvg)
}
period <- 10
onAvg <- seriesAvg(f=f_On, period=period)
offAvg <- seriesAvg(f=f_Off, period=period)
plot_ly() %>%
# Average On fMRI time-series
add_trace(x=c(1:period), y=onAvg, type="scatter", mode="lines", name="On Average") %>%
# On fMRI time-series
add_trace(x=c(1:period), y=offAvg, type="scatter", mode="lines", name="Off Average") %>%
# Repeated measurement break points
layout(title=
"3D fMRI Simulation: On & Off Time-series\n (averaged over all 8 repeats) at Voxel (44,42,33)",
legend = list(orientation='h', y=-0.2))
Next, we’ll define and apply the Laplace Transform (LT) and its inverse (ILT) and use them to show the analytical kime-surface reconstruction.
# create the LT
NuLT = function(datax, datay, inputz, k = 3, fitwarning = FALSE,
mirror = FALSE, range = 2*pi) {
datax = as.numeric(datax)
datay = as.numeric(datay)
n = length(datax)
x1 = n/(n+0.5)*((datax-min(datax))/(max(datax)-min(datax)))*range
if(mirror){
x1 = c(x1,rev(2*range - x1))/2
n = 2*n
datay = c(datay, rev(datay))
plot(x1, datay)
}
#generate the coefficients in indefinite integral of t^n*exp(-zt)
coef = 1;
coefm = as.matrix(coef)
for(i in 1:k){
coefm = cbind(coefm,0)
coef = c(coef*i,1)
coefm = rbind(coefm,coef)
}
# these coefficients ordered by ^0, ^1, ^2, ... in column format
# compute 1, z, z^2...,z^k
zz = cbind(1,inputz)
zt = inputz
for (i in 2:k){
zt = zt*inputz
zz = cbind(zz,zt)
}
zd = zt*inputz
# compute 1, x, x^2...,x^k
tx = x1;
xm = cbind(1,x1)
for (i in 2:k){
tx = tx*x1
xm = cbind(xm,tx)
}
# sum over intervals
result = 0*inputz
ii = 1
while(ii+k <= n) {
A = xm[seq(ii,ii+k),c(0:k+1)]
b = datay[seq(ii,ii+k)]
# polyfit might be faster when using polynomial basis, while matrix inverse, `solve()`,
# is the more general approach that works for any function basis
polyc = as.numeric(solve(A,b))
#ordered by ^0, ^1, ^2, ... in column format
# Enter a new function variable qualityCheck=FALSE
# check fit quality; this step can be skipped for speed/efficiency
# if (qualityCheck) { .... }
if (fitwarning){
xx = seq(A[1,2],A[k+1,2],length.out = 100);
yy = polyval(rev(polyc),xx)
if(max(abs(yy-mean(b)))>2*max(abs(b-mean(b)))){
print(c("Warning: Poor Fit at ",ii,", Largest Deviation is",
max(abs(yy-mean(b)))))
print(c("Spline Polynomial is", polyc),3);
#print(c(polyval(rev(polyc),A[,2]),b))
plot(xx, yy, main="Polynomial fit", ylab="", type="l", col="blue")
lines(A[,2],b, col="red")
legend("topleft",c("fit","data"),fill=c("blue","red"))
print(" ")
}
}
# Use vector/matrix operations to avoid looping,
# some of the expressions look weird
# May need to actually compare the efficiency/speed of
# vector based vs. standard numeric calculations
m1 = t(t(polyc*coefm)*A[1,])
m11 = as.numeric(tapply(m1, col(m1)-row(m1), sum))[0:k+1]
m2 = t(t(polyc*coefm)*A[k+1,])
m22 = as.numeric(tapply(m2, col(m2)-row(m2), sum))[0:k+1]
intgl = (exp(-inputz*A[1,2])*colSums(t(zz)*m11)-
exp(-inputz*A[k+1,2])*colSums(t(zz)*m22))/zd
result = result+intgl
ii=ii+k
}
# Computations over the last interval
if(ii < n){
nk = n-ii;
A = xm[seq(ii,ii+nk),c(0:nk+1)]
b = datay[seq(ii,ii+nk)]
nc = as.numeric(solve(A,b))
nc = c(nc,seq(0,0,length.out = k-nk))
A = xm[seq(ii,ii+nk),]
m1 = t(t(nc*coefm)*A[1,])
m11 = as.numeric(tapply(m1, col(m1)-row(m1), sum))[0:k+1]
m2 = t(t(nc*coefm)*A[nk+1,])
m22 = as.numeric(tapply(m2, col(m2)-row(m2), sum))[0:k+1]
# cc = colSums(coefm*polyc)
intgl = (exp(-inputz*A[1,2])*colSums(t(zz)*m11)-
exp(-inputz*A[nk+1,2])*colSums(t(zz)*m22))/zd
result = result+intgl
}
#offset = 0.05*pi
#result = result+datay[n]*(exp(-2*pi*inputz)-exp(-(2*pi+offset)*inputz))/inputz
return(result)
}
Let’s test the discrete LT using the \(\sin(x)\) function.
datax = seq(1,200)
n = length(datax)
x1 = n/(n+0.5)*((datax-min(datax))/(max(datax)-min(datax)))*(2*pi)
datay = sin(x1) # test function!!!
Lout = 61
range_limit <- 20
x2 <- seq(from = 0.1, to = range_limit, length.out = Lout)
y2 <- seq(from = 0, to = range_limit, length.out = Lout)
x2_g = x2 %o% seq(1,1, length.out = Lout)
y2_g = seq(1,1, length.out = Lout)%o%y2
z2_grid = x2 %o% y2
argz = as.vector(x2_g + 1i*y2_g)
LTz = NuLT(x1, datay, argz)
rec1 = matrix(LTz,nrow = Lout)
recm = abs(rec1)
recr = Re(rec1)
reci = Im(rec1)
ph = reci/recm
surf_color <- atan2(reci,recr)
# colorscale = cbind(seq(-pi, pi, by=1/10, rainbow(11)))
p <- plot_ly(hoverinfo="none", showscale = TRUE) %>%
add_trace(x = x2_g, y = y2_g, z = log(recm),
# F-magnitude or Re(F), # z = Im(z2_grid), # Real or Imaginary part of F(t)
surfacecolor=surf_color, # colorscale=colorscale, #Phase-based color
type = 'surface', opacity=1, visible=T) %>%
layout(title =
"Laplace transform of the Time-Series dataa [5,] NMHC(GT), Color = phase(Z)",
showlegend = TRUE,
scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1, z=1.0),
xaxis=list(title= "Real"), yaxis=list(title = "Imaginary"),
zaxis=list(title= "Z = log(mag(Lf(x+iy)))")) ) # 1:1:1 aspect ratio
p
Next, apply the discrete LT to the average-On
(onAvg
) and average-Off signals
(offAvg
), interpolating from their original size, \(n=10\), to a new supersampled size
\(n=200\), and transforming the time
support from \(t\in[1:10]\) in
increments of \(\Delta t=1\),
to \(t' \in[0,2\pi)\), in
increments of \(\Delta
t'=\frac{n}{n+0.5}\times \frac{1}{(200-1)2\pi}.\)
This numerical longitudinal data preprocessing is done purely to establish some homologies in the structure of the LT domain, i.e., the input space signals (time-series), and the LT Range, i.e., the output space manifold (kime-surface). See the DSPA2 signal interpolation appendix to find out how to regularize either regularly (equally-spaced) sampled longitudinal data or irregularly (unequally-spaced) sampled longitudinal data.
datax = seq(1,200)
n = length(datax)
x_old <- c(1:10)
x1 = n/(n+0.5)*((datax-min(datax))/(max(datax)-min(datax)))*(2*pi)
# longitudinal data series are: onAvg & offAvg
xnew <- x1
spl_onAvg <- spline(x=x_old, y=onAvg, xmin=min(x_old), xmax=max(x_old), n=200)
spl_offAvg <- spline(x=x_old, y=offAvg, xmin=min(x_old), xmax=max(x_old), n=200)
plot_ly(type="scatter", mode="markers") %>%
add_trace(x=x_old, y=onAvg, name="Raw ON-Avg",
marker=list(size=20, color="lightblue", sizemode="area")) %>%
add_trace(x=x_old, y=offAvg, name="Raw OFF-Avg",
marker=list(size=20, color="orange", sizemode="area")) %>%
add_trace(x=spl_onAvg$x, y=spl_onAvg$y, type="scatter", mode="markers+lines",
name="Spline ON-Avg Model", marker=list(size=8, color="blue"),
line = list(color = "blue", width = 4)) %>%
add_trace(x=spl_offAvg$x, y=spl_offAvg$y, type="scatter", mode="markers+lines",
name="Spline OFF-Avg Model", marker=list(size=8, color="red"),
line = list(color = "red", width = 4)) %>%
layout(title="Spline Modeling of 1D ON and OFF fMRI data\n (averaged over repeated samples)",
legend = list(orientation = 'h', y=-0.2)) %>%
hide_colorbar()
# datay = sin(x1) # test function!!!
# Define a time-series to kime-surface plotting function
plotKimeSurface <- function(datay = sin(x1),
title="Laplace transform of the Time-Series, Color = phase(Z)") {
Lout = 61
range_limit <- 20
x2 <- seq(from = 0.1, to = range_limit, length.out = Lout)
y2 <- seq(from = 0, to = range_limit, length.out = Lout)
x2_g = x2 %o% seq(1,1, length.out = Lout)
y2_g = seq(1,1, length.out = Lout)%o%y2
z2_grid = x2 %o% y2
argz = as.vector(x2_g + 1i*y2_g)
LTz = NuLT(x1, datay, argz)
rec1 = matrix(LTz,nrow = Lout)
recm = abs(rec1)
recr = Re(rec1)
reci = Im(rec1)
ph = reci/recm
surf_color <- atan2(reci,recr)
# colorscale = cbind(seq(-pi, pi, by=1/10, rainbow(11)))
p <- plot_ly(hoverinfo="none", showscale = TRUE) %>%
add_trace(x = x2_g, y = y2_g, z = log(recm),
# F-magnitude or Re(F), # z=Im(z2_grid), # Real or Imaginary part of F(t)
surfacecolor=surf_color, # colorscale=colorscale, #Phase-based color
type = 'surface', opacity=1, visible=T) %>%
layout(title = title, showlegend = TRUE,
scene = list(aspectmode="manual", aspectratio=list(x=1, y=1, z=1.0),
xaxis=list(title= "Real"), yaxis=list(title = "Imaginary"),
zaxis=list(title= "Z = log(mag(Lf(x+iy)))")) ) # 1:1:1 aspect ratio
return(p)
}
pOn <- plotKimeSurface(datay=spl_onAvg$y,
title="ON fMRI ((44,42,33): Laplace transform of Avg Time-Series, Color = phase(Z)")
pOff <- plotKimeSurface(datay=spl_offAvg$y,
title="OFF fMRI ((44,42,33): Laplace transform of Avg Time-Series, Color = phase(Z)")
combineWidgets(pOn, pOff)
# synch both 3D scenes in 1 3D plot
# install.packages("manipulateWidget")
# subplot(pOn, pOff, nrows = 1, shareX = TRUE)
# main_plot <- subplot(pOn, pOff, nrows = 2, shareX = TRUE, margin = 0.06) %>%
# layout(title = "Kimesurfaces On (activation) and Off (rest) fMRI voxel (44,42,33)",
# scene = list(domain = list(x = c(0, 0.5), y = c(0.5, 1)), aspectmode = "cube"),
# scene2 = list(domain = list(x = c(0.2, 0.7), y = c(0.5, 1)), aspectmode = "cube"))
#
# main_plot %>%
# htmlwidgets::onRender(
# "function(x, el) {
# x.on('plotly_relayout', function(d) {
# const camera = Object.keys(d).filter((key) => /\\.camera$/.test(key));
# if (camera.length) {
# const scenes = Object.keys(x.layout).filter((key) => /^scene\\d*/.test(key));
# const new_layout = {};
# scenes.forEach(key => {
# new_layout[key] = {...x.layout[key], camera: {...d[camera]}};
# });
# Plotly.relayout(x, new_layout);
# }
# });
# }")
For most of these examples, users may need to review/explore/extend the core functionailty in the TCIU-package source-code (CRAN) to customize the TCIU functionailty (GitHub) for the application specific needs!
The function fmri_image()
is used to create images for
front view, side view, and top view of the fMRI image.
The function fmri_ts_forecast()
fits the ARIMA
models for each voxel (spatial volume element) location.
This function is based on the TSplot_gen()
function from
package TSplotly
.
If there are concrete spatial locations (regions, voxels, coordiantes) that are of specific interest, then one needs to focus on these locations. We’ll demonstrate this using fMRI brain activation (simulated) in the somatosensory (motor) cortex.
The function fmri_simulate_func()
is used to
simulate real-valued fMRI data with specified dimensions
(locations and extents).
The integrated function fmri_stimulus_detect()
is
designed to apply multiple analytical methods for activation detecttion,
in this case in the (human brain) motor area.
Examples of parametric and non-parametric statistical tests already built in include:
fmri_post_hoc()
.# 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)
Summarize the corresponding p-values.
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.8987 1.0000 0.8463 1.0000 1.0000
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000096 0.8598920 1.0000000 0.8579245 1.0000000 1.0000000
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000025 0.1877401 1.0000000 0.7476717 1.0000000 1.0000000
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 1.0000 1.0000 0.9971 1.0000 1.0000
Generate 2D plots of the 2D p-value (projection images) in sagittal, coronal and axial views.
Plot without hemodynamic contours.
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 contours.
The function fmri_pval_comparison_3d()
is used to
visualize two p-value maps showing their differences in detecting
stimulated brain areas in 3D scenes. Since our original fMRI is too big
to use here for different statistical tests, in this example, we compare
the differences of stimulated parts of two different fMRI data using the
same mask.
The function fmri_pval_comparison_2d()
displays the
p-values (generated by different statistical tests on the same fMRI
data) exposing their difference via 2D plots. For simplicity here we
only compare two different 3D arrays of p-values.
Below we demonstrate a three-tier statistical analysis (tri-phase) of regional activation.
Notes:
First, identify large local (regional areas) assocaited with activations/task stimuli.