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_formulas (e.g., GEE, QIF, OrderedModel).

Code Sample, a copy-pastable example if possible

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()