Examples

Sign-flip vs permutation

The class extendedmosaicperm.factor.ExtendMosaicFactorTest adds a sign_flipping flag on top of the original mosaicperm.factor.MosaicFactorTest.

import numpy as np
from mosaicperm import statistics
from extendedmosaicperm.factor import ExtendMosaicFactorTest

T, N, K = 60, 12, 3
rng = np.random.default_rng(1)
F = rng.standard_normal((T, K))
B = rng.standard_normal((K, N))
Y = F @ B + rng.standard_normal((T, N))
exposures = B.T

# permutation-based inference
m_perm = ExtendMosaicFactorTest(
    outcomes=Y,
    exposures=exposures,
    test_stat=statistics.mean_maxcorr_stat,
    sign_flipping=False,
    seed=0,
)
m_perm.fit(nrand=200, verbose=False)

# sign-flip inference
m_sign = ExtendMosaicFactorTest(
    outcomes=Y,
    exposures=exposures,
    test_stat=statistics.mean_maxcorr_stat,
    sign_flipping=True,
    seed=0,
)
m_sign.fit(nrand=200, verbose=False)

print("perm p-value:", m_perm.pval)
print("sign p-value:", m_sign.pval)

Ridge-regularized residuals

You can switch from OLS-based residuals to RidgeCV with a simple ridge_alphas argument.

import numpy as np
from mosaicperm import statistics
from extendedmosaicperm.factor import ExtendMosaicFactorTest

T, N, K = 60, 12, 3
rng = np.random.default_rng(2)
F = rng.standard_normal((T, K))
B = rng.standard_normal((K, N))
Y = F @ B + rng.standard_normal((T, N))
exposures = B.T

alphas = np.logspace(-3, 3, 5)

m_ridge = ExtendMosaicFactorTest(
    outcomes=Y,
    exposures=exposures,
    test_stat=statistics.mean_maxcorr_stat,
    ridge_alphas=alphas,
    sign_flipping=False,
    seed=0,
)
m_ridge.fit(nrand=200, verbose=False)
print("ridge p-value:", m_ridge.pval)

Adaptive tiling

Instead of manually specifying tiles, you can let the test build them from residual covariance via the adaptive_tiling flag.

import numpy as np
from mosaicperm import statistics
from extendedmosaicperm.factor import ExtendMosaicFactorTest

T, N, K = 60, 12, 3
rng = np.random.default_rng(3)
F = rng.standard_normal((T, K))
B = rng.standard_normal((K, N))
Y = F @ B + rng.standard_normal((T, N))
exposures = B.T

m_adapt = ExtendMosaicFactorTest(
    outcomes=Y,
    exposures=exposures,
    adaptive_tiling=True,
    test_stat=statistics.mean_maxcorr_stat,
    sign_flipping=True,
    seed=0,
)
m_adapt.fit(nrand=200, verbose=False)
print("adaptive p-value:", m_adapt.pval)

Simulation experiments

The extendedmosaicperm.experiments submodule contains experiment drivers for Monte Carlo comparisons:

A minimal example:

from extendedmosaicperm.experiments.sign_flip import SignFlipExperiment

exp = SignFlipExperiment(
    n_sims=10,
    nrand=50,
    seed=0,
    violation_strengths=[0.0, 0.1],
)
exp.run()
df = exp.summarize()
print(df.head())