Tests: Detecting abnormal region from either one pre and one post...
In ygu427/iSPREADR: imProved Spatial Regression Analysis for DTI Data in R

Description Usage Arguments Details Value Author(s) References See Also Examples

These two functions provides user level interface to conduct hypothesis testing based on spatial smoothing. Briefly, this procedure can be described as: 1. kernel smoothing; 2. compute the mean diff image; 3. compute the original N-stat; 4. use permutation to generate null distributions; 5. global hypothesis testing based on N-stats; 6. local hypothesis testing based on mean diff image, both FWER and FDR adjusted p-values. Please see the pre-print for more details.

pre.post.test(grids, pre, post, perms=1000, balanced=TRUE,
norm=c("L1","L2","Linf"), method=c("null", "normal"), MTP="BH", ...)

rep.test(grids, pre.list, post.list, perms=50,
rand.comb=0, balanced=TRUE, norm=c("L1","L2","Linf"),
method=c("null", "normal"), MTP="BH", ...)

`grids`	A list of monotonically increasing coordinate grids. For example, `grids <- list(grid.x=seq(0, 100, 2), grid.y=exp(seq(1, 5, .2))).`
`pre`	An m-dim array defined on `grids`, the "pre" scan.
`pre.list`	A list of m-dim array defined on `grids`, the "pre" scans.
`post`	An m-dim array defined on `grids`, the "post" scan.
`post.list`	A list of m-dim array defined on `grids`, the "post" scans. pre/post lists may have different number of scans.
`perms`	Number of spatial permutations to run. See 'Details' for more explanations.
`rand.comb`	The number of random scan permutations. See 'Details' for an explanation. The default value is 0, which means `rep.test()` will exhaust all combinations.
`balanced`	A logical variable used by the built-in multiple testing procedure which controls FWER by simultaneous confidence band. Balanced: I assume the difference is symmetrically distributed around zero, so `abs(diff)` is used in the hypothesis testing. Unbalanced: conduct two separate one-sided test for `diff` to obtain `p.up`, `p.down`; the two sided p-value is then defined as `min(1.0, 2*min(p.up, p.down))`. This p-value is conservative and asymptotically exact. Default: balanced=TRUE.
`norm`	A list of functional norms to be used as summary statistics for pre/post test and as kernels in N-distance for repeated scans test. Implemented norms: L1, L2, and Linf. Since the computational cost is very small, by default these functions compute all three of them.
`method`	Smoothing method used by the generalized bootstrap p-value estimation. Valid methods: gld (genral lambda distribution, depends on package `gld` and is quite slow), sn (skew-normal, depends on package `sn`), normal, and null (the usual nonparametric estimate). Default value: `c("null", "normal")`.
`MTP`	Specify the multiple testing procedure used to compute `fdr.adj.p`. Default to `BH`. See help of `p.adjust()` for more options.
`...`	Other options passed to different smoothers such as a choice of kernel (`kernel=box`) or various options of function `Tps()` implemented in fields.

For repeated scan test (rep.test()), the total number of permutations is the number of spatial permutations (controlled by perms) times the number of pre/post group combinations of scan permutations due to the symmetry of the underlying hypothesis testing procedure. Using this symmetry property can greatly reduce the number of kernel smoothings, which is the most expensive computation.

Let there be n1 pre scans and n2 post scans, there are n1+n2 choose n1 of them, which can be large but not astronomically large since nobody gets scanned for more than 10 times in a study in practice and rep.test() can happily compute the case of 8 versus 8 scans (set perm=1 of course) within minutes. Unless rand.comb is set manually, rep.test() exhausts all the combinations because I employ a computation trick to do this very efficiently (see our paper about this). Consequently you don't need to specify a very large number of spatial permutations.

However, for the sake of writing a robust function, and just in case someone will test this program in a completely unintended environment, parameter rand.comb can be specified so this function will do rand.comb many random scan permutations for each spatial permutation (controlled by perms).

For usual nonparametric p-value estimation I recommend using 1000 total permutations, so you need to set perms=1000 for pre.post.test(), or perms=50 for rep.test() with 3 pre scans and 3 post scans (because 6 choose 3 is 20, 20 times 50 makes 1000). Because of this symmetry property, processing many repeated scans byrep.test() is usually much faster than doing a pre/post test!

To further reduce computation cost, you may want to consider using smoothed (generalized bootstrap) p-value estimation, which reduces the variance (so that you can reduce permutations) at the cost of a small biasedness (so either the type I error rate may go up or the detection power may go down). For these tests I recommend 200 total permutations.

A list with components

`p.global`	A list of global p-values global p-values computed by resampling the specified functional norms of the difference map (`pre.post.test()`) or the N-distances with specified norms (kernels).
`rawp`	A list of arrays of unadjusted per-voxel p-values computed from permutations.
`scb.adj.p`	A list of arrays of per-voxel p-values adjusted for familywise error rate by a simultaneous confidence band with different generalized bootstrap method.
`wy.adj.p`	An array of per-voxel p-values adjusted for familywise error rate by the Westfall-Young stepdown procedure.
`fdr.adj.p`	An array of per-voxel p-values adjusted for false discovery rate by the Benjamini-Hochberg procedure.
`fit.diff`	An array of fitted (smoothed) difference map.
`band.up`	A vector of upper bounds of fitted difference map for each permutation.
`band.down`	A vector of lower bounds of fittted difference map for each permutation.

Xing Qiu

[my own paper], [westfall young], [papers about general bootstrap]

genboot.test, Smooth, Ndist, spatial.perm

## Generate a small data
Nx <- 32; Ny <- 32; Nz <- 3
for (i in 1:5) assign(paste("FA",i,sep=""), array(rnorm(Nx*Ny*Nz), c(Nx,Ny,Nz)))
grid.x <- (1:Nx)*2-1; grid.y <- (1:Ny)*2-1; grid.z <- (1:Nz)*2-1; 
grids <- list(grid.x, grid.y, grid.z)
Xlist <- list(FA1, FA2, FA3); Ylist <- list(FA4, FA5)

## takes about a minute on my PC (4x 2.33GHz), YMMV.
## Not run: r1 <- pre.post.test(grids, FA1, FA2)

## takes about 20 seconds on my PC (4x 2.33GHz), YMMV.
## Not run: r2 <- rep.test(grids, Xlist, Ylist, method="normal")