Development of PowerMap: a Software Package for Statistical Power Calculation in Neuroimaging Studies

Abstract

Although there are a number of statistical software tools for voxel-based massively univariate analysis of neuroimaging data, such as fMRI (functional MRI), PET (positron emission tomography), and VBM (voxel-based morphometry), very few software tools exist for power and sample size calculation for neuroimaging studies. Unlike typical biomedical studies, outcomes from neuroimaging studies are 3D images of correlated voxels, requiring a correction for massive multiple comparisons. Thus, a specialized power calculation tool is needed for planning neuroimaging studies. To facilitate this process, we developed a software tool specifically designed for neuroimaging data. The software tool, called PowerMap, implements theoretical power calculation algorithms based on non-central random field theory. It can also calculate power for statistical analyses with FDR (false discovery rate) corrections. This GUI (graphical user interface)-based tool enables neuroimaging researchers without advanced knowledge in imaging statistics to calculate power and sample size in the form of 3D images. In this paper, we provide an overview of the statistical framework behind the PowerMap tool. Three worked examples are also provided, a regression analysis, an ANOVA (analysis of variance), and a two-sample T-test, in order to demonstrate the study planning process with PowerMap. We envision that PowerMap will be a great aide for future neuroimaging research.

Appendix A

For three independent random variables X, Y, and Z, (means μ _X , μ _Y , and μ _Z , and variances $\sigma_X^2,\sigma_Y^2,\;{\text{and}}\;\sigma_Z^2$, respectively), the variance of the product XYZ is

$$ \matrix{ {Var\left( {XYZ} \right) = \sigma_X^2\sigma_Y^2\sigma_Z^2 + \mu_X^2\sigma_Y^2\sigma_Z^2 + \mu_Y^2\sigma_X^2\sigma_Z^2 + \mu_Z^2\sigma_X^2\sigma_Y^2} \\ { + \mu_X^2\mu_Y^2\sigma_Z^2 + \mu_X^2\mu_Z^2\sigma_Y^2 + \mu_Y^2\mu_Z^2\sigma_X^2} \\ }<!end array> $$

(A1)

Another useful result is that, for a chi-square random variable W with df = v, then the expected value of W ^r for real-valued r is given by

$$ E\left[ {{W^r}} \right] = \frac{{{2^r}\left( {\frac{v}{2} + r} \right)}}{{\Gamma \left( {\frac{v}{2}} \right)}} $$

(A2)

where Γ denotes the gamma function (Johnson et al. 1995). Worsley has shown that the gradient of a 3D T-random field T, ∇T, follows the distribution

$$ \nabla T = {m^{{\frac{1}{2}}}}\left( {1 + \frac{{{T^2}}}{m}} \right){S^{{ - \frac{1}{2}}}}Z $$

where T is the value of a T-random random field with df = m at the location where the gradient is calculated, S is a chi-square random variable with df = m + 1, and Z is a 3 × 3-dimensional normal random variable with mean zero and variance defined by a matrix Λ, with T, S, and Z are all independent (Worsley 1994). The gradient of T depends on the value of T itself. The variance matrix Λ needs to be estimated based on an observed T-statistic image. Since ∇T can be seen as the product of three independent random variables, namely, $X = {m^{{1/2}}}\left( {{{{1 + {T^2}}} \left/ {m} \right.}} \right),Y = {S^{{ \frac{{ - 1}}{2} }}}$, and Z, we can use (A1) to calculate its variance. Using (A2), we can calculate the means and the variances of these variables:

$$ \begin{array}{*{20}{c}} {{{\mu }_{X}} = {{m}^{{\frac{1}{2}}}}\left( {\frac{{m - 1}}{{m - 2}}} \right)\sigma _{X}^{2} = \frac{{2m\left( {m - 1} \right)}}{{{{{\left( {m - 2} \right)}}^{2}}\left( {m - 4} \right)}}} \\ {{{\mu }_{Y}} = \frac{{{{2}^{{\frac{1}{2}}}}\Gamma \left( {m/2} \right)}}{{\Gamma \left( {\frac{{m + 1}}{2}} \right)}}\sigma _{Y}^{2} = \frac{1}{{m - 1}} - \mu _{Y}^{2}} \\ {{{\mu }_{Z}} = \left( {\begin{array}{*{20}{c}} {000} \hfill \\ {000} \hfill \\ {000} \hfill \\ \end{array} } \right)\sigma _{Z}^{2} = \Lambda } \\ \end{array} $$

Since μ _Z is zero, some terms in (A1) cancel and the variance of ∇T can be obtained as

$$ Var\left( {\nabla T} \right) = \left( {\sigma_X^2\sigma_Y^2 + \mu_X^2\sigma_Y^2 + \mu_Y^2\sigma_X^2 + \mu_X^2\mu_Y^2} \right)\Lambda $$

and this can be solved for Λ to yield (4). In order to derive the similar results for ∇F, we focus on the variance of ∇G where G is a random field defined as (m/n)F. In other words, G is a random field defined as the ratio U/V of two chi-square random fields U (with df = m) and V (df = n). Worsley (Worsley 1994) has described that ∇G can be described as

$$ \nabla G = 2{G^{{ - \frac{1}{2}}}}\left( {1 + G} \right){W^{{ - \frac{1}{2}}}}Z $$

where G is the value of a G-random field at the location where the gradient is calculated, with G defined as (m/n)F, F is an F-random variable with df = (m,n), W is a chi-square random variable with df = m + n, and Z is a 3 × 3-dimensional normal random variable as described above. The gradient of G depends on G itself. If we let X = 2G ^1/2 (1 + G), Y = W^-1/2, and Z as it is, then their means and variances are

$$ \begin{array}{*{20}{c}} {{{\mu }_{X}} = \left( {m + n - 2} \right)\left\{ {\Gamma \left( {\frac{{m + 1}}{2}} \right)\Gamma \left( {\frac{{n - 3}}{2}} \right){{{\left( {\Gamma \left( {\frac{m}{2}} \right)\Gamma \left( {\frac{n}{2}} \right)} \right)}}^{{ - 1}}}} \right\}} \\ {\sigma _{X}^{2} = \left( {m + n - 2} \right)\left( {m + n - 4} \right)\frac{m}{{\left( {n - 2} \right)\left( {n - 4} \right)\left( {n - 6} \right)}} - \mu _{X}^{2}} \\ {{{\mu }_{Y}} = \frac{1}{{\sqrt {2} }}\frac{{\Gamma \left( {\frac{{m + n - 1}}{2}} \right)}}{{\Gamma \left( {\frac{{m + n}}{2}} \right)}}\quad \sigma _{Y}^{2} = \frac{1}{{m + n - 2}} - \mu _{Y}^{2}} \\ {{{\mu }_{Z}} = \left( {\begin{array}{*{20}{c}} 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill \\ \end{array} } \right)\quad \sigma _{Z}^{2} = \Lambda } \\ \end{array} $$

Then we can use these results in (A1) to derive the variance of ∇G. As in the T-random field described above, since μ _Z is zero, some terms in (A1) cancel and the variance of ∇G can be obtained as

$$ Var\left( {\nabla G} \right) = \left( {\sigma_X^2\sigma_Y^2 + \mu_X^2\sigma_Y^2 + \mu_Y^2\sigma_X^2 + \mu_X^2\mu_Y^2} \right) $$

Finally, since F = (n/m)G, ∇F = (n/m)∇G. Therefore Var(∇F) = (n/m)² Var(∇G) and this can be solved for Λ as in the T-random field, resulting in (3).