The goal of fmri is to perform an fMRI analysis as described in Tabelow et al. (2006) https://doi.org/10.1016/j.neuroimage.2006.06.029, Polzehl et al. (2010) https://doi.org/10.1016/j.neuroimage.2010.04.241, Tabelow and Polzehl (2011) https://doi.org/10.18637/jss.v044.i11.
You can install the released version of fmri from CRAN with:
install.packages("fmri")
And the development version from GitHub with:
# install.packages("devtools")
::install_github("muschellij2/fmri") devtools
This is a basic example which shows you how to solve a common problem:
library(fmri)
## basic example code
<- function(y,h=1) {
gkernsm <- function(d) {
grid <- d%/%2+1
d0 <- seq(0,1,length=d0)
gd if (2*d0==d+1) gd <- c(gd,-gd[d0:2]) else gd <- c(gd,-gd[(d0-1):2])
gd
}<- dim(y)
dy if (is.null(dy)) dy<-length(y)
<- length(dy)
ldy if (length(h)!=ldy) h <- rep(h[1],ldy)
<- switch(ldy,dnorm(grid(dy),0,2*h/dy),
kern outer(dnorm(grid(dy[1]),0,2*h[1]/dy[1]),
dnorm(grid(dy[2]),0,2*h[2]/dy[2]),"*"),
outer(outer(dnorm(grid(dy[1]),0,2*h[1]/dy[1]),
dnorm(grid(dy[2]),0,2*h[2]/dy[2]),"*"),
dnorm(grid(dy[3]),0,2*h[3]/dy[3]),"*"))
<- kern/sum(kern)
kern <- sum(kern^2)
kernsq list(gkernsm = convolve(y,kern,conj=TRUE),kernsq=kernsq)
}
<- function(){
create.mask <- array(0,dim=c(65,65,26))
mask 5:10,5:10,] <- 1
mask[7:8,7:8,] <- 0
mask[8:10,8:10,] <- 0
mask[14:17,14:17,] <- 1
mask[16:17,16:17,] <- 0
mask[21:23,21:23,] <- 1
mask[22:23,23,] <- 0
mask[23,22,] <- 0
mask[27:28,27:28,] <- 1
mask[28,28,] <- 0
mask[5:7,29:33,] <- 1
mask[7,32:33,] <- 0
mask[14:15,30:33,] <- 1
mask[15,30,] <- 0
mask[21,31:33,] <- 1
mask[22,33,] <- 1
mask[27,32:33,] <- 1
mask[29:33,5:7,] <- 1
mask[32:33,7,] <- 0
mask[30:33,14:15,] <- 1
mask[30,15,] <- 0
mask[31:33,21,] <- 1
mask[33,22,] <- 1
mask[32:33,27,] <- 1
mask[34:65,1:33,] <- mask[32:1,1:33,]
mask[1:33,34:65,] <- mask[1:33,32:1,]
mask[34:65,34:65,] <- mask[32:1,32:1,]
mask[
mask
}
<- function(signal=1.5,efactor=1.2){
create.sig <- array(0,dim=c(65,65,26))
sig 29:37,38:65,] <- signal
sig[38:65,38:65,] <- signal * efactor
sig[38:65,29:37,] <- signal * efactor^2
sig[38:65,1:28,] <- signal * efactor^3
sig[29:37,1:28,] <- signal * efactor^4
sig[1:28,1:28,] <- signal * efactor^5
sig[1:28,29:37,] <- signal * efactor^6
sig[1:28,38:65,] <- signal * efactor^7
sig[* create.mask()
sig
}# some values describing the data
<- 1.5
signal <- 20
noise <- .3
arfactor
# maximaum bandwidth for adaptive smoothing
<- 3.06
hmax
# datacube dimension
<- 65
i <- 65
j <- 26
k <- 107
scans
# define needed arrays
<- array(0,dim=c(i,j,k,scans))
ttt <- array(0,dim=c(i,j,k))
sig
# create the mask for activation
<- create.mask()
mask
# assign amplitudes of signals to activated areas
<- create.sig(signal) sig
# expected BOLD response for some stimulus
<- signal * fmri.stimulus(scans, c(18, 48, 78), 15, 2)
hrf
# create time series
dim(sig) <- c(i*j*k,1)
dim(hrf) <- c(1,scans)
<- sig %*% hrf
sig4 dim(sig) <- c(i,j,k)
dim(sig4) <- c(i,j,k,scans)
RNGversion("3.5.0")
#> Warning in RNGkind("Mersenne-Twister", "Inversion", "Rounding"): non-uniform
#> 'Rounding' sampler used
# create noise with spatial and temporal correlation
set.seed(1)
<- rnorm(i*j*k*scans,0,noise)
noisy4 dim(noisy4) <- c(i,j,k,scans)
for (t in 2:scans) noisy4[,,,t] <- noisy4[,,,t] + arfactor*noisy4[,,,t-1]
for (t in 1:scans) noisy4[,,,t] <- gkernsm(noisy4[,,,t],c(0.8,0.8,0.4))$gkernsm
# finally we got the data
<- sig4 + noisy4
ttt <- list(ttt=writeBin(as.numeric(ttt),raw(),4),
data1 dim=c(i,j,k,scans),weights=c(1,1,2),mask=array(1,c(i,j,k)),
delta = rep(1, 4))
class(data1) <- "fmridata"
# create design matrix and estimate parameters from linear model
<- fmri.stimulus(scans, c(18, 48, 78), 15, 2)
hrf <- fmri.design(hrf)
z <- fmri.lm(data1,z) spm
# adaptively smooth the spm
<- fmri.smooth(spm,hmax=hmax,lkern="Gaussian")
resultaws <- fmri.pvalue(resultaws)
detectaws <- apply(detectaws$pvalue<0.05,c(1,2),sum)
pmask
# smooth non adaptively
<- fmri.smooth(spm,hmax=hmax,adaptation="none",lkern="Gaussian")
resultnonaws <- fmri.pvalue(resultnonaws)
detectnonaws <- apply(detectnonaws$pvalue<0.05,c(1,2),sum)
npmask
# at last show some nice images
par(mfrow=c(1,3),mar=c(2.5,2.5,3,0.5))
image(1:i,1:j,sig[,,1],col=grey(0:255/255),xlab="",ylab="",main="True activation")
image(1:i,1:j,pmask,xlab="",ylab="",main="Detected using AWS")
image(1:i,1:j,npmask,xlab="",ylab="",main="Detected using Gaussian filter")