1""" Testing models module
2"""
3
4import numpy as np
5import pytest
6
7from numpy.testing import assert_array_almost_equal
8from nilearn.glm import OLSModel
9
10
11N = 10
12X = np.c_[np.linspace(- 1, 1, N), np.ones((N,))]
13Y = np.r_[range(5), range(1, 6)]
14MODEL = OLSModel(X)
15RESULTS = MODEL.fit(Y)
16
17""" R script
18
19::
20
21    X = cbind(0:9 * 2/9 -1, 1)
22    Y = as.matrix(c(0:4, 1:5))
23    results = lm(Y ~ X-1)
24    print(results)
25    print(summary(results))
26
27gives::
28
29    Call:
30    lm(formula = Y ~ X - 1)
31
32    Coefficients:
33    X1     X2
34    1.773  2.500
35
36    Residuals:
37        Min      1Q  Median      3Q     Max
38    -1.6970 -0.6667  0.0000  0.6667  1.6970
39
40    Coefficients:
41    Estimate Std. Error t value Pr(>|t|)
42    X1   1.7727     0.5455   3.250   0.0117 *
43    X2   2.5000     0.3482   7.181 9.42e-05 ***
44    ---
45
46    Residual standard error: 1.101 on 8 degrees of freedom
47    Multiple R-squared: 0.8859, Adjusted R-squared: 0.8574
48    F-statistic: 31.06 on 2 and 8 DF,  p-value: 0.0001694
49"""
50
51
52def test_model():
53    # Check basics about the model fit
54    # Check we fit the mean
55    assert_array_almost_equal(RESULTS.theta[1], np.mean(Y))
56    # Check we get the same as R
57    assert_array_almost_equal(RESULTS.theta, [1.773, 2.5], 3)
58    percentile = np.percentile
59    pcts = percentile(RESULTS.resid, [0, 25, 50, 75, 100])
60    assert_array_almost_equal(pcts, [-1.6970, -0.6667, 0, 0.6667, 1.6970], 4)
61
62
63def test_t_contrast():
64    # Test individual t against R
65    assert_array_almost_equal(RESULTS.t(0), 3.25)
66    assert_array_almost_equal(RESULTS.t(1), 7.181, 3)
67    # And contrast
68    assert_array_almost_equal(RESULTS.Tcontrast([1, 0]).t, 3.25)
69    assert_array_almost_equal(RESULTS.Tcontrast([0, 1]).t, 7.181, 3)
70    # Input matrix checked for size
71    with pytest.raises(ValueError):
72        RESULTS.Tcontrast([1])
73    with pytest.raises(ValueError):
74        RESULTS.Tcontrast([1, 0, 0])
75    # And shape
76    with pytest.raises(ValueError):
77        RESULTS.Tcontrast(np.array([1, 0])[:, None])
78
79
80def test_t_output():
81    # Check we get required outputs
82    exp_t = RESULTS.t(0)
83    exp_effect = RESULTS.theta[0]
84    exp_sd = exp_effect / exp_t
85    res = RESULTS.Tcontrast([1, 0])
86    assert_array_almost_equal(res.t, exp_t)
87    assert_array_almost_equal(res.effect, exp_effect)
88    assert_array_almost_equal(res.sd, exp_sd)
89    res = RESULTS.Tcontrast([1, 0], store=('effect',))
90    assert res.t is None
91    assert_array_almost_equal(res.effect, exp_effect)
92    assert res.sd is None
93    res = RESULTS.Tcontrast([1, 0], store=('t',))
94    assert_array_almost_equal(res.t, exp_t)
95    assert res.effect is None
96    assert res.sd is None
97    res = RESULTS.Tcontrast([1, 0], store=('sd',))
98    assert res.t is None
99    assert res.effect is None
100    assert_array_almost_equal(res.sd, exp_sd)
101    res = RESULTS.Tcontrast([1, 0], store=('effect', 'sd'))
102    assert res.t is None
103    assert_array_almost_equal(res.effect, exp_effect)
104    assert_array_almost_equal(res.sd, exp_sd)
105
106
107def test_f_output():
108    # Test f_output
109    res = RESULTS.Fcontrast([1, 0])
110    exp_f = RESULTS.t(0) ** 2
111    assert_array_almost_equal(exp_f, res.F)
112    # Test arrays work as well as lists
113    res = RESULTS.Fcontrast(np.array([1, 0]))
114    assert_array_almost_equal(exp_f, res.F)
115    # Test with matrix against R
116    res = RESULTS.Fcontrast(np.eye(2))
117    assert_array_almost_equal(31.06, res.F, 2)
118    # Input matrix checked for size
119    with pytest.raises(ValueError):
120        RESULTS.Fcontrast([1])
121    with pytest.raises(ValueError):
122        RESULTS.Fcontrast([1, 0, 0])
123    # And shape
124    with pytest.raises(ValueError):
125        RESULTS.Fcontrast(np.array([1, 0])[:, None])
126
127
128def test_f_output_new_api():
129    res = RESULTS.Fcontrast([1, 0])
130    assert_array_almost_equal(res.effect, RESULTS.theta[0])
131    assert_array_almost_equal(res.covariance, RESULTS.vcov()[0][0])
132
133
134def test_conf_int():
135    lower_, upper_ = RESULTS.conf_int()
136    assert (lower_ < upper_).all()
137    assert (lower_ > upper_ - 10).all()
138    lower_, upper_ = RESULTS.conf_int(cols=[1]).T
139    assert lower_ < upper_
140    assert lower_ > upper_ - 10
141