Merge develop

DESCRIPTION

@@ -64,4 +64,4 @@ OS_type: unix
 BugReports: https://github.com/Mouse-Imaging-Centre/RMINC/issues
 URL: https://github.com/Mouse-Imaging-Centre/RMINC,
     https://wiki.mouseimaging.ca/display/MICePub/RMINC
-RoxygenNote: 7.0.2
+RoxygenNote: 7.2.3

NAMESPACE

@@ -117,6 +117,7 @@ export(compare_models)
 export(components_to_map)
 export(connected_components)
 export(create_mesh)
+export(flexLm)
 export(getRMINCTestData)
 export(hanatAnova)
 export(hanatFDR)

R/flexLm.R

@@ -0,0 +1,87 @@
+#' flexLm
+#'
+#' @param formula 
+#' @param data 
+#' @param y 
+#' @param ... 
+#'
+#' @return
+#' @export
+#'
+#' @examples
+flexLm <- function(formula, data, y, ...) {
+  # collect all the datavars that were given in the ... flexible argument list
+  dataVars <- list(...)
+
+  # they have to be named; throw an error if not
+  if(any(names(dataVars)=="") | any(is.null(names(dataVars)))) {
+    stop("All data variables have to be named")
+  }
+  # could improve this by varnames=lapply(substitute(list(...))[-1], deparse) to get variable names used
+
+  # dimensions of y and all dataVars have to be the same
+  if (length(dataVars) > 0) {
+    for (i in 1:length(dataVars)) {
+      if (any(dim(dataVars[[i]]) != dim(y))) {
+        stop("All data variables and the y variable have to be the same dimension")
+      }
+    }
+  }  
+
+  # now we set up the model matrix, using fake data in the first instance
+  fakeData <- data
+  fakeData$y <- y[1,]
+  fakeData2 <- fakeData
+  if (length(dataVars) > 0) {
+    for (i in 1:length(dataVars)) {
+      # two sets of fake data - will let us evaluate whether the formula shortcut is working
+      fakeData[,names(dataVars)[i]] <- rnorm(nrow(fakeData))
+      fakeData2[,names(dataVars)[i]] <- 0
+    }
+  }
+
+  # generate two model matrices, one from each set of fake data
+  mmatrix <- model.matrix(formula(formula), fakeData)
+  mmatrix2 <- model.matrix(formula(formula), fakeData2)
+
+  # find the matrix columns that differ between the model matrices - these
+  # will be the ones that depend on data that varies at every voxel/vertex/whatever
+  mmatrixDiff <- mmatrix
+  mmatrixDiff[,] <- as.integer( !(mmatrix == mmatrix2))
+
+  # now test wheter a simple multiplication of the two matrices produces the correct model matrix
+  # if it does we can be much faster in looping over voxels/vertices/whatever, whereas
+  # if it doesn't we have to evaluate the formula every time.
+  # the shortcut should work if there are no transformations of the voxel data (i.e. taking squares, splines, etc)
+
+  shortcutWorks <- all( mmatrix == (mmatrix2 + (mmatrixDiff * as.vector(fakeData[,names(dataVars)][[1]]) )))
+
+  # create the output matrix
+  out <- matrix(nrow=nrow(y), ncol=ncol(mmatrix)*2)
+
+  # now loop over every voxel/vertex/whatever
+  for (i in 1:nrow(y)) {
+
+    # if the matrix multiplication shortcut works use it
+    if (shortcutWorks) {
+      tmp <- fastLm(mmatrix2 + (mmatrixDiff * dataVars[[1]][i,]), y[i,])
+    } else {
+      # otherwise re-evaluate the formula at every voxel after assign the right
+      # values to the data frame
+      if (length(dataVars) > 0) {
+        for (j in 1:length(dataVars)) {
+          data[,names(dataVars)[j]] <- dataVars[[j]][i,]
+        }
+      }
+      mmatrix <- model.matrix(formula(formula), data)
+      tmp <- fastLm(mmatrix, y[i,])
+    }
+    # keep the coefficients and t statistics
+    out[i,] <- c(tmp$coefficients, tmp$coefficients/tmp$stderr)
+  }
+
+  # assign names to the outputs and return
+  nnames <- c(names(tmp$coefficients), paste("tstatistic", names(tmp$coefficients),  sep=":"))
+  colnames(out) <- nnames
+  return(out)
+}

