Python mass_univariate.permuted_ols函数代码示例

def test_permuted_ols_nocovar(random_state=0):
    rng = check_random_state(random_state)
    # design parameters
    n_samples = 50
    # create design
    target_var = rng.randn(n_samples, 1)
    tested_var = rng.randn(n_samples, 1)
    # compute t-scores with linalg or statsmodels
    ref_score = get_tvalue_with_alternative_library(tested_var, target_var)
    # permuted OLS
    _, own_score, _ = permuted_ols(
        tested_var, target_var, model_intercept=False,
        n_perm=0, random_state=random_state)
    assert_array_almost_equal(ref_score, own_score, decimal=6)

    # test with ravelized tested_var
    _, own_score, _ = permuted_ols(
        np.ravel(tested_var), target_var, model_intercept=False,
        n_perm=0, random_state=random_state)
    assert_array_almost_equal(ref_score, own_score, decimal=6)

    ### Adds intercept (should be equivalent to centering variates)
    # permuted OLS
    _, own_score_intercept, _ = permuted_ols(
        tested_var, target_var, model_intercept=True,
        n_perm=0, random_state=random_state)
    target_var -= target_var.mean(0)
    tested_var -= tested_var.mean(0)
    # compute t-scores with linalg or statsmodels
    ref_score_intercept = get_tvalue_with_alternative_library(
        tested_var, target_var, np.ones((n_samples, 1)))
    assert_array_almost_equal(ref_score_intercept, own_score_intercept,
                              decimal=6)

def test_permuted_ols_intercept_statsmodels_withcovar(random_state=0):
    """

    This test has a statsmodels dependance. There seems to be no simple,
    alternative way to perform a F-test on a linear model including
    covariates.

    """
    try:
        from statsmodels.regression.linear_model import OLS
    except:
        warnings.warn("Statsmodels is required to run this test")
        raise nose.SkipTest

    rng = check_random_state(random_state)
    # design parameters
    n_samples = 50
    # create design
    target_var = rng.randn(n_samples, 1)
    tested_var = np.ones((n_samples, 1))
    confounding_vars = rng.randn(n_samples, 2)
    # statsmodels OLS
    ols = OLS(target_var, np.hstack((tested_var, confounding_vars))).fit()
    fvals = ols.f_test([[1.0, 0.0, 0.0]]).fvalue
    # permuted OLS
    _, orig_scores, _ = permuted_ols(tested_var, target_var, confounding_vars, n_perm=0, random_state=random_state)
    # same thing but with model_intercept=True to check it has no effect
    _, orig_scores_addintercept, _ = permuted_ols(
        tested_var, target_var, confounding_vars, model_intercept=True, n_perm=0, random_state=random_state
    )
    assert_array_almost_equal(fvals, orig_scores, decimal=6)
    assert_array_almost_equal(orig_scores, orig_scores_addintercept, decimal=6)

def test_permuted_ols_withcovar(random_state=0):
    """

    """
    rng = check_random_state(random_state)
    # design parameters
    n_samples = 50
    # create design
    target_var = rng.randn(n_samples, 1)
    tested_var = rng.randn(n_samples, 1)
    confounding_vars = rng.randn(n_samples, 2)
    # compute t-scores with linalg or statsmodels
    ref_score = get_tvalue_with_alternative_library(tested_var, target_var,
                                                    confounding_vars)
    # permuted OLS
    _, own_score, _ = permuted_ols(
        tested_var, target_var, confounding_vars, model_intercept=False,
        n_perm=0, random_state=random_state)
    assert_array_almost_equal(ref_score, own_score, decimal=6)

    ### Adds intercept
    # permuted OLS
    _, own_scores_intercept, _ = permuted_ols(
        tested_var, target_var, confounding_vars, model_intercept=True,
        n_perm=0, random_state=random_state)
    # compute t-scores with linalg or statsmodels
    confounding_vars = np.hstack((confounding_vars, np.ones((n_samples, 1))))
    alt_score_intercept = get_tvalue_with_alternative_library(
        tested_var, target_var, confounding_vars)
    assert_array_almost_equal(alt_score_intercept, own_scores_intercept,
                              decimal=6)

