1 Longitudinal Data (Timeseries/Kimeseries) Example

This timeseries demonstration shows the effects of indexing timeseries (univariate) data only using time and compares the representation of timeseries and kimeseries, which has profound impact on the subsequent data analytics. TO keep this application grounded, we will use real 4D fMRI data (\(x=64 \times y=64\times z=21\times t=180\)), but only focus on one spatial location (\({\bf{x}}=(x,y,z)\in R^3\)). More details are provided in DSPA Chapter 3.

# install EBImage
# source("https://bioconductor.org/biocLite.R")
# biocLite("EBImage")
library(EBImage)
require(brainR)
library(spatstat) 

# 1. download the 4D fMRI data
fMRIURL <- "http://socr.umich.edu/HTML5/BrainViewer/data/fMRI_FilteredData_4D.nii.gz"
fMRIFile <- file.path(tempdir(), "fMRI_FilteredData_4D.nii.gz")
download.file(fMRIURL, dest=fMRIFile, quiet=TRUE)
(fMRIVolume <- readNIfTI(fMRIFile, reorient=FALSE))

## NIfTI-1 format
##   Type            : nifti
##   Data Type       : 4 (INT16)
##   Bits per Pixel  : 16
##   Slice Code      : 0 (Unknown)
##   Intent Code     : 0 (None)
##   Qform Code      : 1 (Scanner_Anat)
##   Sform Code      : 0 (Unknown)
##   Dimension       : 64 x 64 x 21 x 180
##   Pixel Dimension : 4 x 4 x 6 x 3
##   Voxel Units     : mm
##   Time Units      : sec

# dimensions: 64 x 64 x 21 x 180 ; 4mm x 4mm x 6mm x 3 sec 

fMRIVolDims <- dim(fMRIVolume); fMRIVolDims

## [1]  64  64  21 180

# time_dim <- fMRIVolDims[4]; time_dim ## 180

# 2. extract the time-corse of one voxel (25, 25, 12)  # 64  64 21 180
x1 <- c(1:180)
y1 <- loess(fMRIVolume[25, 25, 12, ]~ x1, family = "gaussian")

# windows(width=10, height=8) # For Windows Users
# dev.new(height = 10, width = 8) # For Mac OS Users

plot(fMRIVolume[25, 25, 12, ], type='l', 
     main="Time Series of 3D Voxel \n (x=25, y=25, z=12)", col="blue", 
     xlab = "Time", ylab = "fMRIVolume[25, 25, 12, ] Intensities")
lines(x1, smooth(fMRIVolume[25, 25, 12, ]), col = "red", lwd = 2)
lines(ksmooth(x1, fMRIVolume[25, 25, 12, ], kernel = "normal", bandwidth = 5), col = "green", lwd = 3)
legend("bottomright", legend=c("(raw) fMRI", "smooth(fMRI)", "ksmooth(fMRI"),
       col=c("blue", "red", "green"), lty=1, lwd=4, cex=1.1, y.intersp=0.6,
       x.intersp=1.0, title = "Voxel (25, 25, 12)", bty = "n")

## quartz_off_screen 
##                 2

# 3. FT of 1D time-course
y2 <- fMRIVolume[25, 25, 12, ]
X2 = fft(y2); plot(fftshift1D(log(Re(X2)+2)), main = "log(fftshift1D(Re(FFT(timeseries))))") 

## Warning in log(Re(X2) + 2): NaNs produced

X2_mag <- sqrt(Re(X2)^2+Im(X2)^2); plot(log(fftshift1D(Re(X2_mag))), main = "log(Magnitude(FFT(timeseries)))") 

X2_phase <- atan2(Im(X2), Re(X2)); plot(fftshift1D(X2_phase), main = "Shift(Phase(FFT(timeseries)))")

##  Implicit Automated IFT
hat_X2 = Re(fft(X2, inverse = T)/length(X2))
plot(hat_X2, main = "Re(IFT(FFT(timeseries)))")   # point plot represents the IFT(FT(data))
lines(x1, y2, col = "red", lwd = 2)               # lines represent the raw data time-course (should coincide)

## Manually invert the FT (IFT) using the magnitudes and phases
Real2 = X2_mag * cos(X2_phase)
Imaginary2 = X2_mag * sin(X2_phase)
man_hat_X2 = Re(fft(Real2 + 1i*Imaginary2, inverse = T)/length(X2))
plot(man_hat_X2, type="l", lwd=2, main = "Manual IFT (Magnitude, Phase) Synthesis")
lines(x1, y2, col = "red", lty=2, lwd = 2)  # overlay original data to confirm PERFECT manual reconstruction

# IFT reconstruction (synthesis) using fMRI-Magnitude and Nil-Phase
Real_phase0 = X2_mag * cos(0)
Imaginary_phase0 = X2_mag * sin(0)
ift_NilPhase_X2mag = Re(fft(Real_phase0 + 1i*Imaginary_phase0, inverse = T)/length(X2))
# windows(width=10, height=8)
# dev.new(height = 10, width = 8) # For Mac OS Users
plot(ift_NilPhase_X2mag, col="red", type="l", lty=2, lwd=2, 
     ylim=c(10270, 10500), main=sprintf('Signal Synthesis: IFT (Magnitude=Real, Phase=Nil)\n Correlation(Real, Recon) = %s', format(cor(ift_NilPhase_X2mag, y2), digits=3)), xlab = "Time", 
     ylab = "IFT(NilPhase_Reconstruction)[25, 25, 12, ]", cex=1.3)
