Repository: statsmodels/statsmodels (Python)

Issue:
two-sample ztest may have bug 
#### Describe the bug

The implementation of the two sample ztest does not appear to match documented algorithms

Specifically, the calculation of `var_pooled` below deviates slightly from how it is calculated elsewhere. 

https://github.com/statsmodels/statsmodels/blob/114d812133c0b8360c6d2e6059414e0706007977/statsmodels/stats/weightstats.py#L1543-L1561

For reference consult the following documentation:

https://en.wikipedia.org/wiki/Student%27s_t-test#Independent_two-sample_t-test

The code most closely resembles the "Equal or unequal sample sizes, similar variances" equation, but on line 1552 should subtract `ddof` from the nobs e.g. `var_pooled = (nobs1 - ddof) * x1_var + (nobs2 - ddof) * x2_var`

However we should not assume that the variances are similar, and should be using the more general "Equal or unequal sample sizes, unequal variances" equation.   


Hints:
Hi @gorj-tessella,
It seems there is a confusion between the z-test and the t-test. Your implementation aligns more with a t-test, considering the assumption of unequal variances. However, from my understanding, a z-test always assumes equal variances. Let me know if there is something I am missing.
unequal variance would also be useful in the z-test.
I don't remember why the option was not implemented.

However, the Satterthwaite degrees of freedom correction in the t-test will not apply and is irrelevant for the z-test.

I can implement the case of unequal variance.