def test_permuted_ols_nocovar_multivariate(random_state=0):
    """Test permuted_ols with multiple tested variates and no covariate.

    It is equivalent to fitting several models with only one tested variate.

    """
    rng = check_random_state(random_state)
    # design parameters
    n_samples = 50
    n_targets = 10
    n_regressors = 2
    # create design
    target_vars = rng.randn(n_samples, n_targets)
    tested_var = rng.randn(n_samples, n_regressors)
    # compute t-scores with linalg or statsmodels
    ref_scores = get_tvalue_with_alternative_library(tested_var, target_vars)
    # permuted OLS
    _, own_scores, _ = permuted_ols(
        tested_var, target_vars, model_intercept=False,
        n_perm=0, random_state=random_state)
    assert_array_almost_equal(ref_scores, own_scores, decimal=6)

    ### Adds intercept (should be equivalent to centering variates)
    # permuted OLS
    _, own_scores_intercept, _ = permuted_ols(
        tested_var, target_vars, model_intercept=True,
        n_perm=0, random_state=random_state)
    target_vars -= target_vars.mean(0)
    tested_var -= tested_var.mean(0)
    # compute t-scores with linalg or statsmodels
    ref_scores_intercept = get_tvalue_with_alternative_library(
            tested_var, target_vars, np.ones((n_samples, 1)))
    assert_array_almost_equal(ref_scores_intercept, own_scores_intercept,
                              decimal=6)

def test_permuted_ols_sklearn_nocovar(random_state=0):
    rng = check_random_state(random_state)
    # design parameters
    n_samples = 50
    # create design
    target_var = rng.randn(n_samples, 1)
    tested_var = rng.randn(n_samples, 1)
    # scikit-learn F-score
    fvals, _ = f_regression(target_var, tested_var, center=False)
    # permuted OLS
    _, orig_scores, _ = permuted_ols(tested_var, target_var, model_intercept=False, n_perm=0, random_state=random_state)
    assert_array_almost_equal([fvals], orig_scores, decimal=6)

    # test with ravelized tested_var
    _, orig_scores, _ = permuted_ols(
        np.ravel(tested_var), target_var, model_intercept=False, n_perm=0, random_state=random_state
    )
    assert_array_almost_equal([fvals], orig_scores, decimal=6)

    ### Adds intercept (should be equivalent to centering variates)
    # permuted OLS
    _, orig_scores_addintercept, _ = permuted_ols(
        tested_var, target_var, model_intercept=True, n_perm=0, random_state=random_state
    )
    target_var -= target_var.mean(0)
    tested_var -= tested_var.mean(0)
    # scikit-learn F-score
    fvals_addintercept, _ = f_regression(target_var, tested_var, center=True)
    assert_array_almost_equal([fvals_addintercept], orig_scores_addintercept, decimal=6)

def test_permuted_ols_withcovar_multivariate(random_state=0):
    """Test permuted_ols with multiple tested variates and covariates.

    It is equivalent to fitting several models with only one tested variate.

    """
    rng = check_random_state(random_state)
    # design parameters
    n_samples = 50
    n_targets = 10
    n_covars = 2
    # create design
    target_vars = rng.randn(n_samples, n_targets)
    tested_var = rng.randn(n_samples, 1)
    confounding_vars = rng.randn(n_samples, n_covars)
    # compute t-scores with linalg or statmodels
    ref_scores = get_tvalue_with_alternative_library(tested_var, target_vars,
                                                     confounding_vars)
    # permuted OLS
    _, own_scores, _ = permuted_ols(
        tested_var, target_vars, confounding_vars, model_intercept=False,
        n_perm=0, random_state=random_state)
    assert_almost_equal(ref_scores, own_scores, decimal=6)

    ### Adds intercept
    # permuted OLS
    _, own_scores_intercept, _ = permuted_ols(
        tested_var, target_vars, confounding_vars, model_intercept=True,
        n_perm=0, random_state=random_state)
    # compute t-scores with linalg or statmodels
    confounding_vars = np.hstack((confounding_vars, np.ones((n_samples, 1))))
    ref_scores_intercept = get_tvalue_with_alternative_library(
        tested_var, target_vars, confounding_vars)
    assert_array_almost_equal(ref_scores_intercept,
                              own_scores_intercept, decimal=6)

