Repository: scipy/scipy
Pull request: ENH: allow single observation for equal variance in `stats.ttest_ind`
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->

Closes #18069

#### What does this implement/fix?

Single observation in `ttest_ind` are currently not possible due to a cancellation. This allows it by setting a null variance when only a single observation is available.


Adding warning filtering is causing some issues. Hence skipping the CI for now:

```python
FAILED scipy/stats/tests/test_stats.py::TestMoments::test_skewness - Failed: DID NOT WARN. No warnings of type (<class 'RuntimeWarning'>,) were emitted. The list of emitted warnings is: [].
FAILED scipy/stats/tests/test_stats.py::TestMoments::test_kurtosis - Failed: DID NOT WARN. No warnings of type (<class 'RuntimeWarning'>,) were emitted. The list of emitted warnings is: [].
FAILED scipy/stats/tests/test_stats.py::TestStudentTest::test_onesample - Failed: DID NOT WARN. No warnings of type (<class 'RuntimeWarning'>,) were emitted. The list of emitted warnings is: [].
FAILED scipy/stats/tests/test_stats.py::test_ttest_rel - Failed: DID NOT WARN. No warnings of type (<class 'RuntimeWarning'>,) were emitted. The list of emitted warnings is: [].
FAILED scipy/stats/tests/test_stats.py::test_ttest_ind - ZeroDivisionError: float division by zero
FAILED scipy/stats/tests/test_stats.py::test_gh5686 - ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
FAILED scipy/stats/tests/test_stats.py::TestDescribe::test_describe_scalar - Failed: DID NOT WARN. No warnings of type (<class 'RuntimeWarning'>,) were emitted. The list of emitted warnings is: [].
FAILED scipy/stats/tests/test_stats.py::test_normalitytests - Failed: DID NOT WARN. No warnings of type (<class 'RuntimeWarning'>,) were emitted. The list of emitted warnings is: [].
```



Problem statement:
ENH: stats.ttest_ind is inconsistent with R. It does not allow testing against a single datapoint
### Is your feature request related to a problem? Please describe.

It is frequently useful to perform a two-sample t-test for the case when one of the samples has only one datapoint, and variances are equal. Mathematically, only the total variance is required for the test, so variance of 0 can be assumed for the degenerate dataset.

### Describe the solution you'd like.

I would like for the command

ttest_ind(np.random.normal(0,1,100), [1], equal_var=True)

to produce a test statistic and a p-value, not (NAN, NAN) as it does at the moment. Can be made consistent with R, where one can do the following without problem

t.test(rnorm(100),1,var.equal = TRUE)

### Describe alternatives you've considered.

_No response_

### Additional context (e.g. screenshots, GIFs)

_No response_

Hints from issue discussion:
Often the formula for the t statistic is expressed in terms of the individual sample variances. Can the formula be reduced to show that the NaN you get is just an artifact of the way the formula is expressed? Or is it an assumption, as you wrote, that the sample variance is zero?

If it is an assumption, does the test still produce uniformly distributed p-values under the null hypothesis, or is there a modification to the null distribution?

What I'm after is whether this a correct thing to do, an approximation, or just something that gives non-NaN results.

(At first glance it looks like it's correct that there is an `(n-1)/(n-1)` cancellation - not an assumption - but I can't really look until later.)
I am sorry. The solution was recommended to me by stats SE community, to them it makes common sense. I am trying to find the proof for this case online, but so far I have not succeeded. Here is the discussion on SE, feel free to join in

https://stats.stackexchange.com/questions/606893

Looking at the implementation of the ttest_ind in scipy, I can certainly tell that it attempts to calculate the variance of both arrays in all cases, which should result in a NAN if the array only has one datapoint. Thus, one could extend the functionality of the ttest_ind by allowing exactly one of the arrays to have length 1. This can be done by putting a bunch of IF statements: throwing an exception of both have length 1, otherwise setting the corresponding variance to 0 manually if only one of the arrays have length 1.

So the edit to the function is simple. It still remains to be seen if this is actually correct. I am hoping if somebody with more statistics knowledge could comment
No problem; I'll investigate later today. It looks like it's correct, not an assumption, at least for equal variances.
Yup, I can't find anything wrong with this. The formulas I'm finding present the pooled variance as a weighted average of the sample variances, but symbolically, the `(n-1)` weights cancel with the `(n-1)`s in the denominators of the unbiased sample variances, and we're left with the sum of squared differences of the sample means. This is zero for the sample with only one observation.

