| SOCR ≫ | TCIU Website ≫ | TCIU GitHub ≫ |
Spacekime Analytics (Time Complexity and Inferential Uncertainty)
Basic TCIU Protocol for Predictive Spacekime Analytics using Longitudinal Data
SOCR Team
04/06/2024
This Spacekime TCIU Learning Module presents the core elements of spacekime analytics including:
- Import of repeated measurement longitudinal data,
- Numeric (stitching) and analytic (Laplace) kime-surface reconstruction from time-series data,
- Forward prediction modeling extrapolating the process behavior beyond the observed time-span \([0,T]\),
- Group comparison discrimination between cohorts based on the structure and properties of their corresponding kime-surfaces. For instance, statistically quantify the differences between two or more groups,
- Unsupervised clustering and classification of individuals, traits, and other latent characteristics of cases included in the study,
- Construct low-dimensional visual representations of large repeated measurement data across multiple individuals as pooled kime-surfaces (parameterized 2D manifolds),
- Statistical comparison and quantitative contrasting of kime-surface differences.
1 Preliminary setup
TCIU and other R package dependencies …
2 Longitudinal Data Import
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 brain3 Time-series graphs
3.1 Interactive time-series visualization
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.
4 Kime-series/Kime-surfaces (spacekime analytics protocol)
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.
4.1 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 a large range.
- Generate the
plot_lytext 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). - 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_lyto display in 3D the kime-series as 2D manifolds (kime-surfaces) over the Cartesian domain.
4.2 Function main step: Time-series to Kime-surfaces Mapping
4.2.1 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.
# 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)4.2.2 Structural Data Preprocessing
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.
4.2.3 Intensity Data Preprocessing
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))4.2.4 Generate
plotly labels
Generate 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).
4.2.5 Cartesian representation
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))
}
}4.2.6 Cartesian space interpolation
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)4.3 Numerical Time-series to Kime-Surface Transformation (Stitching)
4.3.1 Generate a long data-frame
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.
4.3.2 Display kime-surfaces
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:
- the On kime-surface
- the Off kime-surface, and
- the difference On - Off kime-surface.
4.4 Analytical Time-series to Kime-Surface Transformation (Laplace)
4.4.1 Data Preparation
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.
4.4.2 Discrete LT
# 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
pNext, 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);
# }
# });
# }")5 Interactive plotly Example
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!
5.1 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.
5.2 Forecasting with time series
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.
6 Motor area detection
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.
6.1 fMRI data simulation
The function fmri_simulate_func() is used to
simulate real-valued fMRI data with specified dimensions
(locations and extents).
6.2 Stimulus detection
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:
- “t-test” and “Wilcoxon-test”, which can be applied to all real*-valued fMRI data,
- “Hotelling’s T2”, “Wilk’s Lambda” and “gLRT” methods, which can be applied to all complex-valued fMRI data,
- “on_off_diff” and “HRF” methods, which 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 the
function
fmri_post_hoc().
6.2.1 Examples
# 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.00007 Motor area visualization
7.1 Visualization and comparison of p-value
7.1.1 3D visualization for p-value
7.1.2 2D visualization for p-value
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.
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))
}
7.2 Comparison of performance of different methods on the same fMRI data
7.2.1 3D p-value comparison
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.
7.2.2 2D p-value comparison
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.
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)
8 Tri-phase ROI-based Spacekime Analytics
Below we demonstrate a three-tier statistical analysis (tri-phase) of regional activation.
Notes:
- For large volumes, these calculations may be expensive (i.e., stats mapping may take a significant time)!
- The following examples are specific to \(4D\) fMRI spacekime analytics. In most situations, these processing steps need to be adapted to the concrete data format and the specific aplication domain. In practice, some of the TCIU package functions may be used as templates to extend this spacekime functionality to alternative data formats and different data structures.
8.1 Phase 1: Detect Potential Activated ROI
First, identify large local (regional areas) assocaited with activations/task stimuli.