Statistical Inference, Multiple Comparisons, Random Field Theory

1
Data Analysis framework in Brain Imaging STAT
992 Image Analysis March 23, 2004 Moo K.
Chung Department of Statistics Department of
Biostatistics and Medical Informatics W.M. Keck
Brain Laboratory University of Wisconsin-Madison
http//www.stat.wisc.edu/mchung
2
Acknowledgement Presentation based on Will
Penny, Jason Lerch and Thomas Nicholss
PowerPoint slides Some images based on Thomas
Hoffmans research
3
Brain Image Analysis

• Brain Image Analysis is a total science.
• Before Images are transformed into a proper
  format for data analysis, it will go through a
  bunch of image processing procedures.
• If we can get extract proper data out of images,
  half of the problems are already solved.

4
image data
parameter estimates
designmatrix
kernel

• General Linear Model
• model fitting
• statistic image

realignment motioncorrection
Random Field Theory
smoothing
fMRI processing steps
normalisation
StatisticalParametric Map
anatomicalreference
corrected p-values
5
MRI processing steps
Montreal Neurological Institute image processing
pipeline
6
Non-uniformity correction
Native
nu_correct native.mnc corrected.mnc
7
Classification/segmentation
classify_clean final.mnc classified.mnc
8
Tissue classification/segmentation Clustering
algorithm based on maximum likelihood mixture
model (Hartigan, 1975)
Automatic skull stripping
9
3 different tissue types
Binary masks 0 or 1
Gaussian kernel smoothing
Gray matter White matter
10
Masking
cortical_surface classified.mnc mask.obj
1.5 surface_mask2 classified.mnc mask.obj
masked.mnc
11
Image registration
12
Corpus callosum (CC) is the white matter brain
substructure that connects hemispheres.
We have 16 autistic subjects and 12 normal
subjects. Quantify the CC shape difference
between two groups?
Group 1
Group 2
13
How do we compare shapes?
Pixel by pixel comparison causes anatomical
mismatching. Solution image registration. The
aim of image registration is to find a smooth
one-to-one mapping that matches homologous
anatomies together.
14
Nonlinear image registration
Estimate a
continuous 3D map that matches two brain images.
15
Thompson, D. W. (1917) On growth and form.
Cambridge University Press, Cambridge.
He introduced deformable grid and deformation of
homologous biological structures
16
Nonlinear image registration based on basis
function expansion
warping
300 subjects average template Warping into
average blurred template reduce the probability
of complete mismatch.
Warped brain
17
Large scale automatic image analysis
. . . .
Subject 1
Subject 2
Subject 498
Subject 3
Subject 499
Subject 500
template
500 MRIs will be warped into a template and
anatomical differences can be compared at a
common reference frame.
18
Estimating Nonlinear Image Registration

• Elastic deformation based
• Fluid dynamics based
• Intensity correlation based
• Bayesian approach

19
Image Registration
Similarity measure
Variational approach
PDE approach
20
3rd order polynomial warping
x' y'
z' -1.212161e02 -1.692117e02
1.239336e00 2.846339e00
9.860215e-01 -4.670216e-03 x
4.541216e-01 3.344188e00 -1.022118e-02
y 2.277959e00 1.849708e00
9.958860e-01 z -9.744798e-03
-4.951447e-03 1.253383e-05 x2
-4.519879e-03 -4.248561e-03 -2.655254e-05
xy -9.122374e-04 -9.371881e-03
4.040382e-05 y2 -1.624103e-02
-3.371953e-03 2.356452e-06 xz
-3.519974e-03 -2.799626e-02 9.228041e-05
yz -1.572948e-02 -1.688950e-03
5.386545e-05 z2 2.495023e-05
4.120123e-06 3.604820e-08 x3
3.232645e-06 1.739698e-05 1.044795e-07
x2y 1.074305e-05 4.357408e-06
-9.302004e-10 xy2 -1.059526e-06
1.699618e-05 1.166377e-07 y3
5.512034e-06 9.330769e-06 -2.219099e-08
x2z 1.275631e-05 -9.233413e-06
1.236940e-07 xyz -5.236010e-07
3.234824e-05 -6.819396e-07 y2z
9.506628e-05 1.214112e-05 -1.238024e-07
xz2 2.016546e-05 1.475354e-04
-1.693465e-08 yz2 3.377913e-06
-7.093638e-05 -2.757074e-07 z3
21
Curve registration by dynamic time warping
algorithm -Thomas Hoffmann, Honors B.Sc. thesis.
22
Image smoothing
23
Kernel smoothing
24
Anisotropic Gaussian kernel smoothing It will
smooth out signals along the eigenvectors. The
amount of smoothing is proportional to the
eigenvalues. So it will basically smooth out
along the principal eigenvector.
Isotropic kernel Anisotropic
kernel
25
Isotropic Gaussian kernel smoothing
Principal eigenvalues gt 0.6
10mm FWHM
20mm FWHM
26
Anisotropic Smoothing Edge Enhancement
27
Before
After diffusion smoothing
28
Inner surface gray/white matter interface
initial mean curvature
diffusion smoothing estimate
29
0.01
Flattened map showing smoothing
0.00
initial mean curvature
20 iterations
100 iterations
30
WHY WE SMOOTH? See next slide Smooth T random
fields
6.5
2.0
-2.0
-6.5
31
AutocorrelationPrecoloring

