Statistical Inference and Random Field Theory

Statistical Inference and RandomField Theory Will Penny SPM short course, Kyoto, Japan, 2002

image data parameter estimates designmatrix kernel • General Linear Model • model fitting • statistic image realignment &motioncorrection Random Field Theory smoothing normalisation StatisticalParametric Map anatomicalreference corrected p-values

Overview …a voxel by voxel hypothesis testing approach • reliably identify regions showing a significant experimental effect • Assessment of statistic images • multiple comparisons • random field theory • smoothness • spatial levels of inference & power

Null hypothesis H test statistic null distributions Hypothesis test control Type I error incorrectly reject H test level or error rate Pr(“reject” H | H)  p –value min a at which Hrejected Pr(T t | H) characterising “surprise” t –distribution, 32 df. F –distribution, 10,32 df. Classical hypothesis testing

Multiple comparisons terminology • Family of hypotheses • Hk k = {1,…,K} • H = H1H2 … HkHK • Familywise Type I error • weak control – omnibus test • Pr(“reject” HH)  • “anything, anywhere”? • strong control – localising test • Pr(“reject” HW HW)   W: W   & HW • “anything, & where”? Activation is zero everywhere eg. Look at average activation over volume eg. Look at maxima of statistical field

The Bonferroni correction Given a family of N independent voxels and a voxel-wise error rate v the Family-Wise Error rate (FWE) or ‘corrected’ error rate is α = 1 – (1-v)N ~ Nv Therefore, to ensure a particular FWE we choose v = α / N A Bonferroni correction is appropriate for independent tests If v=0.05 then over 100 voxels we’ll get 5 voxel-wise type I errors. But we’ll get a much higher α. To ensure α=0.05 we need v=0.0005 !

The Bonferroni correction Independent Voxels Spatially Correlated Voxels Bonferroni is too conservative for brain images

Consider a statistic image as a lattice representation of a continuous random field Use results from continuous random field theory Random Field Theory Lattice representation

Topological measure threshold an image at u excursion set Au (Au)=# blobs - # holes At high u, (Au)=# blobs Reject HΩ if Euler char non-zero α Pr((Au)> 0 ) Expected Euler char p–value (at high u) αE[(Au)] Euler Characteristic (EC)

Example – 2D Gaussian images • α = R (4 ln 2) (2π) -3/2 u exp (-u2/2) Voxel-wise threshold, u Number of Resolution Elements (RESELS), R N=100x100 voxels, Smoothness FWHM=10, gives R=10x10=100

Example – 2D Gaussian images • α = R (4 ln 2) (2π) -3/2 u exp (-u2/2) For R=100 and α=0.05 RFT gives u=3.8 Using R=100 in a Bonferroni correction gives u=3.3 Friston et al. (1991) J. Cer. Bl. Fl. M.

Developments 2D Gaussian fields 3D Gaussian fields 3D t-fields Friston et al. (1991) J. Cer. Bl. Fl. M. (Not EC Method) Worsley et al. (1992) J. Cer. Bl. Fl. M. Worsley et al. (1993) Quant. Brain. Func.

Restricted search regions Box and frame have same number of voxels Box has 16 markers Frame has 32 markers

Au  Unified Theory • General form for expected Euler characteristic • 2, F, & t fields • restricted search regions α = S Rd (W)rd (u) rd (u): EC density; depends on type of field (eg. Gaussian, t) and the threshold, u. Rd (W): RESEL count; depends on the search region – how big, how smooth, what shape ? Worsley et al. (1996), HBM

Au  Unified Theory • General form for expected Euler characteristic • 2, F, & t fields • restricted search regions α = S Rd (W)rd (u) rd (u):d-dimensional EC density – E.g. Gaussian RF: r0(u) = 1- (u) r1(u) = (4 ln2)1/2 exp(-u2/2) / (2p) r2(u) = (4 ln2) exp(-u2/2) / (2p)3/2 r3(u) = (4 ln2)3/2 (u2 -1) exp(-u2/2) / (2p)2 r4(u) = (4 ln2)2 (u3 -3u) exp(-u2/2) / (2p)5/2 Rd (W): RESEL count R0(W) = (W) Euler characteristic of W R1(W) = resel diameter R2(W) = resel surface area R3(W) = resel volume Worsley et al. (1996), HBM

Resel Counts for Brain Structures FWHM=20mm

Functional Imaging Data • The Random Fields are the component fields, Y = Xw +E, e=E/σ • We can only estimate the component fields, using estimates of w and σ • To apply RFT we need the RESEL count which requires smoothness estimates

Y = X  +  ? errors Component fields voxels ? =  + parameters design matrix data matrix scans  variance componentfields

