| Title: | Mixture fMRI Clustering Analysis |
|---|---|
| Description: | Utilizing model-based clustering (unsupervised) for functional magnetic resonance imaging (fMRI) data. The developed methods (Chen and Maitra (2023) <doi:10.1002/hbm.26425>) include 2D and 3D clustering analyses (for p-values with voxel locations) and segmentation analyses (for p-values alone) for fMRI data where p-values indicate significant level of activation responding to stimulate of interesting. The analyses are mainly identifying active voxel/signal associated with normal brain behaviors. Analysis pipelines (R scripts) utilizing this package (see examples in 'inst/workflow/') is also implemented with high performance techniques. |
| Authors: | Wei-Chen Chen [aut, cre], Ranjan Maitra [aut], Dan Nettleton [aut, ctb], Pierre Lafaye De Micheaux [aut, ctb] (Threshold functions from AnalyzeFMRI), Jonathan L Marchini [aut, ctb] (Threshold functions from AnalyzeFMRI) |
| Maintainer: | Wei-Chen Chen <[email protected]> |
| License: | Mozilla Public License 2.0 |
| Version: | 0.1-4 |
| Built: | 2024-10-18 12:31:02 UTC |
| Source: | CRAN |
Utilizing model-based clustering (unsupervised) for fMRI data especially in a distributed manner. The methods includes 2D and 3D clustering analyses and segmentation analyses for fMRI signals where p-values are significant levels of active voxels which respond to stimulate of interesting. The analyses are mainly identifying active voxels/signals from normal brain behaviors. Workflows are also implemented utilizing high performance techniques.
The main function of this package is fclust() that implements
model-based clustering algorithm for fMRI signal data and provides
unsupervised clustering results for the data. Several workflows implemented
with high-performance computing techniques are also built in for automatically
process clustering, hypothesis, cluster merging, and visualizations.
Wei-Chen Chen and Ranjan Maitra.
Chen, W.-C. and Maitra, R. (2023) “A practical model-based segmentation approach for improved activation detection in single-subject functional magnetic resonance imaging studies”, Human Brain Mapping, 44(16), 5309–5335. (doi:10.1002/hbm.26425)
fclust(), set.global().
library(MixfMRI, quietly = TRUE) .rem <- function(){ demo(fclust3d,'MixfMRI',ask=FALSE,echo=FALSE) demo(fclust2d,'MixfMRI',ask=FALSE,echo=FALSE) }library(MixfMRI, quietly = TRUE) .rem <- function(){ demo(fclust3d,'MixfMRI',ask=FALSE,echo=FALSE) demo(fclust2d,'MixfMRI',ask=FALSE,echo=FALSE) }
Main algorithms implemented in fclust.
ecm.step.gbd(PARAM.org) apecma.step.gbd(PARAM.org) em.step.gbd(PARAM.org)ecm.step.gbd(PARAM.org) apecma.step.gbd(PARAM.org) em.step.gbd(PARAM.org)
PARAM.org |
an initialized |
These are main algorithms implemented in fclust().
Return an optimized PARAM.
Wei-Chen Chen and Ranjan Maitra.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
set.global(), fclust(), PARAM,
PARAM.org.
library(MixfMRI, quietly = TRUE) library(EMCluster, quietly = TRUE) # .FC.CT$algorithm <- "em" # .FC.CT$model.X <- "V" # .FC.CT$ignore.X <- TRUE .FC.CT$check.X.unit <- FALSE ### Test toy1. set.seed(1234) X.gbd <- toy1$X.gbd PV.gbd <- toy1$PV.gbd PARAM <- set.global(X.gbd, PV.gbd, K = 2) PARAM.new <- initial.em.gbd(PARAM) PARAM.toy1 <- em.step.gbd(PARAM.new) id.toy1 <- .MixfMRIEnv$CLASS.gbd print(PARAM.toy1$ETA) RRand(toy1$CLASS.gbd, id.toy1) .rem <- function(){ ### Test toy2. set.seed(1234) X.gbd <- toy2$X.gbd PV.gbd <- toy2$PV.gbd PARAM <- set.global(X.gbd, PV.gbd, K = 3) PARAM.new <- initial.em.gbd(PARAM) PARAM.toy2 <- em.step.gbd(PARAM.new) id.toy2 <- .MixfMRIEnv$CLASS.gbd print(PARAM.toy2$ETA) RRand(toy2$CLASS.gbd, id.toy2) }library(MixfMRI, quietly = TRUE) library(EMCluster, quietly = TRUE) # .FC.CT$algorithm <- "em" # .FC.CT$model.X <- "V" # .FC.CT$ignore.X <- TRUE .FC.CT$check.X.unit <- FALSE ### Test toy1. set.seed(1234) X.gbd <- toy1$X.gbd PV.gbd <- toy1$PV.gbd PARAM <- set.global(X.gbd, PV.gbd, K = 2) PARAM.new <- initial.em.gbd(PARAM) PARAM.toy1 <- em.step.gbd(PARAM.new) id.toy1 <- .MixfMRIEnv$CLASS.gbd print(PARAM.toy1$ETA) RRand(toy1$CLASS.gbd, id.toy1) .rem <- function(){ ### Test toy2. set.seed(1234) X.gbd <- toy2$X.gbd PV.gbd <- toy2$PV.gbd PARAM <- set.global(X.gbd, PV.gbd, K = 3) PARAM.new <- initial.em.gbd(PARAM) PARAM.toy2 <- em.step.gbd(PARAM.new) id.toy2 <- .MixfMRIEnv$CLASS.gbd print(PARAM.toy2$ETA) RRand(toy2$CLASS.gbd, id.toy2) }
Calculate contiguous clusters of locations in a 3D array that are above some threshold and with some minimum size.
cluster.threshold(x, nmat = NULL, level.thr = 0.5, size.thr)cluster.threshold(x, nmat = NULL, level.thr = 0.5, size.thr)
x |
A 3D array |
nmat |
A matrix with 3 columns specifying the neighbourhood system. Default is 6 nearest neighbours in 3D. |
level.thr |
The level at which to threshold the array values. Default is 0.5 and is designed to cluster 0-1 arrays. |
size.thr |
The cluster size threshold. |
Note: This function is directly copied from "AnalyzeFMRI".
Returns an array of the same size as x with a 1 at all locations which have a value above level.thr and are in a cluster of similiar locations with size greater than size.thr.
J. L. Marchini
x <- array(0, dim = c(64, 64, 21)) x[10:20, 10:20, 1:5] <- 1 x[30:40, 30:40, 6:7] <- 1 x[50, 50, 8:9] <- 1 a <- cluster.threshold(x, size.thr = 400) sum(x) ## should be 849 sum(a) ## should be 605x <- array(0, dim = c(64, 64, 21)) x[10:20, 10:20, 1:5] <- 1 x[30:40, 30:40, 6:7] <- 1 x[50, 50, 8:9] <- 1 a <- cluster.threshold(x, size.thr = 400) sum(x) ## should be 849 sum(a) ## should be 605
Compute q-values Benjamini and Hochberg's (1995) approach for controlling FDR.
qvalue(p, method = c("BH1995", "BY2001"))qvalue(p, method = c("BH1995", "BY2001"))
p |
a p-value vector. |
method |
using method by either BH1995 or BY2001 |
This function compute q-values using Benjamini and Hochberg's (1995)
approach for controlling FDR. The function bh.fdr is originally
written by Dr. Dan Nettleton.
The Benjamini and Yeekutieli's (2001) approach for controlling FDR using
the function by.fdr is coded by Wei-Chen Chen.
Return corresponding q-values for the input p-values.
Dan Nettleton.
Modified by Wei-Chen Chen.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
library(MixfMRI, quietly = TRUE) set.seed(1234) da <- gendataset(phantom = shepp1fMRI, overlap = 0.01) p <- da$pval[!is.na(da$pval)][1:100] qvalue(p)library(MixfMRI, quietly = TRUE) set.seed(1234) da <- gendataset(phantom = shepp1fMRI, overlap = 0.01) p <- da$pval[!is.na(da$pval)][1:100] qvalue(p)
The function computes statistics for log odds ratio of posterior probability.
logor.stat(x, fcobj, post.z, cov.param = NULL, cov.post.z = NULL, cov.logit.z = NULL, all.x = FALSE, drop.ETA1 = FALSE)logor.stat(x, fcobj, post.z, cov.param = NULL, cov.post.z = NULL, cov.logit.z = NULL, all.x = FALSE, drop.ETA1 = FALSE)
x |
an input list of two elements |
fcobj |
a |
post.z |
a matrix of |
cov.param |
a covariance matrix of |
cov.post.z |
a covariance list of length equal to number of active
voxels, which is also a return of |
cov.logit.z |
a covariance list of length equal to number of active
voxels, which is also a return of |
all.x |
all cov matrices for all observations are returned if TRUE, while for only active observations (those of class ids are greater than 1) if FALSE. |
drop.ETA1 |
if drop the |
For posterior probability, this function compute log odd ratio, cov matrix of log odd ratio, degrees of freedom, and testing statistics.
A list is returned with four elements: log.or, cov.log.or,
df, and test.stat.
Wei-Chen Chen and Ranjan Maitra.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
post.prob(), cov.param(), cov.post.z(),
cov.logit.z().
library(MixfMRI, quietly = TRUE) .FC.CT$model.X <- "I" .FC.CT$CONTROL$debug <- 0 K <- 3 ### Fit toy1. set.seed(1234) X.gbd <- toy1$X.gbd X.range <- apply(X.gbd, 2, range) X.gbd <- t((t(X.gbd) - X.range[1,]) / (X.range[2,] - X.range[1,])) PV.gbd <- toy1$PV.gbd fcobj <- fclust(X.gbd, PV.gbd, K = K, min.1st.prop = 0.5) ### Test log odds ratio. x <- list(X.gbd = X.gbd, PV.gbd = PV.gbd) post.z <- post.prob(x, fcobj) lor <- logor.stat(x, fcobj, post.z) ### Check if 95% CE covers log odd ratio = 1. id <- !is.na(lor$df) id.cover.0 <- which(lor$test.stat[id] < pchisq(0.95, lor$df[id])) ### Get voxels needed for merging. id.active <- which(fcobj$class != 1) id.merge <- id.active[id][id.cover.0] ### Check results. post.z[id.merge,] cbind(toy1$X.gbd[id.merge,], toy1$PV.gbd[id.merge])library(MixfMRI, quietly = TRUE) .FC.CT$model.X <- "I" .FC.CT$CONTROL$debug <- 0 K <- 3 ### Fit toy1. set.seed(1234) X.gbd <- toy1$X.gbd X.range <- apply(X.gbd, 2, range) X.gbd <- t((t(X.gbd) - X.range[1,]) / (X.range[2,] - X.range[1,])) PV.gbd <- toy1$PV.gbd fcobj <- fclust(X.gbd, PV.gbd, K = K, min.1st.prop = 0.5) ### Test log odds ratio. x <- list(X.gbd = X.gbd, PV.gbd = PV.gbd) post.z <- post.prob(x, fcobj) lor <- logor.stat(x, fcobj, post.z) ### Check if 95% CE covers log odd ratio = 1. id <- !is.na(lor$df) id.cover.0 <- which(lor$test.stat[id] < pchisq(0.95, lor$df[id])) ### Get voxels needed for merging. id.active <- which(fcobj$class != 1) id.merge <- id.active[id][id.cover.0] ### Check results. post.z[id.merge,] cbind(toy1$X.gbd[id.merge,], toy1$PV.gbd[id.merge])
These functions compute posterior probabilities, Fisher information with covariance matrix of parameters, covariance matrix of posterior probabilities, and covariance matrix of logit posterior probabilities.
post.prob(x, fcobj) cov.param(x, fcobj, post.z, drop.ETA1 = FALSE) cov.post.z(x, fcobj, post.z, cov.param = NULL, all.x = FALSE, drop.ETA1 = FALSE) cov.logit.z(x, fcobj, post.z, cov.param = NULL, cov.post.z = NULL, all.x = FALSE, drop.ETA1 = FALSE)post.prob(x, fcobj) cov.param(x, fcobj, post.z, drop.ETA1 = FALSE) cov.post.z(x, fcobj, post.z, cov.param = NULL, all.x = FALSE, drop.ETA1 = FALSE) cov.logit.z(x, fcobj, post.z, cov.param = NULL, cov.post.z = NULL, all.x = FALSE, drop.ETA1 = FALSE)
x |
an input list of two elements |
fcobj |
a |
post.z |
a matrix of |
cov.param |
a covariance matrix of |
cov.post.z |
a covariance list of length equal to number of active
voxels, which is also a return of |
all.x |
all cov matrices for all observations are returned if TRUE, while for only active observations (those of class ids are greater than 1) if FALSE. |
drop.ETA1 |
if drop the |
These functions are required to compute covariance matrices of parameters and posterior probabilities.
Use post.prob() to get the posterior probabilities.
Input the returns of post.prob()
to cov.param() to obtain the cov matrix for parameters
(inversed Fisher information obtained from inner product of gradient
of log observed data likelihood). A list is returned with I
for Fisher information, and cov for the covariance matrix which
is inverted by ginv().
Input the returns of post.prob() and cov.param()
to cov.post.z() to obtain the cov matrix for posterior
probabilities by the multivariate delta method on the cov matrix for
parameters.
Input the returns of post.prob(), cov.param(), and
cov.post.z() to cov.logit.z() to obtain cov matrix
for logit posterior probabilities by the multivariate delta method on
cov matrix of posterior probabilities.
A matrix or a list is returned.
The cov.param() will return a list containing two elements
I for the Fisher information, and cov for the covariance matrix
by generalized inversed of the Fisher information. The dimension of both
elements are d * d where d = K * 7 - 4 for 2D data and
d = K * 9 - 4 for 3D data if drop.ETA1 = TRUE, otherwise
they are d = K * 7 - 3 and d = K * 9 -4, respectively.
The cov.post.z() will return a list containing cov matrices of
posterior probabilities for each valid/selected voxel.
The cov.logit.z() will return a list containing cov matrices of
logit posterior probabilities for each valid/selected voxel.
Wei-Chen Chen and Ranjan Maitra.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
EMCluster::lmt(), lmt.I().
library(MixfMRI, quietly = TRUE) library(EMCluster, quietly = TRUE) .FC.CT$model.X <- "I" .FC.CT$CONTROL$debug <- 0 K <- 3 .rem <- function(){ ### Fit toy1. set.seed(1234) X.gbd <- toy1$X.gbd X.range <- apply(X.gbd, 2, range) X.gbd <- t((t(X.gbd) - X.range[1,]) / (X.range[2,] - X.range[1,])) PV.gbd <- toy1$PV.gbd fcobj <- fclust(X.gbd, PV.gbd, K = K, min.1st.prop = 0.5) ### Test cov matrix of posterior z and logit posterior z. x <- list(X.gbd = X.gbd, PV.gbd = PV.gbd) post.z <- post.prob(x, fcobj) cov.param <- cov.param(x, fcobj, post.z = post.z) cov.post.z <- cov.post.z(x, fcobj, post.z = post.z, cov.param = cov.param$cov) cov.logit.z <- cov.logit.z(x, fcobj, post.z = post.z, cov.param = cov.param$cov, cov.post.z = cov.post.z) ### Compute cov matrix of log odds ratio for all k > 1. A <- cbind(rep(-1, K - 1), diag(1, K - 1)) logit.p <- log(post.z[fcobj$class != 1,] / (1 - post.z[fcobj$class != 1,])) log.or <- logit.p %*% t(A) cov.log.or <- lapply(cov.logit.z, function(x) A %*% x %*% t(A)) ### Check if 0 vector covered by 95% confidence ellipsoid. id <- 1 plot(log.or[id,], xlim = log.or[id, 1] + c(-5, 5), ylim = log.or[id, 2] + c(-5, 5), main = "1st observation", xlab = "x", ylab = "y") plotBN(log.or[id,], cov.log.or[[id]]) points(0, 0, col = 2) }library(MixfMRI, quietly = TRUE) library(EMCluster, quietly = TRUE) .FC.CT$model.X <- "I" .FC.CT$CONTROL$debug <- 0 K <- 3 .rem <- function(){ ### Fit toy1. set.seed(1234) X.gbd <- toy1$X.gbd X.range <- apply(X.gbd, 2, range) X.gbd <- t((t(X.gbd) - X.range[1,]) / (X.range[2,] - X.range[1,])) PV.gbd <- toy1$PV.gbd fcobj <- fclust(X.gbd, PV.gbd, K = K, min.1st.prop = 0.5) ### Test cov matrix of posterior z and logit posterior z. x <- list(X.gbd = X.gbd, PV.gbd = PV.gbd) post.z <- post.prob(x, fcobj) cov.param <- cov.param(x, fcobj, post.z = post.z) cov.post.z <- cov.post.z(x, fcobj, post.z = post.z, cov.param = cov.param$cov) cov.logit.z <- cov.logit.z(x, fcobj, post.z = post.z, cov.param = cov.param$cov, cov.post.z = cov.post.z) ### Compute cov matrix of log odds ratio for all k > 1. A <- cbind(rep(-1, K - 1), diag(1, K - 1)) logit.p <- log(post.z[fcobj$class != 1,] / (1 - post.z[fcobj$class != 1,])) log.or <- logit.p %*% t(A) cov.log.or <- lapply(cov.logit.z, function(x) A %*% x %*% t(A)) ### Check if 0 vector covered by 95% confidence ellipsoid. id <- 1 plot(log.or[id,], xlim = log.or[id, 1] + c(-5, 5), ylim = log.or[id, 2] + c(-5, 5), main = "1st observation", xlab = "x", ylab = "y") plotBN(log.or[id,], cov.log.or[[id]]) points(0, 0, col = 2) }
These functions computes covariance matrix of logit ETA.
cov.logit.ETA(x, fcobj, cov.param = NULL)cov.logit.ETA(x, fcobj, cov.param = NULL)
x |
an input list of two elements |
fcobj |
a |
cov.param |
a covariance matrix of |
These functions are required to compute covariance matrices of logit ETA.
Input the returns of cov.param()
to cov.logit.ETA() to obtain the cov matrix for logit ETA
by the multivariate delta method on the cov matrix for parameters.
A matrix.
Wei-Chen Chen and Ranjan Maitra.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
EMCluster::lmt(), lmt.I().
library(MixfMRI, quietly = TRUE) .FC.CT$model.X <- "I" .FC.CT$CONTROL$debug <- 0 K <- 3 .rem <- function(){ ### Fit toy1. set.seed(1234) X.gbd <- toy1$X.gbd X.range <- apply(X.gbd, 2, range) X.gbd <- t((t(X.gbd) - X.range[1,]) / (X.range[2,] - X.range[1,])) PV.gbd <- toy1$PV.gbd fcobj <- fclust(X.gbd, PV.gbd, K = K, min.1st.prop = 0.5) ### Test cov matrix of posterior z. x <- list(X.gbd = X.gbd, PV.gbd = PV.gbd) post.z <- post.prob(x, fcobj) cov.param <- cov.param(x, fcobj, post.z) cov.logit.ETA <- cov.logit.ETA(x, fcobj, cov.param = cov.param$cov) ### Compute cov matrxi of eta_k - eta_1 for all k > 1. A <- cbind(rep(-1, K - 1), diag(1, K - 1)) ETA <- fcobj$param$ETA log.or <- log(ETA / (1 - ETA)) %*% t(A) cov.log.or <- A %*% cov.logit.ETA %*% t(A) }library(MixfMRI, quietly = TRUE) .FC.CT$model.X <- "I" .FC.CT$CONTROL$debug <- 0 K <- 3 .rem <- function(){ ### Fit toy1. set.seed(1234) X.gbd <- toy1$X.gbd X.range <- apply(X.gbd, 2, range) X.gbd <- t((t(X.gbd) - X.range[1,]) / (X.range[2,] - X.range[1,])) PV.gbd <- toy1$PV.gbd fcobj <- fclust(X.gbd, PV.gbd, K = K, min.1st.prop = 0.5) ### Test cov matrix of posterior z. x <- list(X.gbd = X.gbd, PV.gbd = PV.gbd) post.z <- post.prob(x, fcobj) cov.param <- cov.param(x, fcobj, post.z) cov.logit.ETA <- cov.logit.ETA(x, fcobj, cov.param = cov.param$cov) ### Compute cov matrxi of eta_k - eta_1 for all k > 1. A <- cbind(rep(-1, K - 1), diag(1, K - 1)) ETA <- fcobj$param$ETA log.or <- log(ETA / (1 - ETA)) %*% t(A) cov.log.or <- A %*% cov.logit.ETA %*% t(A) }
These functions based on normal assumption and transformation to derive a (mixture) density function of p-values.
dpval(x, mu = 0, log = FALSE) dmixpval(x, eta, mu)dpval(x, mu = 0, log = FALSE) dmixpval(x, eta, mu)
x |
support of p-values which should be between 0 and 1. |
mu |
hypothetical mean of testing statistics (in normal distribution) for producing p-values. |
log |
if return log of density. |
eta |
mixing proportion of |
Note that eta and mu in dmixpval() are of length
K for K component mixtures.
Corresponding density values (to the input x) are returned.
Wei-Chen Chen and Ranjan Maitra.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
gendataset(), qvalue().
library(MixfMRI, quietly = TRUE) set.seed(1234) da <- gendataset(phantom = shepp1fMRI, overlap = 0.01) x <- da$pval[!is.na(da$pval)][1:100] dpval(x) dmixpval(x, mu = da$mu, eta = da$eta)library(MixfMRI, quietly = TRUE) set.seed(1234) da <- gendataset(phantom = shepp1fMRI, overlap = 0.01) x <- da$pval[!is.na(da$pval)][1:100] dpval(x) dmixpval(x, mu = da$mu, eta = da$eta)
Calculates the Expected Euler Characteristic for a 3D Random Field thesholded a level u.
EC.3D(u, sigma, voxdim = c(1, 1, 1), num.vox, type = c("Normal", "t"), df = NULL)EC.3D(u, sigma, voxdim = c(1, 1, 1), num.vox, type = c("Normal", "t"), df = NULL)
u |
The threshold for the field. |
sigma |
The spatial covariance matrix of the field. |
voxdim |
The dimensions of the cuboid 'voxels' upon which the discretized field is observed. |
num.vox |
The number of voxels that make up the field. |
type |
The marginal distribution of the Random Field (only Normal and t at present). |
df |
The degrees of freedom of the t field. |
The Euler Characteristic (Adler, 1981) is a
topological measure that essentially counts the number of isolated
regions of the random field above the threshold minus the
number of 'holes'. As increases the holes disappear and
counts the number of local maxima. So when
becomes close to the maximum of the random field
we have that
where is the null hypothesis that there is no signicant
positive actiavtion/signal present in the field. Thus the Type I error
of the test can be controlled through knowledge of the Expected Euler characteristic.
Note: This function is directly copied from "AnalyzeFMRI".
The value of the expected Euler Characteristic.
J. L. Marchini
Adler, R. (1981) The Geometry of Random Fields.. New York: Wiley.
Worlsey, K. J. (1994) Local maxima and the expected euler
characteristic of excursion sets of , and
fields. Advances in Applied Probability, 26, 13-42.
EC.3D(4.6, sigma = diag(1, 3), voxdim = c(1, 1, 1), num.vox = 10000) EC.3D(4.6, sigma = diag(1, 3), voxdim = c(1, 1, 1), num.vox = 10000, type = "t", df = 100)EC.3D(4.6, sigma = diag(1, 3), voxdim = c(1, 1, 1), num.vox = 10000) EC.3D(4.6, sigma = diag(1, 3), voxdim = c(1, 1, 1), num.vox = 10000, type = "t", df = 100)
These are datasets used to demo examples and workflows in this package.
Objects may contain several information and data.
pstats is a 3D example.
pval.2d.complex and pval.2d.mag are 2D examples.
shepp0fMRI, shepp1fMRI, shepp2fMRI and
sheppAnat are phantoms generated by Dr. Maitra
for simulation studies with different overlap levels for p-values.
toy1 and toy2 are two 3D toy examples.
Wei-Chen Chen and Ranjan Maitra.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
library(MixfMRI, quietly = TRUE) ### Plotting. demo(shepp,'MixfMRI',ask=FALSE,echo=FALSE)library(MixfMRI, quietly = TRUE) ### Plotting. demo(shepp,'MixfMRI',ask=FALSE,echo=FALSE)
Compute q-values Benjamini and Heller's (2007) approach for controlling FDR for spatial signals.
fdr.bh.p1(p, w = rep(1, length(p)), q = 0.05) fdr.bh.p2(p, w = rep(1, length(p)), q = 0.05)fdr.bh.p1(p, w = rep(1, length(p)), q = 0.05) fdr.bh.p2(p, w = rep(1, length(p)), q = 0.05)
p |
a p-value vector. No NA is allowed and all values are in [0, 1]. |
w |
a weight vector for p-values. |
q |
a desired cutoff for adjusting p-values. |
These functions implement first two procedures in Benjamini and Heller (2007) for controlling FDR for spatial signals.
Return the number of rejected hypotheses and all corresponding q-values for the input p-values.
Wei-Chen Chen.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
qvalue().
library(MixfMRI, quietly = TRUE) set.seed(1234) da <- gendataset(phantom = shepp1fMRI, overlap = 0.01) p <- da$pval[!is.na(da$pval)][1:100] fdr.bh.p1(p) fdr.bh.p2(p)library(MixfMRI, quietly = TRUE) set.seed(1234) da <- gendataset(phantom = shepp1fMRI, overlap = 0.01) p <- da$pval[!is.na(da$pval)][1:100] fdr.bh.p1(p) fdr.bh.p2(p)
Find clusters in 2D or 3D based on a generalized CBA method. The CBA method is originally proposed by Heller, et.al. (2006) using the correlation of two time series as the similarity of two spatial locations.
cba.cor(da.ts, da.m = NULL, adj.dist = TRUE, fun.sim = stats::cor) cba.cor.2d(da.ts, da.m = NULL, adj.dist = TRUE, fun.sim = stats::cor) cba.cor.3d(da.ts, da.m = NULL, adj.dist = TRUE, fun.sim = stats::cor)cba.cor(da.ts, da.m = NULL, adj.dist = TRUE, fun.sim = stats::cor) cba.cor.2d(da.ts, da.m = NULL, adj.dist = TRUE, fun.sim = stats::cor) cba.cor.3d(da.ts, da.m = NULL, adj.dist = TRUE, fun.sim = stats::cor)
da.ts |
a time series array of dimensions |
da.m |
a mask determining inside of brain or not. |
adj.dist |
if adjust correlations by distance. |
fun.sim |
a function computing simility of two locations. |
These functions implement the 2D and 3D versions of CBA proposed by Heller, et.al. (2006).
da.ts should have dimensions x * y * z * t for 3D data
and x * y * time for 2D data. Similarly, da.m would have
x * y * z and x * y correspondingly.
da.m has values 0 or 1 indicating outside or inside a brain,
respectively.
fun.sim(a, B) is a function return similarity between a location
a and N neighboring locations B where a is of dimension
t * 1 and B is of dimension t * N.
Ideally, fun.sim()
should return values of similarity which take values
between 0 and 1 where 0 means totally different and
1 means completely identical of two spatial locations.
By default, stats::cor is used.
See the example section next for user defined functions for
fun.sim().
Return the cluster ids for each voxel. NA for outside of brain if
da.m is provided.
Wei-Chen Chen.
Heller, et.al. (2006) “Cluster-based analysis of FMRI data”, NeuroImage, 33(2), 599-608.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
### Simulated data library(MixfMRI, quietly = TRUE) dim <- c(4, 5, 4, 10) set.seed(123) da.ts <- array(rnorm(prod(dim)), dim = dim) id.class <- suppressWarnings(cba.cor(da.ts)) table(id.class) fun.tanh <- function(a, B){ d <- 1 / apply(B, 2, function(b){ dist(rbind(as.vector(a), b)) }) tanh(d) } id.class.tanh <- suppressWarnings(cba.cor(da.ts, fun.sim = fun.tanh)) table(id.class.tanh) fun.logit <- function(a, B){ d <- dist(t(cbind(a, B)))[1:ncol(B)] (1 / (1 + exp(-d))) * 2 - 1 } id.class.logit <- suppressWarnings(cba.cor(da.ts, fun.sim = fun.logit)) table(id.class.logit) .rem <- function(){ ### Real data # library(AnalyzeFMRI, quietly = TRUE) # library(oro.nifti, quietly = TRUE) # fn <- "pb02_volreg_tlrc.nii" # da <- readNIfTI(fn) # da.ts <- [email protected] # fn <- "mask_anat.nii" # da <- readNIfTI(fn) # da.m <- [email protected] # id.class <- suppressWarnings(cba.cor(da.ts, da.m)) # dim(id.class) <- dim(da.m) # length(table(id.class)) }### Simulated data library(MixfMRI, quietly = TRUE) dim <- c(4, 5, 4, 10) set.seed(123) da.ts <- array(rnorm(prod(dim)), dim = dim) id.class <- suppressWarnings(cba.cor(da.ts)) table(id.class) fun.tanh <- function(a, B){ d <- 1 / apply(B, 2, function(b){ dist(rbind(as.vector(a), b)) }) tanh(d) } id.class.tanh <- suppressWarnings(cba.cor(da.ts, fun.sim = fun.tanh)) table(id.class.tanh) fun.logit <- function(a, B){ d <- dist(t(cbind(a, B)))[1:ncol(B)] (1 / (1 + exp(-d))) * 2 - 1 } id.class.logit <- suppressWarnings(cba.cor(da.ts, fun.sim = fun.logit)) table(id.class.logit) .rem <- function(){ ### Real data # library(AnalyzeFMRI, quietly = TRUE) # library(oro.nifti, quietly = TRUE) # fn <- "pb02_volreg_tlrc.nii" # da <- readNIfTI(fn) # da.ts <- [email protected] # fn <- "mask_anat.nii" # da <- readNIfTI(fn) # da.m <- [email protected] # id.class <- suppressWarnings(cba.cor(da.ts, da.m)) # dim(id.class) <- dim(da.m) # length(table(id.class)) }
Main initialization functions.
initial.em.gbd(PARAM) initial.RndEM.gbd(PARAM)initial.em.gbd(PARAM) initial.RndEM.gbd(PARAM)
PARAM |
a list of uninitialized parameters, as usual, the returned
values of |
initial.em.gbd() takes in a template of PARAM (uninitialized),
and usually is available by calling set.global(), then
return an initialized PARAM which is ready for EM runs.
Internally, there are six different initializations implemented for
the function initial.em.gbd() including prob.extend,
prob.simple, qnorm.extend, qnorm.simple, extend,
and simple. These methods are mainly based on transformation of
original space of data (p-values and voxel locations) into more linear
space such that the Euclidean distance more makes sense (fairly) to
classify data in groups.
initial.RndEM.gbd() implement RndEM initialization algorithm
based on repeated calling initial.em.gbd().
Note that all configurations are included in PARAM set by
set.global().
These functions return an initialized PARAM for EM runs based on
pre-stored configuration within the input uninitialized PARAM.
Wei-Chen Chen and Ranjan Maitra.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
set.global(), fclust(), PARAM.
library(MixfMRI, quietly = TRUE) library(EMCluster, quietly = TRUE) # .FC.CT$algorithm <- "em" # .FC.CT$model.X <- "V" # .FC.CT$ignore.X <- TRUE .FC.CT$check.X.unit <- FALSE ### Test toy1. set.seed(1234) X.gbd <- toy1$X.gbd PV.gbd <- toy1$PV.gbd PARAM <- set.global(X.gbd, PV.gbd, K = 2) PARAM.new <- initial.em.gbd(PARAM) PARAM.toy1 <- em.step.gbd(PARAM.new) id.toy1 <- .MixfMRIEnv$CLASS.gbd print(PARAM.toy1$ETA) RRand(toy1$CLASS.gbd, id.toy1) .rem <- function(){ ### Test toy2. set.seed(1234) X.gbd <- toy2$X.gbd PV.gbd <- toy2$PV.gbd PARAM <- set.global(X.gbd, PV.gbd, K = 3) PARAM.new <- initial.em.gbd(PARAM) PARAM.toy2 <- em.step.gbd(PARAM.new) id.toy2 <- .MixfMRIEnv$CLASS.gbd print(PARAM.toy2$ETA) RRand(toy2$CLASS.gbd, id.toy2) }library(MixfMRI, quietly = TRUE) library(EMCluster, quietly = TRUE) # .FC.CT$algorithm <- "em" # .FC.CT$model.X <- "V" # .FC.CT$ignore.X <- TRUE .FC.CT$check.X.unit <- FALSE ### Test toy1. set.seed(1234) X.gbd <- toy1$X.gbd PV.gbd <- toy1$PV.gbd PARAM <- set.global(X.gbd, PV.gbd, K = 2) PARAM.new <- initial.em.gbd(PARAM) PARAM.toy1 <- em.step.gbd(PARAM.new) id.toy1 <- .MixfMRIEnv$CLASS.gbd print(PARAM.toy1$ETA) RRand(toy1$CLASS.gbd, id.toy1) .rem <- function(){ ### Test toy2. set.seed(1234) X.gbd <- toy2$X.gbd PV.gbd <- toy2$PV.gbd PARAM <- set.global(X.gbd, PV.gbd, K = 3) PARAM.new <- initial.em.gbd(PARAM) PARAM.toy2 <- em.step.gbd(PARAM.new) id.toy2 <- .MixfMRIEnv$CLASS.gbd print(PARAM.toy2$ETA) RRand(toy2$CLASS.gbd, id.toy2) }
These functions test two mixture Gaussian fMRI models with diagonal
covariance matrices and different numbers of clusters.
These functions are similar to the EMCluster::lmt(), but is coded
for fMRI models in MixfMRI.
lmt.I(fcobj.0, fcobj.a, X.gbd, PV.gbd, tau = 0.5, n.mc.E.delta = 1000, n.mc.E.chi2 = 1000, verbose = FALSE) lmt.pv(fcobj.0, fcobj.a, X.gbd, PV.gbd, tau = 0.5, n.mc.E.delta = 1000, n.mc.E.chi2 = 1000, verbose = FALSE)lmt.I(fcobj.0, fcobj.a, X.gbd, PV.gbd, tau = 0.5, n.mc.E.delta = 1000, n.mc.E.chi2 = 1000, verbose = FALSE) lmt.pv(fcobj.0, fcobj.a, X.gbd, PV.gbd, tau = 0.5, n.mc.E.delta = 1000, n.mc.E.chi2 = 1000, verbose = FALSE)
fcobj.0 |
a |
fcobj.a |
a |
X.gbd |
a data matrix of |
PV.gbd |
a p-value vector of signals associated with voxels.
|
tau |
proportion of null and alternative hypotheses. |
n.mc.E.delta |
number of Monte Carlo simulations for expectation of delta (difference of logL). |
n.mc.E.chi2 |
number of Monte Carlo simulations for expectation of chisquare statistics. |
verbose |
if verbose. |
This function calls several subroutines to compute information,
likelihood ratio statistics, degrees of freedom, non-centrality
of chi-squared distributions ... etc. Based on Monte Carlo methods
to estimate parameters of likelihood mixture tests, this function
return a p-value for testing H0: fcobj.0 v.s. Ha: fcobj.a.
lmt.pv() only uses PV.gbd.
A list of class lmt.I are returned.
Wei-Chen Chen and Ranjan Maitra.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
EMCluster::lmt().
library(MixfMRI, quietly = TRUE) library(EMCluster, quietly = TRUE) .FC.CT$model.X <- "I" .FC.CT$check.X.unit <- FALSE .FC.CT$CONTROL$debug <- 0 .rem <- function(){ ### Fit toy1. set.seed(1234) X.gbd <- toy1$X.gbd PV.gbd <- toy1$PV.gbd ret.2 <- fclust(X.gbd, PV.gbd, K = 2) ret.3 <- fclust(X.gbd, PV.gbd, K = 3) ret.4 <- fclust(X.gbd, PV.gbd, K = 4) ret.5 <- fclust(X.gbd, PV.gbd, K = 5) ### ARI RRand(toy1$CLASS.gbd, ret.2$class) RRand(toy1$CLASS.gbd, ret.3$class) RRand(toy1$CLASS.gbd, ret.4$class) RRand(toy1$CLASS.gbd, ret.5$class) ### Test toy1. (lmt.23 <- lmt.I(ret.2, ret.3, X.gbd, PV.gbd)) (lmt.24 <- lmt.I(ret.2, ret.4, X.gbd, PV.gbd)) (lmt.25 <- lmt.I(ret.2, ret.5, X.gbd, PV.gbd)) (lmt.34 <- lmt.I(ret.3, ret.4, X.gbd, PV.gbd)) (lmt.35 <- lmt.I(ret.3, ret.5, X.gbd, PV.gbd)) (lmt.45 <- lmt.I(ret.4, ret.5, X.gbd, PV.gbd)) ### Test toy1 using p-values only. (lmt.pv.23 <- lmt.pv(ret.2, ret.3, X.gbd, PV.gbd)) (lmt.pv.24 <- lmt.pv(ret.2, ret.4, X.gbd, PV.gbd)) (lmt.pv.25 <- lmt.pv(ret.2, ret.5, X.gbd, PV.gbd)) (lmt.pv.34 <- lmt.pv(ret.3, ret.4, X.gbd, PV.gbd)) (lmt.pv.35 <- lmt.pv(ret.3, ret.5, X.gbd, PV.gbd)) (lmt.pv.45 <- lmt.pv(ret.4, ret.5, X.gbd, PV.gbd)) }library(MixfMRI, quietly = TRUE) library(EMCluster, quietly = TRUE) .FC.CT$model.X <- "I" .FC.CT$check.X.unit <- FALSE .FC.CT$CONTROL$debug <- 0 .rem <- function(){ ### Fit toy1. set.seed(1234) X.gbd <- toy1$X.gbd PV.gbd <- toy1$PV.gbd ret.2 <- fclust(X.gbd, PV.gbd, K = 2) ret.3 <- fclust(X.gbd, PV.gbd, K = 3) ret.4 <- fclust(X.gbd, PV.gbd, K = 4) ret.5 <- fclust(X.gbd, PV.gbd, K = 5) ### ARI RRand(toy1$CLASS.gbd, ret.2$class) RRand(toy1$CLASS.gbd, ret.3$class) RRand(toy1$CLASS.gbd, ret.4$class) RRand(toy1$CLASS.gbd, ret.5$class) ### Test toy1. (lmt.23 <- lmt.I(ret.2, ret.3, X.gbd, PV.gbd)) (lmt.24 <- lmt.I(ret.2, ret.4, X.gbd, PV.gbd)) (lmt.25 <- lmt.I(ret.2, ret.5, X.gbd, PV.gbd)) (lmt.34 <- lmt.I(ret.3, ret.4, X.gbd, PV.gbd)) (lmt.35 <- lmt.I(ret.3, ret.5, X.gbd, PV.gbd)) (lmt.45 <- lmt.I(ret.4, ret.5, X.gbd, PV.gbd)) ### Test toy1 using p-values only. (lmt.pv.23 <- lmt.pv(ret.2, ret.3, X.gbd, PV.gbd)) (lmt.pv.24 <- lmt.pv(ret.2, ret.4, X.gbd, PV.gbd)) (lmt.pv.25 <- lmt.pv(ret.2, ret.5, X.gbd, PV.gbd)) (lmt.pv.34 <- lmt.pv(ret.3, ret.4, X.gbd, PV.gbd)) (lmt.pv.35 <- lmt.pv(ret.3, ret.5, X.gbd, PV.gbd)) (lmt.pv.45 <- lmt.pv(ret.4, ret.5, X.gbd, PV.gbd)) }
Likelihood ratio tests for merging clusters.
lrt(PV.gbd, CLASS.gbd, K, H0.alpha = .FC.CT$LRT$H0.alpha, H0.beta = .FC.CT$LRT$H0.beta) lrt2(PV.gbd, CLASS.gbd, K, H0.mean = .FC.CT$LRT$H0.mean, upper.beta = .FC.CT$INIT$BETA.beta.max, proc = c("1", "2", "weight")) lrt.betamean(PV.gbd, CLASS.gbd, K, proc = c("1", "2")) lrt.betaab(PV.gbd, CLASS.gbd, K, proc = c("1", "2"))lrt(PV.gbd, CLASS.gbd, K, H0.alpha = .FC.CT$LRT$H0.alpha, H0.beta = .FC.CT$LRT$H0.beta) lrt2(PV.gbd, CLASS.gbd, K, H0.mean = .FC.CT$LRT$H0.mean, upper.beta = .FC.CT$INIT$BETA.beta.max, proc = c("1", "2", "weight")) lrt.betamean(PV.gbd, CLASS.gbd, K, proc = c("1", "2")) lrt.betaab(PV.gbd, CLASS.gbd, K, proc = c("1", "2"))
PV.gbd |
a p-value vector of signals associated with voxels.
|
CLASS.gbd |
a classification vector of signals associated with voxels.
|
K |
number of clusters. |
H0.alpha |
null hypothesis for the alpha parameter of Beta distribution. |
H0.beta |
null hypothesis for the beta parameter of Beta distribution. |
H0.mean |
null hypothesis for the mean of Beta distribution. |
upper.beta |
BETA.beta.max, maximum value of beta parameter of Beta distribution. |
proc |
q-value procedure for adjusting p-values. |
These functions perform likelihood ratio tests for merging clusters. Only p-values coordinates (Beta density) are tested, while voxel location coordinates (multivariate Normal density) are not involved in testing.
lrt.betamean tests if means of any two pairs of mixture
(p-value) component were the same.
The chi-square distribution with 1 degree of freedom is used.
lrt.betaab tests if alpha and beta of any two pairs of mixture
(p-value) components were the same.
The chi-square distribution with 2 degrees of freedom is used.
Procedure to adjust/select plausible p-values,
proc = "1" uses q-value qvalue(),
proc = "2" uses fdr.bh.p2(), and
proc = "weight" uses a weighted version of fdr.bh.p2().
A matrix contains MLEs of parameters of Beta distribution under the null hypothesis and the union of null and alternative hypotheses. The matrix also contains testing statistics and p-values.
Wei-Chen Chen and Ranjan Maitra.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
library(MixfMRI, quietly = TRUE) set.seed(1234) ### Test 2d data. da <- pval.2d.mag id <- !is.na(da) PV.gbd <- da[id] id.loc <- which(id, arr.ind = TRUE) X.gbd <- t(t(id.loc) / dim(da)) ret <- fclust(X.gbd, PV.gbd, K = 2, min.1st.prop = 0.95) # print(ret) ### p-values of rest clusters. ret.lrt <- lrt(PV.gbd, ret$class, K = 2) print(ret.lrt) .rem <- function(){ ret.lrt2 <- lrt2(PV.gbd, ret$class, K = 3) print(ret.lrt2) }library(MixfMRI, quietly = TRUE) set.seed(1234) ### Test 2d data. da <- pval.2d.mag id <- !is.na(da) PV.gbd <- da[id] id.loc <- which(id, arr.ind = TRUE) X.gbd <- t(t(id.loc) / dim(da)) ret <- fclust(X.gbd, PV.gbd, K = 2, min.1st.prop = 0.95) # print(ret) ### p-values of rest clusters. ret.lrt <- lrt(PV.gbd, ret$class, K = 2) print(ret.lrt) .rem <- function(){ ret.lrt2 <- lrt2(PV.gbd, ret$class, K = 3) print(ret.lrt2) }
Main MixfMRI functions.
fclust(X.gbd, PV.gbd, K = 2, PARAM.init = NULL, min.1st.prop = .FC.CT$INIT$min.1st.prop, max.PV = .FC.CT$INIT$max.PV, class.method = .FC.CT$INIT$class.method[1], RndEM.iter = .FC.CT$CONTROL$RndEM.iter, algorithm = .FC.CT$algorithm[1], model.X = .FC.CT$model.X[1], ignore.X = .FC.CT$ignore.X, stop.unstable = TRUE, MPI.gbd = .FC.CT$MPI.gbd, common.gbd = .FC.CT$common.gbd) set.global(X.gbd, PV.gbd, K = 2, min.1st.prop = .FC.CT$INIT$min.1st.prop, max.PV = .FC.CT$INIT$max.PV, class.method = .FC.CT$INIT$class.method[1], RndEM.iter = .FC.CT$CONTROL$RndEM.iter, algorithm = .FC.CT$algorithm[1], model.X = .FC.CT$model.X[1], ignore.X = .FC.CT$ignore.X, check.X.unit = .FC.CT$check.X.unit, MPI.gbd = .FC.CT$MPI.gbd, common.gbd = .FC.CT$common.gbd)fclust(X.gbd, PV.gbd, K = 2, PARAM.init = NULL, min.1st.prop = .FC.CT$INIT$min.1st.prop, max.PV = .FC.CT$INIT$max.PV, class.method = .FC.CT$INIT$class.method[1], RndEM.iter = .FC.CT$CONTROL$RndEM.iter, algorithm = .FC.CT$algorithm[1], model.X = .FC.CT$model.X[1], ignore.X = .FC.CT$ignore.X, stop.unstable = TRUE, MPI.gbd = .FC.CT$MPI.gbd, common.gbd = .FC.CT$common.gbd) set.global(X.gbd, PV.gbd, K = 2, min.1st.prop = .FC.CT$INIT$min.1st.prop, max.PV = .FC.CT$INIT$max.PV, class.method = .FC.CT$INIT$class.method[1], RndEM.iter = .FC.CT$CONTROL$RndEM.iter, algorithm = .FC.CT$algorithm[1], model.X = .FC.CT$model.X[1], ignore.X = .FC.CT$ignore.X, check.X.unit = .FC.CT$check.X.unit, MPI.gbd = .FC.CT$MPI.gbd, common.gbd = .FC.CT$common.gbd)
X.gbd |
a data matrix of |
PV.gbd |
a p-value vector of signals associated with voxels.
|
K |
number of clusters to be estimated. |
PARAM.init |
initial parameters. |
min.1st.prop |
lower bound of mixing proportion (ETA) of the 1st cluster (uniform). |
max.PV |
upper bound of p-values where initializations pick from. |
class.method |
classification method for initializations. |
RndEM.iter |
number of RndEM iterations. |
algorithm |
either “ecm” (ECM), “apecma” (APECMa) or “em” (EM) algorithm. |
model.X |
either “I” or “V” for covariance matrix. |
ignore.X |
if |
check.X.unit |
if |
stop.unstable |
if |
MPI.gbd |
if MPI (“EGM” algorithm) is used. |
common.gbd |
if |
The fclust() contains initialization and EM algorithms for clustering
fMRI signal data which have two parts: X.gbd for voxel information
either 2D or 3D, PV.gbd for p-value of signals associated with
voxels. Each signal is assumed as a mixture distribution with K
components with mixing proportion ETA, and each component has
two independent coordinates with density functions: Beta and multivariate
Normal distributions.
Beta density:
The 1st component is restricted by min.1st.prop and Beta(1, 1)
distribution. The other K - 1 components have Beta(alpha, beta)
distribution with alpha < 1 < beta.
Multivariate Normal density:
model.X = "I" is for diagonal cov matrix of multivariate Normal
distribution, and "V" for unstructured cov matrix.
ignore.X = TRUE is to ignore X.gbd and normal density,
i.e. only Beta density is used.
Currently, APECMa and EM algorithms are implemented with EGM algorithm to speed up convergence if MPI is available. RndEM initialization is also implemented for better chance of good initial values for convergence.
The set.global() has purposes: create a template/storage of
parameters, save configurations, and called by fclust() to initial
the parameters, such as initial.em.gbd() or
initial.RndEM.gbd().
A list with class fclust by fclust() is returned
which can be summarized by print.fclust().
A list PARAM or PARAM.org is returned by set.global():
N.gbd |
number of observations (within the rank), and should be
equal to |
N.all |
numbers of observations (of all ranks
if |
N |
total number of observations ( |
p |
dimension of an observation (3 for 2D signals, 4 for 3D signals), equivalent to total number of coordinates. |
p.X |
dimension of |
K |
number of clusters. |
ETA |
mixing proportion, length |
log.ETA |
|
BETA |
a list of length |
MU |
a matrix of dimension |
SIGMA |
a list of length |
logL |
log likelihood value. |
min.1st.prop |
carried from input. |
max.PV |
carried from input. |
class.method |
classification method of initializations. |
min.N.CLASS |
|
model.X |
carried from input. |
Wei-Chen Chen and Ranjan Maitra.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
print.fclust().
library(MixfMRI, quietly = TRUE) library(EMCluster, quietly = TRUE) # .FC.CT$algorithm <- "em" # .FC.CT$model.X <- "V" # .FC.CT$ignore.X <- TRUE .FC.CT$check.X.unit <- FALSE set.seed(1234) ### Test toy1. X.gbd <- toy1$X.gbd[, -3] PV.gbd <- toy1$PV.gbd PARAM <- fclust(X.gbd, PV.gbd, K = 2) print(PARAM) id.toy1 <- .MixfMRIEnv$CLASS.gbd print(RRand(toy1$CLASS.gbd, id.toy1)) .rem <- function(){ ### Test toy2. X.gbd <- toy2$X.gbd[, -3] PV.gbd <- toy2$PV.gbd PARAM <- fclust(X.gbd, PV.gbd, K = 3) print(PARAM) id.toy2 <- .MixfMRIEnv$CLASS.gbd print(RRand(toy2$CLASS.gbd, id.toy2)) }library(MixfMRI, quietly = TRUE) library(EMCluster, quietly = TRUE) # .FC.CT$algorithm <- "em" # .FC.CT$model.X <- "V" # .FC.CT$ignore.X <- TRUE .FC.CT$check.X.unit <- FALSE set.seed(1234) ### Test toy1. X.gbd <- toy1$X.gbd[, -3] PV.gbd <- toy1$PV.gbd PARAM <- fclust(X.gbd, PV.gbd, K = 2) print(PARAM) id.toy1 <- .MixfMRIEnv$CLASS.gbd print(RRand(toy1$CLASS.gbd, id.toy1)) .rem <- function(){ ### Test toy2. X.gbd <- toy2$X.gbd[, -3] PV.gbd <- toy2$PV.gbd PARAM <- fclust(X.gbd, PV.gbd, K = 3) print(PARAM) id.toy2 <- .MixfMRIEnv$CLASS.gbd print(RRand(toy2$CLASS.gbd, id.toy2)) }
These sets of controls are used to provide default values in this package.
Objects contain several parameters for methods.
The elements of .FC.CT are default values for main controls of
MixfMRI including
| Elements | Default | Usage |
algorithm |
"apecma" | implemented algorithm |
optim.method |
"BFGS" | optimization method |
model.X |
"I" | cov matrix structure |
ignore.X |
FALSE | if using voxel information |
check.X.unit |
TRUE | if checking X in [0, 1] |
CONTROL |
a list | see CONTROL next for details |
INIT |
a list | see INIT next for details |
LRT |
a list | see LRT next for details |
MPI.gbd |
FALSE | if MPI speedup available |
common.gbd |
TRUE | if X in common gbd format |
The elements of CONTROL are default values for optimization controls
of implemented EM algorithm including
| Elements | Default | Usage |
max.iter |
1000 | maximum number of EM iterations |
abs.err |
1e-4 | absolute error of convergence |
rel.err |
1e-6 | relative error of convergence |
debug |
1 | debugging level |
RndEM.iter |
10 | RndEM iterations |
exp.min |
log(.Machine$double.xmin) | minimum exponential power |
exp.max |
log(.Machine$double.xmax) | maximum exponential power |
sigma.ill |
1e-6 | ill condition limit |
DS.max |
1e+4 | maximum chol() cov matrix |
DS.min |
1e-6 | minimum chol() cov matrix |
The elements of INIT are default values or limitations for initial
parameters implemented for EM algorithm including
| Elements | Default | Usage |
min.1st.prop |
0.8 | minimum proportion of 1st cluster |
max.PV |
0.1 | maximum p-value for initialization |
BETA.alpha.min |
0 + 1e-6 | minimum value of alpha parameter of Beta distribution |
BETA.alpha.max |
1 - 1e-6 | maximum value of alpha parameter of Beta distribution |
BETA.beta.min |
1 + 1e-6 | minimum value of beta parameter of Beta distribution |
BETA.beta.max |
1e+6 | maximum value of beta parameter of Beta distribution |
max.try.iter |
10 | maximum retry iterations if result is unstable |
class.method |
"prob.extned" | classification method at initializations |
The elements of LRT are default values or limitations for likelihood
ratio tests including
| Elements | Default | Usage |
H0.alpha |
1 | null hypothesis alpha parameter of Beta distribution |
H0.beta |
1 | null hypothesis beta parameter of Beta distribution |
H0.mean |
0.05 | null hypothesis mean of Beta distribution |
Wei-Chen Chen and Ranjan Maitra.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
set.global(), fclust().
Main plotting function in MixfMRI.
plotfclust(da, posterior, main = NULL, xlim = NULL, ylim = NULL) plotfclustpv(da, posterior, main = NULL, xlim = NULL, ylim = NULL) plotpv(da, posterior, PARAM, zlim = c(0, 0.01), plot.mean = TRUE, xlab = "", ylab = "", main = NULL, xlim = NULL, ylim = NULL, col = my.YlOrRd(), ignore.bg = FALSE) plotpvlegend(zlim = c(0, 0.01), n.level = 20, main = NULL, col = my.YlOrRd())plotfclust(da, posterior, main = NULL, xlim = NULL, ylim = NULL) plotfclustpv(da, posterior, main = NULL, xlim = NULL, ylim = NULL) plotpv(da, posterior, PARAM, zlim = c(0, 0.01), plot.mean = TRUE, xlab = "", ylab = "", main = NULL, xlim = NULL, ylim = NULL, col = my.YlOrRd(), ignore.bg = FALSE) plotpvlegend(zlim = c(0, 0.01), n.level = 20, main = NULL, col = my.YlOrRd())
da |
a data set to be plotted. |
posterior |
a posterior data set to be plotted. |
PARAM |
a returning parameter object from |
main |
title of the plot. |
xlim |
limits of x-axis. |
ylim |
limits of y-axis. |
zlim |
limits of z-axis. |
xlab |
labels of x-axis. |
ylab |
labels of y-axis. |
plot.mean |
if plotting mean values of each cluster. |
col |
colors to be drawn. |
ignore.bg |
if ignoring the background. |
n.level |
number of levels to be plotted. |
These are example functions to plot results, simulations, and datasets.
Return plots.
Wei-Chen Chen and Ranjan Maitra.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
set.global().
library(MixfMRI, quietly = TRUE) set.seed(1234) .rem <- function(){ ### Check 2d data. da <- pval.2d.complex id <- !is.na(da) PV.gbd <- da[id] hist(PV.gbd, nclass = 100, main = "p-value") ### Test 2d data. id.loc <- which(id, arr.ind = TRUE) X.gbd <- t(t(id.loc) / dim(da)) ret <- fclust(X.gbd, PV.gbd, K = 3) print(ret) ### p-values of rest clusters. ret.lrt <- lrt(PV.gbd, ret$class, K = 3) print(ret.lrt) ret.lrt2 <- lrt2(PV.gbd, ret$class, K = 3) print(ret.lrt2) ### Plotting. par(mfrow = c(2, 2), mar = c(0, 0, 2, 0)) plotpv(da, ret$posterior, ret$param, zlim = c(0.005, 0.008), main = "Mean of Beta Distribution") plotpv(da, ret$posterior, ret$param, plot.mean = FALSE, main = "p-value") par(mar = c(5.1, 4.1, 4.1, 2.1)) plotpvlegend(zlim = c(0.005, 0.008), main = "Mean of Beta Distribution") plotpvlegend(zlim = c(0, 0.01), main = "p-value") }library(MixfMRI, quietly = TRUE) set.seed(1234) .rem <- function(){ ### Check 2d data. da <- pval.2d.complex id <- !is.na(da) PV.gbd <- da[id] hist(PV.gbd, nclass = 100, main = "p-value") ### Test 2d data. id.loc <- which(id, arr.ind = TRUE) X.gbd <- t(t(id.loc) / dim(da)) ret <- fclust(X.gbd, PV.gbd, K = 3) print(ret) ### p-values of rest clusters. ret.lrt <- lrt(PV.gbd, ret$class, K = 3) print(ret.lrt) ret.lrt2 <- lrt2(PV.gbd, ret$class, K = 3) print(ret.lrt2) ### Plotting. par(mfrow = c(2, 2), mar = c(0, 0, 2, 0)) plotpv(da, ret$posterior, ret$param, zlim = c(0.005, 0.008), main = "Mean of Beta Distribution") plotpv(da, ret$posterior, ret$param, plot.mean = FALSE, main = "p-value") par(mar = c(5.1, 4.1, 4.1, 2.1)) plotpvlegend(zlim = c(0.005, 0.008), main = "Mean of Beta Distribution") plotpvlegend(zlim = c(0, 0.01), main = "p-value") }
Print flcust related outputs.
## S3 method for class 'fclust' print(x, ...)## S3 method for class 'fclust' print(x, ...)
x |
an object with the class attributes. |
... |
other arguments to the |
x is the return result from fclust().
A summary of fclust object is printed.
Wei-Chen Chen and Ranjan Maitra.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
set.global(), fclust().
library(MixfMRI, quietly = TRUE) set.seed(1234) ### Check 2d data. da <- pval.2d.complex id <- !is.na(da) PV.gbd <- da[id] # hist(PV.gbd, nclass = 100, main = "p-value") ### Test 2d data. id.loc <- which(id, arr.ind = TRUE) X.gbd <- t(t(id.loc) / dim(da)) ret <- fclust(X.gbd, PV.gbd, K = 2) print(ret)library(MixfMRI, quietly = TRUE) set.seed(1234) ### Check 2d data. da <- pval.2d.complex id <- !is.na(da) PV.gbd <- da[id] # hist(PV.gbd, nclass = 100, main = "p-value") ### Test 2d data. id.loc <- which(id, arr.ind = TRUE) X.gbd <- t(t(id.loc) / dim(da)) ret <- fclust(X.gbd, PV.gbd, K = 2) print(ret)
Generate datasets for MixfMRI simulations
gendataset(phantom, overlap, smooth = FALSE)gendataset(phantom, overlap, smooth = FALSE)
phantom |
a phantom dataset. |
overlap |
a desired overlap level. |
smooth |
if |
This is a function to generate simulated fMRI data based on the input
phantom and the desired overlap level for the fMRI p-value.
Return a list contains eta for mixing proportion,
overlap for the desired level, mu for center of p-values,
class.id for the true classifications where p-values belong to,
tval for the testing statistics, and pval for the p-values
of interesting in simulations.
Wei-Chen Chen and Ranjan Maitra.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
set.global().
.rem <- function(){ library(MixfMRI, quietly = TRUE) set.seed(1234) da <- gendataset(phantom = shepp1fMRI, overlap = 0.01)$pval da2 <- gendataset(phantom = shepp2fMRI, overlap = 0.01)$pval par(mfrow = c(2, 2), mar = rep(0.05, 4)) image(shepp1fMRI[50:210, 50:210], axes = FALSE) image(shepp2fMRI[50:210, 50:210], axes = FALSE) image(da[50:210, 50:210], axes = FALSE) image(da2[50:210, 50:210], axes = FALSE) }.rem <- function(){ library(MixfMRI, quietly = TRUE) set.seed(1234) da <- gendataset(phantom = shepp1fMRI, overlap = 0.01)$pval da2 <- gendataset(phantom = shepp2fMRI, overlap = 0.01)$pval par(mfrow = c(2, 2), mar = rep(0.05, 4)) image(shepp1fMRI[50:210, 50:210], axes = FALSE) image(shepp2fMRI[50:210, 50:210], axes = FALSE) image(da[50:210, 50:210], axes = FALSE) image(da2[50:210, 50:210], axes = FALSE) }
Generate datasets with smoothing for MixfMRI simulations
gcv.smooth2d(y, interval)gcv.smooth2d(y, interval)
y |
a set of p-values in 2d phantom |
interval |
an interval for |
The function is used to smooth for Dr. Maitra's 2d phantom simulation. The smoothing method is based on Garcia (2010), CSDA.
Return a list containing two elements im.smooth and
par.val.
Ranjan Maitra.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
Compute summarized overlap on a given overlap (symmetric) matrix.
summarized.overlap(overlap.mat)summarized.overlap(overlap.mat)
overlap.mat |
an overlap (symmetric) matrix. |
overlap.mat is a p * p matrix containing pair wised overlaps
of p experiments. overlap.mat is assumed a symmetric matrix.
This function returns a summarized overlap based on the input
overlap.mat that charactorizes the overall overlap behavior of the
p experiments.
A single value is returned.
Ranjan Maitra.
Chen, W.-C. and Maitra, R. (2021) “A Practical Model-based Segmentation Approach for Accurate Activation Detection in Single-Subject functional Magnetic Resonance Imaging Studies”, arXiv:2102.03639.
library(MixfMRI, quietly = TRUE) set.seed(1234) p <- 10 # 10 experiments. overlap.mat <- diag(1, p) overlap.mat[lower.tri(overlap.mat)] <- runif(p * (p - 1) / 2) overlap.mat[upper.tri(overlap.mat)] <- t(overlap.mat)[upper.tri(overlap.mat)] summarized.overlap(overlap.mat)library(MixfMRI, quietly = TRUE) set.seed(1234) p <- 10 # 10 experiments. overlap.mat <- diag(1, p) overlap.mat[lower.tri(overlap.mat)] <- runif(p * (p - 1) / 2) overlap.mat[upper.tri(overlap.mat)] <- t(overlap.mat)[upper.tri(overlap.mat)] summarized.overlap(overlap.mat)
Calculate the Bonferroni threshold for n iid tests that results in an overall p-value of p.val. The tests can be distributed as Normal, t or F.
Threshold.Bonferroni(p.val, n, type = c("Normal", "t", "F"), df1 = NULL, df2 = NULL)Threshold.Bonferroni(p.val, n, type = c("Normal", "t", "F"), df1 = NULL, df2 = NULL)
p.val |
The required overall p-value. |
n |
The number of tests. |
type |
The distribution of the tests. One of "Normal", "t" or "F" |
df1 |
The degrees of freedom of the t-distribution or the first degrees of freedom parameter for the F distribution. |
df2 |
The second degrees of freedom parameter for the F distribution. |
Note: This function is directly copied from "AnalyzeFMRI".
Returns the Bonferroni threshold.
Pierre Lafaye De Micheaux and J. L. Marchini.
Threshold.Bonferroni(0.05, 1000) Threshold.Bonferroni(0.05, 1000, type = c("t"), df1 = 20) Threshold.Bonferroni(0.05, 1000, type = c("F"), df1 = 3, df2 = 100)Threshold.Bonferroni(0.05, 1000) Threshold.Bonferroni(0.05, 1000, type = c("t"), df1 = 20) Threshold.Bonferroni(0.05, 1000, type = c("F"), df1 = 3, df2 = 100)
Calculates the False Discovery Rate (FDR) threshold for a given vector of statistic values.
Threshold.FDR(x, q, cV.type = 2, type = c("Normal", "t", "F"), df1 = NULL, df2 = NULL)Threshold.FDR(x, q, cV.type = 2, type = c("Normal", "t", "F"), df1 = NULL, df2 = NULL)
x |
A vector of test statistic values. |
q |
The desired False Discovery Rate threshold. |
cV.type |
A flag that specifies the assumptions about the joint distribution of p-values. Choose cV.type = 2 for fMRI data (see Genovese et al (2001) |
type |
The distribution of the statistic values. Either "Normal", "t" or "F". |
df1 |
The degrees of freedom of the t-distribution or the first degrees of freedom parameter for the F distribution. |
df2 |
The second degrees of freedom parameter for the F distribution. |
Note: This function is directly copied from "AnalyzeFMRI".
Returns the FDR threshold.
J. L. Marchini
Genovese et al. (2001) Thresholding of Statistical Maps in Functional NeuroImaging Using the False Discovery Rate.
x <- c(rnorm(1000), rnorm(100, mean = 3)) Threshold.FDR(x = x, q = 0.20, cV.type = 2)x <- c(rnorm(1000), rnorm(100, mean = 3)) Threshold.FDR(x = x, q = 0.20, cV.type = 2)
Calculates the Random Field theory threshold to give that results in a specified p-value.
Threshold.RF(p.val, sigma, voxdim = c(1, 1, 1), num.vox, type = c("Normal", "t"), df = NULL)Threshold.RF(p.val, sigma, voxdim = c(1, 1, 1), num.vox, type = c("Normal", "t"), df = NULL)
p.val |
The required p-value. |
sigma |
The 3D covariance matrix of the random field. |
voxdim |
The dimesnions of a voxel. |
num.vox |
The number of voxels that constitute the random field. |
type |
The type of random field, "Normal" or "t". |
df |
The degrees of the t distributed field. |
Calculates the threshold that produces an expected Euler characteristic equal to the required p-value.
Note: This function is directly copied from "AnalyzeFMRI".
Returns the Random Field threshold.
J. L. Marchini
a <- Threshold.RF(p.val = 0.05, sigma = diag(1, 3), voxdim = c(1, 1, 1), num.vox = 10000) EC.3D(a, sigma = diag(1, 3), voxdim = c(1, 1, 1), num.vox = 10000)a <- Threshold.RF(p.val = 0.05, sigma = diag(1, 3), voxdim = c(1, 1, 1), num.vox = 10000) EC.3D(a, sigma = diag(1, 3), voxdim = c(1, 1, 1), num.vox = 10000)