Repository: statsmodels/statsmodels (Python)

Issue:
Multi-dimensional indexing error in statsmodels.stats.diagnostic.linear_reset when using Pandas
#### Describe the bug

I tried to do a linear RESET test, but I got an error: `ValueError: Multi-dimensional indexing (e.g. `obj[:, None]`) is no longer supported. Convert to a numpy array before indexing instead.`

Complete error message down below.

#### Code Sample, a copy-pastable example if possible


```python
#!/usr/bin/env python3
import numpy as np
import pandas as pd
import statsmodels.stats.diagnostic as smdiag
import statsmodels.api as sm

n = 100

df = pd.DataFrame()

df["y"] = np.random.rand(n)
df["x"] = np.random.rand(n)

OLS_results = sm.OLS(df["y"], df["x"]).fit()

smdiag.linear_reset(OLS_results)
```
<details>

Complete error Traceback:
```
Traceback (most recent call last):
  File "/home/jonathan/RietNAS/Opleiding/EUR/1. Pre-master/FEB22009S Seminar Case Study/bug_report.py", line 16, in <module>
    smdiag.linear_reset(OLS_results)
  File "/home/jonathan/RietNAS/Opleiding/EUR/1. Pre-master/FEB22009S Seminar Case Study/.venv/lib/python3.10/site-packages/pandas/util/_decorators.py", line 210, in wrapper
    return func(*args, **kwargs)
  File "/home/jonathan/RietNAS/Opleiding/EUR/1. Pre-master/FEB22009S Seminar Case Study/.venv/lib/python3.10/site-packages/statsmodels/stats/diagnostic.py", line 1065, in linear_reset
    aug = res.fittedvalues[:, None]
  File "/home/jonathan/RietNAS/Opleiding/EUR/1. Pre-master/FEB22009S Seminar Case Study/.venv/lib/python3.10/site-packages/pandas/core/series.py", line 1033, in __getitem__
    return self._get_with(key)
  File "/home/jonathan/RietNAS/Opleiding/EUR/1. Pre-master/FEB22009S Seminar Case Study/.venv/lib/python3.10/site-packages/pandas/core/series.py", line 1048, in _get_with
    return self._get_values_tuple(key)
  File "/home/jonathan/RietNAS/Opleiding/EUR/1. Pre-master/FEB22009S Seminar Case Study/.venv/lib/python3.10/site-packages/pandas/core/series.py", line 1082, in _get_values_tuple
    disallow_ndim_indexing(result)
  File "/home/jonathan/RietNAS/Opleiding/EUR/1. Pre-master/FEB22009S Seminar Case Study/.venv/lib/python3.10/site-packages/pandas/core/indexers/utils.py", line 343, in disallow_ndim_indexing
    raise ValueError(
ValueError: Multi-dimensional indexing (e.g. `obj[:, None]`) is no longer supported. Convert to a numpy array before indexing instead.
```

The error does not occur if `y` and `x` are passed as NumPy arrays, only when passed as pandas dataframe.

Using `statsmodels.stats.outliers_influence.reset_ramsey` also doesn't give this same problem.


</details>


If the issue has not been resolved, please file it in the issue tracker.

#### Expected Output

I expect to get the results of the test, as `statsmodels.stats.outliers_influence import reset_ramsey` does.

#### Output of ``import statsmodels.api as sm; sm.show_versions()``

<details>

[paste the output of ``import statsmodels.api as sm; sm.show_versions()`` here below this line]
INSTALLED VERSIONS
------------------
Python: 3.10.6.final.0
OS: Linux 6.2.6-76060206-generic #202303130630~1681329778~22.04~d824cd4 SMP PREEMPT_DYNAMIC Wed A x86_64
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8

statsmodels
===========

Installed: 0.14.0 (/home/jonathan/RietNAS/Opleiding/EUR/1. Pre-master/FEB22009S Seminar Case Study/.venv/lib/python3.10/site-packages/statsmodels)

Required Dependencies
=====================