R/minc_anatomy.R

@@ -809,8 +809,14 @@ anatCombineStructures <- function(vols, method = "jacobians",
 #' }
 #' @export
 anatApply <- function(vols, grouping = NULL, method=mean, ...) {
-  if(is.null(grouping))
+  if(is.null(grouping)) {
     grouping <- factor(1)
+  }
+
+  if (!is.factor(grouping)) {
+    warning(paste("Coercing", deparse(substitute(grouping)), "to a factor\n"))
+    grouping <- as.factor(grouping)
+  }
 
   ngroups <- length(levels(grouping))
   output <- matrix(nrow=ncol(vols), ncol=ngroups)

R/minc_interface.R

@@ -1131,7 +1131,7 @@ parseLmFormula <- function(formula,data,mf)
   else {
     for (nTerm in 2:length(formula[[3]])){
       # Skip if it is an interaction term
-      if(length(formula[[3]][[nTerm]] > 1))
+      if(length(formula[[3]][[nTerm]]) > 1)
         next
       rCommand = paste("term <- data$",formula[[3]][[nTerm]],sep="")
       # Skip if it is a formula symbol (i.e. *)
@@ -1344,20 +1344,27 @@ runRMINCTestbed <- function(..., dataPath = getOption("RMINC_DATA_DIR", tempdir(
 #' and \code{wget}
 #' @export
 getRMINCTestData <- function(dataPath = getOption("RMINC_DATA_DIR", tempdir()), method = "libcurl") {
-
-  downloadPath <- file.path(dataPath, "rminctestdata.tar.gz")
+  
+  downloadPath <- file.path(dataPath, "rminctestdata.zip")
   extractedPath <- file.path(dataPath, "rminctestdata/")
 
   if(!file.exists(extractedPath)){
+
+    # Download test data
     if(!file.exists(downloadPath)){
       dir.create(dataPath, showWarnings = FALSE, recursive = TRUE)
-      download.file("https://wiki.mouseimaging.ca/download/attachments/1654/rminctestdata2.tar.gz",
+      download.file("https://github.com/Mouse-Imaging-Centre/RMINC-test-data/archive/refs/heads/main.zip",
                     destfile = downloadPath,
                     method = method) # changed from "wget" to stop freakouts on mac
     }
 
-    untar(downloadPath, exdir = dataPath)
+    # Extract test data
+    utils::unzip(downloadPath, exdir = dataPath)
+
+    # Move test data path one directory up, to play nice with the path specs in the testthat test_*.R scripts
+    file.rename(file.path(dataPath, "RMINC-test-data-main/rminctestdata/"), file.path(dataPath, "rminctestdata/"))
 
+    # Fix hardcoded paths
     rectifyPaths <-
       function(file){
         readLines(file) %>%

R/minc_lmer.R

@@ -467,9 +467,9 @@ ranef_summary <-
                   , grouping = group_name
                   , se = NULL) %>%
           gather("var", "value", c("tvalue", "beta")) %>%
-          unite_("groupingXgroup", c("grouping", "group"), sep = "") %>%
-          unite_("varXeffectXgroupingXgroup", c("var", "effect", "groupingXgroup"), sep = "-") %>%
-          spread_("varXeffectXgroupingXgroup", "value") %>%
+          unite("groupingXgroup", c("grouping", "group"), sep = "") %>%
+          unite("varXeffectXgroupingXgroup", c("var", "effect", "groupingXgroup"), sep = "-") %>%
+          spread("varXeffectXgroupingXgroup", "value") %>%
           as.matrix %>%
           .[1,]
       }, e = eff, group_name = names(eff), SIMPLIFY = FALSE) %>%

R/minc_vertex_statistics.R

@@ -25,36 +25,67 @@ NULL
 #' @export
 vertexMean <- function(filenames, column=1) 
 {
-  vertexData = vertexTable(filenames, column=column)
-  return(rowMeans(vertexData))
-
+  m_orig <- match.call()
+  vertexData <- vertexTable(filenames, column=column)
+  result <- as.matrix(rowMeans(vertexData))
+
+  attr(result, "likeVolume") <- as.character(filenames[1])
+  attr(result, "filenames")  <- as.character(filenames)
+  attr(result, "stat-type")  <- c("mean")
+  attr(result, "call")       <- m_orig
+
+  colnames(result) <- c("mean")
+
+  return(result)
 } 
 
 #' @describeIn vertexSummaries sum
 #' @export
 vertexSum <- function(filenames, column=1) 
 {
+  m_orig <- match.call()
   vertexData = vertexTable(filenames,column=column)
-  return(rowSums(vertexData))
-
+  result <- as.matrix(rowSums(vertexData))
+  attr(result, "likeVolume") <- as.character(filenames[1])
+  attr(result, "filenames")  <- as.character(filenames)
+  attr(result, "stat-type")  <- c("sum")
+  attr(result, "call")       <- m_orig
+
+  colnames(result) <- c("sum")
+
+  return(result)
 } 
 
 #' @describeIn vertexSummaries var
 #' @export
 vertexVar <- function(filenames, column=1) 
 {
+  m_orig <- match.call()
   vertexData = vertexTable(filenames,column=column)
-  return(apply(vertexData,1,var))
-
+  result<-as.matrix(apply(vertexData,1,var))
+  attr(result, "likeVolume") <- as.character(filenames[1])
+  attr(result, "filenames")  <- as.character(filenames)
+  attr(result, "stat-type")  <- c("var")
+  attr(result, "call")       <- m_orig
+
+  colnames(result) <- c("var")
+  return(result)
 } 
 
 #' @describeIn vertexSummaries standard deviation
 #' @export
 vertexSd<- function(filenames,column=1) 
 {
+  m_orig <- match.call()
   vertexData = vertexTable(filenames,column=column)
-  return(apply(vertexData,1,sd))
-
+  result<-as.matrix(apply(vertexData,1,sd))
+  attr(result, "likeVolume") <- as.character(filenames[1])
+  attr(result, "filenames")  <- as.character(filenames)
+  attr(result, "stat-type")  <- c("sd")
+  attr(result, "call")       <- m_orig
+
+  colnames(result) <- c("sd")
+  return(result)
 }
 
 ### Helper function for applying over rows of a potentially 

R/minc_vis2D.R

@@ -4,9 +4,15 @@
 #' with the correct spatial dimensions assigned. While this might in the future happen
 #' at the time those objects are created, for the moment this utility function
 #' works with older style RMINC objects and extract what you need.
+#' 
+#' At times data is computed within a mask; for example mincTable matrices
+#' can be extracted from masked data. In that case if the mask volume
+#' is passed as an argument then the data will will be reinserted into
+#' those parts of the mask volume greater than zero.
 #'
 #' @param volume The input volume (from mincLm, mincGetVolume, etc.)
 #' @param dimIndex The index into a multidimensional object
+#' @param maskVolume The optional mask from which the data was created
 #'
 #' @note R uses Fortran indexing, so dimension assignment is c(dim[3], dim[2], dim[1])
 #' once dimensions are obtained from any libminc functions (which use C indexing)
@@ -22,7 +28,7 @@
 #' vs <- mincLm(jacobians ~ genotype, gf)
 #' tvol <- mincArray(vs, 6)
 #' }
-mincArray <- function(volume, dimIndex=1) {
+mincArray <- function(volume, dimIndex=1, maskVolume=NULL) {
   # 1d file with no dimensions (such as the output of mincGetVolume)
   if (is.null(dim(volume))) {
     outvol <- volume
@@ -38,8 +44,22 @@ mincArray <- function(volume, dimIndex=1) {
   else {
     outvol <- volume
   }
-  # we now have a 1d vector - get it's dimensions from the sizes attribute, or read it from the likeVolume
-  if (! is.null(attr(volume, "sizes"))) {
+  # we now have a 1d vector 
+  # if there's a maskVolume present assume that the values need to be reinserted
+  # into the mask first
+  if (!is.null(maskVolume)) {
+    tmpout <- maskVolume
+    tmpout[maskVolume > 0] <- outvol
+    outvol <- tmpout
+  }
+
+  # get it's dimensions from the sizes attribute, or read it from the likeVolume or maskVolume
+  if (!is.null(maskVolume)) {
+    if (length(dim(maskVolume)) !=3) {
+      stop("maskVolume must have 3 dimensions (i.e. be a mincArray itself)")
+    }
+  }
+  else if (! is.null(attr(volume, "sizes"))) {
     sizes <- attr(volume, "sizes")
     dim(outvol) <- c(sizes[3], sizes[2], sizes[1]) # C to Fortran dim ordering
   }

R/minc_voxel_statistics.R

@@ -38,6 +38,7 @@ mincSummary <- function(filenames, grouping=NULL, mask=NULL, method="mean", mask
   if (is.null(grouping)) {
     grouping <- rep(1, length(filenames))
   }
+
   if (is.null(maskval)) {
     minmask = 1
     maxmask = 99999999
@@ -76,27 +77,59 @@ mincSummary <- function(filenames, grouping=NULL, mask=NULL, method="mean", mask
 #' @describeIn mincSummary mean
 #' @export
 mincMean <- function(filenames, grouping=NULL, mask=NULL, maskval=NULL) {
+
+  if(!is.null(grouping)) {
+    if (!is.factor(grouping)) {
+      warning(paste("Coercing", deparse(substitute(grouping)), "to a factor\n"))
+      grouping <- as.factor(grouping)
+    }
+  }
+
   result <- mincSummary(filenames, grouping, mask, method="mean", maskval=maskval)
   return(result)
 }
 
 #' @describeIn mincSummary Variance
 #' @export
 mincVar <- function(filenames, grouping=NULL, mask=NULL, maskval=NULL) {
+
+  if(!is.null(grouping)) {
+    if (!is.factor(grouping)) {
+      warning(paste("Coercing", deparse(substitute(grouping)), "to a factor\n"))
+      grouping <- as.factor(grouping)
+    }
+  }
+
   result <- mincSummary(filenames, grouping, mask, method="var", maskval=maskval)
   return(result)
 }
 
 #' @describeIn mincSummary Sum
 #' @export
 mincSum <- function(filenames, grouping=NULL, mask=NULL, maskval=NULL) {
+
+  if(!is.null(grouping)) {
+    if (!is.factor(grouping)) {
+      warning(paste("Coercing", deparse(substitute(grouping)), "to a factor\n"))
+      grouping <- as.factor(grouping)
+    }
+  }
+
   result <- mincSummary(filenames, grouping, mask, method="sum", maskval=maskval)
   return(result)
 }
 
 #' @describeIn mincSummary Standard Deviation
 #' @export
 mincSd <- function(filenames, grouping=NULL, mask=NULL, maskval=NULL) {
+
+  if(!is.null(grouping)) {
+    if (!is.factor(grouping)) {
+      warning(paste("Coercing", deparse(substitute(grouping)), "to a factor\n"))
+      grouping <- as.factor(grouping)
+    }
+  }
+
   result <- mincSummary(filenames, grouping, mask, method="var", maskval=maskval)
   result <- sqrt(result)
   return(result)
@@ -779,6 +812,7 @@ mincTtest <- function(filenames, grouping, mask=NULL, maskval=NULL) {
   # the grouping for a t test should only contain 2 groups. Should 
   # also be converted to a factor if it's not.
   if( ! is.factor(grouping) ){
+    warning(paste("Coercing", deparse(substitute(grouping)), "to a factor\n"))
     grouping <- as.factor(grouping)
   }
   if(length(levels(grouping)) != 2 ){
@@ -825,6 +859,7 @@ mincPairedTtest <- function(filenames, grouping, mask=NULL, maskval=NULL) {
   # group 1 is paired with element 1 from group 2 etc.), the groups need to have 
   # the same length.
   if( ! is.factor(grouping) ){
+    warning(paste("Coercing", deparse(substitute(grouping)), "to a factor\n"))
     grouping <- as.factor(grouping)
   }
   if(length(levels(grouping)) != 2 ){
@@ -873,6 +908,10 @@ mincPairedTtest <- function(filenames, grouping, mask=NULL, maskval=NULL) {
 #' }
 #' @export
 mincWilcoxon <- function(filenames, grouping, mask=NULL, maskval=NULL) {
+  if( ! is.factor(grouping) ){
+    warning(paste("Coercing", deparse(substitute(grouping)), "to a factor\n"))
+    grouping <- as.factor(grouping)
+  }
   result <- mincSummary(filenames, grouping, mask, method="wilcoxon", maskval=maskval)
   result <- as.matrix(result);
   attr(result, "likeVolume") <- filenames[1]