# R/connectivityFunctions.R In arf3DS4: Activated Region Fitting, fMRI data analysis (3D).

```#############################################
# arf3DS4 S4 CONNECTIVITY FUNCTIONS			#
# Wouter D. Weeda			                #
# University of Amsterdam					#
#############################################

#[CONTAINS]
#fitConnectivity			[user]
#makeSingleTrialEvents		[user]
#gamma.fmri
#convol.fmri
#correlationTest
#tsCor
#partialCor
#differenceCor				[user]
#roiConnectivity			[user]
#createFuncs
#setFuncFiles
#setFuncTimings				[user]
#saveFunc					[user]

fitConnectivity <-
function(arfmodel,funcfilename='single_events.nii.gz')
#make a connectivity estimate for a model solution
{
regs = 1:.model.regions(arfmodel)

#make model design matrix
X = matrix(NA,n,length(regs))
theta = matrix(.model.estimates(arfmodel),.model.params(arfmodel))

#set theta amplitudes to one (all)
theta[10,]=rep(1,.model.regions(arfmodel))

p=1
for(i in regs) {
thetavec = as.vector(theta[,i])
p=p+1
}

funcvolume = .fmri.data.datavec(funcdata)
dim(funcvolume) = c(.fmri.data.dims(funcdata)[2],.fmri.data.dims(funcdata)[3],.fmri.data.dims(funcdata)[4],.fmri.data.dims(funcdata)[5])

b =  matrix(NA,length(regs),.fmri.data.dims(funcdata)[5])
Xp = solve(t(X)%*%X)%*%t(X)

p=1
for(tm in 1:.fmri.data.dims(funcdata)[5]) {
y = as.vector(funcvolume[,,,tm])
b[,p] = Xp%*%y
p=p+1
}

#get correlations
arfcor <- new('arfcorrelation')
.arfcorrelation.timebyreg(arfcor) <- t(b)
arfcor <- correlationTest(arfcor)

#save correlations
fn = paste(.model.modeldatapath(arfmodel),.Platform\$file.sep,.fmri.data.filename(funcdata),'.Rda',sep='')
save(arfcor,file=fn)

return(arfcor)

}

makeSingleTrialEvents <-
function(subject,condition,sefilename='single_events',hrf.control=list(a1=6,a2=12,b1=0.9,b2=0.9,ce=0.35),experiment=NULL)
#make single trial estimates
{

#open functional file
filelist <- .data.betafiles(arfdata)
path <- .data.funcDir(arfdata)
sp <- .Platform\$file.sep

datavec = numeric(0)
totdim = 0

for(filename in filelist) {
if(file.exists(filename)) {

info <- getFileInfo(filename)
dirname <- .nifti.fileinfo.filename(info)

#get functional information
timings <- .functional.timings(func)
stimlen <- attr(.functional.timings(func),'stimlen')
if(is.null(stimlen)) stimlen = rep(1,length(timings))

#make array of fmri volume
funcvolume = .fmri.data.datavec(fmrivolume)
dim(funcvolume) = c(.fmri.data.dims(fmrivolume)[2],.fmri.data.dims(fmrivolume)[3],.fmri.data.dims(fmrivolume)[4],.fmri.data.dims(fmrivolume)[5])
n = .fmri.data.dims(fmrivolume)[2]*.fmri.data.dims(fmrivolume)[3]*.fmri.data.dims(fmrivolume)[4]

#create timeseries in seconds and create double gamma
tslen = round(.fmri.data.dims(fmrivolume)[5] * .fmri.data.pixdim(fmrivolume)[5])
stick = rep(0,tslen)
hrf <- gamma.fmri(1:tslen,a1=hrf.control\$a1,a2=hrf.control\$a2,b1=hrf.control\$b1,b2=hrf.control\$b2,ce=hrf.control\$ce)
stick[round(timings)]=1
vecnums = which(stick==1)

#define design matrices and outputs
X = matrix(NA,.fmri.data.dims(fmrivolume)[5],length(vecnums))
beta = matrix(0,n,length(vecnums))

#create single-trial columns in design matrix
for(i in 1:length(vecnums)) {
st_protocol = rep(0,tslen)
len = round(stimlen[i])
if(len<=0) len = 1
stvec = vecnums[i]:(vecnums[i]+len-1)
st_protocol[stvec]=1
bold = convol.fmri(hrf,st_protocol)
bold = bold[seq(1,tslen,.fmri.data.pixdim(fmrivolume)[5])]
X[,i] = bold
}

#Solve the LS equation for each voxel (on DEMEANED data)
Xp = solve(t(X)%*%X)%*%t(X)

i=1;
for(z in 1:.fmri.data.dims(fmrivolume)[4]) {
for(y in 1:.fmri.data.dims(fmrivolume)[3]) {
for(x in 1:.fmri.data.dims(fmrivolume)[2]) {
dat = as.vector(funcvolume[x,y,z,])
dat = dat-mean(dat)
beta[i,] = as.vector(Xp%*%dat)
}
i=i+1;
}
}
}

#concatenate datavecs and increase dimensions of volume
datavec = c(datavec,as.vector(beta))
totdim = totdim + length(vecnums)

}
}

#make new file for filetered single events
.fmri.data.filename(fmrivolume) = sefilename
.fmri.data.dims(fmrivolume)[5] = totdim
.fmri.data.pixdim(fmrivolume)[5] = 1
.fmri.data.datavec(fmrivolume) = datavec
.fmri.data.fullpath(fmrivolume) = .data.funcDir(arfdata)
writeData(fmrivolume,.fmri.data.datavec(fmrivolume))

return(fmrivolume)

}

gamma.fmri <-
#double gamma function by Lourens Waldorp
function(t,a1=6,a2=12,b1=0.9,b2=0.9,ce=0.35,...)
{
hrf <- ((t/a1*b1)^a1)*exp(-1*(t - a1*b1)/b1) -
ce*((t/a2*b2)^a2)*exp(-1*(t - a2*b2)/b2);

kernSize <- length(hrf[abs(hrf)>1e-4]);
attr(hrf,"kernSize") <- kernSize;

return(hrf)
}

convol.fmri <-
#convolve fmri timeseries by Lourens Waldorp
function(hrf, prot, ...)
{
if(is.null(attr(hrf,"kernSize"))){
tmp  <- hrf;
hrf  <- prot;
prot <- tmp;
}
window <- 1:length(prot);
win <- length(window);
protExt <- c(rep(0,win),prot);
conv <- numeric(0);
for(i in window){
conv[i] <- sum( protExt[i:(win+i-1)] * hrf[win:1] );
}

convTime <- as.double(conv[window]);
convSample <- rep(NA,length(window));

stimON <- rep(NA,win);
stimON[which(prot[window]!=0)] <- 0;

attr(conv,"time") <- convTime;
attr(conv,"stimON") <- stimON[1:length(window)];
attr(conv,"stim") <- prot;
attr(conv,"hrf")  <- hrf;
class(conv) <- "fmri"

return(conv)
}

correlationTest <-
function(arfcor)
#calculate correlation test of a matrix
{

data = .arfcorrelation.timebyreg(arfcor)

pmat = matrix(0,ncol(data),ncol(data))
cormat = matrix(1,ncol(data),ncol(data))
numcors = (ncol(data)*ncol(data)-ncol(data))/2

.arfcorrelation.num.corr(arfcor) = numcors

i=1;
for(row in 1:(ncol(data)-1)) {
for(col in (row+1):ncol(data)) {
if((row+1)<=ncol(data)) {
ct = cor.test(data[,row],data[,col])
cormat[row,col] = cormat[col,row] = ct\$estimate
pmat[row,col] = pmat[col,row] = ct\$p.value
i = i + 1
}
}
}

pcor = partialCor(cormat,ncol(data))

.arfcorrelation.corr(arfcor) = cormat
.arfcorrelation.corr.pval(arfcor) = pmat
.arfcorrelation.pacorr(arfcor) = pcor\$pcor
.arfcorrelation.pacorr.pval(arfcor) = pcor\$p

return(arfcor)
}

tsCor <-
function(tsmat,regnames=NULL)
#calculate correlation of timeseries matrix
{

R = cor(t(tsmat))
if(!is.null(regnames)) {
rownames(R) = regnames
colnames(R) = regnames
}

return(R)
}

partialCor <-
function(R,n)
#calculate partial correlations
{
Ri=solve(R)
pC=sigPc=matrix(NA,nrow(R),ncol(R))
for(col in 1:ncol(R)) {
for(row in 1:nrow(R)) {
pC[row,col] = (-1*Ri[row,col]) / sqrt((Ri[row,row]) * (Ri[col,col]))
tval = pC[row,col]/sqrt((1-pC[row,col]^2)/(n-2))
sigPc[row,col]=dt(tval,(n-2))

}
}

return(list(pcor=pC,p=sigPc))

}

differenceCor <-
function(c1,c2,n1,n2=n1)
#test if correlations are different
{

diffscore = zval = pval = z1 = z2 = matrix(0,dim(c1)[1],dim(c1)[2])

for(i in 1:dim(c1)[1]) {
for(j in 1:dim(c1)[2]) {

#fischer z transform correlations
zscore1 = log(abs((c1[i,j]+1)/(c1[i,j]-1)))/2
zscore2 = log(abs((c2[i,j]+1)/(c2[i,j]-1)))/2

#reverse the sign if correlations are negative
if(c1[i,j]<0 & zscore1>0) zscore1=zscore1*-1
if(c2[i,j]<0 & zscore2>0) zscore2=zscore2*-1

#calculate the difference between z values
diffz = zscore1-zscore2

#calculate standard error of correlations
SE = sqrt((1/(n1 - 3)) + (1/(n2 - 3)))

#calculate z-value of the difference
zv = diffz/SE

#check if zv=NaN (only when one of the correlations is one)
if(is.nan(zv)) {
z1[i,j]=0
z2[i,j]=0
zval[i,j] = 0
diffscore[i,j] = 0
pval[i,j]=1
} else {
#fill in matrices with zscores/differences/zvaluesdifference and p-values
z1[i,j]=zscore1
z2[i,j]=zscore2
zval[i,j]=zv
diffscore[i,j]=diffz
pval[i,j] = dnorm(zval[i,j])
}
}
}
return(list(z1=z1,z2=z2,dif=diffscore,z=zval,pval=pval))
}

roiConnectivity <-
function(arfmodel,roidata=setIsoContour(arfmodel,95),funcfilename='single_events.nii.gz',type=c('avg','ev'),evmodel=c('spatial','spatiotemporal','eigenvariate'))
{
#match type argument
type = match.arg(type,c('avg','ev'))
evmodel = match.arg(evmodel,c('spatial','spatiotemporal','eigenvariate'))

#make isocontours of model estimates
roidataarray = fmri2array(roidata)
roi = vector('list',.fmri.data.dims(roidata)[5])
for(i in 1:.fmri.data.dims(roidata)[5]) roi[[i]] = which(as.vector(roidataarray[,,,i])>0)

funcvolume = .fmri.data.datavec(funcdata)
dim(funcvolume) = c(.fmri.data.dims(funcdata)[2],.fmri.data.dims(funcdata)[3],.fmri.data.dims(funcdata)[4],.fmri.data.dims(funcdata)[5])

#create time-series for each volume and blob
b =  matrix(NA,.fmri.data.dims(roidata)[5],.fmri.data.dims(funcdata)[5])

#create weigthed average first spatial eigenvector
if(type=='ev') {

if(evmodel=='spatiotemporal') {

eigenvec = vector('list',.fmri.data.dims(roidata)[5])

for(blob in 1:.fmri.data.dims(roidata)[5]) {
timebyvox = matrix(NA,.fmri.data.dims(funcdata)[5],length(roi[[blob]]))

for(tm in 1:.fmri.data.dims(funcdata)[5]) {
timebyvox[tm,] = as.vector(funcvolume[,,,tm])[roi[[blob]]]
}

blobsvd = svd(timebyvox)
eigenvec[[blob]] = (blobsvd\$u[,1]*blobsvd\$d[1])%*%t(blobsvd\$v[,1])

}
}

if(evmodel=='spatial') {

eigenvec = vector('list',.fmri.data.dims(roidata)[5])

for(blob in 1:.fmri.data.dims(roidata)[5]) {
timebyvox = matrix(NA,.fmri.data.dims(funcdata)[5],length(roi[[blob]]))

for(tm in 1:.fmri.data.dims(funcdata)[5]) {
timebyvox[tm,] = as.vector(funcvolume[,,,tm])[roi[[blob]]]
}

#NN = matrix(.fmri.data.datavec(avgdata)[roi[[blob]]],,1)%*%matrix(.fmri.data.datavec(avgdata)[roi[[blob]]],1,)
NN = t(timebyvox)%*%timebyvox
ev1 = eigen(NN)\$vectors[,1]

rev1 = range(ev1)
wrev1 = which.max(abs(rev1))
if(wrev1==1 & rev1[wrev1]<0) ev1=ev1*-1

eigenvec[[blob]] = ev1

}
}

}

#calculate beta time-series estimates
for(blob in 1:.fmri.data.dims(roidata)[5]) {

if(type=='avg') {
for(tm in 1:.fmri.data.dims(funcdata)[5]) b[blob,tm] = mean((as.vector(funcvolume[,,,tm])[roi[[blob]]]))
}

if(type=='ev') {
if(evmodel=='spatiotemporal') {
for(tm in 1:.fmri.data.dims(funcdata)[5]) b[blob,tm] = mean(eigenvec[[blob]][tm,])
}

if(evmodel=='spatial') {
for(tm in 1:.fmri.data.dims(funcdata)[5]) b[blob,tm] = mean((as.vector(funcvolume[,,,tm])[roi[[blob]]])*eigenvec[[blob]])
}

if(evmodel=='eigenvariate') {

timebyvox = matrix(NA,.fmri.data.dims(funcdata)[5],length(roi[[blob]]))

for(tm in 1:.fmri.data.dims(funcdata)[5]) {
timebyvox[tm,] = as.vector(funcvolume[,,,tm])[roi[[blob]]]
}

b[blob,] = svd(timebyvox)\$u[,1]

}

}

}

#get correlations
arfcor <- new('arfcorrelation')
.arfcorrelation.timebyreg(arfcor) <- t(b)
arfcor <- correlationTest(arfcor)

#save correlations
fn = paste(.model.modeldatapath(arfmodel),.Platform\$file.sep,.fmri.data.filename(funcdata),'_',type,'.Rda',sep='')
save(arfcor,file=fn)

return(arfcor)

}

createFuncs <-
function(arfdata)
#createFuncs creates registration files for each
{

#set separator
sp <- .Platform\$file.sep

#make new registration object
functional<- new('functional')

#get betafiles plus path of registration
filelist <- .data.betafiles(arfdata)
path <- .data.funcDir(arfdata)

#check betafile integrity and make paths in regDir
for(filename in filelist) {
if(file.exists(filename)) {

#get info from betafile and set linkedfile
info <- getFileInfo(filename)
dirname <- .nifti.fileinfo.filename(info)

#create dir and create regfilename
if(!file.exists(paste(path,sp,dirname,sep=''))) dir.create(paste(path,sp,dirname,sep=''))

.functional.filename(functional) <- paste(path,sp,dirname,sp,.data.funcRda(arfdata),sep='')

#save objects
save(functional,file=paste(.functional.filename(functional),sep=''))

} else warning('No betafile found to match functional data to')
}
#return(invisble(functional))
}

setFuncFiles <-
function(func_data='filtered_func_data.nii.gz',experiment=NULL)
###
{
#check experiment
if(is.null(experiment)) {
experiment <- try(get('.experiment',envir=.arfInternal),silent=T)
}

#set separator
sp <- .Platform\$file.sep

for(sdirs in 1:.experiment.subject.num(experiment)) {
spath <- paste(.experiment.path(experiment),sp,.experiment.subjectDir(experiment),sp,.experiment.subject.names(experiment)[sdirs],sep='')

for(cdirs in 1:.experiment.condition.num(experiment)) {

createFuncs(dat)

funcdirs <- list.files(cpath)
funcdirs <- funcdirs[which(file.info(list.files(cpath,full.names=T))\$isdir)]

for(trs in 1:length(funcdirs)) {

if(length(func_data)!=length(funcdirs)) func_data = rep(func_data[1],length(funcdirs))

if(length(functional)>0) {
.functional.fullpath(functional) <- paste(spath,sp,.experiment.funcDir(experiment),sp,sep='')
.functional.functionaldata(functional) <- func_data[trs]

#save regfile object
save(functional,file=paste(.functional.filename(functional),sep=''))

}

}

}
}

setFuncTimings <-
function(subject, condition, run, timings, func_data=NULL, experiment=NULL)
#set the timings of a functional volume
{

#set the timings
.functional.timings(functional) = timings

#if func_data is given, also change funcdata
if(!is.null(func_data)) .functional.functionaldata(functional) <- func_data

#save functional.Rda
saveFunc(functional)

#return
return(functional)

}

function(bfslfile)
{

timings = dat[,1]
attr(timings,'stimlen') = dat[,2]

return(timings)
}

function(subject, condition, run, experiment = NULL)
{
#check experiment
if(is.null(experiment)) {
experiment <- try(get('.experiment',envir=.arfInternal),silent=T)
}

#set separator
sp <- .Platform\$file.sep

funcdirs <- list.files(cpath,full.names=T)
funcdirs <- funcdirs[which(file.info(list.files(cpath,full.names=T))\$isdir)]

functional=NULL

#check run

return(functional)

}

saveFunc <-
function(functional)
#save a functional file
{
save(functional,file=.functional.filename(functional))

}

fitConnectivityFromFiles <-
function(filelist,funcfilename='single_events.nii.gz')
#make a connectivity estimate from a series of model files
{
regs = 1:length(filelist)

#make model design matrix
X = matrix(NA,n,length(regs))

#set theta amplitudes to one (all)
p=1
for(i in regs) {
X[,p] = X[,p]/max(abs(X[,p]))
X[,p] = sqrt(X[,p]^2)
p=p+1
}

funcvolume = .fmri.data.datavec(funcdata)
dim(funcvolume) = c(.fmri.data.dims(funcdata)[2],.fmri.data.dims(funcdata)[3],.fmri.data.dims(funcdata)[4],.fmri.data.dims(funcdata)[5])

b =  matrix(NA,length(regs),.fmri.data.dims(funcdata)[5])
Xp = solve(t(X)%*%X)%*%t(X)

p=1
for(tm in 1:.fmri.data.dims(funcdata)[5]) {
y = as.vector(funcvolume[,,,tm])
b[,p] = Xp%*%y
p=p+1
}

#get correlations
arfcor <- new('arfcorrelation')
.arfcorrelation.timebyreg(arfcor) <- t(b)
arfcor <- correlationTest(arfcor)

#save correlations
fn = paste(.fmri.data.filename(funcdata),'.Rda',sep='')
save(arfcor,file=fn)

return(arfcor)

}

makeBetaSeries <-
function(funcdata,timingsfile,method='LS-S',hrf.control=list(a1=6,a2=12,b1=0.9,b2=0.9,ce=0.35)) {
#extract Beta-series regression based on LS-S method of Mumford et al 2012

timings <- tinfo[,1]
stimlen <- tinof[,2]

#make array of fmri volume
funcvolume = .fmri.data.datavec(fmrivolume)
dim(funcvolume) = c(.fmri.data.dims(fmrivolume)[2],.fmri.data.dims(fmrivolume)[3],.fmri.data.dims(fmrivolume)[4],.fmri.data.dims(fmrivolume)[5])
n = .fmri.data.dims(fmrivolume)[2]*.fmri.data.dims(fmrivolume)[3]*.fmri.data.dims(fmrivolume)[4]

#create timeseries in seconds and create double gamma
tslen = round(.fmri.data.dims(fmrivolume)[5] * .fmri.data.pixdim(fmrivolume)[5])
stick = rep(0,tslen)
hrf <- gamma.fmri(1:tslen,a1=hrf.control\$a1,a2=hrf.control\$a2,b1=hrf.control\$b1,b2=hrf.control\$b2,ce=hrf.control\$ce)
stick[round(timings)]=1
vecnums = which(stick==1)

#define design matrices and outputs
X = matrix(NA,.fmri.data.dims(fmrivolume)[5],length(vecnums))
beta = matrix(0,n,length(vecnums))

#create single-trial columns in design matrix
for(i in 1:length(vecnums)) {
st_protocol = rep(0,tslen)
len = round(stimlen[i])
if(len<=0) len = 1
stvec = vecnums[i]:(vecnums[i]+len-1)
st_protocol[stvec]=1
bold = convol.fmri(hrf,st_protocol)
bold = bold[seq(1,tslen,.fmri.data.pixdim(fmrivolume)[5])]
X[,i] = bold
}

#Solve the LS equation for each voxel (on DEMEANED data)
Xp = solve(t(X)%*%X)%*%t(X)

i=1;
for(z in 1:.fmri.data.dims(fmrivolume)[4]) {
for(y in 1:.fmri.data.dims(fmrivolume)[3]) {
for(x in 1:.fmri.data.dims(fmrivolume)[2]) {
dat = as.vector(funcvolume[x,y,z,])
dat = dat-mean(dat)
beta[i,] = as.vector(Xp%*%dat)
}
i=i+1;
}
}
}

#concatenate datavecs and increase dimensions of volume
datavec = c(datavec,as.vector(beta))
totdim = totdim + length(vecnums)

#for each stimulus create two BOLD regressors

#append beta-estimates to functional data

}

extractBetaSeries <-