cython: Not installed
numpy: 1.24.3 (/home/jonathan/RietNAS/Opleiding/EUR/1. Pre-master/FEB22009S Seminar Case Study/.venv/lib/python3.10/site-packages/numpy)
scipy: 1.10.1 (/home/jonathan/RietNAS/Opleiding/EUR/1. Pre-master/FEB22009S Seminar Case Study/.venv/lib/python3.10/site-packages/scipy)
pandas: 2.0.1 (/home/jonathan/RietNAS/Opleiding/EUR/1. Pre-master/FEB22009S Seminar Case Study/.venv/lib/python3.10/site-packages/pandas)
    dateutil: 2.8.2 (/home/jonathan/RietNAS/Opleiding/EUR/1. Pre-master/FEB22009S Seminar Case Study/.venv/lib/python3.10/site-packages/dateutil)
patsy: 0.5.3 (/home/jonathan/RietNAS/Opleiding/EUR/1. Pre-master/FEB22009S Seminar Case Study/.venv/lib/python3.10/site-packages/patsy)

Optional Dependencies
=====================

matplotlib: 3.7.1 (/home/jonathan/RietNAS/Opleiding/EUR/1. Pre-master/FEB22009S Seminar Case Study/.venv/lib/python3.10/site-packages/matplotlib)
    backend: agg 
cvxopt: Not installed
joblib: 1.2.0 (/home/jonathan/RietNAS/Opleiding/EUR/1. Pre-master/FEB22009S Seminar Case Study/.venv/lib/python3.10/site-packages/joblib)

Developer Tools
================

IPython: 8.13.2 (/home/jonathan/RietNAS/Opleiding/EUR/1. Pre-master/FEB22009S Seminar Case Study/.venv/lib/python3.10/site-packages/IPython)
    jinja2: Not installed
sphinx: Not installed
    pygments: 2.15.1 (/home/jonathan/RietNAS/Opleiding/EUR/1. Pre-master/FEB22009S Seminar Case Study/.venv/lib/python3.10/site-packages/pygments)
pytest: Not installed
virtualenv: Not installed

</details>



Hints:
looks like a bug with pandas compatibility
The function should use the unwrapped results instance internally.


`smdiag.linear_reset(OLS_results._results)` should have worked as workaround, however it breaks on the isinstance check.

So the only workaround I see right now (except for monkey patching fittedvalues) is to reestimate the model without pandas:

```
res_np = sm.OLS(OLS_results.model.endog, OLS_results.model.exog).fit()
smdiag.linear_reset(res_np)
<class 'statsmodels.stats.contrast.ContrastResults'>
<Wald test (chi2): statistic=63.544945204592665, p-value=1.5899756467736924e-14, df_denom=2>
```


Yeah that sounds like a fair workaround.

