R/spatialbayesianlm.R

In stnava/ANTsR: ANTs in R: Quantification Tools for Biomedical Images

#' spatially constrained bayesian regression function.
#'
#' Take a standard lm result and use bayesian regression to impose spatial
#' regularity.
#'
#'
#' @param mylm standard lm result of the form mylm<-lm(ymat~.)
#' @param ymat outcome matrix - usually from imaging data
#' @param mask mask with non-zero entries n-columns of ymat
#' @param smth smoothness parameter
#' @param priorWeight weight on the prior
#' @param nhood size of neighborhood
#' @param regweights weights on rows - size of ymat
#' @param smoothcoeffmat prior coefficient matrix
#' @return bayesian regression solution is output as a list of images
#' @author Avants BB
#' @examples
#'
#' # make some simple data
#' \dontrun{
#' if (!exists("fn")) fn <- getANTsRData("pcasl")
#' asl <- antsImageRead(fn)
#' tr <- antsGetSpacing(asl)[4]
#' aslmean <- getAverageOfTimeSeries(asl)
#' aslmask <- getMask(aslmean, lowThresh = mean(aslmean), cleanup = TRUE)
#' pcaslpre <- aslPerfusion(asl,
#'   dorobust = 0, useDenoiser = NULL, skip = 1,
#'   useBayesian = 0, moreaccurate = 0, verbose = T, mask = aslmask
#' )
#' # user might compare to useDenoiser=FALSE
#' pcasl.parameters <- list(sequence = "pcasl", m0 = pcaslpre$m0)
#' aslmat <- timeseries2matrix(asl, aslmask)
#' tc <- as.factor(rep(c("C", "T"), nrow(aslmat) / 2))
#' dv <- computeDVARS(aslmat)
#' perfmodel <- lm(aslmat ~ tc + stats::poly(dv, 4)) # standard model
#' ssp <- spatialbayesianlm(perfmodel, aslmat, aslmask,
#'   priorWeight = 1.e2, smth = 1.6, nhood = rep(2, 3)
#' )
#' plot(ssp[[1]], slices = "2x16x2", axis = 3)
#' }
#'
#' @export spatialbayesianlm
spatialbayesianlm <- function(
    mylm, ymat, mask, smth = 1, priorWeight = 1, nhood = NA,
    regweights = NA, smoothcoeffmat = NA) {
  if (sum(mask == 1) != ncol(ymat)) {
    print("spatialbayesianlm: mask does not match ymat dimensions")
    return(NA)
  }
  if (all(is.na(nhood))) {
    nhood <- rep(1, mask@dimension)
  }
  if (all(is.na(regweights))) {
    regweights <- ymat
    regweights[] <- 1
  }
  if (is.null(dim(regweights))) {
    regweights <- ymat
    regweights[] <- 1
  }
  if (!all(dim(regweights) == dim(ymat))) {
    regweights <- ymat
    regweights[] <- 1
  }
  if (all(is.na(smoothcoeffmat))) {
    smoothcoeffmat <- mylm$coefficients
  }
  nmatimgs <- list()
  for (i in 1:nrow(smoothcoeffmat)) {
    temp <- antsImageClone(mask)
    temp[mask == 1] <- smoothcoeffmat[i, ]
    temp <- smoothImage(temp, smth)
    nmatimgs[[i]] <- getNeighborhoodInMask(temp, mask, nhood, boundary.condition = "mean")
    smoothcoeffmat[i, ] <- temp[mask == 1]
  }
  covmat <- cov(t(smoothcoeffmat))
  invcov <- solve(covmat + diag(ncol(covmat)) * 1e-06)
  betaideal <- rep(0, ncol(ymat))
  blmX <- model.matrix(mylm)
  betamatrix <- ymat[1:(ncol(blmX) - 1), ] * 0
  posteriorI <- rep(0, ncol(ymat))
  for (v in 1:ncol(ymat)) {
    parammat <- nmatimgs[[1]][, v]
    for (k in 2:length(nmatimgs)) {
      parammat <- cbind(parammat, nmatimgs[[k]][
        ,
        v
      ])
    }
    pcov <- cov(parammat)
    locinvcov <- tryCatch(solve(pcov), error = function(e) {
      return(invcov)
    })
    if (typeof(locinvcov) == "character") {
      locinvcov <- invcov
    }
    prior <- (smoothcoeffmat[, v])
    blm <- bayesianlm(blmX, ymat[, v], prior, locinvcov * priorWeight, regweights = regweights[
      ,
      v
    ])
    betamatrix[, v] <- blm$beta
    posteriorI[v] <- blm$posteriorProbability
  }
  mylist <- matrixToImages(betamatrix, mask)
  mylist <- lappend(mylist, makeImage(mask, posteriorI))
  return(mylist)
}