^  Estimated component fields voxels ? ? =  + parameters design matrix errors data matrix scans parameterestimates • estimate   residuals  estimated variance = Each row is an estimated component field estimatedcomponentfields

Roughness || Point Response Function PRF Gaussian PRF fx 0 0  = 0 fy 0 0 0 fz || = (4ln(2))3/2 / (fx fy fz) RESEL COUNT R3() = () / (fx fy fz) α = R3() (4ln(2))3/2 (u 2 -1) exp(-u 2/2) / (2)2 Smoothness Estimation Approximate the peak of the Covariance function with a Gaussian

Multiple comparisons terminology • Family of hypotheses • Hk k = {1,…,K} • H = H1H2 … HkHK • Familywise Type I error • weak control – omnibus test • Pr(“reject” HH)  • “anything, anywhere”? • strong control – localising test • Pr(“reject” HW HW)   W: W   & HW • “anything, & where”? Activation is zero everywhere eg. Look at average activation over volume eg. Look at maxima of statistical field

Cluster and Set-level Inference • We can increase sensitivity by trading off anatomical specificity • Given a voxel level threshold u, we can compute the likelihood (under the null hypothesis) of getting n or more connectedcomponents in the excursion set ie. a cluster containing at least n voxels CLUSTER-LEVEL INFERENCE • Similarly, we can compute the likelihood of getting c clusters each having at least n voxels SET-LEVEL INFERENCE

Primary threshold u examine connected components of excursion set Suprathreshold clusters Reject HW for clusters of voxels W of size S > s Localisation (Strong control) at cluster level increased power esp. high resolutions (f MRI) Thresholds, p –values Pr(Smax > s H )  Nosko, Friston, (Worsley) Poisson occurrence (Adler) Assumme form for Pr(S=s|S>0) Suprathreshold cluster tests 5mm FWHM 10mm FWHM 15mm FWHM (2mm2 pixels)

n=12 n=82 n=32 Levels of inference voxel-level P(c  1 | n > 0, t  4.37) = 0.048 (corrected) At least one cluster with unspecified number of voxels above threshold set-level P(c  3 | n  12, u  3.09) = 0.019 At least 3 clusters above threshold Parameters u - 3.09 k - 12 voxels S - 323 voxels FWHM - 4.7 voxels D - 3 cluster-level P(c  1 | n  82, t  3.09) = 0.029 (corrected) At least one cluster with at least 82 voxels above threshold

SPM results...

Model fit & assumptions valid distributional results Multivariate normality of component images Covariance function of component images must be - Stationary(pre SPM99) - Can be nonstationary (SPM99 onwards) - Twice differentiable Smoothness smoothness » voxel size lattice approximation smoothness estimation practically FWHM 3 VoxDim otherwise conservative RFT Assumptions

Ch5 Ch4 Multiple Comparisons,& Random Field Theory Worsley KJ, Marrett S, Neelin P, Evans AC (1992) “A three-dimensional statistical analysis for CBF activation studies in human brain”Journal of Cerebral Blood Flow and Metabolism12:900-918 Worsley KJ, Marrett S, Neelin P, Vandal AC, Friston KJ, Evans AC (1995) “A unified statistical approach for determining significant signals in images of cerebral activation”Human Brain Mapping 4:58-73 Friston KJ, Worsley KJ, Frackowiak RSJ, Mazziotta JC, Evans AC (1994)“Assessing the Significance of Focal Activations Using their Spatial Extent”Human Brain Mapping 1:214-220 Cao J (1999)“The size of the connected components of excursion sets of 2, t and F fields”Advances in Applied Probability (in press) Worsley KJ, Marrett S, Neelin P, Evans AC (1995)“Searching scale space for activation in PET images”Human Brain Mapping 4:74-90 Worsley KJ, Poline J-B, Vandal AC, Friston KJ (1995)“Tests for distributed, non-focal brain activations”NeuroImage 2:183-194 Friston KJ, Holmes AP, Poline J-B, Price CJ, Frith CD (1996)“Detecting Activations in PET and fMRI: Levels of Inference and Power”Neuroimage 4:223-235

Summary • We should correct for multiple comparisons • We can use Random Field Theory (RFT) • RFT requires (i) a good lattice approximation to underlying multivariate Gaussian fields, (ii) that these fields are continuous with a twice differentiable correlation function • To a first approximation, RFT is a Bonferroni correction using RESELS. • We only need to correct for the volume of interest. • Depending on nature of signal we can trade-off anatomical specificity for signal sensitivity with the use of cluster-level inference.

Statistical Inference and Random Field Theory