Patch:
diff --git a/statsmodels/stats/weightstats.py b/statsmodels/stats/weightstats.py
--- a/statsmodels/stats/weightstats.py
+++ b/statsmodels/stats/weightstats.py
@@ -1510,10 +1510,10 @@ def ztest(
            'larger' :   H1: difference in means larger than value
            'smaller' :  H1: difference in means smaller than value
 
-    usevar : str, 'pooled'
-        Currently, only 'pooled' is implemented.
+    usevar : str, 'pooled' or 'unequal'
         If ``pooled``, then the standard deviation of the samples is assumed to be
-        the same. see CompareMeans.ztest_ind for different options.
+        the same. If ``unequal``, then the standard deviation of the sample is
+        assumed to be different.
     ddof : int
         Degrees of freedom use in the calculation of the variance of the mean
         estimate. In the case of comparing means this is one, however it can
@@ -1528,35 +1528,38 @@ def ztest(
 
     Notes
     -----
-    usevar not implemented, is always pooled in two sample case
-    use CompareMeans instead.
+    usevar can be pooled or unequal in two sample case
 
     """
     # TODO: this should delegate to CompareMeans like ttest_ind
     #       However that does not implement ddof
 
-    # usevar is not used, always pooled
+    # usevar can be pooled or unequal
 
-    if usevar != "pooled":
-        raise NotImplementedError('only usevar="pooled" is implemented')
+    if usevar not in {"pooled", "unequal"}:
+        raise NotImplementedError('usevar can only be "pooled" or "unequal"')
 
     x1 = np.asarray(x1)
     nobs1 = x1.shape[0]
     x1_mean = x1.mean(0)
     x1_var = x1.var(0)
+
     if x2 is not None:
         x2 = np.asarray(x2)
         nobs2 = x2.shape[0]
         x2_mean = x2.mean(0)
         x2_var = x2.var(0)
-        var_pooled = nobs1 * x1_var + nobs2 * x2_var
-        var_pooled /= nobs1 + nobs2 - 2 * ddof
-        var_pooled *= 1.0 / nobs1 + 1.0 / nobs2
+        if usevar == "pooled":
+            var = nobs1 * x1_var + nobs2 * x2_var
+            var /= nobs1 + nobs2 - 2 * ddof
+            var *= 1.0 / nobs1 + 1.0 / nobs2
+        elif usevar == "unequal":
+            var = x1_var / (nobs1 - ddof) + x2_var / (nobs2 - ddof)
     else:
-        var_pooled = x1_var / (nobs1 - ddof)
+        var = x1_var / (nobs1 - ddof)
         x2_mean = 0
 
-    std_diff = np.sqrt(var_pooled)
+    std_diff = np.sqrt(var)
     # stat = x1_mean - x2_mean - value
     return _zstat_generic(x1_mean, x2_mean, std_diff, alternative, diff=value)
 


Test patch:
diff --git a/statsmodels/stats/tests/test_weightstats.py b/statsmodels/stats/tests/test_weightstats.py
--- a/statsmodels/stats/tests/test_weightstats.py
+++ b/statsmodels/stats/tests/test_weightstats.py
@@ -677,6 +677,45 @@ def test_ztest_ztost():
 ztest_larger_mu_1s.method = 'One-sample z-Test'
 ztest_larger_mu_1s.data_name = 'x'
 
+# > zt = z.test(x, sigma.x=0.46436662631627995, y, sigma.y=0.7069805008424409)
+# > cat_items(zt, "ztest_unequal.")
+ztest_unequal = Holder()
+ztest_unequal.statistic = 6.12808151466544
+ztest_unequal.p_value = 8.89450168270109e-10
+ztest_unequal.conf_int = np.array([1.19415646579981, 2.31720717056382])
+ztest_unequal.estimate = np.array([7.01818181818182, 5.2625])
+ztest_unequal.null_value = 0
+ztest_unequal.alternative = 'two.sided'
+ztest_unequal.usevar = 'unequal'
+ztest_unequal.method = 'Two-sample z-Test'
+ztest_unequal.data_name = 'x and y'
+
+# > zt = z.test(x, sigma.x=0.46436662631627995, y, sigma.y=0.7069805008424409, alternative="less")
+# > cat_items(zt, "ztest_smaller_unequal.")
+ztest_smaller_unequal = Holder()
+ztest_smaller_unequal.statistic = 6.12808151466544
+ztest_smaller_unequal.p_value = 0.999999999555275
+ztest_smaller_unequal.conf_int = np.array([np.nan, 2.22692874913371])
+ztest_smaller_unequal.estimate = np.array([7.01818181818182, 5.2625])
+ztest_smaller_unequal.null_value = 0
+ztest_smaller_unequal.alternative = 'less'
+ztest_smaller_unequal.usevar = 'unequal'
+ztest_smaller_unequal.method = 'Two-sample z-Test'
+ztest_smaller_unequal.data_name = 'x and y'
+
+# > zt = z.test(x, sigma.x=0.46436662631627995, y, sigma.y=0.7069805008424409, alternative="greater")
+# > cat_items(zt, "ztest_larger_unequal.")
+ztest_larger_unequal = Holder()
+ztest_larger_unequal.statistic = 6.12808151466544
+ztest_larger_unequal.p_value = 4.44725034576265e-10
+ztest_larger_unequal.conf_int = np.array([1.28443488722992, np.nan])
+ztest_larger_unequal.estimate = np.array([7.01818181818182, 5.2625])
+ztest_larger_unequal.null_value = 0
+ztest_larger_unequal.alternative = 'greater'
+ztest_larger_unequal.usevar = 'unequal'
+ztest_larger_unequal.method = 'Two-sample z-Test'
+ztest_larger_unequal.data_name = 'x and y'
+
 
 alternatives = {'less': 'smaller',
                 'greater': 'larger',
@@ -733,6 +772,14 @@ def test(self):
                           alternative=alternatives[tc.alternative])
             assert_allclose(ci, tc_conf_int - tc.null_value, rtol=1e-10)
 
+        # unequal variances
+        for tc in [ztest_unequal, ztest_smaller_unequal, ztest_larger_unequal]:
+            zstat, pval = ztest(x1, x2, value=tc.null_value,
+                                alternative=alternatives[tc.alternative],
+                                usevar="unequal")
+            assert_allclose(zstat, tc.statistic, rtol=1e-10)
+            assert_allclose(pval, tc.p_value, rtol=1e-10, atol=1e-16)
+
         # 1 sample test copy-paste
         d1 = self.d1
         for tc in [ztest_mu_1s, ztest_smaller_mu_1s, ztest_larger_mu_1s]: