Exploring random field theory for correcting for multiple comparisons
The problem of false positives with multiple statistical tests is an old one. One standard method for dealing with this problem is to use the Bonferroni correction. For the Bonferroni correction, you set your p value threshold for accepting a test as being significant as $alpha$ / (number of tests), where $alpha$ is the false positive rate you are prepared to accept. $alpha$ is often 0.05, or one false positive in 20 repeats of your experiment. Thus, for a statistical map with 200000 voxels, the Bonferroni corrected p value would be 0.05 / 200000 = [equivalent Z] 5.03. We could then threshold our Z map to show us only Z scores higher than 5.03, and be confident that all the remaining Z scores are unlikely to have occurred by chance. For some functional imaging data this is a perfectly reasonable approach, but in most cases the Bonferroni threshold will be considerably too conservative. This is because, for most stastistic maps, the Z scores at each voxel are highly correlated with their neighbours.
An example can show why the Bonferroni correction is too conservative with non-independent tests. Let us first make an example image out of random numbers. We generate 16384 random numbers, and then put them into a 128 by 128 array. This results in a 2D image of spatially independent random numbers.
%pylab inline
import numpy as np
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
# Constants for image simulations etc
shape = [128, 128] # No of pixels in X, Y
n_voxels = np.prod(shape) # Total pixels
s_fwhm = 8 # Smoothing in number of pixels in x, y
seed = 1966 # Seed for random no generator
alpha = 0.05 # Default alpha level
# Image of independent random nos
np.random.seed(seed) # Seed the generator to get same numbers each time
test_img = np.random.standard_normal(shape)
plt.figure()
plt.imshow(test_img)
plt.set_cmap('bone')
plt.xlabel('Pixel position in X')
plt.ylabel('Pixel position in Y')
plt.title('Image 1 - array of independent random numbers')
<matplotlib.text.Text at 0x3ad9f70>
# The inverse normal CDF
import scipy.stats
inv_n_cdf = scipy.stats.norm().ppf
# Bonferroni threshold for this image
bonf_thresh = inv_n_cdf(1 - (alpha / n_voxels))
bonf_thresh
4.5227713755927113
# Divide into FWHM chunks and fill square from mean value
sqmean_img = test_img.copy()
tarr = np.ones((s_fwhm, s_fwhm)) * s_fwhm # mutliply up to unit variance
for i in range(0, shape[0], s_fwhm):
i_slice = slice(i, i+s_fwhm)
for j in range(0, shape[1], s_fwhm):
j_slice = slice(j, j+s_fwhm)
vals = sqmean_img[i_slice, j_slice]
sqmean_img[i_slice, j_slice] = tarr * vals.mean()
plt.figure()
plt.imshow(sqmean_img)
plt.set_cmap('bone')
# axis xy
plt.xlabel('Pixel position in X')
plt.ylabel('Pixel position in Y')
plt.title('Taking means over %s by % elements from image 1' % (s_fwhm, s_fwhm))
<matplotlib.text.Text at 0x4d13270>
# smooth random number image
import scipy.ndimage as spn
sd = s_fwhm / np.sqrt(8.*np.log(2)) # sigma for this FWHM
stest_img = spn.filters.gaussian_filter(test_img, sd, mode='wrap')
# Restore smoothed image to unit variance - this code is tough
# 2d Gaussian kernel
half_dims = [s // 2 for s in shape]
assert np.all(half_dims == np.divide(shape, 2.0)) # check divisible by 2
y_inds = range(-half_dims[1], half_dims[1]+1)
x_inds = range(-half_dims[0], half_dims[0]+1)
[x,y] = np.meshgrid(y_inds, x_inds)
gf = np.exp(-(x*x + y*y) / (2*sd*sd))
gf = gf/np.sum(gf)
# expectation of variance for this kernel
AG = np.fft.fft2(gf)
Pag = AG * np.conj(AG) # Power of the noise
COV = np.real(np.fft.ifft2(Pag))
svar = COV[0, 0]
# smoothed image to unit variance
scf = np.sqrt(1 / svar)
stest_img = stest_img * scf
# display smoothed image
plt.figure()
plt.imshow(stest_img)
plt.set_cmap('bone')
# axis xy
plt.xlabel('Pixel position in X')
plt.ylabel('Pixel position in Y')
plt.title('Image 1 - smoothed with Gaussian kernel of FWHM %s by %s pixels' %
(s_fwhm, s_fwhm))
<matplotlib.text.Text at 0x58d4690>
# No of resels
resels = np.prod(np.divide(shape, s_fwhm))
resels
256
We make a function to show the thresholded image:
def show_threshed(img, th):
thimg = (img > th)
plt.figure()
plt.imshow(thimg)
plt.set_cmap('bone')
# axis xy
plt.xlabel('Pixel position in X')
plt.ylabel('Pixel position in Y')
plt.title('Smoothed image thresholded at Z > %s' % th)
# threshold at 2.75 and display
show_threshed(stest_img, 2.75)
# threshold at 3.5 and display
show_threshed(stest_img, 3.5)
We can predict the Euler characteristic from the smoothness:
# expected EC at various Z thresholds, for two dimensions
plt.figure()
Z = np.arange(0, 5, 0.01)
expEC = (resels* (4 * np.log(2)) * ((2*np.pi)**(-3./2)) * Z) * np.exp((Z ** 2)*(-0.5))
plt.plot(Z, expEC)
plt.xlabel('Z score threshold')
plt.ylabel('Expected EC for thresholded image')
plt.title('Expected EC for smoothed image with %s resels' % resels)
<matplotlib.text.Text at 0x5985a50>
How accurate is this? We can do a simulation:
# simulation to test the RFT threshold.
# find approx threshold from the vector above
tmp = (Z > 3) & (expEC<=alpha)
alphaTH = Z[tmp][0]
repeats = 1000
falsepos = np.zeros((repeats,))
maxes = np.zeros((repeats,))
edgepix = s_fwhm # to add edges to image - see below
big_size = [s + edgepix * 2 for s in shape]
i_slice = slice(edgepix + 1, shape[0] + edgepix + 1)
j_slice = slice(edgepix + 1, shape[1] + edgepix + 1)
for i in range(repeats):
timg = np.random.standard_normal(size=big_size) # gives extra edges to throw away
stimg = spn.filters.gaussian_filter(timg, sd, mode='wrap')
# throw away edges to avoid artefactually high values
# at image edges generated by the smoothing
stimg = stimg[i_slice, j_slice]
# Reset variance using scale factor from calculation above
stimg *= scf
falsepos[i] = np.any(stimg >= alphaTH)
maxes[i] = stimg.max()
print('False positive rate in simulation was %s (%s expected)' %
(sum(falsepos) / float(repeats), alpha))
False positive rate in simulation was 0.04 (0.05 expected)