pasta package

Submodules

pasta.pasta module

pasta.pasta.cov2nef(c_para1, c_para2)[source]

Compute effective sample size from covariance matrices c_para1 and c_para2.

Is a computational efficient implementation equivalent to:

nef=real(1/(trace(B*fc_para1*B*fc_para2)/(trace(B*fc_para1)*trace(B*fc_para2)))+1);
Parameters:

  • c_para1 (ndarray) – Covariance matrix.

  • c_para2 (ndarray) – Covariance matrix.

Returns:

nef – Effective sample size.

Return type:

float

pasta.pasta.covariance_estimation(x, coord=None, D=None, dim=None, M=None, qd=0.7, xparc=None, max_clusters=10, min_cluster_size=500, min_clusters=1, M_cluster=None, nugget=True)[source]

Compute the covariance matrix for a single map x using PaSTA or PaSTA-NS.

This can be particularly useful when pairwise association between a large number of maps needs to be evaluated. Compute the covariance matrix for each data separately and save for later use can avoid repetitive covariance estimation in the effective_sample_size_estimation function.

Statistical significance between two maps can be inferred by loading saved covariance matrices (cov1 and cov2) of two maps (x and y), and compute effective sample size and p-values following steps below:

get submatrix of covariance matrices for points that are valid in both maps - valid = np.isfinite(x) & np.isfinite(y), cov1 = cov1[np.ix_(valid, valid)], cov2 = cov2[np.ix_(valid, valid)], x = x[valid], y = y[valid]

compute nef - cov2nef(cov1, cov2)

compute test statistics such as Pearson correlation coefficient — rX, p_naive = pearsonr(x, y)

compute significance p-value from test statistics and effective sample size: nef2p(rX, nef)

Inputs are same as in effective_sample_size_estimation but with y and yparc removed

Returns:

covmat – Covariance matrix of map x in shape (N, N), where rows and columns corresponding to invalid observations in x (e.g., NaN, Inf) are set to np.nan and need to be removed before computing nef.

Return type:

ndarray (N, N)

pasta.pasta.effective_sample_size_estimation(x, y, coord=None, D=None, dim=None, M=None, qd=0.7, xparc=None, yparc=None, max_clusters=10, min_cluster_size=500, min_clusters=1, M_cluster=None, nugget=True)[source]

Main function that runs PaSTA and PaSTA-NS to compute effective sample size and autocorrelation-corrected p-values.

Leave xparc=None and yparc=None will run PaSTA, while setting to ‘auto’ or user-specified parcellation np int array (N,) will run PaSTA-NS.