• Temporally blur, smooth your data
• This induces more dependence!
• But we exactly know the form of the dependence
  induced
• Assume that intrinsic autocorrelation is
  negligible relative to smoothing
• Then we know autocorrelation exactly
• Correct GLM inferences based on known
  autocorrelation

Friston, et al., To smooth or not to smooth
NI 12196-208 2000
32
AutocorrelationPrewhitening

• Statistically optimal solution
• If know true autocorrelation exactly, canundo
  the dependence
• De-correlate your data, your model
• Then proceed as with independent data
• Problem is obtaining accurate estimates of
  autocorrelation
• Some sort of regularization is required
• Spatial smoothing of some sort

33
fMRI example
34
Basic fMRI Example

• Time series at each voxel.

35
fMRI time series modeling
one stimulus
two stimuli
36
A Linear Model
error

• Linear in parameters ?1 ?2

b1
b2
Time
e
x1
x2
Intensity
37
in matrix form.
N Number of scans, p Number of regressors
38
subject 20
left
right
left
right
39
subject 41
left
right
left
right
40
OSL fitting of subject 20 left amygdala
41
OSL fitting of subject 20 right amygdala
42
Generalized Least Squares (GLS) Estimation
43
HRF snake attacking snake crawling
fish swimming
AFNI result subject 20 right amygdala
44
HRF reconvolved with the initial stimuli
Black HR based on OSL Red HR based on GSL In
this particular example, GSL can get the dip OSL
can not get.
45
(No Transcript)
46
Whitening by GSL correlation
47
Multiple Testing Problem

• Inference on statistic images
• Fit GLM at each voxel
• Create statistic images of effect
• Which of 100,000 voxels are significant?
• ?0.05 ? 5,000 false positives!

48
MCP SolutionsMeasuring False Positives

• Familywise Error Rate (FWER)
• Familywise Error
• Existence of one or more false positives
• FWER is probability of familywise error corrected
  P-value
• False Discovery Rate (FDR)
• R voxels declared active, V falsely so
• Observed false discovery rate V/R
• FDR E(V/R) Q-value
• This is a relative measure.

49
FWER MCP SolutionsRandom Field Theory

• Euler Characteristic ?u
• Topological Measure
• blobs - holes
• At high thresholds,just counts blobs
• FWER P(Max voxel ? u Ho) P(One or more
  blobs Ho) ? P(?u ? 1 Ho) ? E(?u Ho)

Threshold
Random Field
Suprathreshold Sets
50
Example 2D Gaussian images
a R (4 ln 2) (2p) -3/2 u exp (-u2/2)
For R100 and a0.05 RFT gives u3.8
Using R100 in a Bonferroni correction gives
u3.3
Friston et al. (1991) J. Cer. Bl. Fl. M.
51
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.
52
Unified Theory

• General form for expected Euler characteristic
• ?2, F, t fields restricted search regions
• a 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
53
Expected EC for a stationary Gaussian field
Minkowski functionals
EC density
Adler (1981) , Worsley (1996)
54
RFT Assumptions

• 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

55
Random Field Theory

• Closed form results for E(?u)
• Z, t, F, Chi-Squared Continuous RFs
• Results depend only on RESELs
• RESolution ELements
• A volume element of size FWHMx ? FWHMy ? FWHMz
• RESEL count
• Voxel-size-independent measure of volume
• Inferences do not depend on stationarity
• Smoothness can vary
• (It is cluster size results that require
  stationarity)

56
Random Field Theory Limitations

