Skip to content
This repository has been archived by the owner on Apr 2, 2020. It is now read-only.

Commit

Permalink
Removed lots of old, now unused, code
Browse files Browse the repository at this point in the history
  • Loading branch information
@bthirion
bthirion committed Sep 5, 2018
1 parent 5c32ab3 commit 3be4fb7
Show file tree
Hide file tree
Showing 4 changed files with 18 additions and 294 deletions.
16 changes: 7 additions & 9 deletions nistats/contrasts.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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_]
Expand Down
211 changes: 1 addition & 210 deletions nistats/model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
# Inverse t cumulative distribution
inv_t_cdf = t_distribution.ppf


class LikelihoodModelResults(object):
''' Class to contain results from likelihood models '''

Expand Down Expand Up @@ -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 ('<T contrast: effect=%s, sd=%s, t=%s, df_den=%d>' %
(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 '<F contrast: F=%s, df_den=%d, df_num=%d>' % \
(repr(self.F), self.df_den, self.df_num)
12 changes: 10 additions & 2 deletions nistats/tests/test_contrasts.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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():
Expand Down Expand Up @@ -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)
73 changes: 0 additions & 73 deletions nistats/tests/test_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)

0 comments on commit 3be4fb7

Please sign in to comment.