def test_sided_test2(random_state=0):
    """Check that two-sided can actually recover positive and negative effects.
    """
    # create design
    target_var1 = np.arange(0, 10).reshape((-1, 1))  # positive effect
    target_var = np.hstack((target_var1, - target_var1))

    tested_var = np.arange(0, 20, 2)
    # permuted OLS
    # one-sided
    neg_log_pvals_onesided, _, _ = permuted_ols(
        tested_var, target_var, model_intercept=False,
        two_sided_test=False, n_perm=100, random_state=random_state)
    # one-sided (other side)
    neg_log_pvals_onesided2, _, _ = permuted_ols(
        tested_var, -target_var, model_intercept=False,
        two_sided_test=False, n_perm=100, random_state=random_state)
    # two-sdided
    neg_log_pvals_twosided, _, _ = permuted_ols(
        tested_var, target_var, model_intercept=False,
        two_sided_test=True, n_perm=100, random_state=random_state)

    assert_array_almost_equal(neg_log_pvals_onesided[0],
                              neg_log_pvals_onesided2[0][::-1])
    assert_array_almost_equal(neg_log_pvals_onesided + neg_log_pvals_onesided2,
                              neg_log_pvals_twosided)

def test_permuted_ols_statsmodels_withcovar_multivariate(random_state=0):
    """Test permuted_ols with multiple tested variates and covariates.

    It is equivalent to fitting several models with only one tested variate.

    This test has a statsmodels dependance. There seems to be no simple,
    alternative way to perform a F-test on a linear model including
    covariates.

    """
    try:
        from statsmodels.regression.linear_model import OLS
    except:
        warnings.warn("Statsmodels is required to run this test")
        raise nose.SkipTest

    rng = check_random_state(random_state)
    # design parameters
    n_samples = 50
    n_targets = 10
    n_covars = 2
    # create design
    target_vars = rng.randn(n_samples, n_targets)
    tested_var = rng.randn(n_samples, 1)
    confounding_vars = rng.randn(n_samples, n_covars)
    # statsmodels OLS
    fvals = np.empty((n_targets, 1))
    test_matrix = np.array([[1.0] + [0.0] * n_covars])
    for i in range(n_targets):
        ols = OLS(target_vars[:, i], np.hstack((tested_var, confounding_vars)))
        fvals[i] = ols.fit().f_test(test_matrix).fvalue[0][0]
    # permuted OLS
    _, orig_scores, _ = permuted_ols(
        tested_var, target_vars, confounding_vars, model_intercept=False, n_perm=0, random_state=random_state
    )
    assert_almost_equal(fvals, orig_scores, decimal=6)

    ### Adds intercept
    # permuted OLS
    _, orig_scores_addintercept, _ = permuted_ols(
        tested_var, target_vars, confounding_vars, model_intercept=True, n_perm=0, random_state=random_state
    )
    # statsmodels OLS
    confounding_vars = np.hstack((confounding_vars, np.ones((n_samples, 1))))
    fvals_addintercept = np.empty((n_targets, 1))
    test_matrix = np.array([[1.0] + [0.0] * (n_covars + 1)])
    for i in range(n_targets):
        ols = OLS(target_vars[:, i], np.hstack((tested_var, confounding_vars)))
        fvals_addintercept[i] = ols.fit().f_test(test_matrix).fvalue[0][0]
    assert_array_almost_equal(fvals_addintercept, orig_scores_addintercept, decimal=6)