After making this fix, as a sanity check:
```python3
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
rng = np.random.default_rng(246834602926842)

x = rng.normal(size=(10000, 2))
y = rng.normal(size=(10000, 1))
res = stats.ttest_ind(x, y, equal_var=True, axis=-1)
plt.hist(res.pvalue)
plt.xlabel('p-value')
plt.ylabel('frequency')
plt.title('p-values produced by ttest_ind under the null hypothesis')
```
![image](https://user-images.githubusercontent.com/6570539/222023862-ecf9b4a1-496d-4fcf-9a62-5076b87025f6.png)

Looks quite uniform to me.  Also `stats.ttest_ind([2, 3, 5], [1.5], equal_var=True)` gives the same result as `t.test(c(2, 3, 5), c(1.5), equal.var=TRUE)` in R.

There are a few ways to get the same numerical result, including setting the sample variance to zero (although this _looks_ a bit suspicious, so it would deserve an explanatory comment). **PR welcome**, but please keep in mind that we only want to change the `equal_var=True` case. I haven't looked at the unequal variance math lately, but I doubt it works the same way, and R produces a `not enough 'y' observations` error in this case. 
@mdhaber I suppose that would help support the issue we've seen on Dunnett's test with the test case having a single obs in the last group.
Sure. If you can open a PR, we can probably merge it today.
👍 for visibility also putting this here. I propose to add the following in `_equal_var_ttest_denom`.

```python
    if n1 == 1:
        v1 = 0
    if n2 == 1:
        v2 = 0
```

Then the only thing is that there is the `RuntimeWarning`. Shall I just filter it out in the root function?

Yes, that and a comment would be good. I would also be OK with one-liners in this case, e.g. `v1 = 0 if n1 == 1 else v1`.

Preferably, the "Precision loss occurred in moment calculation due to " "catastrophic cancellation." would not be emitted by `_moment`  at all if there is only one observation.

The others due to division by zero in `var *= np.divide(n, n-ddof)` can safely be silenced with `np.errstate`, I think.

Commit message: TST: add single observation test

[skip ci]

Patch to review:
diff --git a/scipy/stats/_stats_py.py b/scipy/stats/_stats_py.py
index c78e7b71c18a..6e3cb5d4b998 100644
--- a/scipy/stats/_stats_py.py
+++ b/scipy/stats/_stats_py.py
@@ -1131,7 +1131,8 @@ def _moment(a, moment, axis, *, mean=None):
                               keepdims=True) / np.abs(mean)
         with np.errstate(invalid='ignore'):
             precision_loss = np.any(rel_diff < eps)
-        if precision_loss:
+        n = a.shape[axis] if axis is not None else a.size
+        if precision_loss and n > 1:
             message = ("Precision loss occurred in moment calculation due to "
                        "catastrophic cancellation. This occurs when the data "
                        "are nearly identical. Results may be unreliable.")
@@ -1155,7 +1156,8 @@ def _var(x, axis=0, ddof=0, mean=None):
     var = _moment(x, 2, axis, mean=mean)
     if ddof != 0:
         n = x.shape[axis] if axis is not None else x.size
-        var *= np.divide(n, n-ddof)  # to avoid error on division by zero
+        with np.errstate(divide='ignore', invalid='ignore'):
+            var *= np.divide(n, n-ddof)  # to avoid error on division by zero
     return var
 
 
@@ -6963,6 +6965,10 @@ def _unequal_var_ttest_denom(v1, n1, v2, n2):
 
 
 def _equal_var_ttest_denom(v1, n1, v2, n2):
+    # null variance in case there is a single observation
+    v1 = 0 if n1 == 1 else v1
+    v2 = 0 if n2 == 1 else v2
+
     df = n1 + n2 - 2.0
     svar = ((n1 - 1) * v1 + (n2 - 1) * v2) / df
     denom = np.sqrt(svar * (1.0 / n1 + 1.0 / n2))
diff --git a/scipy/stats/tests/test_stats.py b/scipy/stats/tests/test_stats.py
index 9ff9664c5542..5355e200ab4b 100644
--- a/scipy/stats/tests/test_stats.py
+++ b/scipy/stats/tests/test_stats.py
@@ -5221,6 +5221,23 @@ def test_ttest_ind_from_stats_inputs_zero():
     assert_equal(result, [np.nan, np.nan])
 
 
+def test_ttest_single_observation():
+    rng = np.random.default_rng(246834602926842)
+    x = rng.normal(size=(10000, 2))
+    y = rng.normal(size=(10000, 1))
+    res = stats.ttest_ind(x, y, equal_var=True, axis=-1)
+
+    assert_allclose(np.percentile(res.pvalue, 25), 0.25, atol=1e-2)
+    assert_allclose(np.percentile(res.pvalue, 50), 0.5, atol=1e-2)
+    assert_allclose(np.percentile(res.pvalue, 75), 0.75, atol=1e-2)
+
+    res = stats.ttest_ind(y, x, equal_var=True, axis=-1)
+
+    assert_allclose(np.percentile(res.pvalue, 25), 0.25, atol=1e-2)
+    assert_allclose(np.percentile(res.pvalue, 50), 0.5, atol=1e-2)
+    assert_allclose(np.percentile(res.pvalue, 75), 0.75, atol=1e-2)
+
+
 def test_ttest_1samp_new():
     n1, n2, n3 = (10,15,20)
     rvn1 = stats.norm.rvs(loc=5,scale=10,size=(n1,n2,n3))


Code review comments:
diff_hunk: @@ -1155,7 +1156,8 @@ def _var(x, axis=0, ddof=0, mean=None):
     var = _moment(x, 2, axis, mean=mean)
     if ddof != 0:
         n = x.shape[axis] if axis is not None else x.size
-        var *= np.divide(n, n-ddof)  # to avoid error on division by zero
+        with np.errstate(divide='ignore', invalid='ignore'):
path: scipy/stats/_stats_py.py
text: @user1:
My fault. For now, go in `ttest_ind` and preferably only when `equal_var` is True. We shouldn't always silence the warning right now because usually this division by zero results in NaNs that the user will see, but in the case of `ttest_ind` with `equal_var=True` it doesn't.

The rest looks good! Although some of the errors will probably persist after this change. Can't investigate those right now but maybe sometimes `n` is an array in `_equal_var_ttest_denom`, like with `nan_policy='omit'`?

@author:
I moved the silencing to the main function as you suggested.

diff_hunk: @@ -6963,6 +6965,10 @@ def _unequal_var_ttest_denom(v1, n1, v2, n2):
 
 
 def _equal_var_ttest_denom(v1, n1, v2, n2):
+    # null variance in case there is a single observation
+    v1 = 0 if n1 == 1 else v1
path: scipy/stats/_stats_py.py
text: @user1:
Not sure but maybe this being 0 is instead of NaN is causing the divide by 0 error. That can be fixed by using `np.divide` instead or by making an argument of the division a NumPy 0.

@author:
I added `np.divide` and it did solve an issue.