R/iModel.R
In Tokazama/iClass: More intuitive syntax for manipulating imaging data
#' Class iModel
#'
#' Object containing a fitted statistical model for image data.
#'
#' @slot iData iData reference class object for fitted model.
#' @slot Y Pointer to image matrix stored as h5 dataset.
#' @slot X Information about the design:
#' \itemize{
#' \item{X} {Design matrix.}
#' \item{K} {Filter information specific to each image.}
#' \item{W} {weights}
#' \item{KWX} {Weighted design matrix information.}
#' \item{pKX} {Pseudoinverse of X.}
#' \item{V} {Non-sphericity matrix.}
#' \item{betaCov} {Covariance of correlation coefficients.}
#' \item{trRV} {Trace of RV.}
#' \item{trRVRV} {Trace of RVRV.}
#' \item{rdf} {Residual degrees of freedom.}
#' }
#' @slot beta Pointer to beta matrix stored as h5 dataset.
#' @slot res Pointer to residual matrix stored as h5 dataset.
#' @slot mrss Pointer to mean residual sum of squares matrix stored as h5 dataset.
#' @slot xVi Information about intrinsic temporal non-sphericity:
#' \itemize{
#' \item{Vi} {List of non-sphericity components.}
#' \item{V} {Non-sphericity matrix (Cov(e) = sigma^2*V)}
#' \item{h} {Hyperparameters.}
#' \item{Cy} {Covariance of response matrix.}
#' }
#' @slot dims Model dimensions:
#' \itemize{
#' \item{npred} {Number of predictors.}
#' \item{nimg} {Number of images.}
#' \item{nvox} {Number of voxels.}
#' \item{fwhm} {full-width at half-maxima}
#' \item{resels} {Resolution elements.}
#' \item{rpvImage} {Resels per voxel image.}
#' }
#' @slot C list of contrasts:
#' \itemize{
#' \item{name} {Name of the contrast.}
#' \item{c} {Contrast weights.}
#' \item{X1} {Reamaining design space (orthogonal to X0)}
#' \item{X0} {Reduced design matrix.}
#' \item{iX0} {Indicates how contrat was specified.}
#' \item{dims} {Dimensions for contrast}
#' \itemize{
#' \item{trRV} {Trace of RV.}
#' \item{trRVRV} {Trace of RVRV.}
#' \item{idf} {Degrees of interest.}
#' }
#' \item{fieldType} {Type of statyistical field being fitted.}
#' \item{contrastImage} {Object of class antsImage representing contrast.}
#' \item{clusterImage} {Object of class antsImage representing the thresholded contrastImage.}
#' \item{results} {Either a list or data.frame containing the results of a given contrast.}
#' \item{sthresh} {Statistical threshold.}
#' \item{cthresh} {Cluster threshold.}
#' }
#' @slot control Control values for fitting the model (see \code{\link{iControl}}).
#'
#' @author Zachary P. Christensen
#' @seealso \code{\link{ilm}}, \code{\link{iglm}}
#'
#' @rdname iModel-class
iModel <- setClass("iModel",
slots = list(
file = "H5File",
beta = "DataSet",
res = "DataSet",
mrss = "DataSet",
C = "list",
location = "character"),
contains = "iDesign")
#' @export
#' @rdname iModel-class
ilModel <- setClass("ilModel", contains = "iModel")
#' @export
#' @rdname iModel-class
iglModel <- setClass("iglModel", contains = "iModel")
#' @export
setMethod("show", "iModel", function(object) {
cat("iModel: \n\n")
cat(" Call = ", object@method[1], "\n")
if (length(object@method) > 1)
cat(" Optimization = ", object@method[2], "\n")
cat(" Predictors = ", paste(colnames(object@X$X), collapse = ", "), "\n")
cat(" Images = ", object@dims$nimg, "\n")
cat(" Residual degrees of freedom = ", object@X$rdf, "\n")
cat(" Voxels = ", object@dims$nvox, "\n")
cat(" Location = ", object@location, "\n")
if (!is.null(object@dims$resels)) {
cat(" FWHM = ", round(object@dims$fwhm, 2), "\n")
cat(" Resels = ", round(object@dims$resels), "\n")
}
cat("--- \n\n")
ncon <- length(object@C)
if (ncon == 0)
cat("Contrasts not set.\n")
else {
for (i in seq_len(ncon)) {
cat("Contrast", object@C[[i]]$name, "\n")
cat(" Contrast weights: \n")
c <- t(object@C[[i]]$c)
colnames(c) <- colnames(object@X$X)
print(c)
cat("\n")
if (object@control$rft) {
if (class(object@C[[i]]$results) == "character") {
cat(object@C[[i]]$results, "\n\n")
} else {
cat(" Set-level: \n")
cat(" Clusters = ", ncol(object@C[[i]]$results$clusterLevel), "\n")
cat(" p-value = ", object@C[[i]]$results$setLevel, "\n\n")
cat(" Cluster-Level: \n")
print(round(object@C[[i]]$results$clusterLevel, 3))
cat("\n")
cat(" Peak-level: \n")
print(round(object@C[[i]]$results$peakLevel, 3))
cat("\n")
}
} else {
if (class(object@C[[i]]$results) == "character") {
cat(object@C[[i]]$results, "\n\n")
} else
print(object@C[[i]]$results)
}
cat("Interest degrees of freedom = ", object@C[[i]]$dims$idf, "\n")
cat(" Statistical threshold = ", round(object@C[[i]]$sthresh, 2), "\n")
cat(" Cluster threshold = ", object@C[[i]]$cthresh, "\n")
cat(" Threshold type = ", object@C[[i]]$threshType, "\n")
cat("---\n\n")
}
}
})
#' Fit Images to Linear Models
#'
#' @param formula Object of class \code{formula} (or one that can be coerced to that class): a symbolic description of the model to be fitted. The details of model specification are given under ‘Details’.
#' @param iData Object of class \code{\link{iData}} containing data represented in the provided formula.
#' @param weights An optional vector of ‘prior weights’ to be used in the fitting process. Should be NULL or a numeric vector.
#' @param na.action If \code{"na.omit"} then all images with NA values in
#' corresponding variables are omitted. If a function is specified
#' @param control A list of parameters for controlling the fitting process. See \code{\link{iControl}} for details.
#' @param verbose Enables verbose output. (default = \code{TRUE}).
#'
#'
#' @export ilm
ilm <- function(formula, iData, weights, na.action = "na.omit", control = list(), verbose = TRUE) {
resp <- iFormula(formula, iData, subset, weights, na.action, offset, control)
if (class(resp) == "list") {
out <- list()
for (i in seq_len(length(resp))) {
out[[i]] <- iModelMake(res[[i]], call = "ilm")
out[[i]] <- iModelSolve(out[[i]])
}
} else {
out <- iModelMake(iResp, filename, call = "ilm")
out <- iModelSolve(out)
if (out@control$rft) {
if (verbose)
cat("Estimating FWHM/Resels. \n")
tmp <- iModelRFT(object, verbose)
out@dims$fwhm <- out$fwhm
out@dims$rpvImage <- out$rpvImage
out@dims$resels <- out$resels
}
}
if (class(out) == "list")
class(out) <- "multilModel"
return(out)
}
#' Fit Images to Generalized Linear Models
#'
#' @param formula Object of class \code{formula} (or one that can be coerced to that class): a symbolic description of the model to be fitted. The details of model specification are given under ‘Details’.
#' @param iData Object of class \code{\link{iData}} containing data represented in the provided formula.
#' @param weights An optional vector of ‘prior weights’ to be used in the fitting process. Should be NULL or a numeric vector.
#' @param na.action If \code{"na.omit"} then all images with NA values in
#' corresponding variables are omitted. If a function is specified
#' @param method Optimization method "none", "IWLS", or "REML" (default = \code{"none"}).
#' @param its Number of iterations for optimization to undergo (if method = "IWLS" defaualt = \code{200}, if method = "REML", default = \code{32}).
#' @param control A list of parameters for controlling the fitting process. See \code{\link{iControl}} for details.
#' @param verbose Enables verbose output. (default = \code{TRUE}).
#'
#' @export iglm
iglm <- function(formula, iData, subset = NULL, weights = NULL, na.action = "na.omit", method, its, control = list(), verbose = TRUE) {
resp <- iFormula(formula, iData, subset, weights, na.action, control)
if (class(resp) == "list") {
out <- list()
for (i in seq_len(length(resp))) {
if (method == "REML") {
if (missing(its))
its <- 32
out[[i]] <- .iwls_helper(resp[[i]], its, verbose)
} else {
out[[i]] <- iModelMake(resp[[i]], call = c("iglm", "REML"))
out[[i]]@xVi <- iModelOptim(out[[i]], its)
out[[i]] <- iModelUpdate(out[[i]])
out[[i]] <- iModelSolve(out[[i]])
}
if (out@control$rft) {
if (verbose)
cat("Estimating FWHM/Resels. \n")
tmp <- iModelRFT(object, verbose)
out[[i]]@dims$fwhm <- tmp$fwhm
out[[i]]@dims$rpvImage <- tmp$rpvImage
out[[i]]@dims$resels <- tmp$resels
}
}
} else {
if (method == "IWLS")
out <- .iwls_helper(resp, verbose)
else {
out <- iModelMake(resp, call = c("iglm", "REML"))
if (missing(its))
its <- 32
out@xVi <- iModelOptim(out, its)
out <- iModelUpdate(out)
out <- iModelSolve(out)
}
if (out@control$rft) {
if (verbose)
cat("Estimating FWHM/Resels. \n")
tmp <- iModelRFT(object, verbose)
out@dims$fwhm <- out$fwhm
out@dims$rpvImage <- out$rpvImage
out@dims$resels <- out$resels
}
}
if (class(out) == "list")
class(out) <- "multiglModel"
return(out)
}
.iwls_helper <- function(iResp, its, verbose) {
H <- diag(out@X$X %*% tcrossprod(MASS::ginv(crossprod(out@X$X), out@X$X)))
ores <- 1
nres <- 10
n <- 0
W <- matrix(1, out@dims$nimg, out@dims$nvox)
while(max(abs(ores - nres)) > (1e-4)) {
ores <- nres
n <- n + 1
if (n > its) {
warning("ilm could not converge. Maximal number of iterations exceeded.");
break
}
nres <- 0
for (j in seq_len(out@dims$nvox)) {
if (i == 1)
W[, j][is.na(out@Y[, j])] <- 0
out@beta[, j] <- MASS::ginv(crossprod(out@X$X, diag(W[, j])) %*% out@X$X) %*% crossprod(out@X$X, diag(W[, j])) %*% out@Y[, j]
out@res[, j] <- out@Y[, j] - out@X$X %*% out@beta[, j]
out@mrss[1, j] <- sum(out@res[, j]^2) / out@X$trRV
restmp <- out@res[, j]
mad <- mean(abs(restmp - mean(restmp)))
restmp <- restmp / mad
# need to figure this one out
restmp <- restmp * H
restmp <- abs(restmp) - control$os
restmp[restmp < 0] <- 0
wtmp <- (abs(restmp) < 1) * ((1 - restmp^2)^2)
wtmp[is.na(out@Y[, j])] <- 0
wtmp[out@Y[, j] == 0] <- 0
W[, j] <- wtmp
nres <- nres + sum(restmp[!is.na(restmp)]^2)
out@res[, j] <- out@Y[, j] - out@X$X %*% out@beta[, j]
}
if (verbose)
cat(paste("Iterative reweighting finished after ", n, " iterations.", sep = ""))
}
if (out@control$rft) {
if (verbose)
cat("Estimating FWHM/Resels. \n")
smooth <- estSmooth(out@res[], out@mask, out@X$rdf, scaleResid = FALSE, sample = out@control$sar, verbose = verbose)
out@dims$fwhm <- smooth$fwhm
out@dims$rpvImage <- smooth[[2]]
out@dims$resels <- resels(out@mask, smooth$fwhm)
}
}
#' @export
#' @docType methods
#' @details \strong{iModelMake} Make iModel object.
#' @rdname iModel-method
iModelMake <- function(iDesign, call) {
if (!usePkg("h5"))
stop("Please install package h5 in order to use this function.")
if (call[1] == "ilm")
out <- ilModel()
else if (call[1] == "iglm")
out <- iglModel()
out@Y <- iDesign@Y
out@X <- iDesign@X
out@xVi <- iDesign@xVi
out@mask <- antsImageClone(iDesign@mask)
out@dims <- iDesign@dims
out@method <- call
out@control <- iDesign@control
if (!usePkg("h5"))
stop("Please install package h5 in order to use this function.")
out@location <- tempfile(fileext = ".h5")
out@file <- h5file(out@location)
chunksize <- out@Y@chunksize[2]
nchunks <- floor(out@dims$nvox / chunksize)
chunkseq <- seq_len(chunksize)
for (i in seq_len(nchunks)) {
if (i == 1) {
out@file["beta"] <- data.matrix(matrix(0, out@dims$npred, chunksize))
beta <- out@file["beta"]
out@file["res"] <- data.matrix(matrix(0, out@dims$npred, chunksize))
res <- out@file["res"]
out@file["mrss"] <- data.matrix(matrix(0, out@dims$npred, chunksize))
mrss <- out@file["mrss"]
} else {
beta <- cbind(beta, data.matrix(matrix(0, out@dims$npred, chunksize)))
res <- cbind(res, data.matrix(matrix(0, out@dims$npred, chunksize)))
mrss <- cbind(mrss, data.matrix(matrix(0, out@dims$npred, chunksize)))
}
}
chunkvox <- chunksize * nchunks
if (out@dims$nvox > chunkvox) {
chunksize <- out@dims$nvox - chunkvox
beta <- cbind(beta, data.matrix(matrix(0, out@dims$npred, chunksize)))
res <- cbind(res, data.matrix(matrix(0, out@dims$npred, chunksize)))
mrss <- cbind(mrss, data.matrix(matrix(0, out@dims$npred, chunksize)))
}
out@beta <- beta
out@res <- res
out@mrss <- mrss
return(out)
}
#' @export
#' @docType methods
#' @details \strong{iModelSolve} Solve already initialized iModel for
#' coefficients, residuals, and mean residual sum of squares.
#' @rdname iModel-methods
iModelSolve <- function(object, verbose = TRUE) {
chunksize <- object@Y@chunksize[2]
nchunks <- floor(object@dims$nvox / chunksize)
vrange <- seq_len(chunksize)
if (verbose)
progress <- txtProgressBar(min = 0, max = nchunks, style = 3)
for (j in seq_len(nchunks)) {
if (vrange[chunksize] > object@dims$nvox)
vrange <- vrange[1]:object@dims$nvox
KWY <- .filter(object@X$K, object@X$W %*% object@Y[, vrange])
object@beta[, vrange] <- object@X$pKX %*% KWY
object@res[, vrange] <- .res(object@X$KWX, KWY)
object@mrss[, vrange] <- colSums(object@res[, vrange]^2) / object@X$trRV
if (object@control$scr)
object@res[, vrange] <- t(t(object@res[, vrange]) * (1 / as.numeric(object@mrss[, vrange])))
vrange <- vrange + chunksize
if (verbose)
setTxtProgressBar(progress, j)
}
if (verbose)
close(progress)
return(object)
}
#' @export
#' @docType methods
#' @details \strong{iModelUpdate} Primarily used to update X slot after
#' optimization.
#' @rdname iModel-methods
iModelUpdate <- function(object, ...) {
ll <- c(...)
if (!is.null(ll$weights)) {
if (class(weights) == "numeric" && length(weights) == dims$nimg)
object@X$W <- diag(weights)
else if (all(dim(weights) == dims$nimg))
object@X$W <- weights
else
stop("weights must be a matrix of nimg x nimg or a vector of length nimg")
}
# set weighted design matrix
object@X$KWX <- .setx(.filter(object@X$K, object@X$W %*% object@X$X))
# pseudoinverse of X
object@X$pKX <- .pinvx(object@X$KWX)
object@X$V <- .filter(object@X$K, .filter(object@X$K, object@X$W %*% tcrossprod(object@xVi$V, object@X$W)))
object@X$betaCov <- object@X$pKX %*% tcrossprod(object@X$V, object@X$pKX)
out <- .trRV(object@X$KWX, object@X$V)
object@X$trRV <- out$trRV
object@X$trRVRV <- out$trRVRV
object@X$rdf <- out$rdf
return(object)
}
#' @export
#' @docType methods
#' @details \strong{iModelRFT} Solve for images in resel space.
#' @rdname iModel-methods
iModelRFT <- function(object, verbose = TRUE) {
dims <- dim(object@mask)
D <- length(dims)
dimx <- seq_len(dims[1])
dimx1 <- dimx + 1
if (D > 1) {
dimy <- seq_len(dims[2])
dimy1 <- dimy + 1
}
if (D > 2) {
dimz <- seq_len(dims[3])
dimz1 <- dimz + 1
}
sar <- object@control$sar
# check sample
if (sar > nrow(object@Y) | sar == 0)
sar <- seq_len(nrow(object@Y))
else {
sar <- sample(seq_len(nrow(object@Y)), sar)
scale <- (nrow(object@Y) / object@X$rdf) * (1 / length(sar))
}
# set up for loop----
if (D == 1) {
d1 <- m1 <- matrix(0, dims[1] + 1)
maskar <- as.numeric(object@mask)
m1[dimx1] <- maskar
m3 <- ((m1[dimx1] * m1[dimx]))
Vxx <- matrix(0, dims[1])
} else if (D == 2) {
d1 <- m1 <- matrix(0, dims[1] + 1, dims[2] + 1)
maskar <- as.matrix(object@mask)
m1[dimx1, dimy1, dimz1] <- maskar
m3 <- ((m1[dimx1, dimy1] * m1[dimx, dimy1])) *
((m1[dimx1, dimy1] * m1[dimx1, dimy]))
Vxx <- Vyy <- Vxy <- matrix(0, dims[1], dims[2])
} else if (D == 3) {
d1 <- m1 <- array(0, dim = dims + 1)
maskar <- as.array(object@mask)
m1[dimx1, dimy1, dimz1] <- maskar
m3 <- ((m1[dimx1, dimy1, dimz1] * m1[dimx, dimy1, dimz1])) *
((m1[dimx1, dimy1, dimz1] * m1[dimx1, dimy, dimz1])) *
((m1[dimx1, dimy1, dimz1] * m1[dimx1, dimy1, dimz])) # mask to eliminate voxels without neighbors
Vxx <- Vyy <- Vzz <- Vxy <- Vxz <- Vyz <- array(0, dim = dims)
}
# partial derivatives of each image----
if (verbose)
progress <- txtProgressBar(min = 0, max = length(sar), style = 3)
for (i in sar) {
if (D == 1)
d1[dimx1] <- makeImage(object@mask, object@res[i,] / mrss[])[dimx]
else if (D == 2)
d1[dimx1, dimy1] <- makeImage(object@mask, object@res[i,] / mrss[])[dimx, dimy]
else if (D == 3)
d1[dimx1, dimy1, dimz1] <- makeImage(object@mask, object@res[i,] / mrss[])[dimx, dimy, dimz]
if (D == 1) {
dx <- (d1[dimx1] - d1[dimx]) * m3
} else if (D == 2) {
dx <- (d1[dimx1, dimy1] - d1[dimx, dimy1]) * m3
dy <- (d1[dimx1, dimy1] - d1[dimx1, dimy]) * m3
} else if (D == 3) {
dx <- (d1[dimx1, dimy1, dimz1] - d1[dimx, dimy1, dimz1]) * m3
dy <- (d1[dimx1, dimy1, dimz1] - d1[dimx1, dimy, dimz1]) * m3
dz <- (d1[dimx1, dimy1, dimz1] - d1[dimx1, dimy1, dimz]) * m3
}
Vxx <- Vxx + (dx * dx)
if (D > 1) {
Vyy <- Vyy + (dy * dy)
Vxy <- Vxy + (dx * dy)
}
if (D > 2) {
Vzz <- Vzz + (dz * dz)
Vxz <- Vxz + (dx * dz)
Vyz <- Vyz + (dy * dz)
}
if (verbose)
setTxtProgressBar(progress, i)
}
if (verbose)
close(progress)
# scale variances/covariances----
Vxx <- Vxx * scale
if (D > 1) {
Vyy <- Vyy * scale
Vxy <- Vxy * scale
}
if (D > 2) {
Vzz <- Vzz * scale
Vxz <- Vxz * scale
Vyz <- Vyz * scale
}
if (D == 1) {
xyz <- Vxx * m3
} else if (D == 2) {
xyz <- cbind(matrix(Vxx * m3, ncol = 1), matrix(Vyy * m3, ncol = 1))
rpv <- (Vxx * Vyy ) + (Vxy * 2) # this needs to be checked
} else if (D == 3) {
xyz <- cbind(Vxx * m3, Vyy * m3, Vzz * m3)
rpv <- (Vxx * Vyy * Vzz) +
(Vxy * Vyz * Vxz * 2) -
(Vyz * Vyz * Vxx) -
(Vxy * Vxy * Vzz) -
(Vxz * Vxz * Vyy)
}
out <- list()
# make RPV Image----
rpv[rpv < 0] <- 0
rpv <- sqrt(rpv / (4 * log(2))^D)
out$rpvImage <- as.antsImage(rpv * maskar)
# estimate fwhm----
xyz <- sqrt((xyz) / (4 * log(2)))
nvox <- sum(m3)
rpv <- sum(rpv) / nvox
xyz <- colSums(xyz) / nvox
resels <- rpv^(1 / D) * (xyz / prod(xyz)^(1 / D))
out$fwhm <- 1 / resels
out$resels <- resels(object@mask, out$fwhm)
return(out)
}
#' @export
#' @docType methods
#' @details \strong{getImages} Returns antsImages of specified contrasts within model iModel.
#' @rdname iModel-methods
getImages <- function(object, contrast, rpv = FALSE, mask = FALSE) {
out <- list()
for (i in seq_len(length(contrast))) {
out[[i]]$clusterImage <- antsImageClone(object@C[[contrast]]$clusterImage)
out[[i]]$contrastImage <- antsImageClone(object@C[[contrast]]$clusterImage)
}
names(out) <- contrast
if (rpv)
out$rpvImage <- antsImageClone(object@dims$rpvImage)
if (mask)
out$mask <- antsImageClone(object@iData@iList[[object@y]]@mask)
return(out)
}
#' iModel Methods
#'
#' V
#'
#' @param object Object of class iModel.
#' @param filename h5 file to save iModel to.
#' @param iData_dirname Directory for iData component of iModel.
#' @param dirname Directory to place \code{report} output.
#' @param contrast Name of contrast.
#' @param rpv Return resels per voxel image?
#' @param mask Return mask used for iModel?
#' @param verbose Enables verbose output. (default = \code{TRUE}).
#' @param ... Additional named arguments passed to \code{iModelUpdate}.
#'
#' @author Zachary P. Christensen
#'
#' @name iModel-methods
NULL
#' @export
#' @docType methods
#' @details \strong{coef} Retrieve fitted coefficients from iModel object.
#' @rdname iModel-methods
coef.iModel <- function(object) {
object@beta[]
}
#' @export
#' @docType methods
#' @details \strong{fitted} Retrieve fitted values from iModel object.
#' @rdname iModel-methods
fitted.iModel <- function(object) {
object@X$X %*% object@beta[]
}
#' @export
#' @docType methods
#' @details \strong{resid} Retrieve residuals from iModel object.
#' @rdname iModel-methods
resid.iModel <- function(object) {
object@res[]
}
#' @export
#' @docType methods
#' @details \strong{iModelRead} Read/load iModel object.
#' @rdname iModel-methods
iModelRead <- function(filename) {
if (!file.exists(filename))
stop("filename does not exist.")
x <- iModel()
file <- h5file(filename)
# load big matrices
x@mrss <- file["mrss"]
x@beta <- file["beta"]
x@res <- file["res"]
x@Y <- file["Y"]
x@location <- filename
x@y <- file["y"][]
x@method <- file["method"][]
# read control
x@control$cf <- file["control/cf"][]
x@control$scr <- file["control/scr"][]
x@control$sar <- file["control/sar"][]
x@control$n <- file["control/n"][]
x@control$iso <- file["control/iso"][]
x@control$os <- file["control/os"][]
x@control$rft <- file["control/rft"][]
# read dims
x@dims$npred <- file["dims/npred"][]
x@dims$nvox <- file["dims/nvox"][]
x@dims$nimg <- file["dims/nimg"][]
if (x@control$rft) {
x$dims$fwhm <- file["dims/fwhm"][]
x@dims$resels <- file["dims/resels"][]
x@dims$rpvImage <- as.antsImage(file["dims/rpvImage"][])
k = antsSetSpacing(x@dims$rpvImage, h5attr(file["dims/rpvImage"] , "spacing"))
k = antsSetDirection(x@dims$rpvImage, h5attr(file["dims/rpvImage"] , "direction"))
k = antsSetOrigin(x@dims$rpvImage, h5attr(file["dims/rpvImage"] , "origin"))
}
# read X
x@X$X <- file["X/X"][]
colnames(x@X$X) <- h5attr(file["X/X"], "colnames")[]
x@X$W <- file["X/W"][]
x@X$pKX <- file["X/pKX"][]
x@X$V <- file["X/V"][]
x@X$betaCov <- file["X/betaCov"][]
x@X$trRV <- file["X/trRV"][]
x@X$trRVRV <- file["X/trRVRV"][]
x@X$rdf <- file["X/rdf"][]
x@X$KWX$X <- file["X/KWX/X"][]
x@X$KWX$v <- file["X/KWX/v"][]
x@X$KWX$u <- file["X/KWX/u"][]
x@X$KWX$d <- file["X/KWX/d"][]
x@X$KWX$tol <- file["X/KWX/tol"][]
x@X$KWX$rk <- file["X/KWX/rk"][]
K <- data.frame(Filters = x@file[file.path("K", "Filters")][],
HParam = x@file[file.path("K", "HParam")][],
RT = x@file[file.path("K", "RT")][])
x@X$K <- .filter(K)
# read xVi
if (file["xVi/Vi"][] != "null")
x@xVi$Vi <- file["xVi/Vi"][]
if (file["xVi/V"][] != "null")
x@xVi$V <- file["xVi/V"][]
if (file["xVi/h"][] != "null")
x@xVi$h <- file["xVi/h"][]
if (file["xVi/Cy"][] != "null")
x@xVi$Cy <- file["xVi/Cy"][]
# read contrasts
ncon <- file["C/count"][]
for (i in seq_len(x@C)) {
cname <- paste("C/C", i, sep = "")
x@C[[i]]$name <- file[file.path(cname, "name")][]
x@C[[i]]$c <- file[file.path(cname, "c")][]
x@C[[i]]$X1 <- file[file.path(cname, "X1")][]
x@C[[i]]$X0 <- file[file.path(cname, "X0")][]
x@C[[i]]$iX0file[file.path(cname, "iX0")][]
x@C[[i]]$fieldType <- file[file.path(cname, "fieldType")][]
x@C[[i]]$results <- file[file.path(cname, "results")][]
x@C[[i]]$sthresh <- file[file.path(cname, "sthresh")][]
x@C[[i]]$cthresh <- file[file.path(cname, "cthresh")][]
x@C[[i]]$contrastImage <- as.antsImage(file[file.path(cname, "contrastImage")][])
x@C[[i]]$contrastImage <- antsSetSpacingh5attr(file[file.path(cname, "contrastImage")], "spacing")
x@C[[i]]$contrastImage <- antsSetDirection(h5attr(file[file.path(cname, "contrastImage")], "direction"))
x@C[[i]]$contrastImage <- antsSetOrigin(h5attr(file[file.path(cname, "contrastImage")], "origin"))
x@C[[i]]$clusterImage <- as.antsImage(file[file.path(cname, "clusterImage")][])
x@C[[i]]$clusterImage <- antsSetSpacing(h5attr(file[file.path(cname, "clusterImage")], "spacing"))
x@C[[i]]$clusterImage <- antsSetDirection(h5attr(file[file.path(cname, "clusterImage")], "direction"))
x@C[[i]]$clusterImage <- antsSetOrigin(h5attr(file[file.path(cname, "clusterImage")], "origin"))
x@C[[i]]$dims$trMV <- file[file.path(cname, "dims", "trMV")][]
x@C[[i]]$dims$trMVMV <- file[file.path(cname, "dims", "trMVMV")][]
x@C[[i]]$dims$idf <- file[file.path(cname, "dims", "idf")][]
}
h5close(file)
return(x)
}
#' @export
#' @docType methods
#' @details \strong{iModelWrite} Read/load iModel object.
#' @rdname iModel-methods
# @describeIn iModel write/save iModel object
iModelWrite <- function(object, filename) {
if (class(object) != "iModel")
stop("object must be of class iModel.")
if (file.exists(filename)) {
if (object@location != filename)
stop("filename already exists.")
}
if (!usePkg("h5"))
stop("Please install package h5 in order to use this function.")
file <- h5file(filename)
#### write big matrices
chunksize <- object@beta@chunksize[2]
chunkseq <- seq_len(chunksize)
nvox <- object@Y@dim[2]
nchunk <- floor(nvox / chunksize)
for (i in seq_len(nchunk)) {
if (i == 1) {
# initialize files
file["beta"] <- object@beta[, chunkseq]
file["res"] <- object@res[, chunkseq]
file["mrss"] <- object@mrss[, chunkseq]
file["Y"] <- object@Y[, chunkseq]
beta <- file["beta"]
res <- file["res"]
mrss <- file["mrss"]
Y <- file["Y"]
} else {
beta <- cbind(beta, object@beta[, chunkseq])
res <- cbind(res, object@res[, chunkseq])
mrss <- cbind(mrss, object@mrss[, chunkseq])
Y <- cbind(Y, object@Y[, chunkseq])
}
chunkseq <- chunkseq + chunksize
}
# fill in any remaining chunkage
if (nvox >= chunkseq[1]) {
chunkseq <- chunkseq[1]:nvox
beta <- cbind(beta, object@beta[, chunkseq])
res <- cbind(res, object@res[, chunkseq])
mrss <- cbind(mrss, object@mrss[, chunkseq])
Y <- cbind(Y, object@Y[, chunkseq])
}
#### end big matrices
file["y"] <- object@y
file["method"] <- object@method
# write control
file["control/cf"] <- object@control$cf
file["control/scr"] <- object@control$scr
file["control/sar"] <- object@control$sar
file["control/n"] <- object@control$n
file["control/iso"] <- object@control$iso
file["control/os"] <- object@control$os
file["control/rft"] <- object@control$rft
# write dims
file["dims/npred"] <- object@dims$npred
file["dims/nvox"] <- object@dims$nvox
file["dims/nimg"] <- object@dims$nimg
if (object@control$rft) {
file["dims/fwhm"] <- object@dims$fwhm
file["dims/resels"] <- object@dims$resels
file["dims/rpvImage"] <- as.array(object@dims$rpvImage)
h5attr(file["dims/rpvImage"] , "spacing") <- antsGetSpacing(object@dims$rpvImage)
h5attr(file["dims/rpvImage"] , "direction") <- antsGetDirection(object@dims$rpvImage)
h5attr(file["dims/rpvImage"] , "origin") <- antsGetOrigin(object@dims$rpvImage)
}
# write X
file[file.path("K", "Filters")] <- object@iData@iList[[object@y]]@K$Filters
file[file.path("K", "HParam")] <- object@iData@iList[[object@y]]@K$HParam
file[file.path("K", "RT")] <- object@iData@iList[[object@y]]@K$RT
file["X/X"] <- data.matrix(object@X$X)
h5attr(file["X/X"], "colnames") <- colnames(object@X$X)
file["X/W"] <- data.matrix(object@X$W)
file["X/pKX"] <- data.matrix(object@X$pKX)
file["X/V"] <- data.matrix(object@X$V)
file["X/betaCov"] <- data.matrix(object@X$betaCov)
file["X/trRV"] <- object@X$trRV
file["X/trRVRV"] <- object@X$trRVRV
file["X/rdf"] <- object@X$rdf
file["X/KWX/X"] <- object@X$KWX$X
file["X/KWX/v"] <- object@X$KWX$v
file["X/KWX/u"] <- object@X$KWX$u
file["X/KWX/d"] <- object@X$KWX$d
file["X/KWX/tol"] <- object@X$KWX$tol
file["X/KWX/rk"] <- object@X$KWX$rk
# write xVi
file["xVi/Vi"] <- if (is.null(object@xVi$Vi)) "null" else object@xVi$Vi
file["xVi/V"] <- if (is.null(object@xVi$V)) "null" else object@xVi$V
file["xVi/h"] <- if (is.null(object@xVi$h)) "null" else object@xVi$h
file["xVi/Cy"] <- if (is.null(object@xVi$Cy)) "null" else object@xVi$Cy
# write contrasts
file["C/count"] <- length(object@C)
for (i in seq_len(object@C)) {
cname <- paste("C/C", i, sep = "")
file[file.path(cname, "name")] <- object@C[[i]]$name
file[file.path(cname, "c")] <- object@C[[i]]$c
file[file.path(cname, "X1")] <- object@C[[i]]$X1
file[file.path(cname, "X0")] <- object@C[[i]]$X0
file[file.path(cname, "iX0")] <- object@C[[i]]$iX0
file[file.path(cname, "fieldType")] <- object@C[[i]]$fieldType
file[file.path(cname, "results")] <- object@C[[i]]$results
file[file.path(cname, "sthresh")] <- object@C[[i]]$sthresh
file[file.path(cname, "cthresh")] <- object@C[[i]]$cthresh
file[file.path(cname, "contrastImage")] <- as.array(object@C[[i]]$contrastImage)
h5attr(file[file.path(cname, "contrastImage")], "spacing") <- antsGetSpacing(object@C[[i]]$contrastImage)
h5attr(file[file.path(cname, "contrastImage")], "direction") <- antsGetDirection(object@C[[i]]$contrastImage)
h5attr(file[file.path(cname, "contrastImage")], "origin") <- antsGetOrigin(object@C[[i]]$contrastImage)
file[file.path(cname, "clusterImage")] <- as.array(object@C[[i]]$clusterImage)
h5attr(file[file.path(cname, "clusterImage")], "spacing") <- antsGetSpacing(object@C[[i]]$clusterImage)
h5attr(file[file.path(cname, "clusterImage")], "direction") <- antsGetDirection(object@C[[i]]$clusterImage)
h5attr(file[file.path(cname, "clusterImage")], "origin") <- antsGetOrigin(object@C[[i]]$clusterImage)
file[file.path(cname, "dims", "trMV")] <- object@C[[i]]$dims$trMV
file[file.path(cname, "dims", "trMVMV")] <- object@C[[i]]$dims$trMVMV
file[file.path(cname, "dims", "idf")] <- object@C[[i]]$dims$idf
}
h5close(file)
return(TRUE)
}
#' @export
#' @docType methods
#' @details \strong{iModelSolve} Solve already initialized iModel for
#' coefficients, residuals, and mean residual sum of squares.
#' @rdname iModel-methods
iModelSolve <- function(object, verbose = TRUE) {
chunksize <- object@Y@chunksize[2]
nchunks <- floor(object@dims$nvox / chunksize)
vrange <- seq_len(chunksize)
if (verbose)
progress <- txtProgressBar(min = 0, max = nchunks, style = 3)
for (j in seq_len(nchunks)) {
if (vrange[chunksize] > object@dims$nvox)
vrange <- vrange[1]:object@dims$nvox
KWY <- .filter(object@X$K, object@X$W %*% object@iData@iList[[object@y]]@iMatrix[, vrange])
object@beta[, vrange] <- object@X$pKX %*% KWY
object@res[, vrange] <- .res(object@X$KWX, KWY)
object@mrss[, vrange] <- colSums(object@res[, vrange]^2) / object@X$trRV
if (object@control$scr)
object@res[, vrange] <- t(t(object@res[, vrange]) * (1 / as.numeric(object@mrss[, vrange])))
vrange <- vrange + chunksize
if (verbose)
setTxtProgressBar(progress, j)
}
if (verbose)
close(progress)
return(object)
}
#' @export
#' @docType methods
#' @details \strong{iModelUpdate} Primarily used to update X slot after optimization
#' @rdname iModel-methods
iModelUpdate <- function(object, ...) {
ll <- c(...)
if (!is.null(ll$weights)) {
if (class(weights) == "numeric" && length(weights) == dims$nimg)
object@X$W <- diag(weights)
else if (all(dim(weights) == dims$nimg))
object@X$W <- weights
else
stop("weights must be a matrix of nimg x nimg or a vector of length nimg")
}
# set weighted design matrix
object@X$KWX <- .setx(.filter(object@X$K, object@X$W %*% object@X$X))
# pseudoinverse of X
object@X$pKX <- .pinvx(object@X$KWX)
object@X$V <- .filter(object@X$K, .filter(object@X$K, object@X$W %*% tcrossprod(object@xVi$V, object@X$W)))
object@X$betaCov <- object@X$pKX %*% tcrossprod(object@X$V, object@X$pKX)
out <- .trRV(object@X$KWX, object@X$V)
object@X$trRV <- out$trRV
object@X$trRVRV <- out$trRVRV
object@X$rdf <- out$rdf
return(object)
}
#' @export
#' @docType methods
#' @details \strong{model.matrix} Retrieve design matrix from iModel object.
#' @rdname iModel-methods
model.matrix.iModel <- function(object) {
return(object@X$X)
}
#' @export
#' @docType methods
#' @details \strong{report} Creates a .pdf and .html report of for iModel
#' objects (requires package rmarkdown).
#' @rdname iModel-methods
report <- function(object, dirname, surfimg) {
if (!usePkg("rmarkdown"))
stop("Please install package rmarkdown in order to use this function.")
ncon <- length(object@C)
if (missing(surfimg))
surfimg <- object@mask
dirname <- normalizePath(dirname)
if (!file.exists(dirname))
dir.create(dirname)
## render images----
for (i in seq_len(ncon)) {
if (class(object@C[[i]]$results) != "character") {
brain <- renderSurfaceFunction(surfimg = list(surfimg), funcimg = list(object@C[[i]]$clusterImage), alphasurf = 0.1, smoothsval = 1.5)
tmp <- make3ViewPNG(fnprefix = paste(dirname, "/", object@C[[i]]$name, sep = ""))
}
}
md <- file.path(dirname, "out.Rmd")
zz <- file(md, open = "wt")
sink(zz)
sink(zz, type = "message")
# describe the model----
cat("## iModel object fit by call to ", object@method[1])
if (object@control$rft)
cat(" using random field theory. \n")
else
cat(". \n\n")
cat(" Predictors = ", paste(colnames(object@X$X), collapse = ", "), "\n\n")
cat(" Images = ", object@dims$nimg, "\n\n")
cat(" Degrees of freedom = ", object@X$rdf, "\n\n")
cat(" Voxels = ", object@dims$nvox, "\n\n")
cat(" Optimization = ", object@control$opt, "\n\n")
if (object@control$rft && length(object@dims$fwhm > 0)) {
cat(" FWHM = ", round(object@dims$fwhm, 2), "\n\n")
cat(" Resels = ", round(object@dims$resels), "\n\n")
}
cat("---- \n\n")
# contrast results----
cat("## Contrast Results \n\n")
for (i in seq_len(ncon)) {
cat("### ", object@C[[i]]$name, "\n\n")
if (class(object@C[[i]]$results) != "character")
cat(paste("![", object@C[[i]]$name, "](", file.path(dirname, paste(object@C[[i]]$name, ".png", sep = "")), ") \n\n", sep = ""))
## render results----
cat("Contrast weights: \n\n")
c <- t(object@C[[i]]$c)
colnames(c) <- colnames(object@X$X)
print(c)
cat("\n\n")
if (object@control$rft) {
if (class(object@C[[i]]$results) == "character") {
cat(object@C[[i]]$results, "\n\n")
} else {
cat("#### Set-level \n\n")
cat(" Clusters = ", ncol(object@C[[i]]$results$clusterLevel), "\n\n")
cat(" p-value = ", object@C[[i]]$results$setLevel, "\n\n")
cat("#### Cluster-Level: \n\n")
tmp <- knitr::kable(round(object@C[[i]]$results$clusterLevel, 3))
print(tmp)
cat("\n\n")
cat("#### Peak-level: \n\n")
tmp <- knitr::kable(round(object@C[[i]]$results$peakLevel, 3))
print(tmp)
cat("\n\n")
}
} else {
if (class(object@C[[i]]$results) == "character") {
cat(object@C[[i]]$results, "\n\n")
} else {
knitr::kable(object@C[[i]]$results)
}
}
cat("Interest degrees of freedom = ", object@C[[i]]$dims$idf, "\n\n")
cat(" Statistical threshold = ", round(object@C[[i]]$sthresh, 2), "\n\n")
cat(" Cluster threshold = ", object@C[[i]]$cthresh, "\n\n")
cat(" Threshold type = ", object@C[[i]]$threshType, "\n\n")
cat("---- \n\n")
}
sink(type = "message")
sink()
rmarkdown::render(md)
rmarkdown::render(md, pdf_document())
return(TRUE)
}
#' Control parameters for RFT based analyses
#'
#' Auxillary function for controlling \code{\link{iModel}} fitting.
#'
#' @param cf Critical F-threshold for selecting voxels over which the non-sphericity is estimated (default = \code{0.001}).
#' @param scr Logical. scale residuals? (default = \code{TRUE}).
#' @param sar Number of residual images to sample for estimating the FWHM (default = \code{64}).
#' @param n images in conjunction (default = \code{1}).
#' @param iso logical. should images be assumed to be isotropic? (default = \code{TRUE}).
#' @param os offset weighting for iteratively reweighted least squares (default = \code{3}).
#' @param rft logical. should voxels be estimated in resel space for random field theory analysis (default = \code{TRUE}).
#' @return
#'
#' @export iControl
iControl <- function(cf = 0.05, scr = TRUE, sar = 64, n = 1, iso = TRUE, os = 3, rft = TRUE) {
list(cf = cf,
scr = scr,
sar = sar,
n = n,
iso = iso,
os = os,
rft = rft)
}
Tokazama/iClass documentation built on May 9, 2019, 4:51 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.
#' Class iModel
#'
#' Object containing a fitted statistical model for image data.
#'
#' @slot iData iData reference class object for fitted model.
#' @slot Y Pointer to image matrix stored as h5 dataset.
#' @slot X Information about the design:
#' \itemize{
#' \item{X} {Design matrix.}
#' \item{K} {Filter information specific to each image.}
#' \item{W} {weights}
#' \item{KWX} {Weighted design matrix information.}
#' \item{pKX} {Pseudoinverse of X.}
#' \item{V} {Non-sphericity matrix.}
#' \item{betaCov} {Covariance of correlation coefficients.}
#' \item{trRV} {Trace of RV.}
#' \item{trRVRV} {Trace of RVRV.}
#' \item{rdf} {Residual degrees of freedom.}
#' }
#' @slot beta Pointer to beta matrix stored as h5 dataset.
#' @slot res Pointer to residual matrix stored as h5 dataset.
#' @slot mrss Pointer to mean residual sum of squares matrix stored as h5 dataset.
#' @slot xVi Information about intrinsic temporal non-sphericity:
#' \itemize{
#' \item{Vi} {List of non-sphericity components.}
#' \item{V} {Non-sphericity matrix (Cov(e) = sigma^2*V)}
#' \item{h} {Hyperparameters.}
#' \item{Cy} {Covariance of response matrix.}
#' }
#' @slot dims Model dimensions:
#' \itemize{
#' \item{npred} {Number of predictors.}
#' \item{nimg} {Number of images.}
#' \item{nvox} {Number of voxels.}
#' \item{fwhm} {full-width at half-maxima}
#' \item{resels} {Resolution elements.}
#' \item{rpvImage} {Resels per voxel image.}
#' }
#' @slot C list of contrasts:
#' \itemize{
#' \item{name} {Name of the contrast.}
#' \item{c} {Contrast weights.}
#' \item{X1} {Reamaining design space (orthogonal to X0)}
#' \item{X0} {Reduced design matrix.}
#' \item{iX0} {Indicates how contrat was specified.}
#' \item{dims} {Dimensions for contrast}
#' \itemize{
#' \item{trRV} {Trace of RV.}
#' \item{trRVRV} {Trace of RVRV.}
#' \item{idf} {Degrees of interest.}
#' }
#' \item{fieldType} {Type of statyistical field being fitted.}
#' \item{contrastImage} {Object of class antsImage representing contrast.}
#' \item{clusterImage} {Object of class antsImage representing the thresholded contrastImage.}
#' \item{results} {Either a list or data.frame containing the results of a given contrast.}
#' \item{sthresh} {Statistical threshold.}
#' \item{cthresh} {Cluster threshold.}
#' }
#' @slot control Control values for fitting the model (see \code{\link{iControl}}).
#'
#' @author Zachary P. Christensen
#' @seealso \code{\link{ilm}}, \code{\link{iglm}}
#'
#' @rdname iModel-class
iModel <- setClass("iModel",
slots = list(
file = "H5File",
beta = "DataSet",
res = "DataSet",
mrss = "DataSet",
C = "list",
location = "character"),
contains = "iDesign")
#' @export
#' @rdname iModel-class
ilModel <- setClass("ilModel", contains = "iModel")
#' @export
#' @rdname iModel-class
iglModel <- setClass("iglModel", contains = "iModel")
#' @export
setMethod("show", "iModel", function(object) {
cat("iModel: \n\n")
cat(" Call = ", object@method[1], "\n")
if (length(object@method) > 1)
cat(" Optimization = ", object@method[2], "\n")
cat(" Predictors = ", paste(colnames(object@X$X), collapse = ", "), "\n")
cat(" Images = ", object@dims$nimg, "\n")
cat(" Residual degrees of freedom = ", object@X$rdf, "\n")
cat(" Voxels = ", object@dims$nvox, "\n")
cat(" Location = ", object@location, "\n")
if (!is.null(object@dims$resels)) {
cat(" FWHM = ", round(object@dims$fwhm, 2), "\n")
cat(" Resels = ", round(object@dims$resels), "\n")
}
cat("--- \n\n")
ncon <- length(object@C)
if (ncon == 0)
cat("Contrasts not set.\n")
else {
for (i in seq_len(ncon)) {
cat("Contrast", object@C[[i]]$name, "\n")
cat(" Contrast weights: \n")
c <- t(object@C[[i]]$c)
colnames(c) <- colnames(object@X$X)
print(c)
cat("\n")
if (object@control$rft) {
if (class(object@C[[i]]$results) == "character") {
cat(object@C[[i]]$results, "\n\n")
} else {
cat(" Set-level: \n")
cat(" Clusters = ", ncol(object@C[[i]]$results$clusterLevel), "\n")
cat(" p-value = ", object@C[[i]]$results$setLevel, "\n\n")
cat(" Cluster-Level: \n")
print(round(object@C[[i]]$results$clusterLevel, 3))
cat("\n")
cat(" Peak-level: \n")
print(round(object@C[[i]]$results$peakLevel, 3))
cat("\n")
}
} else {
if (class(object@C[[i]]$results) == "character") {
cat(object@C[[i]]$results, "\n\n")
} else
print(object@C[[i]]$results)
}
cat("Interest degrees of freedom = ", object@C[[i]]$dims$idf, "\n")
cat(" Statistical threshold = ", round(object@C[[i]]$sthresh, 2), "\n")
cat(" Cluster threshold = ", object@C[[i]]$cthresh, "\n")
cat(" Threshold type = ", object@C[[i]]$threshType, "\n")
cat("---\n\n")
}
}
})
#' Fit Images to Linear Models
#'
#' @param formula Object of class \code{formula} (or one that can be coerced to that class): a symbolic description of the model to be fitted. The details of model specification are given under ‘Details’.
#' @param iData Object of class \code{\link{iData}} containing data represented in the provided formula.
#' @param weights An optional vector of ‘prior weights’ to be used in the fitting process. Should be NULL or a numeric vector.
#' @param na.action If \code{"na.omit"} then all images with NA values in
#' corresponding variables are omitted. If a function is specified
#' @param control A list of parameters for controlling the fitting process. See \code{\link{iControl}} for details.
#' @param verbose Enables verbose output. (default = \code{TRUE}).
#'
#'
#' @export ilm
ilm <- function(formula, iData, weights, na.action = "na.omit", control = list(), verbose = TRUE) {
resp <- iFormula(formula, iData, subset, weights, na.action, offset, control)
if (class(resp) == "list") {
out <- list()
for (i in seq_len(length(resp))) {
out[[i]] <- iModelMake(res[[i]], call = "ilm")
out[[i]] <- iModelSolve(out[[i]])
}
} else {
out <- iModelMake(iResp, filename, call = "ilm")
out <- iModelSolve(out)
if (out@control$rft) {
if (verbose)
cat("Estimating FWHM/Resels. \n")
tmp <- iModelRFT(object, verbose)
out@dims$fwhm <- out$fwhm
out@dims$rpvImage <- out$rpvImage
out@dims$resels <- out$resels
}
}
if (class(out) == "list")
class(out) <- "multilModel"
return(out)
}
#' Fit Images to Generalized Linear Models
#'
#' @param formula Object of class \code{formula} (or one that can be coerced to that class): a symbolic description of the model to be fitted. The details of model specification are given under ‘Details’.
#' @param iData Object of class \code{\link{iData}} containing data represented in the provided formula.
#' @param weights An optional vector of ‘prior weights’ to be used in the fitting process. Should be NULL or a numeric vector.
#' @param na.action If \code{"na.omit"} then all images with NA values in
#' corresponding variables are omitted. If a function is specified
#' @param method Optimization method "none", "IWLS", or "REML" (default = \code{"none"}).
#' @param its Number of iterations for optimization to undergo (if method = "IWLS" defaualt = \code{200}, if method = "REML", default = \code{32}).
#' @param control A list of parameters for controlling the fitting process. See \code{\link{iControl}} for details.
#' @param verbose Enables verbose output. (default = \code{TRUE}).
#'
#' @export iglm
iglm <- function(formula, iData, subset = NULL, weights = NULL, na.action = "na.omit", method, its, control = list(), verbose = TRUE) {
resp <- iFormula(formula, iData, subset, weights, na.action, control)
if (class(resp) == "list") {
out <- list()
for (i in seq_len(length(resp))) {
if (method == "REML") {
if (missing(its))
its <- 32
out[[i]] <- .iwls_helper(resp[[i]], its, verbose)
} else {
out[[i]] <- iModelMake(resp[[i]], call = c("iglm", "REML"))
out[[i]]@xVi <- iModelOptim(out[[i]], its)
out[[i]] <- iModelUpdate(out[[i]])
out[[i]] <- iModelSolve(out[[i]])
}
if (out@control$rft) {
if (verbose)
cat("Estimating FWHM/Resels. \n")
tmp <- iModelRFT(object, verbose)
out[[i]]@dims$fwhm <- tmp$fwhm
out[[i]]@dims$rpvImage <- tmp$rpvImage
out[[i]]@dims$resels <- tmp$resels
}
}
} else {
if (method == "IWLS")
out <- .iwls_helper(resp, verbose)
else {
out <- iModelMake(resp, call = c("iglm", "REML"))
if (missing(its))
its <- 32
out@xVi <- iModelOptim(out, its)
out <- iModelUpdate(out)
out <- iModelSolve(out)
}
if (out@control$rft) {
if (verbose)
cat("Estimating FWHM/Resels. \n")
tmp <- iModelRFT(object, verbose)
out@dims$fwhm <- out$fwhm
out@dims$rpvImage <- out$rpvImage
out@dims$resels <- out$resels
}
}
if (class(out) == "list")
class(out) <- "multiglModel"
return(out)
}
.iwls_helper <- function(iResp, its, verbose) {
H <- diag(out@X$X %*% tcrossprod(MASS::ginv(crossprod(out@X$X), out@X$X)))
ores <- 1
nres <- 10
n <- 0
W <- matrix(1, out@dims$nimg, out@dims$nvox)
while(max(abs(ores - nres)) > (1e-4)) {
ores <- nres
n <- n + 1
if (n > its) {
warning("ilm could not converge. Maximal number of iterations exceeded.");
break
}
nres <- 0
for (j in seq_len(out@dims$nvox)) {
if (i == 1)
W[, j][is.na(out@Y[, j])] <- 0
out@beta[, j] <- MASS::ginv(crossprod(out@X$X, diag(W[, j])) %*% out@X$X) %*% crossprod(out@X$X, diag(W[, j])) %*% out@Y[, j]
out@res[, j] <- out@Y[, j] - out@X$X %*% out@beta[, j]
out@mrss[1, j] <- sum(out@res[, j]^2) / out@X$trRV
restmp <- out@res[, j]
mad <- mean(abs(restmp - mean(restmp)))
restmp <- restmp / mad
# need to figure this one out
restmp <- restmp * H
restmp <- abs(restmp) - control$os
restmp[restmp < 0] <- 0
wtmp <- (abs(restmp) < 1) * ((1 - restmp^2)^2)
wtmp[is.na(out@Y[, j])] <- 0
wtmp[out@Y[, j] == 0] <- 0
W[, j] <- wtmp
nres <- nres + sum(restmp[!is.na(restmp)]^2)
out@res[, j] <- out@Y[, j] - out@X$X %*% out@beta[, j]
}
if (verbose)
cat(paste("Iterative reweighting finished after ", n, " iterations.", sep = ""))
}
if (out@control$rft) {
if (verbose)
cat("Estimating FWHM/Resels. \n")
smooth <- estSmooth(out@res[], out@mask, out@X$rdf, scaleResid = FALSE, sample = out@control$sar, verbose = verbose)
out@dims$fwhm <- smooth$fwhm
out@dims$rpvImage <- smooth[[2]]
out@dims$resels <- resels(out@mask, smooth$fwhm)
}
}
#' @export
#' @docType methods
#' @details \strong{iModelMake} Make iModel object.
#' @rdname iModel-method
iModelMake <- function(iDesign, call) {
if (!usePkg("h5"))
stop("Please install package h5 in order to use this function.")
if (call[1] == "ilm")
out <- ilModel()
else if (call[1] == "iglm")
out <- iglModel()
out@Y <- iDesign@Y
out@X <- iDesign@X
out@xVi <- iDesign@xVi
out@mask <- antsImageClone(iDesign@mask)
out@dims <- iDesign@dims
out@method <- call
out@control <- iDesign@control
if (!usePkg("h5"))
stop("Please install package h5 in order to use this function.")
out@location <- tempfile(fileext = ".h5")
out@file <- h5file(out@location)
chunksize <- out@Y@chunksize[2]
nchunks <- floor(out@dims$nvox / chunksize)
chunkseq <- seq_len(chunksize)
for (i in seq_len(nchunks)) {
if (i == 1) {
out@file["beta"] <- data.matrix(matrix(0, out@dims$npred, chunksize))
beta <- out@file["beta"]
out@file["res"] <- data.matrix(matrix(0, out@dims$npred, chunksize))
res <- out@file["res"]
out@file["mrss"] <- data.matrix(matrix(0, out@dims$npred, chunksize))
mrss <- out@file["mrss"]
} else {
beta <- cbind(beta, data.matrix(matrix(0, out@dims$npred, chunksize)))
res <- cbind(res, data.matrix(matrix(0, out@dims$npred, chunksize)))
mrss <- cbind(mrss, data.matrix(matrix(0, out@dims$npred, chunksize)))
}
}
chunkvox <- chunksize * nchunks
if (out@dims$nvox > chunkvox) {
chunksize <- out@dims$nvox - chunkvox
beta <- cbind(beta, data.matrix(matrix(0, out@dims$npred, chunksize)))
res <- cbind(res, data.matrix(matrix(0, out@dims$npred, chunksize)))
mrss <- cbind(mrss, data.matrix(matrix(0, out@dims$npred, chunksize)))
}
out@beta <- beta
out@res <- res
out@mrss <- mrss
return(out)
}
#' @export
#' @docType methods
#' @details \strong{iModelSolve} Solve already initialized iModel for
#' coefficients, residuals, and mean residual sum of squares.
#' @rdname iModel-methods
iModelSolve <- function(object, verbose = TRUE) {
chunksize <- object@Y@chunksize[2]
nchunks <- floor(object@dims$nvox / chunksize)
vrange <- seq_len(chunksize)
if (verbose)
progress <- txtProgressBar(min = 0, max = nchunks, style = 3)
for (j in seq_len(nchunks)) {
if (vrange[chunksize] > object@dims$nvox)
vrange <- vrange[1]:object@dims$nvox
KWY <- .filter(object@X$K, object@X$W %*% object@Y[, vrange])
object@beta[, vrange] <- object@X$pKX %*% KWY
object@res[, vrange] <- .res(object@X$KWX, KWY)
object@mrss[, vrange] <- colSums(object@res[, vrange]^2) / object@X$trRV
if (object@control$scr)
object@res[, vrange] <- t(t(object@res[, vrange]) * (1 / as.numeric(object@mrss[, vrange])))
vrange <- vrange + chunksize
if (verbose)
setTxtProgressBar(progress, j)
}
if (verbose)
close(progress)
return(object)
}
#' @export
#' @docType methods
#' @details \strong{iModelUpdate} Primarily used to update X slot after
#' optimization.
#' @rdname iModel-methods
iModelUpdate <- function(object, ...) {
ll <- c(...)
if (!is.null(ll$weights)) {
if (class(weights) == "numeric" && length(weights) == dims$nimg)
object@X$W <- diag(weights)
else if (all(dim(weights) == dims$nimg))
object@X$W <- weights
else
stop("weights must be a matrix of nimg x nimg or a vector of length nimg")
}
# set weighted design matrix
object@X$KWX <- .setx(.filter(object@X$K, object@X$W %*% object@X$X))
# pseudoinverse of X
object@X$pKX <- .pinvx(object@X$KWX)
object@X$V <- .filter(object@X$K, .filter(object@X$K, object@X$W %*% tcrossprod(object@xVi$V, object@X$W)))
object@X$betaCov <- object@X$pKX %*% tcrossprod(object@X$V, object@X$pKX)
out <- .trRV(object@X$KWX, object@X$V)
object@X$trRV <- out$trRV
object@X$trRVRV <- out$trRVRV
object@X$rdf <- out$rdf
return(object)
}
#' @export
#' @docType methods
#' @details \strong{iModelRFT} Solve for images in resel space.
#' @rdname iModel-methods
iModelRFT <- function(object, verbose = TRUE) {
dims <- dim(object@mask)
D <- length(dims)
dimx <- seq_len(dims[1])
dimx1 <- dimx + 1
if (D > 1) {
dimy <- seq_len(dims[2])
dimy1 <- dimy + 1
}
if (D > 2) {
dimz <- seq_len(dims[3])
dimz1 <- dimz + 1
}
sar <- object@control$sar
# check sample
if (sar > nrow(object@Y) | sar == 0)
sar <- seq_len(nrow(object@Y))
else {
sar <- sample(seq_len(nrow(object@Y)), sar)
scale <- (nrow(object@Y) / object@X$rdf) * (1 / length(sar))
}
# set up for loop----
if (D == 1) {
d1 <- m1 <- matrix(0, dims[1] + 1)
maskar <- as.numeric(object@mask)
m1[dimx1] <- maskar
m3 <- ((m1[dimx1] * m1[dimx]))
Vxx <- matrix(0, dims[1])
} else if (D == 2) {
d1 <- m1 <- matrix(0, dims[1] + 1, dims[2] + 1)
maskar <- as.matrix(object@mask)
m1[dimx1, dimy1, dimz1] <- maskar
m3 <- ((m1[dimx1, dimy1] * m1[dimx, dimy1])) *
((m1[dimx1, dimy1] * m1[dimx1, dimy]))
Vxx <- Vyy <- Vxy <- matrix(0, dims[1], dims[2])
} else if (D == 3) {
d1 <- m1 <- array(0, dim = dims + 1)
maskar <- as.array(object@mask)
m1[dimx1, dimy1, dimz1] <- maskar
m3 <- ((m1[dimx1, dimy1, dimz1] * m1[dimx, dimy1, dimz1])) *
((m1[dimx1, dimy1, dimz1] * m1[dimx1, dimy, dimz1])) *
((m1[dimx1, dimy1, dimz1] * m1[dimx1, dimy1, dimz])) # mask to eliminate voxels without neighbors
Vxx <- Vyy <- Vzz <- Vxy <- Vxz <- Vyz <- array(0, dim = dims)
}
# partial derivatives of each image----
if (verbose)
progress <- txtProgressBar(min = 0, max = length(sar), style = 3)
for (i in sar) {
if (D == 1)
d1[dimx1] <- makeImage(object@mask, object@res[i,] / mrss[])[dimx]
else if (D == 2)
d1[dimx1, dimy1] <- makeImage(object@mask, object@res[i,] / mrss[])[dimx, dimy]
else if (D == 3)
d1[dimx1, dimy1, dimz1] <- makeImage(object@mask, object@res[i,] / mrss[])[dimx, dimy, dimz]
if (D == 1) {
dx <- (d1[dimx1] - d1[dimx]) * m3
} else if (D == 2) {
dx <- (d1[dimx1, dimy1] - d1[dimx, dimy1]) * m3
dy <- (d1[dimx1, dimy1] - d1[dimx1, dimy]) * m3
} else if (D == 3) {
dx <- (d1[dimx1, dimy1, dimz1] - d1[dimx, dimy1, dimz1]) * m3
dy <- (d1[dimx1, dimy1, dimz1] - d1[dimx1, dimy, dimz1]) * m3
dz <- (d1[dimx1, dimy1, dimz1] - d1[dimx1, dimy1, dimz]) * m3
}
Vxx <- Vxx + (dx * dx)
if (D > 1) {
Vyy <- Vyy + (dy * dy)
Vxy <- Vxy + (dx * dy)
}
if (D > 2) {
Vzz <- Vzz + (dz * dz)
Vxz <- Vxz + (dx * dz)
Vyz <- Vyz + (dy * dz)
}
if (verbose)
setTxtProgressBar(progress, i)
}
if (verbose)
close(progress)
# scale variances/covariances----
Vxx <- Vxx * scale
if (D > 1) {
Vyy <- Vyy * scale
Vxy <- Vxy * scale
}
if (D > 2) {
Vzz <- Vzz * scale
Vxz <- Vxz * scale
Vyz <- Vyz * scale
}
if (D == 1) {
xyz <- Vxx * m3
} else if (D == 2) {
xyz <- cbind(matrix(Vxx * m3, ncol = 1), matrix(Vyy * m3, ncol = 1))
rpv <- (Vxx * Vyy ) + (Vxy * 2) # this needs to be checked
} else if (D == 3) {
xyz <- cbind(Vxx * m3, Vyy * m3, Vzz * m3)
rpv <- (Vxx * Vyy * Vzz) +
(Vxy * Vyz * Vxz * 2) -
(Vyz * Vyz * Vxx) -
(Vxy * Vxy * Vzz) -
(Vxz * Vxz * Vyy)
}
out <- list()
# make RPV Image----
rpv[rpv < 0] <- 0
rpv <- sqrt(rpv / (4 * log(2))^D)
out$rpvImage <- as.antsImage(rpv * maskar)
# estimate fwhm----
xyz <- sqrt((xyz) / (4 * log(2)))
nvox <- sum(m3)
rpv <- sum(rpv) / nvox
xyz <- colSums(xyz) / nvox
resels <- rpv^(1 / D) * (xyz / prod(xyz)^(1 / D))
out$fwhm <- 1 / resels
out$resels <- resels(object@mask, out$fwhm)
return(out)
}
#' @export
#' @docType methods
#' @details \strong{getImages} Returns antsImages of specified contrasts within model iModel.
#' @rdname iModel-methods
getImages <- function(object, contrast, rpv = FALSE, mask = FALSE) {
out <- list()
for (i in seq_len(length(contrast))) {
out[[i]]$clusterImage <- antsImageClone(object@C[[contrast]]$clusterImage)
out[[i]]$contrastImage <- antsImageClone(object@C[[contrast]]$clusterImage)
}
names(out) <- contrast
if (rpv)
out$rpvImage <- antsImageClone(object@dims$rpvImage)
if (mask)
out$mask <- antsImageClone(object@iData@iList[[object@y]]@mask)
return(out)
}
#' iModel Methods
#'
#' V
#'
#' @param object Object of class iModel.
#' @param filename h5 file to save iModel to.
#' @param iData_dirname Directory for iData component of iModel.
#' @param dirname Directory to place \code{report} output.
#' @param contrast Name of contrast.
#' @param rpv Return resels per voxel image?
#' @param mask Return mask used for iModel?
#' @param verbose Enables verbose output. (default = \code{TRUE}).
#' @param ... Additional named arguments passed to \code{iModelUpdate}.
#'
#' @author Zachary P. Christensen
#'
#' @name iModel-methods
NULL
#' @export
#' @docType methods
#' @details \strong{coef} Retrieve fitted coefficients from iModel object.
#' @rdname iModel-methods
coef.iModel <- function(object) {
object@beta[]
}
#' @export
#' @docType methods
#' @details \strong{fitted} Retrieve fitted values from iModel object.
#' @rdname iModel-methods
fitted.iModel <- function(object) {
object@X$X %*% object@beta[]
}
#' @export
#' @docType methods
#' @details \strong{resid} Retrieve residuals from iModel object.
#' @rdname iModel-methods
resid.iModel <- function(object) {
object@res[]
}
#' @export
#' @docType methods
#' @details \strong{iModelRead} Read/load iModel object.
#' @rdname iModel-methods
iModelRead <- function(filename) {
if (!file.exists(filename))
stop("filename does not exist.")
x <- iModel()
file <- h5file(filename)
# load big matrices
x@mrss <- file["mrss"]
x@beta <- file["beta"]
x@res <- file["res"]
x@Y <- file["Y"]
x@location <- filename
x@y <- file["y"][]
x@method <- file["method"][]
# read control
x@control$cf <- file["control/cf"][]
x@control$scr <- file["control/scr"][]
x@control$sar <- file["control/sar"][]
x@control$n <- file["control/n"][]
x@control$iso <- file["control/iso"][]
x@control$os <- file["control/os"][]
x@control$rft <- file["control/rft"][]
# read dims
x@dims$npred <- file["dims/npred"][]
x@dims$nvox <- file["dims/nvox"][]
x@dims$nimg <- file["dims/nimg"][]
if (x@control$rft) {
x$dims$fwhm <- file["dims/fwhm"][]
x@dims$resels <- file["dims/resels"][]
x@dims$rpvImage <- as.antsImage(file["dims/rpvImage"][])
k = antsSetSpacing(x@dims$rpvImage, h5attr(file["dims/rpvImage"] , "spacing"))
k = antsSetDirection(x@dims$rpvImage, h5attr(file["dims/rpvImage"] , "direction"))
k = antsSetOrigin(x@dims$rpvImage, h5attr(file["dims/rpvImage"] , "origin"))
}
# read X
x@X$X <- file["X/X"][]
colnames(x@X$X) <- h5attr(file["X/X"], "colnames")[]
x@X$W <- file["X/W"][]
x@X$pKX <- file["X/pKX"][]
x@X$V <- file["X/V"][]
x@X$betaCov <- file["X/betaCov"][]
x@X$trRV <- file["X/trRV"][]
x@X$trRVRV <- file["X/trRVRV"][]
x@X$rdf <- file["X/rdf"][]
x@X$KWX$X <- file["X/KWX/X"][]
x@X$KWX$v <- file["X/KWX/v"][]
x@X$KWX$u <- file["X/KWX/u"][]
x@X$KWX$d <- file["X/KWX/d"][]
x@X$KWX$tol <- file["X/KWX/tol"][]
x@X$KWX$rk <- file["X/KWX/rk"][]
K <- data.frame(Filters = x@file[file.path("K", "Filters")][],
HParam = x@file[file.path("K", "HParam")][],
RT = x@file[file.path("K", "RT")][])
x@X$K <- .filter(K)
# read xVi
if (file["xVi/Vi"][] != "null")
x@xVi$Vi <- file["xVi/Vi"][]
if (file["xVi/V"][] != "null")
x@xVi$V <- file["xVi/V"][]
if (file["xVi/h"][] != "null")
x@xVi$h <- file["xVi/h"][]
if (file["xVi/Cy"][] != "null")
x@xVi$Cy <- file["xVi/Cy"][]
# read contrasts
ncon <- file["C/count"][]
for (i in seq_len(x@C)) {
cname <- paste("C/C", i, sep = "")
x@C[[i]]$name <- file[file.path(cname, "name")][]
x@C[[i]]$c <- file[file.path(cname, "c")][]
x@C[[i]]$X1 <- file[file.path(cname, "X1")][]
x@C[[i]]$X0 <- file[file.path(cname, "X0")][]
x@C[[i]]$iX0file[file.path(cname, "iX0")][]
x@C[[i]]$fieldType <- file[file.path(cname, "fieldType")][]
x@C[[i]]$results <- file[file.path(cname, "results")][]
x@C[[i]]$sthresh <- file[file.path(cname, "sthresh")][]
x@C[[i]]$cthresh <- file[file.path(cname, "cthresh")][]
x@C[[i]]$contrastImage <- as.antsImage(file[file.path(cname, "contrastImage")][])
x@C[[i]]$contrastImage <- antsSetSpacingh5attr(file[file.path(cname, "contrastImage")], "spacing")
x@C[[i]]$contrastImage <- antsSetDirection(h5attr(file[file.path(cname, "contrastImage")], "direction"))
x@C[[i]]$contrastImage <- antsSetOrigin(h5attr(file[file.path(cname, "contrastImage")], "origin"))
x@C[[i]]$clusterImage <- as.antsImage(file[file.path(cname, "clusterImage")][])
x@C[[i]]$clusterImage <- antsSetSpacing(h5attr(file[file.path(cname, "clusterImage")], "spacing"))
x@C[[i]]$clusterImage <- antsSetDirection(h5attr(file[file.path(cname, "clusterImage")], "direction"))
x@C[[i]]$clusterImage <- antsSetOrigin(h5attr(file[file.path(cname, "clusterImage")], "origin"))
x@C[[i]]$dims$trMV <- file[file.path(cname, "dims", "trMV")][]
x@C[[i]]$dims$trMVMV <- file[file.path(cname, "dims", "trMVMV")][]
x@C[[i]]$dims$idf <- file[file.path(cname, "dims", "idf")][]
}
h5close(file)
return(x)
}
#' @export
#' @docType methods
#' @details \strong{iModelWrite} Read/load iModel object.
#' @rdname iModel-methods
# @describeIn iModel write/save iModel object
iModelWrite <- function(object, filename) {
if (class(object) != "iModel")
stop("object must be of class iModel.")
if (file.exists(filename)) {
if (object@location != filename)
stop("filename already exists.")
}
if (!usePkg("h5"))
stop("Please install package h5 in order to use this function.")
file <- h5file(filename)
#### write big matrices
chunksize <- object@beta@chunksize[2]
chunkseq <- seq_len(chunksize)
nvox <- object@Y@dim[2]
nchunk <- floor(nvox / chunksize)
for (i in seq_len(nchunk)) {
if (i == 1) {
# initialize files
file["beta"] <- object@beta[, chunkseq]
file["res"] <- object@res[, chunkseq]
file["mrss"] <- object@mrss[, chunkseq]
file["Y"] <- object@Y[, chunkseq]
beta <- file["beta"]
res <- file["res"]
mrss <- file["mrss"]
Y <- file["Y"]
} else {
beta <- cbind(beta, object@beta[, chunkseq])
res <- cbind(res, object@res[, chunkseq])
mrss <- cbind(mrss, object@mrss[, chunkseq])
Y <- cbind(Y, object@Y[, chunkseq])
}
chunkseq <- chunkseq + chunksize
}
# fill in any remaining chunkage
if (nvox >= chunkseq[1]) {
chunkseq <- chunkseq[1]:nvox
beta <- cbind(beta, object@beta[, chunkseq])
res <- cbind(res, object@res[, chunkseq])
mrss <- cbind(mrss, object@mrss[, chunkseq])
Y <- cbind(Y, object@Y[, chunkseq])
}
#### end big matrices
file["y"] <- object@y
file["method"] <- object@method
# write control
file["control/cf"] <- object@control$cf
file["control/scr"] <- object@control$scr
file["control/sar"] <- object@control$sar
file["control/n"] <- object@control$n
file["control/iso"] <- object@control$iso
file["control/os"] <- object@control$os
file["control/rft"] <- object@control$rft
# write dims
file["dims/npred"] <- object@dims$npred
file["dims/nvox"] <- object@dims$nvox
file["dims/nimg"] <- object@dims$nimg
if (object@control$rft) {
file["dims/fwhm"] <- object@dims$fwhm
file["dims/resels"] <- object@dims$resels
file["dims/rpvImage"] <- as.array(object@dims$rpvImage)
h5attr(file["dims/rpvImage"] , "spacing") <- antsGetSpacing(object@dims$rpvImage)
h5attr(file["dims/rpvImage"] , "direction") <- antsGetDirection(object@dims$rpvImage)
h5attr(file["dims/rpvImage"] , "origin") <- antsGetOrigin(object@dims$rpvImage)
}
# write X
file[file.path("K", "Filters")] <- object@iData@iList[[object@y]]@K$Filters
file[file.path("K", "HParam")] <- object@iData@iList[[object@y]]@K$HParam
file[file.path("K", "RT")] <- object@iData@iList[[object@y]]@K$RT
file["X/X"] <- data.matrix(object@X$X)
h5attr(file["X/X"], "colnames") <- colnames(object@X$X)
file["X/W"] <- data.matrix(object@X$W)
file["X/pKX"] <- data.matrix(object@X$pKX)
file["X/V"] <- data.matrix(object@X$V)
file["X/betaCov"] <- data.matrix(object@X$betaCov)
file["X/trRV"] <- object@X$trRV
file["X/trRVRV"] <- object@X$trRVRV
file["X/rdf"] <- object@X$rdf
file["X/KWX/X"] <- object@X$KWX$X
file["X/KWX/v"] <- object@X$KWX$v
file["X/KWX/u"] <- object@X$KWX$u
file["X/KWX/d"] <- object@X$KWX$d
file["X/KWX/tol"] <- object@X$KWX$tol
file["X/KWX/rk"] <- object@X$KWX$rk
# write xVi
file["xVi/Vi"] <- if (is.null(object@xVi$Vi)) "null" else object@xVi$Vi
file["xVi/V"] <- if (is.null(object@xVi$V)) "null" else object@xVi$V
file["xVi/h"] <- if (is.null(object@xVi$h)) "null" else object@xVi$h
file["xVi/Cy"] <- if (is.null(object@xVi$Cy)) "null" else object@xVi$Cy
# write contrasts
file["C/count"] <- length(object@C)
for (i in seq_len(object@C)) {
cname <- paste("C/C", i, sep = "")
file[file.path(cname, "name")] <- object@C[[i]]$name
file[file.path(cname, "c")] <- object@C[[i]]$c
file[file.path(cname, "X1")] <- object@C[[i]]$X1
file[file.path(cname, "X0")] <- object@C[[i]]$X0
file[file.path(cname, "iX0")] <- object@C[[i]]$iX0
file[file.path(cname, "fieldType")] <- object@C[[i]]$fieldType
file[file.path(cname, "results")] <- object@C[[i]]$results
file[file.path(cname, "sthresh")] <- object@C[[i]]$sthresh
file[file.path(cname, "cthresh")] <- object@C[[i]]$cthresh
file[file.path(cname, "contrastImage")] <- as.array(object@C[[i]]$contrastImage)
h5attr(file[file.path(cname, "contrastImage")], "spacing") <- antsGetSpacing(object@C[[i]]$contrastImage)
h5attr(file[file.path(cname, "contrastImage")], "direction") <- antsGetDirection(object@C[[i]]$contrastImage)
h5attr(file[file.path(cname, "contrastImage")], "origin") <- antsGetOrigin(object@C[[i]]$contrastImage)
file[file.path(cname, "clusterImage")] <- as.array(object@C[[i]]$clusterImage)
h5attr(file[file.path(cname, "clusterImage")], "spacing") <- antsGetSpacing(object@C[[i]]$clusterImage)
h5attr(file[file.path(cname, "clusterImage")], "direction") <- antsGetDirection(object@C[[i]]$clusterImage)
h5attr(file[file.path(cname, "clusterImage")], "origin") <- antsGetOrigin(object@C[[i]]$clusterImage)
file[file.path(cname, "dims", "trMV")] <- object@C[[i]]$dims$trMV
file[file.path(cname, "dims", "trMVMV")] <- object@C[[i]]$dims$trMVMV
file[file.path(cname, "dims", "idf")] <- object@C[[i]]$dims$idf
}
h5close(file)
return(TRUE)
}
#' @export
#' @docType methods
#' @details \strong{iModelSolve} Solve already initialized iModel for
#' coefficients, residuals, and mean residual sum of squares.
#' @rdname iModel-methods
iModelSolve <- function(object, verbose = TRUE) {
chunksize <- object@Y@chunksize[2]
nchunks <- floor(object@dims$nvox / chunksize)
vrange <- seq_len(chunksize)
if (verbose)
progress <- txtProgressBar(min = 0, max = nchunks, style = 3)
for (j in seq_len(nchunks)) {
if (vrange[chunksize] > object@dims$nvox)
vrange <- vrange[1]:object@dims$nvox
KWY <- .filter(object@X$K, object@X$W %*% object@iData@iList[[object@y]]@iMatrix[, vrange])
object@beta[, vrange] <- object@X$pKX %*% KWY
object@res[, vrange] <- .res(object@X$KWX, KWY)
object@mrss[, vrange] <- colSums(object@res[, vrange]^2) / object@X$trRV
if (object@control$scr)
object@res[, vrange] <- t(t(object@res[, vrange]) * (1 / as.numeric(object@mrss[, vrange])))
vrange <- vrange + chunksize
if (verbose)
setTxtProgressBar(progress, j)
}
if (verbose)
close(progress)
return(object)
}
#' @export
#' @docType methods
#' @details \strong{iModelUpdate} Primarily used to update X slot after optimization
#' @rdname iModel-methods
iModelUpdate <- function(object, ...) {
ll <- c(...)
if (!is.null(ll$weights)) {
if (class(weights) == "numeric" && length(weights) == dims$nimg)
object@X$W <- diag(weights)
else if (all(dim(weights) == dims$nimg))
object@X$W <- weights
else
stop("weights must be a matrix of nimg x nimg or a vector of length nimg")
}
# set weighted design matrix
object@X$KWX <- .setx(.filter(object@X$K, object@X$W %*% object@X$X))
# pseudoinverse of X
object@X$pKX <- .pinvx(object@X$KWX)
object@X$V <- .filter(object@X$K, .filter(object@X$K, object@X$W %*% tcrossprod(object@xVi$V, object@X$W)))
object@X$betaCov <- object@X$pKX %*% tcrossprod(object@X$V, object@X$pKX)
out <- .trRV(object@X$KWX, object@X$V)
object@X$trRV <- out$trRV
object@X$trRVRV <- out$trRVRV
object@X$rdf <- out$rdf
return(object)
}
#' @export
#' @docType methods
#' @details \strong{model.matrix} Retrieve design matrix from iModel object.
#' @rdname iModel-methods
model.matrix.iModel <- function(object) {
return(object@X$X)
}
#' @export
#' @docType methods
#' @details \strong{report} Creates a .pdf and .html report of for iModel
#' objects (requires package rmarkdown).
#' @rdname iModel-methods
report <- function(object, dirname, surfimg) {
if (!usePkg("rmarkdown"))
stop("Please install package rmarkdown in order to use this function.")
ncon <- length(object@C)
if (missing(surfimg))
surfimg <- object@mask
dirname <- normalizePath(dirname)
if (!file.exists(dirname))
dir.create(dirname)
## render images----
for (i in seq_len(ncon)) {
if (class(object@C[[i]]$results) != "character") {
brain <- renderSurfaceFunction(surfimg = list(surfimg), funcimg = list(object@C[[i]]$clusterImage), alphasurf = 0.1, smoothsval = 1.5)
tmp <- make3ViewPNG(fnprefix = paste(dirname, "/", object@C[[i]]$name, sep = ""))
}
}
md <- file.path(dirname, "out.Rmd")
zz <- file(md, open = "wt")
sink(zz)
sink(zz, type = "message")
# describe the model----
cat("## iModel object fit by call to ", object@method[1])
if (object@control$rft)
cat(" using random field theory. \n")
else
cat(". \n\n")
cat(" Predictors = ", paste(colnames(object@X$X), collapse = ", "), "\n\n")
cat(" Images = ", object@dims$nimg, "\n\n")
cat(" Degrees of freedom = ", object@X$rdf, "\n\n")
cat(" Voxels = ", object@dims$nvox, "\n\n")
cat(" Optimization = ", object@control$opt, "\n\n")
if (object@control$rft && length(object@dims$fwhm > 0)) {
cat(" FWHM = ", round(object@dims$fwhm, 2), "\n\n")
cat(" Resels = ", round(object@dims$resels), "\n\n")
}
cat("---- \n\n")
# contrast results----
cat("## Contrast Results \n\n")
for (i in seq_len(ncon)) {
cat("### ", object@C[[i]]$name, "\n\n")
if (class(object@C[[i]]$results) != "character")
cat(paste("![", object@C[[i]]$name, "](", file.path(dirname, paste(object@C[[i]]$name, ".png", sep = "")), ") \n\n", sep = ""))
## render results----
cat("Contrast weights: \n\n")
c <- t(object@C[[i]]$c)
colnames(c) <- colnames(object@X$X)
print(c)
cat("\n\n")
if (object@control$rft) {
if (class(object@C[[i]]$results) == "character") {
cat(object@C[[i]]$results, "\n\n")
} else {
cat("#### Set-level \n\n")
cat(" Clusters = ", ncol(object@C[[i]]$results$clusterLevel), "\n\n")
cat(" p-value = ", object@C[[i]]$results$setLevel, "\n\n")
cat("#### Cluster-Level: \n\n")
tmp <- knitr::kable(round(object@C[[i]]$results$clusterLevel, 3))
print(tmp)
cat("\n\n")
cat("#### Peak-level: \n\n")
tmp <- knitr::kable(round(object@C[[i]]$results$peakLevel, 3))
print(tmp)
cat("\n\n")
}
} else {
if (class(object@C[[i]]$results) == "character") {
cat(object@C[[i]]$results, "\n\n")
} else {
knitr::kable(object@C[[i]]$results)
}
}
cat("Interest degrees of freedom = ", object@C[[i]]$dims$idf, "\n\n")
cat(" Statistical threshold = ", round(object@C[[i]]$sthresh, 2), "\n\n")
cat(" Cluster threshold = ", object@C[[i]]$cthresh, "\n\n")
cat(" Threshold type = ", object@C[[i]]$threshType, "\n\n")
cat("---- \n\n")
}
sink(type = "message")
sink()
rmarkdown::render(md)
rmarkdown::render(md, pdf_document())
return(TRUE)
}
#' Control parameters for RFT based analyses
#'
#' Auxillary function for controlling \code{\link{iModel}} fitting.
#'
#' @param cf Critical F-threshold for selecting voxels over which the non-sphericity is estimated (default = \code{0.001}).
#' @param scr Logical. scale residuals? (default = \code{TRUE}).
#' @param sar Number of residual images to sample for estimating the FWHM (default = \code{64}).
#' @param n images in conjunction (default = \code{1}).
#' @param iso logical. should images be assumed to be isotropic? (default = \code{TRUE}).
#' @param os offset weighting for iteratively reweighted least squares (default = \code{3}).
#' @param rft logical. should voxels be estimated in resel space for random field theory analysis (default = \code{TRUE}).
#' @return
#'
#' @export iControl
iControl <- function(cf = 0.05, scr = TRUE, sar = 64, n = 1, iso = TRUE, os = 3, rft = TRUE) {
list(cf = cf,
scr = scr,
sar = sar,
n = n,
iso = iso,
os = os,
rft = rft)
}
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.