def test_permuted_ols_check_h0_noeffect_signswap(random_state=0):
    rng = check_random_state(random_state)
    # design parameters
    n_samples = 100
    # create dummy design with no effect
    target_var = rng.randn(n_samples, 1)
    tested_var = np.ones((n_samples, 1))
    # permuted OLS
    # We check that h0 is close to the theoretical distribution, which is
    # known for this simple design (= t(n_samples - dof)).
    perm_ranges = [10, 100, 1000]  # test various number of permutations
    all_kstest_pvals = []
    # we compute the Mean Squared Error between cumulative Density Function
    # as a proof of consistency of the permutation algorithm
    all_mse = []
    for i, n_perm in enumerate(np.repeat(perm_ranges, 10)):
        pval, orig_scores, h0 = permuted_ols(
            tested_var, target_var, model_intercept=False,
            n_perm=n_perm, two_sided_test=False, random_state=i)
        assert_equal(h0.size, n_perm)
        # Kolmogorov-Smirnov test
        kstest_pval = stats.kstest(h0, stats.t(n_samples).cdf)[1]
        all_kstest_pvals.append(kstest_pval)
        mse = np.mean(
            (stats.t(n_samples).cdf(np.sort(h0))
             - np.linspace(0, 1, h0.size + 1)[1:]) ** 2)
        all_mse.append(mse)
    all_kstest_pvals = np.array(all_kstest_pvals).reshape(
        (len(perm_ranges), -1))
    all_mse = np.array(all_mse).reshape((len(perm_ranges), -1))
    # check that a difference between distributions is not rejected by KS test
    assert_array_less(0.01 / (len(perm_ranges) * 10.), all_kstest_pvals)
    # consistency of the algorithm: the more permutations, the less the MSE
    assert_array_less(np.diff(all_mse.mean(1)), 0)

def test_sided_test(random_state=0):
    rng = check_random_state(random_state)
    # design parameters
    n_samples = 50
    # create design
    target_var = rng.randn(n_samples, 100)
    tested_var = rng.randn(n_samples, 1)
    # permuted OLS
    # one-sided
    neg_log_pvals_onesided, _, _ = permuted_ols(
        tested_var, target_var, model_intercept=False,
        two_sided_test=False, n_perm=100, random_state=random_state)
    # two-sdided
    neg_log_pvals_twosided, _, _ = permuted_ols(
        tested_var, target_var, model_intercept=False,
        two_sided_test=True, n_perm=100, random_state=random_state)
    assert_equal(
        np.sum(neg_log_pvals_twosided - neg_log_pvals_onesided > 0), 0)

def test_permuted_ols_intercept_sklearn_nocovar(random_state=0):
    rng = check_random_state(random_state)
    # design parameters
    n_samples = 50
    # create design
    target_var = rng.randn(n_samples, 1)
    tested_var = np.ones((n_samples, 1))
    # scikit-learn F-score
    fvals, _ = f_regression(target_var, tested_var, center=False)
    # permuted OLS
    neg_log_pvals, orig_scores, _ = permuted_ols(
        tested_var, target_var, confounding_vars=None, n_perm=10, random_state=random_state
    )
    assert_array_less(neg_log_pvals, 1.0)  # ensure sign swap is correctly done
    # same thing but with model_intercept=True to check it has no effect
    _, orig_scores_addintercept, _ = permuted_ols(
        tested_var, target_var, confounding_vars=None, model_intercept=True, n_perm=0, random_state=random_state
    )
    assert_array_almost_equal([fvals], orig_scores, decimal=6)
    assert_array_almost_equal(orig_scores, orig_scores_addintercept, decimal=6)