• Sufficient smoothness
• FWHM smoothness 3-4 times voxel size
• Smoothness estimation
• Estimate is biased when images not so smooth
• Estimate is an estimatesds on corr. p-values
• Multivariate normality
• Virtually impossible to check
• Several layers of approximations
• E.g., t field results conservative for low df

57
Smoothness Estimation

• 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)
• a R3(?) (4ln(2))3/2 (u 2 -1) exp(-u 2/2) /
  (2?)2

Approximate the peak of the Covariance
function with a Gaussian
58
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 connected components 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

59
Suprathreshold cluster tests

• Primary threshold u
• examine connected components of excursion set
• Suprathreshold clusters
• Reject HW for clusters of voxels W of size S gt s?
• Localisation (Strong control)
• at cluster level
• increased power
• esp. high resolutions (f MRI)
• Thresholds, p values
• Pr(S?max gt s? ? H? ) ? ?
• Nosko, Friston, (Worsley)
• Poisson occurrence (Adler)
• Assumme form for Pr(SsSgt0)

5mm FWHM
10mm FWHM
15mm FWHM
(2mm2 pixels)
60
Poisson Clumping Heuristic
Expected number of clusters
pcluster volume gt k
Expected cluster volume
EC density (??
Smoothness
Search volume (R)
61
Multiple Comparisons, Random Field Theory
Ch5
Ch4

• Worsley KJ, Marrett S, Neelin P, Evans AC (1992)
• A three-dimensional statistical analysis for
  CBF activation studies in human brainJournal of
  Cerebral Blood Flow and Metabolism 12900-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
  activationHuman Brain Mapping 458-73
• Friston KJ, Worsley KJ, Frackowiak RSJ, Mazziotta
  JC, Evans AC (1994)Assessing the Significance
  of Focal Activations Using their Spatial
  ExtentHuman Brain Mapping 1214-220
• Cao J (1999)The size of the connected
  components of excursion sets of ?2, t and F
  fieldsAdvances in Applied Probability (in
  press)
• Worsley KJ, Marrett S, Neelin P, Evans AC
  (1995)Searching scale space for activation in
  PET imagesHuman Brain Mapping 474-90
• Worsley KJ, Poline J-B, Vandal AC, Friston KJ
  (1995)Tests for distributed, non-focal brain
  activationsNeuroImage 2183-194
• Friston KJ, Holmes AP, Poline J-B, Price CJ,
  Frith CD (1996)Detecting Activations in PET and
  fMRI Levels of Inference and PowerNeuroimage
  4223-235

62
Controlling FWER Permutation Test

• Parametric methods
• Assume distribution ofmax statistic under
  nullhypothesis
• Nonparametric methods
• Use data to find distribution of max
  statisticunder null hypothesis
• Any max statistic!
• Due to Fisher

63
Measuring False PositivesFWER vs FDR
Noise
SignalNoise
64
Control of Per Comparison Rate at 10
Percentage of Null Pixels that are False Positives
Control of Familywise Error Rate at 10
FWE
Occurrence of Familywise Error
Control of False Discovery Rate at 10
Percentage of Activated Pixels that are False
Positives
65
Controlling FDRBenjamini Hochberg

• Select desired limit q on E(FDR)
• Order p-values, p(1) ? p(2) ? ... ? p(V)
• Let r be largest i such that
• Reject all hypotheses corresponding to p(1),
  ... , p(r).

1
p(i)
p-value
i/V ? q
0
0
1
i/V
66
Benjamini HochbergVarying Signal Extent
p 0.019274 z 2.07
7
67
FDR Software for SPM
http//www.sph.umich.edu/nichols/FDR
68
References

• KM Petersson, TE Nichols, J-B Poline, and AP
  Holmes.Statistical limitations in functional
  neuroimaging I. Non-inferential methods and
  statistical models. Statistical limitations in
  functional neuroimaging II. Signal detection and
  statistical inference. Philosophical
  Transactions of the Royal Society Biological
  Sciences, 3541239-1281, 1999.
• CR Genovese, N Lazar and TE Nichols.Thresholding
  of Statistical Maps in Functional Neuroimaging
  Using the False Discovery Rate. NeuroImage,
  15870-878, 2002.

http//www.sph.umich.edu/nichols
69
Application Surface data
70
Corrected P-value
Morphological descriptor
Cortical surface area
2D Euler characteristic density of t random field
71
Rejection regions.
Red Tissue growth Blue Tissue loss Yellow
Structure displacement
72
Yellow Hotellings T2 field
73
(No Transcript)
74
Cortical thickness change t-map
75
Curvature change t-map