Repository: statsmodels/statsmodels
Pull request: Fix formula eval depth in select models
- [X] closes #9047
- [X] closes #9037
- [X] tests added / passed. 
- [X] code/documentation is well formatted.  
- [X] properly formatted commit message. See 
      [NumPy's guide](https://docs.scipy.org/doc/numpy-1.15.1/dev/gitwash/development_workflow.html#writing-the-commit-message). 

<details>


**Notes**:

* It is essential that you add a test when making code changes. Tests are not 
  needed for doc changes.
* When adding a new function, test values should usually be verified in another package (e.g., R/SAS/Stata).
* When fixing a bug, you must add a test that would produce the bug in main and
  then show that it is fixed with the new code.
* New code additions must be well formatted. Changes should pass flake8. If on Linux or OSX, you can
  verify you changes are well formatted by running 
  ```
  git diff upstream/main -u -- "*.py" | flake8 --diff --isolated
  ```
  assuming `flake8` is installed. This command is also available on Windows 
  using the Windows System for Linux once `flake8` is installed in the 
  local Linux environment. While passing this test is not required, it is good practice and it help 
  improve code quality in `statsmodels`.
* Docstring additions must render correctly, including escapes and LaTeX.

</details>


Problem statement:
BUG: Cox Model from formula doesn't grab the correct environment
#### Describe the bug

Constructing a proportional hazard model from a formula with `smf.phreg` (aka `PHReg.from_formula`) doesn't use the proper environment.

Specifically, `PHReg` overrides the `from_formula` class method but doesn't adjust or set the `eval_env` keyword. It calls `Model.from_formula` to parse the formula, but `Model.from_formula` doesn't know that it was called indirectly and therefore underestimates its depth by one. In the default case, `Model.from_formula` pulls the environment from where it was called, but this is `hazard_regression.py` rather than the user's environment (where `PHReg.from_formula` was called). User-defined symbols aren't found, though internal ones like `np` are available.

This problem should also show up for all the other subclass implementations of `from_formula`s (e.g., `GEE`, `QIF`, `OrderedModel`).

#### Code Sample, a copy-pastable example if possible


```python
import pandas as pd
import statsmodels.formula.api as smf

data = pd.DataFrame(
	[
		{"y": 2, "x": 2},
		{"y": 2, "x": 4},
		{"y": 4, "x": 6},
		{"y": 4, "x": 8},
	]
)

# Works
smf.phreg("y ~ x", data=data).fit()

def plus_one(x):
	return x + 1

# Works
smf.ols("y ~ plus_one(x)", data=data).fit()

# PatsyError from NameError: name 'plus_one' is not defined
smf.phreg("y ~ plus_one(x)", data=data).fit()
```

#### Expected Output

A clear and concise description of what you expected to happen.

#### Output of ``import statsmodels.api as sm; sm.show_versions()``

<details>

INSTALLED VERSIONS
------------------
Python: 3.11.5.final.0
OS: Darwin
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8

statsmodels
===========

Installed: 0.14.0

Required Dependencies
=====================

cython: Not installed
numpy: 1.26.1
scipy: 1.11.3
pandas: 2.1.1
    dateutil: 2.8.2
patsy: 0.5.3

Optional Dependencies
=====================

matplotlib: Not installed
cvxopt: Not installed
joblib: Not installed

Developer Tools
================

IPython: Not installed
    jinja2: Not installed
sphinx: Not installed
    pygments: Not installed
pytest: Not installed
virtualenv: Not installed

</details>

Hints from issue discussion:

Commit message: TST: Add tests for eval depth

Patch to review:
diff --git a/statsmodels/genmod/tests/test_gee.py b/statsmodels/genmod/tests/test_gee.py
index 1965335879c..8ddc9151590 100644
--- a/statsmodels/genmod/tests/test_gee.py
+++ b/statsmodels/genmod/tests/test_gee.py
@@ -1245,6 +1245,46 @@ def test_formulas(self):
 
         check_wrapper(rslt2)
 
+    def test_formula_environment(self):
+        """Test that GEE uses the right environment for formulas."""
+
+        n = 100
+        rng = np.random.default_rng(34234)
+        X1 = rng.normal(size=n)
+        Y = X1 + rng.normal(size=n)
+        Time = rng.uniform(size=n)
+        groups = np.kron(lrange(20), np.ones(5))
+
+        data = pd.DataFrame({"Y": Y, "X1": X1, "Time": Time, "groups": groups})
+
+        va = cov_struct.Autoregressive(grid=False)
+        family = families.Gaussian()
+
+        def times_two(x):
+            return 2 * x
+
+        mat = np.concatenate((np.ones((n, 1)), times_two(X1[:, None])), axis=1)
+        result_direct = gee.GEE(
+            Y, mat, groups, time=Time, family=family, cov_struct=va
+        ).fit()
+        assert result_direct is not None
+
+        result_formula = gee.GEE.from_formula(
+            "Y ~ times_two(X1)",
+            groups,
+            data,
+            time=Time,
+            family=family,
+            cov_struct=va,
+        ).fit()
+        assert result_formula is not None
+
+        assert_almost_equal(
+            result_direct.params,
+            result_formula.params,
+            decimal=8,
+        )
+
     def test_compare_logit(self):
 
         vs = cov_struct.Independence()
diff --git a/statsmodels/genmod/tests/test_qif.py b/statsmodels/genmod/tests/test_qif.py
index 6e247b6c459..c0b663f08d6 100644
--- a/statsmodels/genmod/tests/test_qif.py
+++ b/statsmodels/genmod/tests/test_qif.py
@@ -111,3 +111,37 @@ def test_formula(cov_struct):
     if not isinstance(cov_struct, QIFIndependence):
         _ = result2.bic
         _ = result2.aic
+
+
+def test_formula_environment():
+    """Test that QIF uses the right environment for formulas."""
+
+    rng = np.random.default_rng(3423)
+
+    x1 = rng.normal(size=100)
+    y = x1 + rng.normal(size=100)
+    groups = np.kron(np.arange(25), np.ones(4))
+
+    def times_two(x):
+        return 2 * x
+
+    cov_struct = QIFIndependence()
+
+    result_direct = QIF(
+        y,
+        times_two(x1).reshape(-1, 1),
+        groups=groups,
+        cov_struct=cov_struct
+    ).fit()
+
+    df = pd.DataFrame({"y": y, "x1": x1, "groups": groups})
+
+    result_formula = QIF.from_formula(
+        "y ~ 0 + times_two(x1)",
+        groups="groups",
+        cov_struct=cov_struct,
+        data=df
+    ).fit()
+
+    assert_allclose(result_direct.params, result_formula.params)
+    assert_allclose(result_direct.bse, result_formula.bse)


Code review comments:
diff_hunk: @@ -74,6 +74,12 @@
         model = "I(food/income) ~ income + persons"
         cls.income_fit = BetaModel.from_formula(model, income).fit()
 
+        def times_two(x):
path: statsmodels/othermod/tests/test_beta.py
text: @user1:
## Unused local variable

Variable times_two is not used.

[Show more details](https://github.com/statsmodels/statsmodels/security/code-scanning/2802)

diff_hunk: @@ -305,6 +305,25 @@
                     distr="probit",
                 )
 
+    def test_formula_eval_env(self):
+        def times_two(x):
path: statsmodels/miscmodels/tests/test_ordinal_model.py
text: @user1:
## Unused local variable

Variable times_two is not used.

[Show more details](https://github.com/statsmodels/statsmodels/security/code-scanning/2804)

diff_hunk: @@ -325,6 +325,21 @@
             data=df)
         result1 = model1.fit()
 
+        def times_two(x):
path: statsmodels/regression/tests/test_lme.py
text: @user1:
## Unused local variable

Variable times_two is not used.

[Show more details](https://github.com/statsmodels/statsmodels/security/code-scanning/2803)