def test_permuted_ols_intercept_nocovar_multivariate(random_state=0):
    rng = check_random_state(random_state)
    # design parameters
    n_samples = 50
    n_targets = 10
    # create design
    target_vars = rng.randn(n_samples, n_targets)
    tested_vars = np.ones((n_samples, 1))
    # compute t-scores with nilearn routine
    ref_scores = get_tvalue_with_alternative_library(tested_vars, target_vars)
    # permuted OLS
    _, own_scores, _ = permuted_ols(
        tested_vars, target_vars, confounding_vars=None, n_perm=0,
        random_state=random_state)
    # same thing but with model_intercept=True to check it has no effect
    _, own_scores_intercept, _ = permuted_ols(
        tested_vars, target_vars, confounding_vars=None, model_intercept=True,
        n_perm=0, random_state=random_state)
    assert_array_almost_equal(ref_scores, own_scores, decimal=6)
    assert_array_almost_equal(own_scores, own_scores_intercept, decimal=6)

def test_permuted_ols_intercept_nocovar(random_state=0):
    rng = check_random_state(random_state)
    # design parameters
    n_samples = 50
    # create design
    target_var = rng.randn(n_samples, 1)
    tested_var = np.ones((n_samples, 1))
    # compute t-scores with linalg or statmodels
    t_val_ref = get_tvalue_with_alternative_library(tested_var, target_var)
    # permuted OLS
    neg_log_pvals, orig_scores, _ = permuted_ols(
        tested_var, target_var, confounding_vars=None, n_perm=10,
        random_state=random_state)
    assert_array_less(neg_log_pvals, 1.)  # ensure sign swap is correctly done
    # same thing but with model_intercept=True to check it has no effect
    _, orig_scores_addintercept, _ = permuted_ols(
        tested_var, target_var, confounding_vars=None, model_intercept=True,
        n_perm=0, random_state=random_state)
    assert_array_almost_equal(t_val_ref, orig_scores, decimal=6)
    assert_array_almost_equal(orig_scores, orig_scores_addintercept, decimal=6)

def test_permuted_ols_intercept_statsmodels_withcovar(random_state=0):
    rng = check_random_state(random_state)
    # design parameters
    n_samples = 50
    # create design
    target_var = rng.randn(n_samples, 1)
    tested_var = np.ones((n_samples, 1))
    confounding_vars = rng.randn(n_samples, 2)
    # compute t-scores with linalg or statmodels
    ref_scores = get_tvalue_with_alternative_library(tested_var, target_var,
                                                     confounding_vars)
    # permuted OLS
    _, own_scores, _ = permuted_ols(
        tested_var, target_var, confounding_vars, n_perm=0,
        random_state=random_state)
    # same thing but with model_intercept=True to check it has no effect
    _, own_scores_intercept, _ = permuted_ols(
        tested_var, target_var, confounding_vars, model_intercept=True,
        n_perm=0, random_state=random_state)
    assert_array_almost_equal(ref_scores, own_scores, decimal=6)
    assert_array_almost_equal(ref_scores, own_scores_intercept, decimal=6)

def test_permuted_ols_intercept_sklearn_nocovar_multivariate(random_state=0):
    rng = check_random_state(random_state)
    # design parameters
    n_samples = 50
    n_targets = 10
    # create design
    target_vars = rng.randn(n_samples, n_targets)
    tested_var = np.ones((n_samples, 1))
    # scikit-learn F-scores
    fvals = np.empty((n_targets, 1))
    for i in range(n_targets):
        fvals[i], _ = f_regression(target_vars[:, i], tested_var, center=False)
    # permuted OLS
    _, orig_scores, _ = permuted_ols(
        tested_var, target_vars, confounding_vars=None, n_perm=0, random_state=random_state
    )
    # same thing but with model_intercept=True to check it has no effect
    _, orig_scores_addintercept, _ = permuted_ols(
        tested_var, target_vars, confounding_vars=None, model_intercept=True, n_perm=0, random_state=random_state
    )
    assert_array_almost_equal(fvals, orig_scores, decimal=6)
    assert_array_almost_equal(orig_scores, orig_scores_addintercept, decimal=6)