diff --git a/nistats/contrasts.py b/nistats/contrasts.py index fb5f80cd..bdb1a396 100644 --- a/nistats/contrasts.py +++ b/nistats/contrasts.py @@ -8,8 +8,9 @@ from warnings import warn import numpy as np import scipy.stats as sps +from scipy.linalg import sqrtm -from .utils import multiple_mahalanobis, z_score +from .utils import z_score DEF_TINY = 1e-50 @@ -54,18 +55,15 @@ def compute_contrast(labels, regression_result, con_val, contrast_type=None): '"{0}" is not a known contrast type. Allowed types are {1}'. format(contrast_type, acceptable_contrast_types)) + var_ = np.zeros((1, 1, labels.size)) + effect_ = np.zeros((dim, labels.size)) if contrast_type == 't': - effect_ = np.zeros((1, labels.size)) - var_ = np.zeros(labels.size) for label_ in regression_result: label_mask = labels == label_ - resl = regression_result[label_].Tcontrast(con_val) - effect_[:, label_mask] = resl.effect.T - var_[label_mask] = (resl.sd ** 2).T + reg = regression_result[label_] + effect_[:, label_mask] = con_val.dot(reg.theta) + var_[:, :, label_mask] = reg.vcov(matrix=np.atleast_2d(con_val)) elif contrast_type == 'F': - from scipy.linalg import sqrtm - effect_ = np.zeros((dim, labels.size)) - var_ = np.zeros(labels.size) for label_ in regression_result: label_mask = labels == label_ reg = regression_result[label_] diff --git a/nistats/model.py b/nistats/model.py index 007eead0..411ac1f9 100644 --- a/nistats/model.py +++ b/nistats/model.py @@ -17,6 +17,7 @@ # Inverse t cumulative distribution inv_t_cdf = t_distribution.ppf + class LikelihoodModelResults(object): ''' Class to contain results from likelihood models ''' @@ -156,213 +157,3 @@ def vcov(self, matrix=None, column=None, dispersion=None, other=None): return tmp[:, :, np.newaxis] * dispersion if matrix is None and column is None: return self.cov * dispersion - - def Tcontrast(self, matrix, store=('t', 'effect', 'sd'), dispersion=None): - """ Compute a Tcontrast for a row vector `matrix` - - To get the t-statistic for a single column, use the 't' method. - - Parameters - ---------- - matrix : 1D array-like - contrast matrix - - store : sequence, optional - components of t to store in results output object. Defaults to all - components ('t', 'effect', 'sd'). - - dispersion : None or float, optional - - Returns - ------- - res : ``TContrastResults`` object - """ - matrix = np.asarray(matrix) - # 1D vectors assumed to be row vector - if matrix.ndim == 1: - matrix = matrix[None] - if matrix.shape[0] != 1: - raise ValueError("t contrasts should have only one row") - if matrix.shape[1] != self.theta.shape[0]: - raise ValueError("t contrasts should be length P=%d, " - "but this is length %d" % (self.theta.shape[0], - matrix.shape[1])) - store = set(store) - if not store.issubset(('t', 'effect', 'sd')): - raise ValueError('Unexpected store request in %s' % store) - st_t = st_effect = st_sd = effect = sd = None - if 't' in store or 'effect' in store: - effect = np.dot(matrix, self.theta) - if 'effect' in store: - st_effect = np.squeeze(effect) - if 't' in store or 'sd' in store: - sd = np.sqrt(self.vcov(matrix=matrix, dispersion=dispersion)) - if 'sd' in store: - st_sd = np.squeeze(sd) - if 't' in store: - st_t = np.squeeze(effect * pos_recipr(sd)) - return TContrastResults(effect=st_effect, t=st_t, sd=st_sd, - df_den=self.df_resid) - - def Fcontrast(self, matrix, dispersion=None, invcov=None): - """ Compute an Fcontrast for a contrast matrix `matrix`. - - Here, `matrix` M is assumed to be non-singular. More precisely - - .. math:: - - M pX pX' M' - - is assumed invertible. Here, :math:`pX` is the generalized inverse of - the design matrix of the model. There can be problems in non-OLS models - where the rank of the covariance of the noise is not full. - - See the contrast module to see how to specify contrasts. In particular, - the matrices from these contrasts will always be non-singular in the - sense above. - - Parameters - ---------- - matrix : 1D array-like - contrast matrix - - dispersion : None or float, optional - If None, use ``self.dispersion`` - - invcov : None or array, optional - Known inverse of variance covariance matrix. - If None, calculate this matrix. - - Returns - ------- - f_res : ``FContrastResults`` instance - with attributes F, df_den, df_num - - Notes - ----- - For F contrasts, we now specify an effect and covariance - """ - matrix = np.asarray(matrix) - # 1D vectors assumed to be row vector - if matrix.ndim == 1: - matrix = matrix[None] - if matrix.shape[1] != self.theta.shape[0]: - raise ValueError("F contrasts should have shape[1] P=%d, " - "but this has shape[1] %d" % (self.theta.shape[0], - matrix.shape[1])) - ctheta = np.dot(matrix, self.theta) - if matrix.ndim == 1: - matrix = matrix.reshape((1, matrix.shape[0])) - if dispersion is None: - dispersion = self.dispersion - q = matrix.shape[0] - if invcov is None: - invcov = inv(self.vcov(matrix=matrix, dispersion=1.0)) - F = np.add.reduce(np.dot(invcov, ctheta) * ctheta, 0) *\ - pos_recipr((q * dispersion)) - F = np.squeeze(F) - return FContrastResults( - effect=ctheta, covariance=self.vcov( - matrix=matrix, dispersion=dispersion[np.newaxis]), - F=F, df_den=self.df_resid, df_num=invcov.shape[0]) - - def conf_int(self, alpha=.05, cols=None, dispersion=None): - ''' The confidence interval of the specified theta estimates. - - Parameters - ---------- - alpha : float, optional - The `alpha` level for the confidence interval. - ie., `alpha` = .05 returns a 95% confidence interval. - - cols : tuple, optional - `cols` specifies which confidence intervals to return - - dispersion : None or scalar - scale factor for the variance / covariance (see class docstring and - ``vcov`` method docstring) - - Returns - ------- - cis : ndarray - `cis` is shape ``(len(cols), 2)`` where each row contains [lower, - upper] for the given entry in `cols` - - Examples - -------- - >>> from numpy.random import standard_normal as stan - >>> from nistats.regression import OLSModel - >>> x = np.hstack((stan((30,1)),stan((30,1)),stan((30,1)))) - >>> beta=np.array([3.25, 1.5, 7.0]) - >>> y = np.dot(x,beta) + stan((30)) - >>> model = OLSModel(x).fit(y) - >>> confidence_intervals = model.conf_int(cols=(1,2)) - - Notes - ----- - Confidence intervals are two-tailed. - TODO: - tails : string, optional - `tails` can be "two", "upper", or "lower" - ''' - if cols is None: - lower = self.theta - inv_t_cdf(1 - alpha / 2, self.df_resid) *\ - np.sqrt(np.diag(self.vcov(dispersion=dispersion))) - upper = self.theta + inv_t_cdf(1 - alpha / 2, self.df_resid) *\ - np.sqrt(np.diag(self.vcov(dispersion=dispersion))) - else: - lower, upper = [], [] - for i in cols: - lower.append( - self.theta[i] - inv_t_cdf(1 - alpha / 2, self.df_resid) * - np.sqrt(self.vcov(column=i, dispersion=dispersion))) - upper.append( - self.theta[i] + inv_t_cdf(1 - alpha / 2, self.df_resid) * - np.sqrt(self.vcov(column=i, dispersion=dispersion))) - return np.asarray(list(zip(lower, upper))) - - -class TContrastResults(object): - """ Results from a t contrast of coefficients in a parametric model. - - The class does nothing, it is a container for the results from T contrasts, - and returns the T-statistics when np.asarray is called. - """ - - def __init__(self, t, sd, effect, df_den=None): - if df_den is None: - df_den = np.inf - self.t = t - self.sd = sd - self.effect = effect - self.df_den = df_den - - def __array__(self): - return np.asarray(self.t) - - def __str__(self): - return ('' % - (self.effect, self.sd, self.t, self.df_den)) - - -class FContrastResults(object): - """ Results from an F contrast of coefficients in a parametric model. - - The class does nothing, it is a container for the results from F contrasts, - and returns the F-statistics when np.asarray is called. - """ - - def __init__(self, effect, covariance, F, df_num, df_den=None): - if df_den is None: - df_den = np.inf - self.effect = effect - self.covariance = covariance - self.F = F - self.df_den = df_den - self.df_num = df_num - def __array__(self): - return np.asarray(self.F) - - def __str__(self): - return '' % \ - (repr(self.F), self.df_den, self.df_num) diff --git a/nistats/tests/test_contrasts.py b/nistats/tests/test_contrasts.py index 1bdc3234..f091561b 100644 --- a/nistats/tests/test_contrasts.py +++ b/nistats/tests/test_contrasts.py @@ -6,6 +6,7 @@ from nistats.contrasts import compute_contrast, _fixed_effect_contrast from numpy.testing import assert_almost_equal +from scipy.linalg import sqrtm def test_Tcontrast(): @@ -97,12 +98,19 @@ def test_contrast_values(): # t test cval = np.eye(q)[0] con = compute_contrast(lab, res, cval) - t_ref = list(res.values())[0].Tcontrast(cval).t + summary_stats = res.values()[0] + t_ref = cval.dot(summary_stats.theta) /\ + np.ravel(np.sqrt(summary_stats.vcov(matrix=np.atleast_2d(cval)))) assert_almost_equal(np.ravel(con.stat()), t_ref) # F test cval = np.eye(q)[:3] con = compute_contrast(lab, res, cval) - F_ref = list(res.values())[0].Fcontrast(cval).F + cbeta = np.atleast_2d(np.dot(cval, summary_stats.theta)) + invcov = np.linalg.inv( + summary_stats.vcov(matrix=cval, dispersion=1.)) + wcbeta = np.dot(sqrtm(invcov), cbeta) + F_ref = np.sum(wcbeta ** 2, 0) / 3 / summary_stats.dispersion + # F_ref = list(res.values())[0].Fcontrast(cval).F # Note that the values are not strictly equal, # this seems to be related to a bug in Mahalanobis assert_almost_equal(np.ravel(con.stat()), F_ref, 3) diff --git a/nistats/tests/test_model.py b/nistats/tests/test_model.py index d08aa3ea..5fd816fb 100644 --- a/nistats/tests/test_model.py +++ b/nistats/tests/test_model.py @@ -63,76 +63,3 @@ def test_model(): pcts = percentile(RESULTS.resid, [0, 25, 50, 75, 100]) assert_array_almost_equal(pcts, [-1.6970, -0.6667, 0, 0.6667, 1.6970], 4) - -def test_t_contrast(): - # Test indivudual t against R - assert_array_almost_equal(RESULTS.t(0), 3.25) - assert_array_almost_equal(RESULTS.t(1), 7.181, 3) - # And contrast - assert_array_almost_equal(RESULTS.Tcontrast([1, 0]).t, 3.25) - assert_array_almost_equal(RESULTS.Tcontrast([0, 1]).t, 7.181, 3) - # Input matrix checked for size - assert_raises(ValueError, RESULTS.Tcontrast, [1]) - assert_raises(ValueError, RESULTS.Tcontrast, [1, 0, 0]) - # And shape - assert_raises(ValueError, RESULTS.Tcontrast, np.array([1, 0])[:, None]) - - -def test_t_output(): - # Check we get required outputs - exp_t = RESULTS.t(0) - exp_effect = RESULTS.theta[0] - exp_sd = exp_effect / exp_t - res = RESULTS.Tcontrast([1, 0]) - assert_array_almost_equal(res.t, exp_t) - assert_array_almost_equal(res.effect, exp_effect) - assert_array_almost_equal(res.sd, exp_sd) - res = RESULTS.Tcontrast([1, 0], store=('effect',)) - assert_equal(res.t, None) - assert_array_almost_equal(res.effect, exp_effect) - assert_equal(res.sd, None) - res = RESULTS.Tcontrast([1, 0], store=('t',)) - assert_array_almost_equal(res.t, exp_t) - assert_equal(res.effect, None) - assert_equal(res.sd, None) - res = RESULTS.Tcontrast([1, 0], store=('sd',)) - assert_equal(res.t, None) - assert_equal(res.effect, None) - assert_array_almost_equal(res.sd, exp_sd) - res = RESULTS.Tcontrast([1, 0], store=('effect', 'sd')) - assert_equal(res.t, None) - assert_array_almost_equal(res.effect, exp_effect) - assert_array_almost_equal(res.sd, exp_sd) - - -def test_f_output(): - # Test f_output - res = RESULTS.Fcontrast([1, 0]) - exp_f = RESULTS.t(0) ** 2 - assert_array_almost_equal(exp_f, res.F) - # Test arrays work as well as lists - res = RESULTS.Fcontrast(np.array([1, 0])) - assert_array_almost_equal(exp_f, res.F) - # Test with matrix against R - res = RESULTS.Fcontrast(np.eye(2)) - assert_array_almost_equal(31.06, res.F, 2) - # Input matrix checked for size - assert_raises(ValueError, RESULTS.Fcontrast, [1]) - assert_raises(ValueError, RESULTS.Fcontrast, [1, 0, 0]) - # And shape - assert_raises(ValueError, RESULTS.Fcontrast, np.array([1, 0])[:, None]) - - -def test_f_output_new_api(): - res = RESULTS.Fcontrast([1, 0]) - assert_array_almost_equal(res.effect, RESULTS.theta[0]) - assert_array_almost_equal(res.covariance, RESULTS.vcov()[0][0]) - - -def test_conf_int(): - lower_, upper_ = RESULTS.conf_int() - assert_true((lower_ < upper_).all()) - assert_true((lower_ > upper_ - 10).all()) - lower_, upper_ = RESULTS.conf_int(cols=[1]).T - assert_true(lower_ < upper_) - assert_true(lower_ > upper_ - 10)