Simulation 1:
- Three-group permutation test
- n=120 (40 per group)
- effect size = effect size (0, 0, 8)
- Number of voxels: 10 x 10 x 4
- Significant area size: 2 x 2 x 1
- Multiple comparison correction (p.adjust, method=BH)
- Report false discovery rates.
- Random samples: 10,20,100
- Permuting numbers: 1000 vs 10,000
img.dim <- c(10,10,4)
n <- 120
set.seed(1)
voxel_data <- array(rnorm(prod(img.dim)*n),dim=c(img.dim,n))
#Create dataset with effect size (0, 0, 0.8)
voxel_data_set1 <- voxel_data
voxel_data_set1[4:5,4:5,3,81:120] <- voxel_data_set1[4:5,4:5,3,81:120] + 0.8
#Show the means have actually changed
lapply(1:80, function(x) mean(voxel_data_set1[4:5,4:5,3,x])) %>% unlist() %>% mean
## [1] 0.03646248
lapply(81:120, function(x) mean(voxel_data_set1[4:5,4:5,3,x])) %>% unlist() %>% mean
## [1] 0.7858647
#Create the nifti file
simulation_nifti <- nifti(voxel_data_set1, dim = c(10, 10, 4, 120))
#writeNIfTI(simulation_nifti, "simulation1_data")
#Look-up FSL implentation in R
cov_mat_randomise <- rep(1:3, rep(40,3)) %>% as.factor
cov_mat_parperm <- cbind(rep(c(0,1,0),c(40,40,40)), rep(c(0,0,1),c(40,40,40)))
xtx <- xmat_trans(cov_mat_parperm)
#Dataset1 -- show functions from parPerm work
ymat_set1 <- ymat_trans(simulation_nifti)
cores <- detectCores()
cl <- makeCluster(cores[1]-1)
registerDoParallel(cl)
result_1000_set1 <- run_permtest_multicat(x.mat = xtx, y = ymat_set1, columns = 1:2, split = 2, num_perms = 1000)
result_10000_set1 <- run_permtest_multicat(x.mat = xtx, y = ymat_set1, columns = 1:2, split = 2, num_perms = 10000)
system.time(run_permtest_multicat(x.mat = xtx, y = ymat_set1, columns = 1:2, split = 2, num_perms = 10000))
## user system elapsed
## 0.028 0.001 7.340
quick_df <- cbind(result_1000_set1, result_10000_set1) %>% as.data.frame()
names(quick_df) <- c("Set1000_B0", "Set1000_B1", "Set10000_B0", "Set10000_B1")
quick_df %>%
mutate_if(is.numeric, function(x) p.adjust(p = x, method = "BH")) %>%
arrange(Set10000_B1) %>%
head() %>%
kable()
| 0.9974937 |
0.0000000 |
0.9988 |
0.0000000 |
| 0.9974937 |
0.0000000 |
0.9988 |
0.0200000 |
| 0.9980000 |
0.0000000 |
0.9988 |
0.0266667 |
| 0.9974937 |
0.0000000 |
0.9988 |
0.0400000 |
| 0.7500000 |
0.0800000 |
0.7400 |
0.2560000 |
| 0.9974937 |
0.4666667 |
0.9988 |
0.4866667 |