I'm not sure what the internals of this looks like, but a call `y = np.array(y)` somewhere in `linear_reset` probably fixes it, but I'll leave it to the devs to devise a better solution (:
I have the same issue for `.compare_lm_test()`.

```
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

n = 100

df = pd.DataFrame()

for v in ["y","x","z","k"]:
    df[v] = np.random.rand(n)

res0 = smf.ols("y ~ x + z + k", data=df).fit()
res1 = smf.ols("y ~ x", data=df).fit()

res0.compare_lm_test(res1)
```

The above code generates the following error

```
ValueError: Multi-dimensional indexing (e.g. `obj[:, None]`) is no longer supported. Convert to a numpy array before indexing instead.
```

Python: 3.10.11
pandas: 2.0.1
statsmodels   0.14.0
You can downgrade your pandas as well to get around this.  

Should be fixed in next point release.


Patch:
diff --git a/statsmodels/stats/_adnorm.py b/statsmodels/stats/_adnorm.py
--- a/statsmodels/stats/_adnorm.py
+++ b/statsmodels/stats/_adnorm.py
@@ -5,6 +5,8 @@
 Author: Josef Perktold and Scipy developers
 License : BSD-3
 """
+import warnings
+
 import numpy as np
 from scipy import stats
 
@@ -64,8 +66,12 @@ def anderson_statistic(x, dist='norm', fit=True, params=(), axis=0):
     sl2 = [slice(None)] * x.ndim
     sl2[axis] = slice(None, None, -1)
     sl2 = tuple(sl2)
-    s = np.sum((2 * i[sl1] - 1.0) / nobs * (np.log(z) + np.log1p(-z[sl2])),
-               axis=axis)
+    with warnings.catch_warnings():
+        warnings.filterwarnings(
+            "ignore", message="divide by zero encountered in log1p"
+        )
+        ad_values = (2 * i[sl1] - 1.0) / nobs * (np.log(z) + np.log1p(-z[sl2]))
+        s = np.sum(ad_values, axis=axis)
     a2 = -nobs - s
     return a2
 
diff --git a/statsmodels/stats/diagnostic.py b/statsmodels/stats/diagnostic.py
--- a/statsmodels/stats/diagnostic.py
+++ b/statsmodels/stats/diagnostic.py
@@ -1072,7 +1072,7 @@ def linear_reset(res, power=3, test_type="fitted", use_f=False,
             raise ValueError("power must contains distinct integers all >= 2")
     exog = res.model.exog
     if test_type == "fitted":
-        aug = res.fittedvalues[:, None]
+        aug = np.asarray(res.fittedvalues)[:, None]
     elif test_type == "exog":
         # Remove constant and binary
         aug = res.model.exog
@@ -1293,7 +1293,7 @@ def linear_lm(resid, exog, func=None):
     if func is None:
         def func(x):
             return np.power(x, 2)
-
+    exog = np.asarray(exog)
     exog_aux = np.column_stack((exog, func(exog[:, 1:])))
 
     nobs, k_vars = exog.shape


Test patch:
diff --git a/statsmodels/stats/tests/test_diagnostic.py b/statsmodels/stats/tests/test_diagnostic.py
--- a/statsmodels/stats/tests/test_diagnostic.py
+++ b/statsmodels/stats/tests/test_diagnostic.py
@@ -1957,3 +1957,43 @@ def test_small_skip(reset_randomstate):
 # M2 + fit(M1)-exp(fit(M2)) < 2.22e-16 ***
 # ---
 # Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+
+@pytest.mark.smoke
+def test_diagnostics_pandas(reset_randomstate):
+    # GH 8879
+    n = 100
+    df = pd.DataFrame(
+        {
+            "y": np.random.rand(n),
+            "x": np.random.rand(n),
+            "z": np.random.rand(n)}
+    )
+    y, x = df["y"], add_constant(df["x"])
+
+    res = OLS(df["y"], add_constant(df[["x"]])).fit()
+    res_large = OLS(df["y"], add_constant(df[["x", "z"]])).fit()
+    res_other = OLS(df["y"], add_constant(df[["z"]])).fit()
+    smsdia.linear_reset(res_large)
+    smsdia.linear_reset(res_large, test_type="fitted")
+    smsdia.linear_reset(res_large, test_type="exog")
+    smsdia.linear_reset(res_large, test_type="princomp")
+    smsdia.het_goldfeldquandt(y, x)
+    smsdia.het_breuschpagan(res.resid, x)
+    smsdia.het_white(res.resid, x)
+    smsdia.het_arch(res.resid)
+    smsdia.acorr_breusch_godfrey(res)
+    smsdia.acorr_ljungbox(y)
+    smsdia.linear_rainbow(res)
+    smsdia.linear_lm(res.resid, x)
+    smsdia.linear_harvey_collier(res)
+    smsdia.acorr_lm(res.resid)
+    smsdia.breaks_cusumolsresid(res.resid)
+    smsdia.breaks_hansen(res)
+    smsdia.compare_cox(res, res_other)
+    smsdia.compare_encompassing(res, res_other)
+    smsdia.compare_j(res, res_other)
+    smsdia.recursive_olsresiduals(res)
+    smsdia.recursive_olsresiduals(
+        res, order_by=np.arange(y.shape[0] - 1, 0 - 1, -1)
+    )
+    smsdia.spec_white(res.resid, x)