Parameters:

  • x (ndarray (N,)) – Spatial map data to evaluate association. Can contain missing values such as NaN and Inf.

  • y (ndarray (N,)) – Spatial map data to evaluate association. Can contain missing values such as NaN and Inf.

  • coord (ndarray (N, 3) or None) – Spatial coordinates for observations. When unknown and left as None, the function requires D to run PaSTA, and D and dim to run PaSTA-NS.

  • D (ndarray (N, N) or None) – Distance matrix. When left as None, computed from coord.

  • dim (int or None) – Spatial dimension of data. When left as None, computed as coord.shape[1] if needed.

  • M (int or None) – Number of lag distances to evaluate when estimating variogram — important hyperparameter that determines the quality of variogram estimation, large values preferred. When set to None, use 3*sqrt(N) as default.

  • qd (float (0, 1]) – Determine the coverage of lag distances evaluated in variogram, with maximum distance evaluated being qd*np.max(D) — important hyperparameter that determines the quality of variogram estimation, large values preferred. Default 0.7.

  • xparc (None, 'auto', or ndarray (N,)) – Parcellation setting for map x. If ndarray, index should begin from 0, i.e., 0 to Np - 1 if Np parcels specified.

  • yparc (None, 'auto', or ndarray (N,)) – Parcellation setting for map y. If ndarray, index should begin from 0, i.e., 0 to Np - 1 if Np parcels specified.

  • max_clusters (int) – Maximum number of parcellations allowed in PaSTA-NS.

  • min_clusters (int) – Minimum number of parcellations allowed in PaSTA-NS.

  • min_cluster_size (int) – Minimum size of parcellations (# observations per parcel).

  • M_cluster (int or None) – Number of lag distances to evaluate in PaSTA-NS parcels when estimating their variograms. When set to None, default to 3*sqrt(Np), where Np is the number of observations in each parcel.

  • nugget (bool) – Indicator of whether use nugget in variogram models or not. Default True because nugget helps with discontinuity at short distances, and setting to False can result in problems especially when data are nonstationary.

Returns:

  • pef (float) – Significance p-values based on PaSTA/PaSTA-NS.

  • rX (float) – Pearson correlation coefficient between x and y.

  • nef (float) – Effective sample size estimated.

  • run_status (int) – 1 indicates successful run, and 0 indicates unsuccessful run such as when nef < 2 and data are too smooth to infer significance.

  • n_parc (ndarray (2,)) – [xn_parc, yn_parc] that indicates the number of parcels for each map in PaSTA-NS.

  • p_naive (float) – Significance with independence assumption and without controlling for autocorrelation.

  • fc_para1 (ndarray) – Covariance matrix for map x, with Nvx indicating the number of valid (finite value) observations in map x.

  • fc_para2 (ndarray) – Covariance matrix for map y, with Nvy indicating the number of valid (finite value) observations in map y.

pasta.pasta.estimate_variogram(D, data, M: int, qd: float)[source]

Estimate the empirical variogram from distance matrix between vertices, and data value at each vertex. Estimation performed in M bins, ranging from min_distance to qd * max_distance, where min_distance and max_distance are the min and max distance in the distance matrix.

Parameters:

  • D (ndarray (N, N)) – Distance matrix between all vertices.

  • data (ndarray (N,)) – Data value at each vertex.

  • M (int) – Number of bins to estimate variogram.

  • qd (float) – Determine the maximum distance to evaluate variogram.

Returns:

  • v (ndarray (M,)) – Estimated variogram values, i.e., semivariance.

  • h (ndarray (M,)) – Lag distances.

Notes

This is similar to variogram estimation in BrainSMASH but determining the max distance evaluated in a different way.

pasta.pasta.fit_covariance_blocks(x, D, n_clusters, point_cluster_idx, M_cluster, qd, nugget, exponent)[source]

Fit variogram model for each parcel and compute the diagonal blocks of nonstationary covariance matrix.

Parameters:

  • x (ndarray (N,)) – Spatial map data to evaluate association. All values are valid.

  • D (ndarray (N, N)) – Distance matrix.

  • n_clusters (int) – Number of parcels.

  • point_cluster_idx (ndarray (N,)) – Int array specifying parcellation settings for map x, ranging from 0 to NP-1 if NP parcels.

  • M_cluster (int or None) – Number of lag distances to evaluate in parcel when estimating their variograms. When set to None, default to 3*sqrt(Np), where Np is the number of observations in each parcel.

  • qd (float (0, 1])

  • nugget (bool) – Indicator of whether use nugget in variogram models or not.

  • exponent (float) – Shape parameter estimated using global stationary variogram model. This will be kept the same across parcels to obtain valid nonstationary covariance expression (i.e., PSD matrix).

Returns:

  • c_para (ndarray (N, N)) – Covariance matrix for map x, where within parcel covariance are estimated but cross-parcel elements are set to 0.

  • b (ndarray (n_clusters, 4)) – Stable variogram model parameters, each row corresponds a parcel.

pasta.pasta.fit_variogram(h, v, D, PrecomputedVariance=None, nugget: bool = True)[source]

Fit a stable variogram model to an empirical variogram.

Parameters:

  • h ((M,) ndarray) – Empirical lag distances.

  • v ((M,) ndarray) – Empirical semivariance evaluated at lag distances h.

  • D ((N, N) ndarray) – Pairwise distance matrix between spatial locations.

  • PrecomputedVariance (float or None, optional) – Precomputed sill (total variance) as initial guess for optimization. If None, the sill is estimated as the maximum value of v.

  • nugget (bool, default=True) – Whether to include a nugget term in the fitted model.

Returns:

  • c_para ((N, N) ndarray) – Fitted covariance matrix derived from the stable variogram model.

  • b ((4,) ndarray) – Estimated stable model parameters in the order (sill, range, exponent, nugget).

  • f (callable) – Variogram function. f(h) returns the semivariance at lag distance h.

  • fcov (callable) – Covariance function. fcov(h) returns the covariance at lag distance h.

Notes

The fitted model follows a stable variogram parameterization. The nugget term should be included for better fit.

pasta.pasta.fit_variogram_fixed_exponent(h, v, D, exponent, PrecomputedVariance=None, nugget: bool = True)[source]

Same as fit_variogram, but for stable model with predetermined range parameter.

Parameters:

  • h (ndarray)

  • v (ndarray)

  • D (ndarray)

  • exponent (float)

  • PrecomputedVariance (float or None)

  • nugget (bool)

Returns:

  • c_para (ndarray)

  • b (ndarray)

  • f (callable)

  • fcov (callable)

pasta.pasta.nef2p(rX, nef)[source]

Infer statistical significance p-value from test statistics rX and effective sample size nef.

Parameters:

  • rX (float) – Test statistic.

  • nef (float) – Effective sample size.

Returns:

p – Statistical significance p-value.

Return type:

float

pasta.pasta.parc_data(parc, c_para, b, D, coord, max_clusters, min_clusters, min_cluster_size, map_idx)[source]

parcellate data depending on setting to account for nonstationarity

3 scenarios: 1. parc is None: do not parcellate and return covariance matrix c_para as is (new variable name fc_para) 2. parc is string ‘auto’: determine the number of parcels based on estiamted range and shape parameter from stable variogram model (i.e., b[1] and b[2]), and parcellate data using spatial clustering 3. parc is user specified np int array with shape (M,) with each int indicating a unique parcel: return parc as is (new variable name parc_out), raise warning if risk of over-parcellation (compared to ‘auto’)

Parameters:

  • parc (either None, 'auto', or (N,)) – specifying the setting of parcellation

  • c_para (covariance matrix estimated from PaSTA, i.e., without parcel)

  • b (stable variogram model parameters estimated from PaSTA)

  • D (distance matrix of data (N, N))

  • coord ((N, 3) spatial coordinates of data or None) – If None, spatial clustering is conducted based on the distance matrix D with KMedoids

  • max_clusters (maximum number of parcellation, set to avoid over-parcellation at weak autocorrelation, e.g., spatial independence)

  • min_clusters (minimum number of parcellation, set to 1 will allow PaSTA-NS to collapse to PaSTA. This was used to mandate parcellation and test difference between PaSTA and PaSTA-NS in the manuscript.)

  • min_cluster_size (set to avoid too small parcellations, we set to 500 in fsaverage5 mesh with 10k vertices)

  • map_idx (used to identify the map evaluated when raise warning)

Returns:

  • parc_out (None when there is no subdivision of parcels, or (N,) of int where each unique int indicate a parcel)

  • n_parc (number of parcels)

  • unique_parcs (index of unique parcels)

  • fc_para (covariance matrix, either as c_para (when no parcellation) or zeros (parcellated))

pasta.pasta.process_convolution_crossblocks(c_para, b, x, D, n_clusters, point_cluster_idx, dim, exponent)[source]

Process convolution to infer the cross-parcel covariance of nonstationary covariance matrix.

Parameters:

  • c_para (ndarray (N, N)) – Covariance matrix output from fit_covariance_blocks, where covariances are estimated for within parcel pairs but not cross-parcel.

  • b (ndarray (n_clusters, 4)) – Fitted stable variogram model parameters, each row per parcel.

  • x (ndarray (N,)) – Map data.

  • D (ndarray (N, N)) – Distance matrix.

  • n_clusters (int) – Number of parcels.

  • point_cluster_idx (ndarray (N,)) – Int array starting from 0 indicating the membership of each point to parcels.

  • dim (int) – Spatial dimension of the data.

  • exponent (float) – Shape parameter fitted using the global stationary variogram, kept the same for valid PSD covariance matrix.

Returns:

c_para – Nonstationary covariance matrix by process convolution.

Return type:

ndarray (N, N)

pasta.pasta.stable_covariance_func(h, b)[source]

Covariance function based on stable variogram model for observations with distance h.

Equivalent to: (sill + nugget) - (sill * (1 - exp(-(h / range)**shape)) + nugget) = sill * exp(-(h / range)**shape) for h > 0.

When h == 0, set to sill + nugget, which is not computed here.

Parameters:

  • h (float or ndarray) – Lag distance at which to compute covariance.

  • b (ndarray) – Parameters for stable models.

Returns:

Covariance at distance h.

Return type:

float or ndarray

pasta.pasta.stable_variogram(h, b1, b2, b3, b4)[source]

Stable variogram model without nugget, defined as: semivariance = sill * (1 - exp(-(h / range)**shape))

Parameters:

  • h (float or ndarray) – Lag distance to be evaluated.

  • b1 (float) – Sill.

  • b2 (float) – Range parameter.

  • b3 (float) – Shape.

Returns:

Semivariance at lag distance h.

Return type:

float or ndarray

pasta.pasta.stable_variogram_fixed_exp(h, b1, b2, b3, fixed_exp)[source]

Stable variogram with prespecified shape parameter, with nugget.

Parameters:

  • h (float or ndarray)

  • b1 (float)

  • b2 (float)

  • b3 (float)

  • fixed_exp (float)

Return type:

float or ndarray

pasta.pasta.stable_variogram_fixed_exp_no_nugget(h, b1, b2, fixed_exp)[source]

Stable variogram with prespecified shape parameter, without nugget (i.e., nugget set to 0).

Parameters:

  • h (float or ndarray)

  • b1 (float)

  • b2 (float)

  • fixed_exp (float)

Return type:

float or ndarray

pasta.pasta.stable_variogram_no_nugget(h, b1, b2, b3)[source]

stable variogram model without nugget, defined as: semivariance = sill * (1-exp(-(h/range)**shape))

Parameters:

  • h (float or ndarray) – lag distance to be evaluated

  • b1 (float) – sill

  • b2 (float) – range parameter

  • b3 (float) – shape

Returns:

senuvaruabce at lag distance h

Return type:

float or ndarray

Module contents