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

```#############################################
# arf3DS4 S4 VAR/COVAR FUNCTIONS	  		#
# Wouter D. Weeda							#
# University of Amsterdam					#
#############################################

#[CONTAINS]
#makeDerivs			[user]
#makeResiduals
#makeWeightsBin
#varcov				[user]
#BIC				[user]
#RMSEA				[user]
#wald				[user]
#mcpCorrect			[user]
#approxHessianR
#getlocsquare
#makeBWlocations

makeDerivs <-
function(arfmodel,method=c('viaR','direct'))
#make derivs calculates analytical first order derivatives based on the estimated model parameters
{

method=match.arg(method[1],c('viaR','direct'))

if(.model.valid(arfmodel)) {

#set the filename
fn <- paste(.model.modeldatapath(arfmodel),.Platform\$file.sep,.model.derivativeFile(arfmodel),sep='')

#calculate and write the derivatives
if(method=='viaR') {
con = file(fn,'wb')
writeBin(derivs,con,double(),endian=.Platform\$endian)
close(con)
}

if(method=='direct') {
}

return(invisible(TRUE))

} else return(invisible(FALSE))

}

makeResiduals <-
function(arfmodel)
#make residuals makes residuals for each run
{
if(.model.valid(arfmodel)) {
#set separator
sp <- .Platform\$file.sep

#open a binary connection to the residualfile
con <- file(paste(.model.modeldatapath(arfmodel),sp,.model.residualFile(arfmodel),sep=''),'wb')

if(length(nonbrain)>0) model = model[-nonbrain]

#write the residuals per run
for(bfile in .model.betafiles(arfmodel)) {

if(length(nonbrain)>0) data = data[-nonbrain]
writeBin(data-model,con,double(),endian=.Platform\$endian)

}
#close the connection
close(con)

return(invisible(TRUE))

} else return(invisible(FALSE))

}

makeWeightsBin <-
function(arfmodel)
#make binary weightsfile
{

weightdata = .fmri.data.datavec(weights)

if(length(rem)>0) weightdata = weightdata[-rem]

con <- file(paste(.model.modeldatapath(arfmodel),.Platform\$file.sep,.model.weightFile(arfmodel),sep=''),'wb')
writeBin(weightdata,con,double(),endian=.Platform\$endian)
close(con)

}

varcov <-
function(arfmodel)
#calculate sandwich estimate
{
#set separator
sp <- .Platform\$file.sep

if(.model.valid(arfmodel)) {

#check for deriv and resid
if(!file.exists(paste(.model.modeldatapath(arfmodel),sp,.model.residualFile(arfmodel),sep=''))) makeResiduals(arfmodel)
if(!file.exists(paste(.model.modeldatapath(arfmodel),sp,.model.derivativeFile(arfmodel),sep=''))) makeDerivs(arfmodel)
if(!file.exists(paste(.model.modeldatapath(arfmodel),sp,.model.weightFile(arfmodel),sep=''))) makeWeightsBin(arfmodel)

st_time <- Sys.time()

#try to invert hessian
hessian <- try(solve(.model.hessian(arfmodel)),silent=T)

#if Hessian is not ok, try pseudoinverse
if(!is.null(attr(hessian,'class')))  {
#library(corpcor)
hessian <- try(pseudoinverse(.model.hessian(arfmodel)),silent=T)
.model.warnings(arfmodel) <- c(.model.warnings(arfmodel),paste('[varcov] Using Pseudoinverse (hessian probably singular)',sep=''))
pseudo = T
} else { pseudo = F}

#check if hessian is good
if(is.null(attr(hessian,'class')))  {
n = .model.n(arfmodel)

#perform the inner_sandwich procedure
if(.model.sandwichmethod(arfmodel)[1]=='diag') B <- try(.C('innerSWdiag',as.integer(n),as.integer(.model.regions(arfmodel)*.model.params(arfmodel)),as.integer(.model.runs(arfmodel)),as.character(paste(.model.modeldatapath(arfmodel),sp,.model.derivativeFile(arfmodel),sep='')),as.character(paste(.model.modeldatapath(arfmodel),sp,.model.residualFile(arfmodel),sep='')),as.character(paste(.model.modeldatapath(arfmodel),sp,.model.weightFile(arfmodel),sep='')),as.double(numeric((.model.regions(arfmodel)*.model.params(arfmodel))^2)))[[7]],silent=T)
if(.model.sandwichmethod(arfmodel)[1]=='fast') B <- try(.C('innerSWfast',as.integer(n),as.integer(.model.regions(arfmodel)*.model.params(arfmodel)),as.integer(.model.runs(arfmodel)),as.character(paste(.model.modeldatapath(arfmodel),sp,.model.derivativeFile(arfmodel),sep='')),as.character(paste(.model.modeldatapath(arfmodel),sp,.model.residualFile(arfmodel),sep='')),as.character(paste(.model.modeldatapath(arfmodel),sp,.model.weightFile(arfmodel),sep='')),as.double(numeric((.model.regions(arfmodel)*.model.params(arfmodel))^2)))[[7]],silent=T)
if(.model.sandwichmethod(arfmodel)[1]=='full') B <- try(.C('innerSWfull',as.integer(n),as.integer(.model.regions(arfmodel)*.model.params(arfmodel)),as.integer(.model.runs(arfmodel)),as.character(paste(.model.modeldatapath(arfmodel),sp,.model.derivativeFile(arfmodel),sep='')),as.character(paste(.model.modeldatapath(arfmodel),sp,.model.residualFile(arfmodel),sep='')),as.character(paste(.model.modeldatapath(arfmodel),sp,.model.weightFile(arfmodel),sep='')),as.character(paste(.model.modeldatapath(arfmodel),sp,'mean',.model.residualFile(arfmodel),sep='')),as.double(numeric((.model.regions(arfmodel)*.model.params(arfmodel))^2)))[[8]],silent=T)
if(.model.sandwichmethod(arfmodel)[1]=='bw') {
#save original model values
old_n = .model.n(arfmodel)

#set mask and n to full
.model.n(arfmodel) <- length(.fmri.data.datavec(weights))

#remake Resids,Derivs,Weights
makeResiduals(arfmodel)
makeDerivs(arfmodel)
makeWeightsBin(arfmodel)

#make bandwidth and locations matrices
bw = as.numeric(.model.sandwichmethod(arfmodel)[2])
escapevar = .model.n(arfmodel)+1
Lv = makeBWLocations(arfmodel,escapevar)

#run BandWidth SW
B <- try(.C('innerSWbw',as.integer(.model.n(arfmodel)),as.integer(.model.regions(arfmodel)*.model.params(arfmodel)),as.integer(.model.runs(arfmodel)),as.integer(((bw*2)+1)^3),as.integer(escapevar),as.integer(Lv),as.character(paste(.model.modeldatapath(arfmodel),sp,.model.derivativeFile(arfmodel),sep='')),as.character(paste(.model.modeldatapath(arfmodel),sp,.model.residualFile(arfmodel),sep='')),as.character(paste(.model.modeldatapath(arfmodel),sp,.model.weightFile(arfmodel),sep='')),as.character(paste(.model.modeldatapath(arfmodel),sp,'mean',.model.residualFile(arfmodel),sep='')),as.double(numeric((.model.regions(arfmodel)*.model.params(arfmodel))^2)))[[11]],silent=T)

#reset original values
.model.n(arfmodel) = old_n
}

if(.model.sandwichmethod(arfmodel)[1]=='bwf') {
#save original model values
old_n = .model.n(arfmodel)

#set mask and n to full
.model.n(arfmodel) <- length(.fmri.data.datavec(weights))

#remake Resids,Derivs,Weights
makeResiduals(arfmodel)
makeDerivs(arfmodel)
makeWeightsBin(arfmodel)

#make bandwidth and locations matrices
bw = as.numeric(.model.sandwichmethod(arfmodel)[2])
escapevar = .model.n(arfmodel)+1
Lv = makeBWLocations(arfmodel,escapevar)

#run BandWidth SW
B <- try(.C('innerSWbwfast',as.integer(.model.n(arfmodel)),as.integer(.model.regions(arfmodel)*.model.params(arfmodel)),as.integer(.model.runs(arfmodel)),as.integer(((bw*2)+1)^3),as.integer(escapevar),as.integer(Lv),as.character(paste(.model.modeldatapath(arfmodel),sp,.model.derivativeFile(arfmodel),sep='')),as.character(paste(.model.modeldatapath(arfmodel),sp,.model.residualFile(arfmodel),sep='')),as.character(paste(.model.modeldatapath(arfmodel),sp,.model.weightFile(arfmodel),sep='')),as.character(paste(.model.modeldatapath(arfmodel),sp,'mean',.model.residualFile(arfmodel),sep='')),as.double(numeric((.model.regions(arfmodel)*.model.params(arfmodel))^2)))[[11]],silent=T)

#reset original values
.model.n(arfmodel) = old_n
}

#check if innersandwich works
if(is.null(attr(B,'class'))) {
#set B to be a matrix
dim(B) <- c(.model.regions(arfmodel)*.model.params(arfmodel),.model.regions(arfmodel)*.model.params(arfmodel))

if(!pseudo) SW <- try(solve(.5*.model.hessian(arfmodel))%*%B%*%solve(.5*.model.hessian(arfmodel)),silent=T)
if(pseudo) SW <- try(pseudoinverse(.5*.model.hessian(arfmodel))%*%B%*%pseudoinverse(.5*.model.hessian(arfmodel)),silent=T)

#check if outersandwich works
if(is.null(attr(SW,'class'))) {
.model.varcov(arfmodel) <- SW
.model.valid(arfmodel) <- TRUE
} else { #outersandiwch not good
.model.warnings(arfmodel) <- c(.model.warnings(arfmodel),paste('[varcov] outerSandwich did not compute: ',gsub('\n','',SW),sep=''))
.model.valid(arfmodel) <- FALSE
}
} else { #innersandiwch not good
.model.warnings(arfmodel) <- c(.model.warnings(arfmodel),paste('[varcov] innerSandwich did not compute: ',gsub('\n','',B),sep=''))
.model.valid(arfmodel) <- FALSE
}
} else { #hessian not good
.model.warnings(arfmodel) <- c(.model.warnings(arfmodel),paste('[varcov] Hessian is singular: ',gsub('\n','',hessian),sep=''))
.model.valid(arfmodel) <- FALSE
}

en_time <- Sys.time()
.model.proctime(arfmodel)[1,2] <- as.numeric(difftime(en_time,st_time,units='sec'))

} else { #model not valid
.model.warnings(arfmodel) <- c(.model.warnings(arfmodel),'[varcov] No valid model. var/cov not calculated')
.model.valid(arfmodel) <- FALSE
}

#remove derivatives and residuals
#fn <- paste(.model.modeldatapath(arfmodel),.Platform\$file.sep,.model.derivativeFile(arfmodel),sep='')
#if(file.exists(fn)) file.remove(fn)
#fn <- paste(.model.modeldatapath(arfmodel),.Platform\$file.sep,.model.derivativeFile(arfmodel),sep='')
#if(file.exists(fn)) file.remove(fn)
#fn <- paste(.model.modeldatapath(arfmodel),.Platform\$file.sep,.model.residualFile(arfmodel),sep='')
#if(file.exists(fn)) file.remove(fn)
#fn <- paste(.model.modeldatapath(arfmodel),.Platform\$file.sep,'mean',.model.residualFile(arfmodel),sep='')
#if(file.exists(fn)) file.remove(fn)

#save the modelInfo
saveModel(arfmodel)

#return the varcov
return(invisible(arfmodel))
}

BIC <-
#calculate the BIC
{
if(.model.valid(arfmodel)) {

#read in weights, used in constant
n <- .model.n(arfmodel)
W <- .fmri.data.datavec(Wdata)[1:(.fmri.data.dims(Wdata)[2]*.fmri.data.dims(Wdata)[3]*.fmri.data.dims(Wdata)[4])]

#calculate the determinant of the weights
dtm <- prod(W)

#check if determinant is valid
if(!is.na(dtm) & !is.nan(dtm)) {
if(is.numeric(try(log(n))) & is.numeric(try(log(dtm))) & is.numeric(log(.model.minimum(arfmodel)))) {
cons <- try((2*(((n/2)*log(2*pi))+((1/2)*log(dtm))+((1/2)*(.model.minimum(arfmodel))))),silent=T)

if(cons!=0) {
if(log(dtm)==-Inf) {
dtm=1e-323
#.model.warnings(arfmodel) <- c(.model.warnings(arfmodel),'Determinant (in BIC) was -Inf, set to minimum value 1e-323.')
}
if(log(dtm)==Inf) {
dtm=1e308
#.model.warnings(arfmodel) <- c(.model.warnings(arfmodel),'Determinant (in BIC) was Inf, set to maximum value 1e308.')
}
cons <- (2*(((n/2)*log(2*pi))+((1/2)*log(dtm))+((1/2)*(.model.minimum(arfmodel)))))
}

} else { #logs not valid
.model.warnings(arfmodel) <- c(.model.warnings(arfmodel),'Error calculating logs. BIC not calculated')
.model.valid(arfmodel) <- FALSE
}
} else { #determinant of W is not valid
.model.warnings(arfmodel) <- c(.model.warnings(arfmodel),'Invalid determinant. BIC not calculated')
.model.valid(arfmodel) <- FALSE
}

#check if constant is a number and calculate BIC
if(is.numeric(cons)) {
.model.fit(arfmodel)[1,1]  <- cons + (((.model.regions(arfmodel)*.model.params(arfmodel)))*log(n))

} else { #constant is not a number
.model.warnings(arfmodel) <- c(.model.warnings(arfmodel),'Constant invalid. BIC not calculated')
.model.valid(arfmodel) <- FALSE
}

} else	{
.model.warnings(arfmodel) <- c(.model.warnings(arfmodel),'No valid model. BIC not calculated')
.model.valid(arfmodel) <- FALSE
}

#save the modelInfo
saveModel(arfmodel)

return(invisible(arfmodel))

}

RMSEA <-
#RMSEA calculates root mean square errors of models
{

#check model validity
if(.model.valid(arfmodel)) {

#set number of voxels
n <- .model.n(arfmodel)

#Hotellings T
HTs = .model.runs(arfmodel)*.model.minimum(arfmodel)

#noncentrality parameter
ncp = max(c((HTs-((n-(.model.regions(arfmodel)*.model.params(arfmodel)))/.model.runs(arfmodel))),0))

#check if ncp < 0
if(ncp<0) .model.warnings(arfmodel) <- c(.model.warnings(arfmodel),'Noncentrality parameter smaller than zero, ncp is set to zero')

#calculate RMSEA
eps = sqrt(ncp/(n-(.model.regions(arfmodel)*.model.params(arfmodel))))

if(checkVersion(.model.version(arfmodel),1,2,7)) .model.fit(arfmodel)[1,2] = eps else .model.warnings(arfmodel) <- c(.model.warnings(arfmodel),'RMSEA not calculated for models prior to v1.2-7')

} else {
#if not good warn
.model.warnings(arfmodel) <- c(.model.warnings(arfmodel),'No valid model. RMSEA not calculated')
.model.valid(arfmodel) <- FALSE
}

#save the modelInfo
saveModel(arfmodel)

return(invisible(arfmodel))

}

function(theta)
#get determinant of sigma with derivatives
{

#set relevant parameters
theta3=theta[1]
theta4=theta[2]
theta5=theta[3]
theta6=theta[4]
theta7=theta[5]
theta8=theta[6]

#expression for det(sigma)
det_sigma <- expression((theta3^2)*(theta4^2)*(theta5^2)-(theta3^2)*(theta8*theta4*theta5)*(theta8*theta4*theta5)-(theta6*theta4*theta3)*(theta6*theta4*theta3)*(theta5^2)+(theta6*theta4*theta3)*(theta7*theta5*theta3)*(theta8*theta4*theta5)+(theta7*theta5*theta3)*(theta6*theta4*theta3)*(theta8*theta4*theta5)-(theta7*theta5*theta3)*(theta4^2)*(theta7*theta5*theta3))

#return derivative and value

}

wald <-
#calculate Wald statistics
{

if(.model.modeltype(arfmodel)!='gauss') stop('wald statistics can only be calculated for the full model')

if(length(.model.varcov(arfmodel))==0) stop('(co)variance matrix not yet calculated, cannot compute Wald statistics!')

#define function to calculate Wald statistic
W <- function(a,A,C) t(a)%*%solve(A%*%C%*%t(A))%*%a

if(.model.valid(arfmodel)) {

#if no design matrix is specified in the waldobject, make the default matrix (zero-filled)
if(dim(.wald.consts(waldobject))[1]==0) .wald.consts(waldobject) <- matrix(0,.model.regions(arfmodel),5)

# get dimensions (set number of voxels)
n <- .model.n(arfmodel)

#set relevant matrix sizes and dfs
.wald.stats(waldobject) <- matrix(0,.model.regions(arfmodel),5)
.wald.pvalues(waldobject) <- matrix(0,.model.regions(arfmodel),5)
.wald.df1(waldobject) <- rep(n,5)
.wald.df2(waldobject) <- .wald.df1(waldobject)-rep(.model.regions(arfmodel)*.model.params(arfmodel),5)

#perform hypothesis tests for each region and for locations, extent and amplitude
for(region in 1:.model.regions(arfmodel)) {

#select the 10*10 vcov matrix and estimates for each region, with determinant and deriv
theta <- .model.estimates(arfmodel)[((1+(region-1)*.model.params(arfmodel)):(region*.model.params(arfmodel)))]
C <- .model.varcov(arfmodel)[((1+(region-1)*.model.params(arfmodel)):(region*.model.params(arfmodel))),((1+(region-1)*.model.params(arfmodel)):(region*.model.params(arfmodel)))]

#define the a matrix (containing hypotheses), uses info from the designmatrix
a <- c(theta[1]-.wald.consts(waldobject)[region,1],theta[2]-.wald.consts(waldobject)[region,2],theta[3]-.wald.consts(waldobject)[region,3],sigma\$value-.wald.consts(waldobject)[region,4],theta[10]-.wald.consts(waldobject)[region,5])

#define the A matrix (containing the derivatives of a)
A <- matrix(0,5,10)
A[1,1] <- 1
A[2,2] <- 1
A[3,3] <- 1
A[5,10] <- 1

#perform tests for locations (calc stats and p-values)
.wald.stats(waldobject)[region,1] <- W(a[1],A[1,1],C[1,1])
.wald.pvalues(waldobject)[region,1] <- 1-pf(.wald.stats(waldobject)[region,1],.wald.df1(waldobject)[1],.wald.df2(waldobject)[1])
.wald.stats(waldobject)[region,2] <- W(a[2],A[2,2],C[2,2])
.wald.pvalues(waldobject)[region,2] <- 1-pf(.wald.stats(waldobject)[region,2],.wald.df1(waldobject)[2],.wald.df2(waldobject)[2])
.wald.stats(waldobject)[region,3] <- W(a[3],A[3,3],C[3,3])
.wald.pvalues(waldobject)[region,3] <- 1-pf(.wald.stats(waldobject)[region,3],.wald.df1(waldobject)[3],.wald.df2(waldobject)[3])

#perform tests for spatial extent (calc stats and p-values)
.wald.stats(waldobject)[region,4] <- W(a[4],matrix(A[4,(4:9)],1,6),C[(4:9),(4:9)])
.wald.pvalues(waldobject)[region,4] <- 1-pf(.wald.stats(waldobject)[region,4],.wald.df1(waldobject)[4],.wald.df2(waldobject)[4])

#perform tests for amplitude (calc stats and p-values)
.wald.stats(waldobject)[region,5] <- W(a[5],A[5,10],C[10,10])
.wald.pvalues(waldobject)[region,5] <- 1-pf(.wald.stats(waldobject)[region,5],.wald.df1(waldobject)[5],.wald.df2(waldobject)[5])
}

.model.wald(arfmodel) <- waldobject

#save the model file
#saveModel(arfmodel)

} else	warning('No valid model. wald statistics not calculated.')

#save the modelInfo
saveModel(arfmodel)

return(invisible(arfmodel))

}

mcpCorrect <-
#calulate the multiple comparison correction on fmri data (uncorrected,bonferroni or FDR)
{
veclen <- length(.fmri.data.datavec(fmridata))
if(adj.n) n <- length(.fmri.data.datavec(fmridata)[.fmri.data.datavec(fmridata)!=0]) else n <- length(.fmri.data.datavec(fmridata))

pseq = seq(sig.steps,1,-1)

which <- match.arg(type)

if(which=='bonferroni') {
sigvec = numeric(veclen)
for(i in 1:sig.steps) 	sigvec[pvec[i]>(1-pt(abs(.fmri.data.datavec(fmridata)),df))]=1/pseq[i]

}

if(which=='FDR') {
fdr.value = sort(1-pt(.fmri.data.datavec(fmridata),df))
i = 1
while(fdr.value[i]<=((i*q)/(n*cv))) {
i = i + 1
if(i==length(fdr.value)) break()
}

i = max(i-1,i)
fdr=fdr.value[i]
pvec = fdr / seq(1,sig.steps)
sigvec = numeric(veclen)
for(i in 1:sig.steps) 	sigvec[pvec[i]>(1-pt(abs(.fmri.data.datavec(fmridata)),df))]=1/pseq[i]
thres.t = abs(qt(fdr,df))
thres.p = fdr

}

if(which=='uncorrected') {
sigvec = numeric(veclen)
for(i in 1:sig.steps) 	sigvec[pvec[i]>(1-pt(abs(.fmri.data.datavec(fmridata)),df))]=1/pseq[i]

}

.fmri.data.datavec(overlay_fmridata) <- sigvec

}

function(arfmodel)
{
sp = .Platform\$file.sep

fn = paste(.model.modeldatapath(arfmodel),sp,.model.fullmodelDataFile(arfmodel),sep='')

rm(dat)
p = .model.regions(arfmodel)*.model.params(arfmodel)
fn = paste(.model.modeldatapath(arfmodel),sp,.model.derivativeFile(arfmodel),sep='')

if(file.exists(fn)) {
con <- file(fn,open='rb')
dim(dfvec) = c(n,p)
close(con)
return(dfvec)
}
return(NULL)
}

approxHessianR <-
function(arfmodel)
#approxHessian calculates the approximate hessian (used only to check C code)
{

hessian = 2 * (t(df)%*%(df*(1/W)))

return(hessian)

}

getlocsquare <-
#return voxel-vector locations from a box with width bw*2 + 1
{
fw = (bw*2)+1
lv = numeric(fw^3)

xv = (x-bw):(x+bw)
yv = (y-bw):(y+bw)
zv = (z-bw):(z+bw)

i=1;
for(zz in zv) {
for(yy in yv) {
for(xx in xv) {
if(xx<1 | xx>dx | yy<1 | yy>dy | zz<1 | zz>dz) lv[i]=escapevar else lv[i] = (xx + (yy-1)*(dx) + (zz-1)*(dx)*(dy))-1
#if(lv[i]!=escapevar) {
#if(mask[(xx + (yy-1)*(dx) + (zz-1)*(dx)*(dy))]==0) lv[i]=escapevar;
#}
i=i+1
}
}
}

#returns location in c-style (starting at zero)
return(lv)

}

makeBWLocations <-
function(arfmodel,escapevar)
#make Bandwidth location matrix
{

bw = as.numeric(.model.sandwichmethod(arfmodel)[2])

dx = .fmri.data.dims(dat)[2]
dy = .fmri.data.dims(dat)[3]
dz = .fmri.data.dims(dat)[4]

#make and fill locations matrix
i=j=1
for(z in 1:dz) {
for(y in 1:dy) {
for(x in 1:dx) {
i=i+1
}
j=j+1
}
}
}

return(as.vector(lv))
}

```

## Try the arf3DS4 package in your browser

Any scripts or data that you put into this service are public.

arf3DS4 documentation built on May 2, 2019, 5:16 p.m.