#' Convenience wrapper for 2-view eigenanatomy decomposition.
#'
#' Decomposes two matrices into paired sparse eigenevectors to maximize
#' canonical correlation. Note: we do not scale the matrices internally.
#' We leave scaling choices to the user.
#'
#' @param inmatrix input as inmatrix=list(mat1,mat2). n by p input matrix and n
#' by q input matrix , spatial variable lies along columns.
#' @param inmask optional pair of antsImage masks
#' @param sparseness a c(.,.) pair of values e.g c(0.01,0.1) enforces an
#' unsigned 99 percent and 90 percent sparse solution for each respective view
#' @param nvecs number of eigenvector pairs
#' @param its number of iterations, 10 or 20 usually sufficient
#' @param cthresh cluster threshold pair
#' @param statdir temporary directory if you want to look at full output
#' @param perms number of permutations. settings permutations greater than 0
#' will estimate significance per vector empirically. For small datasets, these
#' may be conservative. p-values depend on how one scales the input matrices.
#' @param uselong enforce solutions of both views to be the same - requires
#' matrices to be the same size
#' @param z subject space (low-dimensional space) sparseness value
#' @param smooth smooth the data (only available when mask is used)
#' @param robust rank transform input matrices
#' @param mycoption enforce 1 - spatial orthogonality, 2 - low-dimensional
#' orthogonality or 0 - both
#' @param initializationList initialization for first view
#' @param initializationList2 initialization for 2nd view
#' @param ell1 gradient descent parameter, if negative then l0 otherwise use l1
#' @param priorWeight Scalar value weight on prior between 0 (prior is weak)
#' and 1 (prior is strong). Only engaged if initialization is used
#' @param verbose activates verbose output to screen
#' @param rejector rejects small correlation solutions
#' @param maxBased boolean that chooses max-based thresholding
#' @return outputs a decomposition of a pair of matrices
#' @author Avants BB
#' @examples
#'
#' mat <- replicate(100, rnorm(20))
#' mat2 <- replicate(100, rnorm(20))
#' mat <- scale(mat)
#' mat2 <- scale(mat2)
#' mydecom <- sparseDecom2(
#' inmatrix = list(mat, mat2),
#' sparseness = c(0.1, 0.3),
#' nvecs = 3, its = 3, perms = 0
#' )
#' wt <- 0.666
#' mat3 <- mat * wt + mat2 * (1 - wt)
#' mydecom <- sparseDecom2(
#' inmatrix = list(mat, mat3),
#' sparseness = c(0.2, 0.2), nvecs = 5, its = 10, perms = 5
#' )
#'
#' \dontrun{
#' # a masked example
#' im <- antsImageRead(getANTsRData("r64"))
#' mask <- thresholdImage(im, 250, Inf)
#' dd <- sum(mask == 1)
#' mat1 <- matrix(rnorm(dd * 10), nrow = 10)
#' mat2 <- matrix(rnorm(dd * 10), nrow = 10)
#' initlist <- list()
#' for (nvecs in 1:2) {
#' init1 <- antsImageClone(mask)
#' init1[mask == 1] <- rnorm(dd)
#' initlist <- lappend(initlist, init1)
#' }
#' ff <- sparseDecom2(
#' inmatrix = list(mat1, mat2), inmask = list(mask, mask),
#' sparseness = c(0.1, 0.1), nvecs = length(initlist), smooth = 1,
#' cthresh = c(0, 0), initializationList = initlist, ell1 = 11
#' )
#' ### now SNPs ###
#' rf <- usePkg("randomForest")
#' bg <- usePkg("BGLR")
#' if (bg & rf) {
#' data(mice)
#' snps <- mice.X
#' numericalpheno <- as.matrix(mice.pheno[, c(4, 5, 13, 15)])
#' numericalpheno <- residuals(lm(numericalpheno ~
#' as.factor(mice.pheno$Litter)))
#' nfolds <- 6
#' train <- sample(rep(c(1:nfolds), 1800 / nfolds))
#' train <- (train < 4)
#' snpd <- sparseDecom2(
#' inmatrix = list(
#' (as.matrix(snps[train, ])),
#' numericalpheno[train, ]
#' ), nvecs = 20, sparseness = c(0.001, -0.5),
#' its = 3, ell1 = 0.1, z = -1
#' )
#' for (j in 3:3) {
#' traindf <- data.frame(
#' bmi = numericalpheno[train, j],
#' snpse = as.matrix(snps[train, ]) %*% as.matrix(snpd$eig1)
#' )
#' testdf <- data.frame(
#' bmi = numericalpheno[!train, j],
#' snpse = as.matrix(snps[!train, ]) %*% as.matrix(snpd$eig1)
#' )
#' myrf <- randomForest(bmi ~ ., data = traindf)
#' preddf <- predict(myrf, newdata = testdf)
#' print(cor.test(preddf, testdf$bmi))
#' plot(preddf, testdf$bmi)
#' }
#' } # check bg and rf
#' }
#'
#' @export sparseDecom2
sparseDecom2 <- function(
inmatrix,
inmask = list(NULL, NULL),
sparseness = c(0.01, 0.01),
nvecs = 3,
its = 20,
cthresh = c(0, 0),
statdir = NA,
perms = 0,
uselong = 0,
z = 0,
smooth = 0,
robust = 0,
mycoption = 0,
initializationList = list(),
initializationList2 = list(),
ell1 = 10,
priorWeight = 0,
verbose = FALSE,
rejector = 0,
maxBased = FALSE) {
idim <- 3
# safety 1 & 2
if (!is.null(inmask[[1]])) {
if (!is.antsImage(inmask[[1]]) && !is.na(inmask[[1]])) {
if (ncol(inmatrix[[1]]) != sum(inmask[[1]])) {
stop("Number of columns in matrix X must equal the size of the maskx")
}
}
}
if (!is.null(inmask[[2]])) {
if (!is.antsImage(inmask[[2]]) && !is.na(inmask[[2]])) {
if (ncol(inmatrix[[2]]) != sum(inmask[[2]])) {
stop("Number of columns in matrix Y must equal the size of the masky")
}
}
}
if (is.antsImage(inmask[[1]])) idim <- inmask[[1]]@dimension
if (is.antsImage(inmask[[2]])) idim <- inmask[[2]]@dimension
if (!is.antsImage(inmask[[1]])) {
maskx <- new("antsImage", "float", idim)
} else {
maskx <- antsImageClone(inmask[[1]])
}
if (!is.antsImage(inmask[[2]])) {
masky <- new("antsImage", "float", idim)
} else {
masky <- antsImageClone(inmask[[2]])
}
if (nrow(inmatrix[[1]]) != nrow(inmatrix[[2]])) {
stop("Matrices must have same number of rows")
}
inmask <- list(maskx, masky)
verbose <- as.numeric(verbose)
if (robust > 0) {
inputMatrices <- list(
robustMatrixTransform(inmatrix[[1]]),
robustMatrixTransform(inmatrix[[2]])
)
} else {
inputMatrices <- inmatrix
}
# helper function allows easier R-based permutation
sccaner <- .sparseDecom2helper2(
inputMatrices,
inmask,
sparseness,
nvecs,
its,
cthresh,
statdir,
uselong,
z,
smooth,
robust,
mycoption,
initializationList,
initializationList2,
ell1,
priorWeight,
verbose,
maxBased
)
ccasummary <- data.frame(
corrs = sccaner$corrs,
pvalues = rep(NA, nvecs)
)
if (perms > 0) {
ccasummary$pvalues <- rep(0, nvecs)
nsubs <- nrow(inputMatrices[[1]])
for (permer in 1:perms)
{
m1 <- data.matrix(inputMatrices[[1]][sample(1:nsubs), ])
m2 <- data.matrix(inputMatrices[[2]][sample(1:nsubs), ])
permmatrix <- list(m1, m2)
sccanerp <- .sparseDecom2helper2(
permmatrix,
inmask,
sparseness,
nvecs,
its,
cthresh,
statdir,
uselong,
z,
smooth,
robust,
mycoption,
initializationList,
initializationList2,
ell1,
priorWeight,
verbose,
maxBased
)
counter <- as.numeric(abs(ccasummary$corrs) < abs(sccanerp$corrs))
ccasummary$pvalues <- ccasummary$pvalues + counter
}
ccasummary$pvalues <- ccasummary$pvalues / perms
}
mynames <- paste(0:(nvecs - 1), sep = "")
if (length(mynames) <= 10) mynames <- paste("00", mynames, sep = "")
if (length(mynames) > 10) {
mynames[1:10] <- paste("00", mynames[1:10], sep = "")
tt <- nvecs
if (tt > 99) tt <- 99
mynames[11:tt] <- paste("0", mynames[11:tt], sep = "")
}
mynames <- paste("Variate", mynames, sep = "")
if (nvecs > 1) {
colnames(sccaner$eig1) <- mynames[1:nvecs]
colnames(sccaner$eig2) <- mynames[1:nvecs]
colnames(sccaner$projections) <- mynames[1:nvecs]
colnames(sccaner$projections2) <- mynames[1:nvecs]
}
if (rejector > 0) {
selector <- abs(ccasummary$corrs) > rejector
if (sum(selector) > 1) {
sccaner$eig1 <- sccaner$eig1[, selector]
sccaner$eig2 <- sccaner$eig2[, selector]
ccasummary <- ccasummary[selector, ]
}
}
return(
list(
projections = sccaner$projections,
projections2 = sccaner$projections2,
eig1 = sccaner$eig1,
eig2 = sccaner$eig2,
ccasummary = ccasummary,
sparseness = sparseness
)
)
}
.sparseDecom2helper2 <- function(
inputMatrices,
inmask,
sparseness,
nvecs,
its,
cthresh,
statdir,
uselong,
z,
smooth,
robust,
mycoption,
initializationList,
initializationList2,
ell1,
priorWeight,
verbose,
maxBased) {
outval <- ANTsRCore::sccanCpp(
inputMatrices[[1]],
inputMatrices[[2]],
inmask[[1]],
inmask[[2]],
sparseness[1],
sparseness[2],
nvecs,
its,
cthresh[1],
cthresh[2],
z,
smooth,
initializationList,
initializationList2,
mycoption,
ell1,
verbose,
priorWeight,
maxBased
)
p1 <- inputMatrices[[1]] %*% t(outval$eig1)
p2 <- inputMatrices[[2]] %*% t(outval$eig2)
outcorrs <- diag(cor(p1, p2))
if (priorWeight < 1.e-10) {
myord <- rev(order(abs(outcorrs)))
outcorrs <- outcorrs[myord]
p1 <- p1[, myord]
p2 <- p2[, myord]
outval$eig1 <- outval$eig1[myord, ]
outval$eig2 <- outval$eig2[myord, ]
}
return(
list(
projections = p1,
projections2 = p2,
eig1 = t(outval$eig1),
eig2 = t(outval$eig2),
corrs = outcorrs
)
)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.