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(c, k, STAT, Z, df, R, n=1, Q=None, expectations_only=False, version='spm12')

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) 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. >>> P,p,Ec,Ek,EN = rft1d.prob.rft(1, 0, 'T', 2.1, [1,8], [1,10]) 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. 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')

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) An instance of the RFTCalculator class. expected — access to RFT expectations p — access to RFT probabilities isf — inverse survival function sf — survival function >>> 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)

Parameters: alpha – upper tail probability (float; 0 < alpha < 1) Quantile corresponding to upper-tail probability alpha. Equivalently: critical threshold at a Type I error rate of alpha. >>> 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) The probability of exceeding the specified heights. >>> 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: >>> calc = rft1d.prob.RFTCalculator('T', (1,8), 101, 15.0) >>> calc.expected.nodes_per_upcrossing(2.7)

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: >>> calc = rft1d.prob.RFTCalculator('T', (1,8), 101, 15.0) >>> calc.expected.number_of_upcrossings(2.7)
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: >>> calc = rft1d.prob.RFTCalculator('T', (1,8), 101, 15.0) >>> calc.expected.number_of_suprathreshold_nodes(2.8)

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: >>> calc = rft1d.prob.RFTCalculator('T', (1,8), 101, 15.0) >>> calc.expected.number_of_suprathreshold_resels(2.9)
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: >>> calc = rft1d.prob.RFTCalculator('T', (1,8), 101, 15.0) >>> calc.expected.resels_per_upcrossing(3.0)

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 Cluster-specific probability value. >>> calc = rft1d.prob.RFTCalculator('T', (1,8), 101, 15.0) >>> calc.p.cluster(0.1, 3.0)
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 Set-specific probability value. >>> calc = rft1d.prob.RFTCalculator('T', (1,8), 101, 15.0) >>> calc.p.set(2, 0.1, 2.7)
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) The probability of exceeding the specified heights. >>> calc = rft1d.prob.RFTCalculator('T', (1,8), 101, 15.0) >>> calc.sf(3.5)

Bonferroni correction¶

prob.p_bonferroni(STAT, 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 The probability of exceeding the specified height. >>> rft1d.prob.p_bonferroni('Z', 3.1, None, 101) #yields 0.098