lines(x1, y2, col = "black", lwd = 2)  # overlay original data to show Nil-Phase effects on reconstruction
legend("top", bty="n", legend=c("(raw) fMRI", "Nil-Phase (Time-only) Reconstruction"),
       col=c("black", "red"), lty=c(1,2), lwd=c(2,2), cex=1.3, x.intersp=0.8)

# cor(ift_NilPhase_X2mag, y2)

## quartz_off_screen 
##                 2

# IFT reconstruction (synthesis) using fMRI-Magnitude and Random-Phase
set.seed(1234)
rand_Phase <- runif(length(X2_mag), -pi, pi)
#  rand_Phase <- runif(1, -pi, pi)
Real_phaseRand = X2_mag * cos(rand_Phase)
Imaginary_phaseRand = X2_mag * sin(rand_Phase)
ift_RandPhase_X2mag = Re(fft(Real_phaseRand + 1i*Imaginary_phaseRand, inverse = T)/length(X2))
plot(ift_RandPhase_X2mag-mean(ift_RandPhase_X2mag) + mean(y2), col="red", type="l", lty=2, lwd=2, 
     ylim=c(10300, 10500), main=sprintf('Signal Synthesis: IFT (Magnitude=Real, Phase=Random)\n Correlation(Real, Recon) = %s', format(cor(ift_RandPhase_X2mag, y2), digits=3)))
lines(x1, y2, col = "black", lwd = 2)  # overlay original data to show Rand-Phase effects on reconstruction
legend("top", bty="n", legend=c("(raw) fMRI", "Random-Phase Reconstruction"),
       col=c("black", "red"), lty=c(1,2), lwd=c(2,2), cex=0.9, x.intersp=0.5)

# IFT reconstruction (synthesis) using fMRI-Magnitude and Phase from a different voxel location (28, 22, 10)
y3 <- fMRIVolume[28, 22, 10, ]
X3 = fft(y3)  # ;  plot(fftshift1D(log(Re(X3)+2)), main = "log(fftshift1D(Re(FFT(timeseries))))") 
# X3_mag <- sqrt(Re(X3)^2+Im(X3)^2); plot(log(fftshift1D(Re(X3_mag))), main = "log(Magnitude(FFT(timeseries)))") 
neighbor_Phase <- atan2(Im(X3), Re(X3)) # ; plot(fftshift1D(X3_phase), main = "Shift(Phase(FFT(timeseries)))")

Real_phaseNeighbor = X2_mag * cos(neighbor_Phase)
Imaginary_phaseNeighbor = X2_mag * sin(neighbor_Phase)
ift_NeighborPhase_X2mag = Re(fft(Real_phaseNeighbor + 1i*Imaginary_phaseNeighbor, inverse = T)/length(X2))
plot(ift_NeighborPhase_X2mag, col="red", type="l", lty=2, lwd=2, 
     ylim=c(10300, 10500), main=sprintf('Signal Synthesis: IFT (Magnitude=Real(25, 25, 12), Phase=Voxel Neighbor(28, 22, 10))\n Correlation(Real, Recon) = %s', format(cor(ift_NeighborPhase_X2mag, y2), digits=3)))
lines(x1, y2, col = "black", lwd = 2)  # overlay original data to show Rand-Phase effects on reconstruction
legend("top", bty="n", legend=c("(raw) fMRI", "Neighbor-Derived-Phase Reconstruction"),
       col=c("black", "red"), lty=c(1,2), lwd=c(2,2), cex=0.9, x.intersp=0.5)

# IFT reconstruction (synthesis) using fMRI-Magnitude and Phase from a highly correlated voxel location
set.seed(1234)
y4 <- y2 + rnorm(n=length(y2), 0, 40) #; plot(y2, y4); cor(y2, y4)
X4 = fft(y4)  # ;  plot(fftshift1D(log(Re(X3)+2)), main = "log(fftshift1D(Re(FFT(timeseries))))") 
# X4_mag <- sqrt(Re(X4)^2+Im(X4)^2); plot(log(fftshift1D(Re(X4_mag))), main = "log(Magnitude(FFT(timeseries)))") 
corr_Phase <- atan2(Im(X4), Re(X4)) # ; plot(fftshift1D(X4_phase), main = "Shift(Phase(FFT(timeseries)))")

Real_phaseCorr = X2_mag * cos(corr_Phase)
Imaginary_phaseCorr = X2_mag * sin(corr_Phase)
ift_CorrPhase_X2mag = Re(fft(Real_phaseCorr + 1i*Imaginary_phaseCorr, inverse = T)/length(X2))

# windows(width=10, height=8)
#dev.new(height = 10, width = 8) # For Mac OS Users

plot(ift_CorrPhase_X2mag, col="red", type="l", lty=2, lwd=2, 
     ylim=c(10250, 10510), main=sprintf('Signal Synthesis: IFT (Magnitude=Real, Phase=Highly-Correlated (%s) Voxel\n Correlation(Real, Recon) = %s', format(cor(y4, y2), digits=3),
                                        format(cor(ift_CorrPhase_X2mag, y2), digits=3)),
     xlab = "Time", ylab = "IFT(Highly-Corr_Reconstruction)[25, 25, 12, ] Intensities")
lines(x1, y2, col = "black", lwd = 2)  # overlay original data to show Rand-Phase effects on reconstruction
legend("top", bty="n", legend=c("(raw) fMRI", "Correlated-Voxel-Phase Reconstruction"),
       col=c("black", "red"), lty=c(1,2), lwd=c(2,2), cex=1.2, x.intersp=0.8)