sanssouci-org / sanssouci.python

Post hoc inference via multiple testing
GNU General Public License v3.0
6 stars 3 forks source link

Example consistent with literature #3

Open bthirion opened 3 years ago

bthirion commented 3 years ago

I suggest that the examples should:

pneuvial commented 3 years ago
* It is a bit frustrating that the example shows null data (no effect)

Indeed. TODO: a simple case with Gaussian test statistics, under equi-correlation.

pneuvial commented 3 years ago

R code to do this:

p <- 50000         # number of voxels
pi0 <- 0.99        # proportion of noise (true null hypotheses)
p0 <- round(pi0*p) # number of true null hypotheses (non active voxels)
n <- 100           # number of subjects
rho <- 0.2         # equi-correlation
SNR <- 2           # signal to noise ratio

# generate noise
Z <- matrix(rnorm(n*p), ncol = p) # p independent N_n(0,1)
w <- rnorm(n)                     # common factor N_n(0,1)
W <- matrix(w, ncol = p, nrow = n, byrow = FALSE)
Y <- sqrt(1-rho)*Z + sqrt(rho)*W

# sanity check: correlations between pairs of variables
co <- cor(Y[, 1:10])
boxplot(co[upper.tri(co)])
mean(co[upper.tri(co)])

# add signal (for false null hypotheses)
categ <- rbinom(n, 1, prob = 0.5) # two balanced categories
X <- Y
colnames(X) <- rep(c("H0", "H1"), times = c(p0, p-p0))
w1 <- which(categ==1)
X[w1, (p0+1):p] <- Y[w1, (p0+1):p] + SNR

# calculate p-values
pval <- sansSouci::rowWelchTests(t(X), categ)$p.value

# sanity check: histogram of p-values
hist(pval)
pneuvial commented 3 years ago

e47cfb3 addresses the following points

(dependency to matplotlib not added)

pneuvial commented 2 years ago

TODO: reproduce the first figures of @alexblnn 's paper to illustrate the methods