rft1d.prob¶
Random Field Theory expectations and probabilities. The core RFT computations are conducted inside prob.rft, and the RFTCalculator class serves as a high-level interface to prob.rft
Core RFT computations¶
- prob.rft(k, STAT, Z, df, R, n=1, Q=None, expectations_only=False, version='spm12', _0d_check=True)¶
Random Field Theory probabilities and expectations using unified Euler Characteristic (EC) theory. This code is based on “spm_P_RF.m” and “spm_P.m” from the spm8 and spm12 Matlab packages which are available from: http://www.fil.ion.ucl.ac.uk/spm/
- Parameters:
c — number of clusters
k — cluster extent (resels)
STAT — test statistic (one of: “Z”, “T”, “F”, “X2”, “T2”)
Z — field height
df — degrees of freedom [df{interest} df{error}]
R — resel counts (0D counts, 1D counts) defining search volume
n — number of test statistic fields in conjunction
Q — number of field nodes (used for Bonferroni comparison)
expectations_only — if True only expectations will be returned
version — “spm8” or “spm12” (see below)
- Returns:
P — corrected P value
p — uncorrected P value
Ec — expected number of upcrossings {c}
Ek — expected resels per upcrossing {k}
EN — expected excursion set resels
NOTE! If expectations_only==True, then only (Ec,Ek,EN) are returned.
- Examples:
>>> P,p,Ec,Ek,EN = rft1d.prob.rft(1, 0, 'T', 2.1, [1,8], [1,10])
- Notes:
1. The spm8 and spm12 Matlab functions on which this code is based were developed by K.Friston and other members of the Wellcome Trust Centre for Neuroimaging. This function makes minor modifications to those procedures to take advantage of the simplicity of the 1D case.
2. Results for the spm8 and spm12 versions can be obtained via the keyword “version”. When expected ECs approach zero, the spm8 and spm12 results will diverge slightly, due to a minor modification in spm12: In spm8: “EC = EC + eps”. In spm12: “EC = np.array([max(ec,eps) for ec in EC])”
3. Setting c and k in particular manners will yield important probabilities. Consider these three cases:
rft(1, 0, STAT, Z, R) — field maximum
rft(1, k, STAT, u, R) — cluster-based inference
rft(c, k, STAT, u, R) — set-based inference
(a) is the probability that Gaussian fields will produce 1 upcrossing with an extent of 0. Thus this pertains to the maximum of height of Gaussian fields, and can be used, for example, for critical threshold computations.
(b) is the probability that Gaussian fields, when thresholded at u, will produce 1 upcrossing with an extent of k. This is used for cluster-level inference (i.e. p values for individual upcrossings).
(c) is the probability that Gaussian fields, when thresholded at u, will produce c upcrossings with a minimum extent of k. This is used for set-level inference (i.e. p values for the entire result).
Warning
Set-based inference (c) is more powerful than cluster-based inference (b), but unlike (b) it has no localizing information; it is a global p value pertaining to the entire excursion set en masse. It will thus always be lower than (b).
4. If Q==None, then no Bonferroni check is made. If Q!=None, the RFT correction will be compared to Bonferroni correction, and the less severe correction will be returned. This will only have an effect for very rough fields, for example: when then second resel count approaches 0.5*Q.
- References:
Hasofer AM (1978) Upcrossings of random fields. Suppl Adv Appl Prob 10:14-21.
Friston KJ et al (1994) Assessing the significance of focal activations using their spatial extent. Human Brain Mapping 1: 210-220.
Worsley KJ et al (1996) A unified statistical approach for determining significant signals in images of cerebral activation. Human Brain Mapping 4:58-73.
High-level interface to RFT computations¶
- class rft1d.prob.RFTCalculator(STAT='Z', df=None, nodes=101, FWHM=10.0, n=1, withBonf=False, version='spm12')¶
A convenience class for high-level access to RFT probabilities.
- Parameters:
STAT — test statistic (one of: “Z”, “T”, “F”, “X2”, “T2”)
df — degrees of freedom [df{interest} df{error}]
nodes — number of field nodes (int) OR a binary field (boolean array)
FWHM — field smoothness (float)
n — number of test statistic fields in conjunction
withBonf — use a Bonferroni correction if less severe than the RFT correction
version — “spm8” or “spm12” (see below)
- Returns:
An instance of the RFTCalculator class.
- Attributes:
expected — access to RFT expectations
p — access to RFT probabilities
- Methods:
isf — inverse survival function
sf — survival function
- Examples:
>>> calc = rft1d.prob.RFTCalculator('T', (1,8), 101, 15.0) >>> calc.expected.number_of_upcrossings(1.0) #yields 1.343 >>> calc.expected.number_of_upcrossings(4.5) #yields 0.0223
- isf(alpha)¶
Inverse survival function. (see also the survival function: RFTCalculator.sf)
- Parameters:
alpha – upper tail probability (float; 0 < alpha < 1)
- Returns:
Quantile corresponding to upper-tail probability alpha. Equivalently: critical threshold at a Type I error rate of alpha.
- Examples:
>>> calc = rft1d.prob.RFTCalculator('T', (1,8), 101, 15.0) >>> calc.isf(0.05)
- sf(u)¶
Survival function. (Equivalent to RFTCalculator.p.upcrossing)
Probability that 1D Gaussian fields with a smoothness FWHM would produce a 1D statistic field whose maximum exceeds u.
- Parameters:
u – threshold (int, float, or sequence of int or float)
- Returns:
The probability of exceeding the specified heights.
- Examples:
>>> calc = rft1d.prob.RFTCalculator('T', (1,8), 101, 15.0) >>> calc.sf(3.5)
Expectations (RFTCalculator.expected)¶
- expected.nodes_per_upcrossing(u)¶
Number of nodes expected for each uprcrossing at threshold u.
- Example:
Warning
This is a node count, so is equivalent to: (FWHM x resels_per_upcrossing) + 1
- expected.number_of_upcrossings(u)¶
Number of upcrossings expected for threshold u.
- Example:
- expected.number_of_suprathreshold_nodes(u)¶
Number of nodes expected in the entire excursion set at threshold u. These nodes can come from multiple upcrossings.
- Example:
Warning
This is a node count, so is equivalent to: (FWHM x number_of_suprathreshold_resels) + number_of_upcrossings
- expected.number_of_suprathreshold_resels(u)¶
Number of resels expected in the entire excursion set at threshold u. These resels can come from multiple upcrossings. One resel contains (1 x FWHM) nodes. Thus this is equivalent to: (FWHM x number_of_suprathreshold_nodes)
- Example:
- expected.resels_per_upcrossing(u)¶
Number of nodes expected for each uprcrossing at threshold u. One resel contains (1 x FWHM) nodes. Thus this is equivalent to: (FWHM x nodes_per_upcrossing)
- Example:
Probabilities (RFTCalculator.p)¶
- p.cluster(k, u)¶
Cluster-level inference.
Probability that 1D Gaussian fields would produce an upcrossing of extent k when thresholded at u.
Warning
The threshold u should generally be chosen objectively. One possibility is to calculate the alpha-based critical threshold using the inverse survival function: RFTCalculator.isf
- Parameters:
k – cluster extent (resels)
u – threshold
- Returns:
Cluster-specific probability value.
- Examples:
- p.set(c, k, u)¶
Set-level inference.
Probability that 1D Gaussian fields would produce at least c upcrossings with a minimum extent of k when thresholded at u. This probability pertains to the entire excursion set.
Warning
The threshold u should generally be chosen objectively. One possibility is to calculate the alpha-based critical threshold using the inverse survival function: RFTCalculator.isf
- Parameters:
c – number of upcrossings
k – minimum cluster extent (resels)
u – threshold
- Returns:
Set-specific probability value.
- Examples:
- p.upcrossing(u)¶
Survival function (equivalent to RFTCalculator.sf)
Probability that 1D Gaussian fields would produce a 1D statistic field whose maximum exceeds u.
- Parameters:
u – threshold (int, float, or sequence of int or float)
- Returns:
The probability of exceeding the specified heights.
- Examples:
Bonferroni correction¶
- prob.p_bonferroni(z, df, Q, n=1)¶
Bonferroni correction.
When fields are very rough a Bonferroni correction might be less severe than the RFT threshold. This function yields Bonferroni-corrected p values based on the number of field nodes Q.
- Parameters:
STAT — test statistic (one of: “Z”, “T”, “F”, “X2”, “T2”)
z — field height
df — degrees of freedom [df{interest} df{error}]
Q — number of field nodes (used for Bonferroni comparison)
n — number of test statistic fields in conjunction
- Returns:
The probability of exceeding the specified height.
- Example: