Nothing
R/sandwichFunctions.R
In arfS4: Activated Region Fitting, fMRI data analysis.
# Sandwich Functions (Derivatives, Residuals, SandwichCalc)
# Wouter D. Weeda, University of Amsterdam
###############################################################################
meanResidualFull <- function(arfmodel) {
## meanResidual calculates the mean residual matrix used in the sandwich estimation
## input is an object of class model
## output is a vector of mean residuals
if(.model.valid(arfmodel)) {
#open data for first trial to get dimensions and trial data
trialdata <- readData(.model.datafiles(arfmodel)[1])
#warn if data is more than 2D
if(.fmri.data.dims(trialdata)[1]>2) warning('residualCalc only works on 2D data, only first two dimensions are used.')
#set dimensions and data
dimx <- .fmri.data.dims(trialdata)[2]
dimy <- .fmri.data.dims(trialdata)[3]
data <- .fmri.data.datavec(trialdata)[1:(dimx*dimy)]
n <- dimx^2*dimy^2
rm(trialdata)
#calculate model based on parameter estimates
model <- .C('gauss',as.double(.model.estimates(arfmodel)),as.integer(.model.regions(arfmodel)*6),as.integer(dimx),as.integer(dimy),as.double(numeric(dimx*dimy)))[[5]]
#caluclate outerproduct of residuals
resids <- .C('outerprod',as.integer(dimx*dimy),as.double(data-model),as.double(vector('numeric',n)))[[3]]
#calculate resids for second and further trials and sum
for(i in 2:.model.trials(arfmodel)) {
data <- .fmri.data.datavec(readData(.model.datafiles(arfmodel)[i]))[1:(dimx*dimy)]
resids <- resids + .C('outerprod',as.integer(dimx*dimy),as.double(data-model),as.double(vector('numeric',n)))[[3]]
}
#divide resids by number of trials squared
resids <- resids/.model.trials(arfmodel)^2
} else warning('No valid model. Residuals not calculated.')
return(invisible(resids))
}
meanResidualFile <- function(arfmodel) {
## meanResidualFile calculates the mean residual matrix used in the sandwich estimation
## input is an object of class model
## output is a vector of mean residuals
if(.model.valid(arfmodel)) {
#open data for first trial to get dimensions and trial data
trialdata <- readData(.model.datafiles(arfmodel)[1])
#warn if data is more than 2D
if(.fmri.data.dims(trialdata)[1]>2) warning('residualCalc only works on 2D data, only first two dimensions are used.')
#set dimensions and data
dimx <- .fmri.data.dims(trialdata)[2]
dimy <- .fmri.data.dims(trialdata)[3]
data <- .fmri.data.datavec(trialdata)[1:(dimx*dimy)]
n <- dimx*dimy
rm(trialdata)
#calculate model based on parameter estimates
model <- .C('gauss',as.double(.model.estimates(arfmodel)),as.integer(.model.regions(arfmodel)*6),as.integer(dimx),as.integer(dimy),as.double(numeric(dimx*dimy)))[[5]]
#caluclate outerproduct of residuals
curwd=getwd()
setwd(paste(.model.fullpath(arfmodel),'/stats',sep=''))
invisible(.C('newResidualFile',as.integer(n)))
invisible(.C('addResidualFile',as.integer(n),as.double(.model.trials(arfmodel)),as.double(data-model)))
file.copy('trialresiduals.arf','residuals.arf',over=T)
#calculate resids for all trials
for(i in 2:.model.trials(arfmodel)) {
data <- .fmri.data.datavec(readData(.model.datafiles(arfmodel)[i]))[1:(dimx*dimy)]
invisible(.C('addResidualFile',as.integer(n),as.double(.model.trials(arfmodel)),as.double(data-model)))
file.copy('trialresiduals.arf','residuals.arf',over=T)
}
file.remove('residuals.arf')
setwd(curwd)
} else warning('No valid model. Residuals not calculated.')
return(TRUE)
}
meanResidualTriangular <- function(arfmodel) {
if(.model.valid(arfmodel)) {
#open data for first trial to get dimensions and trial data
trialdata <- readData(.model.datafiles(arfmodel)[1])
headinf <- readHeader(getFileInfo(.model.datafiles(arfmodel)[1]))
#warn if data is more than 2D
if(.fmri.data.dims(trialdata)[1]>2) warning('residualCalc only works on 2D data, only first two dimensions are used.')
#set dimensions and data
dimx <- .fmri.data.dims(trialdata)[2]
dimy <- .fmri.data.dims(trialdata)[3]
data <- .fmri.data.datavec(trialdata)[1:(dimx*dimy)]
n <- dimx*dimy
rm(trialdata)
#calculate model based on parameter estimates
model <- .C('gauss',as.double(.model.estimates(arfmodel)),as.integer(.model.regions(arfmodel)*6),as.integer(dimx),as.integer(dimy),as.double(numeric(dimx*dimy)))[[5]]
#caluclate outerproduct of residuals
resids <- .C('outerprodtriangular',as.integer(dimx*dimy),as.double(data-model),as.double(vector('numeric',n)),as.double(vector('numeric',(n^2-n)/2)))
residdiag <- resids[[3]]
residtria <- resids[[4]]
#calculate resids for second and further trials and sum
for(i in 2:.model.trials(arfmodel)) {
data <- .fmri.data.datavec(readData(.model.datafiles(arfmodel)[i]))[1:(dimx*dimy)]
resids <- .C('outerprodtriangular',as.integer(dimx*dimy),as.double(data-model),as.double(vector('numeric',n)),as.double(vector('numeric',(n^2-n)/2)))
residdiag <- residdiag + resids[[3]]
residtria <- residtria + resids[[4]]
}
#divide resids by number of trials squared
residdiag <- residdiag/.model.trials(arfmodel)^2
residtria <- residtria/.model.trials(arfmodel)^2
} else warning('No valid model. Residuals not calculated.')
return(invisible(list(diag=residdiag,tria=residtria)))
}
meanResidualDiagonal <- function(arfmodel) {
if(.model.valid(arfmodel)) {
#open data for first trial to get dimensions and trial data
trialdata <- readData(.model.datafiles(arfmodel)[1])
headinf <- readHeader(getFileInfo(.model.datafiles(arfmodel)[1]))
#warn if data is more than 2D
if(.fmri.data.dims(trialdata)[1]>2) warning('residualCalc only works on 2D data, only first two dimensions are used.')
#set dimensions and data
dimx <- .fmri.data.dims(trialdata)[2]
dimy <- .fmri.data.dims(trialdata)[3]
data <- .fmri.data.datavec(trialdata)[1:(dimx*dimy)]
n <- dimx*dimy
rm(trialdata)
#calculate model based on parameter estimates
model <- .C('gauss',as.double(.model.estimates(arfmodel)),as.integer(.model.regions(arfmodel)*6),as.integer(dimx),as.integer(dimy),as.double(numeric(dimx*dimy)))[[5]]
#caluclate outerproduct of residuals
resids <- .C('outerproddiagonal',as.integer(dimx*dimy),as.double(data-model),as.double(vector('numeric',n)))[[3]]
#calculate resids for second and further trials and sum
for(i in 2:.model.trials(arfmodel)) {
data <- .fmri.data.datavec(readData(.model.datafiles(arfmodel)[i]))[1:(dimx*dimy)]
resids <- resids + .C('outerproddiagonal',as.integer(dimx*dimy),as.double(data-model),as.double(vector('numeric',n)))[[3]]
}
#divide resids by number of trials squared
resids <- resids/.model.trials(arfmodel)^2
} else warning('No valid model. Residuals not calculated.')
return(invisible(resids))
}
sandwichFull <- function(arfmodel) {
## sandwich calculates the sandwich variance/covariance estimator
## input is an object of class arfmodel
## output is an object of class arfmodel
if(.model.valid(arfmodel)) {
#calculate residuals, derivatives and the number of parameters
R <- meanResidualFull(arfmodel)
F <- derivatives(arfmodel)
p <- .model.regions(arfmodel)*6
#read in average weights (with dims)
weights <- readData(.model.avgWfile(arfmodel))
if(.fmri.data.dims(weights)[1]>2) warning('sandwich only works on 2D data, only first two dimensions are used.')
n <- (.fmri.data.dims(weights)[2])*(.fmri.data.dims(weights)[3])
W <- .fmri.data.datavec(weights)[1:n]
rm(weights)
#calculate the inner sandwich part B (in A-1BA-1)
B <- matrix(.C('inner_sandwich',as.integer(n),as.integer(p),as.double(F),as.double(W),as.double(R),as.double(rep(0,p*p)))[[6]],p,p)
#calculate the sandwich estimator (using the Hessian returned by nlm)
SW <- try(solve(.5*.model.hessian(arfmodel))%*%B%*%solve(.5*.model.hessian(arfmodel)),silen=T)
#check if alll went well and add to the arfmodel object
if(is.null(attr(SW,'class'))) {
.model.varcov(arfmodel) <- SW
} else {
warning('Failed to calculate Sandwich estimate.')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'no var/covar matrix calculated.\n')
.model.valid(arfmodel) <- FALSE
}
} else {
warning('No valid model. sandwich not calculated')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'no var/covar matrix calculated.\n')
}
return(invisible(arfmodel))
}
sandwichFile <- function(arfmodel) {
## sandwich calculates the sandwich variance/covariance estimator
## input is an object of class arfmodel
## output is an object of class arfmodel
if(.model.valid(arfmodel)) {
#calculate residuals, derivatives and the number of parameters
meanResidualFile(arfmodel)
F <- derivatives(arfmodel)
p <- .model.regions(arfmodel)*6
#read in average weights (with dims)
weights <- readData(.model.avgWfile(arfmodel))
if(.fmri.data.dims(weights)[1]>2) warning('sandwich only works on 2D data, only first two dimensions are used.')
n <- (.fmri.data.dims(weights)[2])*(.fmri.data.dims(weights)[3])
W <- .fmri.data.datavec(weights)[1:n]
rm(weights)
#calculate the inner sandwich part B (in A-1BA-1)
curwd=getwd()
setwd(paste(.model.fullpath(arfmodel),'/stats',sep=''))
B <- matrix(.C('inner_sandwichFile',as.integer(n),as.integer(p),as.double(F),as.double(W),as.double(rep(0,p*p)))[[5]],p,p)
setwd(curwd)
#calculate the sandwich estimator (using the Hessian returned by nlm)
SW <- try(solve(.5*.model.hessian(arfmodel))%*%B%*%solve(.5*.model.hessian(arfmodel)),silen=T)
#check if alll went well and add to the arfmodel object
if(is.null(attr(SW,'class'))) {
.model.varcov(arfmodel) <- SW
} else {
warning('Failed to calculate Sandwich estimate.')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'no var/covar matrix calculated.\n')
.model.valid(arfmodel) <- FALSE
}
} else {
warning('No valid model. sandwich not calculated')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'no var/covar matrix calculated.\n')
}
return(invisible(arfmodel))
}
sandwichTriangular <- function(arfmodel) {
## sandwich calculates the sandwich variance/covariance estimator for large n values
## input is an object of class arfmodel
## output is an object of class arfmodel
if(.model.valid(arfmodel)) {
p <- .model.regions(arfmodel)*6
F <- matrix(derivatives(arfmodel),,p)
SW <- matrix(NA,p,p)
#read in average weights (with dims)
weights <- readData(.model.avgWfile(arfmodel))
if(.fmri.data.dims(weights)[1]>2) warning('sandwich only works on 2D data, only first two dimensions are used.')
n <- (.fmri.data.dims(weights)[2])*(.fmri.data.dims(weights)[3])
W <- .fmri.data.datavec(weights)[1:n]
rm(weights)
leftPart <- solve(.5*.model.hessian(arfmodel))%*%(t(F)*(1/W))
rightPart <- (F*(1/W))%*%solve(.5*.model.hessian(arfmodel))
rT <- meanResidualTriangular(arfmodel)
#check if alll went well and add to the arfmodel object
if(is.null(attr(SW,'class'))) {
.model.varcov(arfmodel) <- SW
} else {
warning('Failed to calculate Sandwich estimate.')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'no var/covar matrix calculated.\n')
.model.valid(arfmodel) <- FALSE
}
} else {
warning('No valid model. sandwich not calculated')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'no var/covar matrix calculated.\n')
}
#NOT WORKING ALWAYS RETURN ERROR
warning('Failed to calculate Sandwich estimate.')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'sandiwchTria not yet implemented!!!!\n')
.model.valid(arfmodel) <- FALSE
return(invisible(arfmodel))
}
sandwichDiagonal <- function(arfmodel) {
## sandwichcalculates the sandwich variance/covariance estimator for large n values
## input is an object of class arfmodel
## output is an object of class arfmodel
if(.model.valid(arfmodel)) {
R <- meanResidualDiagonal(arfmodel)
p <- .model.regions(arfmodel)*6
F <- matrix(derivatives(arfmodel),,p)
#read in average weights (with dims)
weights <- readData(.model.avgWfile(arfmodel))
if(.fmri.data.dims(weights)[1]>2) warning('sandwich only works on 2D data, only first two dimensions are used.')
n <- (.fmri.data.dims(weights)[2])*(.fmri.data.dims(weights)[3])
W <- .fmri.data.datavec(weights)[1:n]
rm(weights)
#calculate the inner sandwich part B (in A-1BA-1)
B <- ((t(F)*(1/W))*R)%*%(F*(1/W))
#calculate the sandwich estimator (using the Hessian returned by nlm)
SW <- try(solve(.5*.model.hessian(arfmodel))%*%B%*%solve(.5*.model.hessian(arfmodel)),silen=T)
#check if alll went well and add to the arfmodel object
if(is.null(attr(SW,'class'))) {
.model.varcov(arfmodel) <- SW
} else {
warning('Failed to calculate Sandwich estimate.')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'no var/covar matrix calculated.\n')
.model.valid(arfmodel) <- FALSE
}
} else {
warning('No valid model. sandwich not calculated')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'no var/covar matrix calculated.\n')
}
return(invisible(arfmodel))
}
derivatives <- function(arfmodel) {
## derivatives calculates the first order derivatives of the parameters to the objective function
## input is an object of class model
## output is a nxp matrix of derivatives
if(.model.valid(arfmodel)) {
#read in average data
fmridata <- readData(.model.avgdatfile(arfmodel))
#warn if data has more than 2 dimensions
if(.fmri.data.dims(fmridata)[1]>2) warning('derivatives only works on 2D data, only first two dimensions are used.')
dimx <- .fmri.data.dims(fmridata)[2]
dimy <- .fmri.data.dims(fmridata)[3]
#calculate the derivatives
deriv <- .C('fderiv',as.integer(dimx*dimy),as.integer(.model.regions(arfmodel)*6),as.integer(dimx),as.integer(dimy),as.double(.model.estimates(arfmodel)),as.double(numeric(.model.regions(arfmodel)*6*dimx*dimy)))[[6]]
} else warning('No valid model. Derivatives not calculated.')
return(invisible(deriv))
}
Try the arfS4 package in your browser
Any scripts or data that you put into this service are public.
arfS4 documentation built on May 2, 2019, 6:14 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# Sandwich Functions (Derivatives, Residuals, SandwichCalc)
# Wouter D. Weeda, University of Amsterdam
###############################################################################
meanResidualFull <- function(arfmodel) {
## meanResidual calculates the mean residual matrix used in the sandwich estimation
## input is an object of class model
## output is a vector of mean residuals
if(.model.valid(arfmodel)) {
#open data for first trial to get dimensions and trial data
trialdata <- readData(.model.datafiles(arfmodel)[1])
#warn if data is more than 2D
if(.fmri.data.dims(trialdata)[1]>2) warning('residualCalc only works on 2D data, only first two dimensions are used.')
#set dimensions and data
dimx <- .fmri.data.dims(trialdata)[2]
dimy <- .fmri.data.dims(trialdata)[3]
data <- .fmri.data.datavec(trialdata)[1:(dimx*dimy)]
n <- dimx^2*dimy^2
rm(trialdata)
#calculate model based on parameter estimates
model <- .C('gauss',as.double(.model.estimates(arfmodel)),as.integer(.model.regions(arfmodel)*6),as.integer(dimx),as.integer(dimy),as.double(numeric(dimx*dimy)))[[5]]
#caluclate outerproduct of residuals
resids <- .C('outerprod',as.integer(dimx*dimy),as.double(data-model),as.double(vector('numeric',n)))[[3]]
#calculate resids for second and further trials and sum
for(i in 2:.model.trials(arfmodel)) {
data <- .fmri.data.datavec(readData(.model.datafiles(arfmodel)[i]))[1:(dimx*dimy)]
resids <- resids + .C('outerprod',as.integer(dimx*dimy),as.double(data-model),as.double(vector('numeric',n)))[[3]]
}
#divide resids by number of trials squared
resids <- resids/.model.trials(arfmodel)^2
} else warning('No valid model. Residuals not calculated.')
return(invisible(resids))
}
meanResidualFile <- function(arfmodel) {
## meanResidualFile calculates the mean residual matrix used in the sandwich estimation
## input is an object of class model
## output is a vector of mean residuals
if(.model.valid(arfmodel)) {
#open data for first trial to get dimensions and trial data
trialdata <- readData(.model.datafiles(arfmodel)[1])
#warn if data is more than 2D
if(.fmri.data.dims(trialdata)[1]>2) warning('residualCalc only works on 2D data, only first two dimensions are used.')
#set dimensions and data
dimx <- .fmri.data.dims(trialdata)[2]
dimy <- .fmri.data.dims(trialdata)[3]
data <- .fmri.data.datavec(trialdata)[1:(dimx*dimy)]
n <- dimx*dimy
rm(trialdata)
#calculate model based on parameter estimates
model <- .C('gauss',as.double(.model.estimates(arfmodel)),as.integer(.model.regions(arfmodel)*6),as.integer(dimx),as.integer(dimy),as.double(numeric(dimx*dimy)))[[5]]
#caluclate outerproduct of residuals
curwd=getwd()
setwd(paste(.model.fullpath(arfmodel),'/stats',sep=''))
invisible(.C('newResidualFile',as.integer(n)))
invisible(.C('addResidualFile',as.integer(n),as.double(.model.trials(arfmodel)),as.double(data-model)))
file.copy('trialresiduals.arf','residuals.arf',over=T)
#calculate resids for all trials
for(i in 2:.model.trials(arfmodel)) {
data <- .fmri.data.datavec(readData(.model.datafiles(arfmodel)[i]))[1:(dimx*dimy)]
invisible(.C('addResidualFile',as.integer(n),as.double(.model.trials(arfmodel)),as.double(data-model)))
file.copy('trialresiduals.arf','residuals.arf',over=T)
}
file.remove('residuals.arf')
setwd(curwd)
} else warning('No valid model. Residuals not calculated.')
return(TRUE)
}
meanResidualTriangular <- function(arfmodel) {
if(.model.valid(arfmodel)) {
#open data for first trial to get dimensions and trial data
trialdata <- readData(.model.datafiles(arfmodel)[1])
headinf <- readHeader(getFileInfo(.model.datafiles(arfmodel)[1]))
#warn if data is more than 2D
if(.fmri.data.dims(trialdata)[1]>2) warning('residualCalc only works on 2D data, only first two dimensions are used.')
#set dimensions and data
dimx <- .fmri.data.dims(trialdata)[2]
dimy <- .fmri.data.dims(trialdata)[3]
data <- .fmri.data.datavec(trialdata)[1:(dimx*dimy)]
n <- dimx*dimy
rm(trialdata)
#calculate model based on parameter estimates
model <- .C('gauss',as.double(.model.estimates(arfmodel)),as.integer(.model.regions(arfmodel)*6),as.integer(dimx),as.integer(dimy),as.double(numeric(dimx*dimy)))[[5]]
#caluclate outerproduct of residuals
resids <- .C('outerprodtriangular',as.integer(dimx*dimy),as.double(data-model),as.double(vector('numeric',n)),as.double(vector('numeric',(n^2-n)/2)))
residdiag <- resids[[3]]
residtria <- resids[[4]]
#calculate resids for second and further trials and sum
for(i in 2:.model.trials(arfmodel)) {
data <- .fmri.data.datavec(readData(.model.datafiles(arfmodel)[i]))[1:(dimx*dimy)]
resids <- .C('outerprodtriangular',as.integer(dimx*dimy),as.double(data-model),as.double(vector('numeric',n)),as.double(vector('numeric',(n^2-n)/2)))
residdiag <- residdiag + resids[[3]]
residtria <- residtria + resids[[4]]
}
#divide resids by number of trials squared
residdiag <- residdiag/.model.trials(arfmodel)^2
residtria <- residtria/.model.trials(arfmodel)^2
} else warning('No valid model. Residuals not calculated.')
return(invisible(list(diag=residdiag,tria=residtria)))
}
meanResidualDiagonal <- function(arfmodel) {
if(.model.valid(arfmodel)) {
#open data for first trial to get dimensions and trial data
trialdata <- readData(.model.datafiles(arfmodel)[1])
headinf <- readHeader(getFileInfo(.model.datafiles(arfmodel)[1]))
#warn if data is more than 2D
if(.fmri.data.dims(trialdata)[1]>2) warning('residualCalc only works on 2D data, only first two dimensions are used.')
#set dimensions and data
dimx <- .fmri.data.dims(trialdata)[2]
dimy <- .fmri.data.dims(trialdata)[3]
data <- .fmri.data.datavec(trialdata)[1:(dimx*dimy)]
n <- dimx*dimy
rm(trialdata)
#calculate model based on parameter estimates
model <- .C('gauss',as.double(.model.estimates(arfmodel)),as.integer(.model.regions(arfmodel)*6),as.integer(dimx),as.integer(dimy),as.double(numeric(dimx*dimy)))[[5]]
#caluclate outerproduct of residuals
resids <- .C('outerproddiagonal',as.integer(dimx*dimy),as.double(data-model),as.double(vector('numeric',n)))[[3]]
#calculate resids for second and further trials and sum
for(i in 2:.model.trials(arfmodel)) {
data <- .fmri.data.datavec(readData(.model.datafiles(arfmodel)[i]))[1:(dimx*dimy)]
resids <- resids + .C('outerproddiagonal',as.integer(dimx*dimy),as.double(data-model),as.double(vector('numeric',n)))[[3]]
}
#divide resids by number of trials squared
resids <- resids/.model.trials(arfmodel)^2
} else warning('No valid model. Residuals not calculated.')
return(invisible(resids))
}
sandwichFull <- function(arfmodel) {
## sandwich calculates the sandwich variance/covariance estimator
## input is an object of class arfmodel
## output is an object of class arfmodel
if(.model.valid(arfmodel)) {
#calculate residuals, derivatives and the number of parameters
R <- meanResidualFull(arfmodel)
F <- derivatives(arfmodel)
p <- .model.regions(arfmodel)*6
#read in average weights (with dims)
weights <- readData(.model.avgWfile(arfmodel))
if(.fmri.data.dims(weights)[1]>2) warning('sandwich only works on 2D data, only first two dimensions are used.')
n <- (.fmri.data.dims(weights)[2])*(.fmri.data.dims(weights)[3])
W <- .fmri.data.datavec(weights)[1:n]
rm(weights)
#calculate the inner sandwich part B (in A-1BA-1)
B <- matrix(.C('inner_sandwich',as.integer(n),as.integer(p),as.double(F),as.double(W),as.double(R),as.double(rep(0,p*p)))[[6]],p,p)
#calculate the sandwich estimator (using the Hessian returned by nlm)
SW <- try(solve(.5*.model.hessian(arfmodel))%*%B%*%solve(.5*.model.hessian(arfmodel)),silen=T)
#check if alll went well and add to the arfmodel object
if(is.null(attr(SW,'class'))) {
.model.varcov(arfmodel) <- SW
} else {
warning('Failed to calculate Sandwich estimate.')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'no var/covar matrix calculated.\n')
.model.valid(arfmodel) <- FALSE
}
} else {
warning('No valid model. sandwich not calculated')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'no var/covar matrix calculated.\n')
}
return(invisible(arfmodel))
}
sandwichFile <- function(arfmodel) {
## sandwich calculates the sandwich variance/covariance estimator
## input is an object of class arfmodel
## output is an object of class arfmodel
if(.model.valid(arfmodel)) {
#calculate residuals, derivatives and the number of parameters
meanResidualFile(arfmodel)
F <- derivatives(arfmodel)
p <- .model.regions(arfmodel)*6
#read in average weights (with dims)
weights <- readData(.model.avgWfile(arfmodel))
if(.fmri.data.dims(weights)[1]>2) warning('sandwich only works on 2D data, only first two dimensions are used.')
n <- (.fmri.data.dims(weights)[2])*(.fmri.data.dims(weights)[3])
W <- .fmri.data.datavec(weights)[1:n]
rm(weights)
#calculate the inner sandwich part B (in A-1BA-1)
curwd=getwd()
setwd(paste(.model.fullpath(arfmodel),'/stats',sep=''))
B <- matrix(.C('inner_sandwichFile',as.integer(n),as.integer(p),as.double(F),as.double(W),as.double(rep(0,p*p)))[[5]],p,p)
setwd(curwd)
#calculate the sandwich estimator (using the Hessian returned by nlm)
SW <- try(solve(.5*.model.hessian(arfmodel))%*%B%*%solve(.5*.model.hessian(arfmodel)),silen=T)
#check if alll went well and add to the arfmodel object
if(is.null(attr(SW,'class'))) {
.model.varcov(arfmodel) <- SW
} else {
warning('Failed to calculate Sandwich estimate.')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'no var/covar matrix calculated.\n')
.model.valid(arfmodel) <- FALSE
}
} else {
warning('No valid model. sandwich not calculated')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'no var/covar matrix calculated.\n')
}
return(invisible(arfmodel))
}
sandwichTriangular <- function(arfmodel) {
## sandwich calculates the sandwich variance/covariance estimator for large n values
## input is an object of class arfmodel
## output is an object of class arfmodel
if(.model.valid(arfmodel)) {
p <- .model.regions(arfmodel)*6
F <- matrix(derivatives(arfmodel),,p)
SW <- matrix(NA,p,p)
#read in average weights (with dims)
weights <- readData(.model.avgWfile(arfmodel))
if(.fmri.data.dims(weights)[1]>2) warning('sandwich only works on 2D data, only first two dimensions are used.')
n <- (.fmri.data.dims(weights)[2])*(.fmri.data.dims(weights)[3])
W <- .fmri.data.datavec(weights)[1:n]
rm(weights)
leftPart <- solve(.5*.model.hessian(arfmodel))%*%(t(F)*(1/W))
rightPart <- (F*(1/W))%*%solve(.5*.model.hessian(arfmodel))
rT <- meanResidualTriangular(arfmodel)
#check if alll went well and add to the arfmodel object
if(is.null(attr(SW,'class'))) {
.model.varcov(arfmodel) <- SW
} else {
warning('Failed to calculate Sandwich estimate.')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'no var/covar matrix calculated.\n')
.model.valid(arfmodel) <- FALSE
}
} else {
warning('No valid model. sandwich not calculated')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'no var/covar matrix calculated.\n')
}
#NOT WORKING ALWAYS RETURN ERROR
warning('Failed to calculate Sandwich estimate.')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'sandiwchTria not yet implemented!!!!\n')
.model.valid(arfmodel) <- FALSE
return(invisible(arfmodel))
}
sandwichDiagonal <- function(arfmodel) {
## sandwichcalculates the sandwich variance/covariance estimator for large n values
## input is an object of class arfmodel
## output is an object of class arfmodel
if(.model.valid(arfmodel)) {
R <- meanResidualDiagonal(arfmodel)
p <- .model.regions(arfmodel)*6
F <- matrix(derivatives(arfmodel),,p)
#read in average weights (with dims)
weights <- readData(.model.avgWfile(arfmodel))
if(.fmri.data.dims(weights)[1]>2) warning('sandwich only works on 2D data, only first two dimensions are used.')
n <- (.fmri.data.dims(weights)[2])*(.fmri.data.dims(weights)[3])
W <- .fmri.data.datavec(weights)[1:n]
rm(weights)
#calculate the inner sandwich part B (in A-1BA-1)
B <- ((t(F)*(1/W))*R)%*%(F*(1/W))
#calculate the sandwich estimator (using the Hessian returned by nlm)
SW <- try(solve(.5*.model.hessian(arfmodel))%*%B%*%solve(.5*.model.hessian(arfmodel)),silen=T)
#check if alll went well and add to the arfmodel object
if(is.null(attr(SW,'class'))) {
.model.varcov(arfmodel) <- SW
} else {
warning('Failed to calculate Sandwich estimate.')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'no var/covar matrix calculated.\n')
.model.valid(arfmodel) <- FALSE
}
} else {
warning('No valid model. sandwich not calculated')
.model.warnings(arfmodel) <- paste(.model.warnings(arfmodel),'no var/covar matrix calculated.\n')
}
return(invisible(arfmodel))
}
derivatives <- function(arfmodel) {
## derivatives calculates the first order derivatives of the parameters to the objective function
## input is an object of class model
## output is a nxp matrix of derivatives
if(.model.valid(arfmodel)) {
#read in average data
fmridata <- readData(.model.avgdatfile(arfmodel))
#warn if data has more than 2 dimensions
if(.fmri.data.dims(fmridata)[1]>2) warning('derivatives only works on 2D data, only first two dimensions are used.')
dimx <- .fmri.data.dims(fmridata)[2]
dimy <- .fmri.data.dims(fmridata)[3]
#calculate the derivatives
deriv <- .C('fderiv',as.integer(dimx*dimy),as.integer(.model.regions(arfmodel)*6),as.integer(dimx),as.integer(dimy),as.double(.model.estimates(arfmodel)),as.double(numeric(.model.regions(arfmodel)*6*dimx*dimy)))[[6]]
} else warning('No valid model. Derivatives not calculated.')
return(invisible(deriv))
}
Try the arfS4 package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.