MixedLM.hessian() returns a tuple instead of just the hessian matrix, causing MixedLM.fit to fail when minimizing with methods that require the hessian

[pcuestas]

Describe the bug

MixedLM.hessian() returns a tuple:

statsmodels/statsmodels/regression/mixed_linear_model.py

Lines 1855 to 1878 in 23faea3

	def hessian(self, params):
	"""
	Returns the model's Hessian matrix.
	Calculates the Hessian matrix for the linear mixed effects
	model with respect to the parameterization in which the
	covariance matrix is represented directly (without square-root
	transformation).
	Parameters
	----------
	params : MixedLMParams or array_like
	The model parameters at which the Hessian is calculated.
	If array-like, must contain the packed parameters in a
	form that is compatible with this model instance.
	Returns
	-------
	hess : 2d ndarray
	The Hessian matrix, evaluated at `params`.
	sing : boolean
	If True, the covariance matrix is singular and a
	pseudo-inverse is returned.
	"""

statsmodels/statsmodels/regression/mixed_linear_model.py

Line 2029 in 23faea3

return hess, sing

However, it is supposed to return just the hessian matrix, as stated in LikelyhoodModel.hessian():

statsmodels/statsmodels/base/model.py

Lines 332 to 345 in 23faea3

	def hessian(self, params):
	"""
	The Hessian matrix of the model.
	Parameters
	----------
	params : ndarray
	The parameters to use when evaluating the Hessian.
	Returns
	-------
	ndarray
	The hessian evaluated at the parameters.
	"""

The fact that MixedLM.hessian() returns a tuple causes optimization to fail when the hessian matrix is used in the minimization algorithm, when the following function is called (within MixedLM.fit()):

statsmodels/statsmodels/base/model.py

Lines 547 to 548 in 23faea3

	def hess(params, *args):
	return -self.hessian(params, *args) / nobs

(The exception TypeError: bad operand type for unary -: 'tuple' is thrown).

[kshedden]

As I recall that Hessian function is intended to be used for standard errors, not for optimization. The optimization uses a different parameterization that is not compatible with the Hessian function you are referring to. Optimization should use a first order or zeroth order approach that do not rely on Hessian matrices.