library(oro.nifti)
## oro.nifti 0.11.4
library(fmri) # for NIFTI file input
library(plotly) # for interactive plots
## Loading required package: ggplot2
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:oro.nifti':
##
## slice
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
fMRIURL <- "https://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
# Extract the 1D time course of a specific 3D voxel (say the one at x=30, y=30, z=10):
y_30_30_10 <- fMRIVolume[30, 30, 10,]
# Smooth the noisy fMRI time-series
x1 <- c(1:180)
smoothed <- smooth(y_30_30_10)
# Plot the time series and smoothed versions
plot_ly(x = ~x1) %>%
add_trace(y = ~y_30_30_10, type = "scatter", mode = "lines", name = "Raw fMRI") %>%
add_trace(y = ~smoothed, type = "scatter", mode = "lines", name = "Smoothed fMRI") %>%
layout(xaxis = list(title = "Time"), yaxis = list(title = "Intensity"), legend = list(orientation = 'h'))
# Define a pair of fMRI timeseries at 2 separate voxel locations
y_20_30_15 <- fMRIVolume[20, 30, 15,]
# Plot y_30_30_10 (x-axis) against y_20_30_15 (y-axis)
plot_ly(x = ~y_30_30_10, y = ~y_20_30_15, type = "scatter", mode = "markers")
# Compute their correlation
correlation <- cor(y_30_30_10, y_20_30_15)
cat("Correlation between y_30_30_10 and y_20_30_15:", correlation, "\n")
## Correlation between y_30_30_10 and y_20_30_15: 0.01868463
print(correlation)
## [1] 0.01868463
# Compute the Auto-correlation Function for the first time-series (y_30_30_10)
acf_values <- acf(y_30_30_10, lag.max = 120, plot = TRUE)
# Print ACF values
print(acf_values$acf)
## , , 1
##
## [,1]
## [1,] 1.0000000000
## [2,] -0.1077614222
## [3,] -0.2432317464
## [4,] 0.0378221226
## [5,] -0.0403623874
## [6,] 0.0058705217
## [7,] -0.0619374496
## [8,] -0.0181982623
## [9,] -0.0298783998
## [10,] 0.0279140367
## [11,] 0.0259744190
## [12,] 0.0053631483
## [13,] -0.1034859610
## [14,] -0.0325388803
## [15,] 0.0768761988
## [16,] 0.0216664397
## [17,] 0.0209782452
## [18,] -0.0471600707
## [19,] -0.0117043383
## [20,] -0.0353201399
## [21,] -0.0708648466
## [22,] 0.0853542370
## [23,] -0.0622161632
## [24,] 0.0274522489
## [25,] 0.0533502916
## [26,] -0.0468415918
## [27,] -0.0842554914
## [28,] 0.0636309657
## [29,] 0.0847311539
## [30,] -0.0495634476
## [31,] 0.0986313875
## [32,] -0.0038299888
## [33,] -0.0710433397
## [34,] -0.0018382144
## [35,] 0.0534300850
## [36,] -0.0275014392
## [37,] -0.1180313758
## [38,] 0.0136572111
## [39,] 0.0538813964
## [40,] 0.1273172986
## [41,] -0.0998392905
## [42,] -0.0389591878
## [43,] 0.1238171099
## [44,] -0.0327841183
## [45,] -0.0350535737
## [46,] -0.0055284990
## [47,] 0.0691522364
## [48,] -0.0006798772
## [49,] -0.1193278025
## [50,] 0.0575230352
## [51,] -0.0072369787
## [52,] -0.0285498304
## [53,] 0.0659605384
## [54,] -0.0519117951
## [55,] -0.0295689704
## [56,] 0.0533286441
## [57,] -0.1262660911
## [58,] 0.0330785089
## [59,] 0.0114936087
## [60,] -0.0541712623
## [61,] 0.0582701454
## [62,] -0.0011790083
## [63,] 0.1054611759
## [64,] -0.0189891556
## [65,] -0.0596794656
## [66,] 0.0050608534
## [67,] -0.0129723773
## [68,] 0.0767462201
## [69,] -0.0015217632
## [70,] -0.0101481983
## [71,] -0.0158983540
## [72,] 0.0154257894
## [73,] 0.0065134498
## [74,] 0.0611219048
## [75,] -0.0140573984
## [76,] -0.1131191821
## [77,] 0.0221636558
## [78,] 0.0128150248
## [79,] -0.0312696397
## [80,] 0.0149544827
## [81,] 0.0068783966
## [82,] -0.0304302321
## [83,] -0.1005918454
## [84,] 0.0061514132
## [85,] 0.0610954510
## [86,] -0.0054151361
## [87,] -0.0228060958
## [88,] 0.0480647775
## [89,] 0.0794179844
## [90,] -0.0218380987
## [91,] -0.0202669387
## [92,] -0.0076702099
## [93,] 0.0193763506
## [94,] -0.0438835101
## [95,] 0.0728278766
## [96,] -0.0015143283
## [97,] -0.0641925797
## [98,] 0.0429503029
## [99,] -0.0537821710
## [100,] 0.0072756926
## [101,] 0.0347928464
## [102,] 0.0490644479
## [103,] 0.0099263351
## [104,] -0.0218872515
## [105,] 0.0011975392
## [106,] -0.0627463292
## [107,] 0.0337946030
## [108,] 0.0248675627
## [109,] -0.0053554694
## [110,] -0.0463453708
## [111,] -0.0667653626
## [112,] 0.0334059810
## [113,] 0.0102425235
## [114,] -0.0172150763
## [115,] -0.0086629148
## [116,] 0.0164145890
## [117,] 0.0598018030
## [118,] 0.0278330981
## [119,] -0.0246544862
## [120,] -0.0715586548
## [121,] 0.0135107856
# Fit a linear model
lm_model <- lm(y_30_30_10 ~ y_20_30_15)
summary(lm_model)
##
## Call:
## lm(formula = y_30_30_10 ~ y_20_30_15)
##
## Residuals:
## Min 1Q Median 3Q Max
## -202.850 -35.762 0.349 39.126 132.700
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.392e+04 4.223e+02 32.961 <2e-16 ***
## y_20_30_15 1.665e-02 6.678e-02 0.249 0.803
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 60.77 on 178 degrees of freedom
## Multiple R-squared: 0.0003491, Adjusted R-squared: -0.005267
## F-statistic: 0.06216 on 1 and 178 DF, p-value: 0.8034
# Plot the linear model
plot_ly(x = ~y_20_30_15, y = ~y_30_30_10, type = "scatter", mode = "markers") %>%
add_trace(y = ~(lm_model$coefficients[1] + lm_model$coefficients[2]*y_20_30_15), mode = "lines", name = "Linear Model") %>%
layout(legend = list(orientation = 'h'))
##It is important to determine whether a time-series is auto-correlated and whether the time series at two voxels are correlated because it provides valuable insights into the underlying structure and relationships within the data. Auto-correlation analysis helps identify patterns or dependencies within a single time series, aiding in understanding its behavior and potential predictive capabilities. On the other hand, examining correlations between time series at different voxels enables the exploration of spatial relationships, crucial in fields like neuroscience or environmental science, where understanding how different regions or locations interact over time is essential.