from __future__ import annotations

from statsmodels.compat.python import lzip

from collections import defaultdict
from functools import reduce
import warnings

import numpy as np
import pandas as pd
from scipy import stats

from statsmodels.base.data import handle_data
from statsmodels.base.optimizer import Optimizer
import statsmodels.base.wrapper as wrap
from statsmodels.formula import handle_formula_data
from statsmodels.formula._manager import FormulaManager
from statsmodels.stats.contrast import (
    ContrastResults,
    WaldTestResults,
    t_test_pairwise,
)
from statsmodels.tools.data import _is_using_pandas
from statsmodels.tools.decorators import (
    cache_readonly,
    cached_data,
    cached_value,
)
from statsmodels.tools.numdiff import approx_fprime
from statsmodels.tools.sm_exceptions import (
    HessianInversionWarning,
    ValueWarning,
)
from statsmodels.tools.tools import nan_dot, recipr
from statsmodels.tools.validation import bool_like

ERROR_INIT_KWARGS = False

_model_params_doc = """Parameters
    ----------
    endog : array_like
        A 1-d endogenous response variable. The dependent variable.
    exog : array_like
        A nobs x k array where `nobs` is the number of observations and `k`
        is the number of regressors. An intercept is not included by default
        and should be added by the user. See
        :func:`statsmodels.tools.add_constant`."""

_missing_param_doc = """\
missing : str
        Available options are 'none', 'drop', and 'raise'. If 'none', no nan
        checking is done. If 'drop', any observations with nans are dropped.
        If 'raise', an error is raised. Default is 'none'."""
_extra_param_doc = """
    hasconst : None or bool
        Indicates whether the RHS includes a user-supplied constant. If True,
        a constant is not checked for and k_constant is set to 1 and all
        result statistics are calculated as if a constant is present. If
        False, a constant is not checked for and k_constant is set to 0.
    **kwargs
        Extra arguments that are used to set model properties when using the
        formula interface."""


class Model:
    __doc__ = """
    A (predictive) statistical model. Intended to be subclassed not used.

    {params_doc}
    {extra_params_doc}

    Attributes
    ----------
    exog_names
    endog_names

    Notes
    -----
    `endog` and `exog` are references to any data provided.  So if the data is
    already stored in numpy arrays and it is changed then `endog` and `exog`
    will change as well.
    """.format(
        params_doc=_model_params_doc,
        extra_params_doc=_missing_param_doc + _extra_param_doc,
    )

    # Maximum number of endogenous variables when using a formula
    # Default is 1, which is more common. Override in models when needed
    # Set to None to skip check
    _formula_max_endog = 1
    # kwargs that are generically allowed, maybe not supported in all models
    _kwargs_allowed = [
        "missing",
        "missing_idx",
        "formula",
        "model_spec",
        "hasconst",
    ]

    def __init__(self, endog, exog=None, **kwargs):
        missing = kwargs.pop("missing", "none")
        hasconst = kwargs.pop("hasconst", None)
        self.data = self._handle_data(endog, exog, missing, hasconst, **kwargs)
        self.k_constant = self.data.k_constant
        self.exog = self.data.exog
        self.endog = self.data.endog
        self._data_attr = []
        self._data_attr.extend(["exog", "endog", "data.exog", "data.endog"])
        if "formula" not in kwargs:  # will not be able to unpickle without these
            self._data_attr.extend(["data.orig_endog", "data.orig_exog"])
        # store keys for extras if we need to recreate model instance
        # we do not need 'missing', maybe we need 'hasconst'
        self._init_keys = list(kwargs.keys())
        if hasconst is not None:
            self._init_keys.append("hasconst")

    def _get_init_kwds(self):
        """return dictionary with extra keys used in model.__init__"""
        kwds = {key: getattr(self, key, None) for key in self._init_keys}

        return kwds

    def _check_kwargs(self, kwargs, keys_extra=None, error=ERROR_INIT_KWARGS):

        kwargs_allowed = [
            "missing",
            "missing_idx",
            "formula",
            "model_spec",
            "hasconst",
        ]
        if keys_extra:
            kwargs_allowed.extend(keys_extra)

        kwargs_invalid = [i for i in kwargs if i not in kwargs_allowed]
        if kwargs_invalid:
            msg = "unknown kwargs " + repr(kwargs_invalid)
            if error is False:
                warnings.warn(msg, ValueWarning, stacklevel=2)
            else:
                raise ValueError(msg)

    def _handle_data(self, endog, exog, missing, hasconst, **kwargs):
        data = handle_data(endog, exog, missing, hasconst, **kwargs)
        # kwargs arrays could have changed, easier to just attach here
        for key in kwargs:
            if key in ["model_spec", "formula"]:  # leave attached to data
                continue
            # pop so we do not start keeping all these twice or references
            try:
                setattr(self, key, data.__dict__.pop(key))
            except KeyError:  # panel already pops keys in data handling
                pass
        return data

    @classmethod
    def from_formula(cls, formula, data, subset=None, drop_cols=None, *args, **kwargs):
        """
        Create a Model from a formula and dataframe.

        Parameters
        ----------
        formula : str or generic Formula object
            The formula specifying the model.
        data : array_like
            The data for the model. See Notes.
        subset : array_like
            An array-like object of booleans, integers, or index values that
            indicate the subset of df to use in the model. Assumes df is a
            `pandas.DataFrame`.
        drop_cols : array_like
            Columns to drop from the design matrix.  Cannot be used to
            drop terms involving categoricals.
        *args
            Additional positional argument that are passed to the model.
        **kwargs
            These are passed to the model with one exception. The
            ``eval_env`` keyword is passed to patsy. It can be either a
            :class:`patsy:patsy.EvalEnvironment` object or an integer
            indicating the depth of the namespace to use. For example, the
            default ``eval_env=0`` uses the calling namespace. If you wish
            to use a "clean" environment set ``eval_env=-1``.

        Returns
        -------
        model
            The model instance.

        Notes
        -----
        data must define __getitem__ with the keys in the formula terms
        args and kwargs are passed on to the model instantiation. E.g.,
        a numpy structured or rec array, a dictionary, or a pandas DataFrame.
        """
        # TODO: provide a docs template for args/kwargs from child models
        # TODO: subset could use syntax. issue #469.
        mgr = FormulaManager()
        if subset is not None:
            data = data.loc[subset]
        eval_env = kwargs.pop("eval_env", None)
        if eval_env is None:
            eval_env = 2
        elif eval_env == -1:
            eval_env = mgr.get_empty_eval_env()
        elif isinstance(eval_env, int):
            eval_env += 1  # we're going down the stack again
        missing = kwargs.get("missing", "drop")
        if missing == "none":  # with patsy it's drop or raise. let's raise.
            missing = "raise"

        tmp = handle_formula_data(data, None, formula, depth=eval_env, missing=missing)
        ((endog, exog), missing_idx, model_spec) = tmp
        max_endog = cls._formula_max_endog
        if max_endog is not None and endog.ndim > 1 and endog.shape[1] > max_endog:
            raise ValueError(
                "endog has evaluated to an array with multiple "
                "columns that has shape {}. This occurs when "
                "the variable converted to endog is non-numeric"
                " (e.g., bool or str).".format(endog.shape)
            )
        if drop_cols is not None and len(drop_cols) > 0:
            cols = [x for x in exog.columns if x not in drop_cols]
            if len(cols) < len(exog.columns):
                exog = exog[cols]
                spec_cols = list(mgr.get_term_names(model_spec))
                for col in drop_cols:
                    try:
                        if mgr.engine == "formulaic" and col == "Intercept":
                            col = "1"
                        spec_cols.remove(col)
                    except ValueError:
                        pass  # OK if not present
                # TODO: Patsy migration, need to add method to handle
                model_spec = model_spec.subset(spec_cols)

        kwargs.update(
            {
                "missing_idx": missing_idx,
                "missing": missing,
                "formula": formula,  # attach formula for unpckling
                "model_spec": model_spec,
            }
        )
        mod = cls(endog, exog, *args, **kwargs)
        mod.formula = formula
        # since we got a dataframe, attach the original
        mod.data.frame = data
        return mod

    @property
    def endog_names(self):
        """
        Names of endogenous variables.
        """
        return self.data.ynames

    @property
    def exog_names(self) -> list[str] | None:
        """
        Names of exogenous variables.
        """
        return self.data.xnames

    def fit(self):
        """
        Fit a model to data.
        """
        raise NotImplementedError

    def predict(self, params, exog=None, *args, **kwargs):
        """
        After a model has been fit predict returns the fitted values.

        This is a placeholder intended to be overwritten by individual models.
        """
        raise NotImplementedError


class LikelihoodModel(Model):
    """
    Likelihood model is a subclass of Model.
    """

    def __init__(self, endog, exog=None, **kwargs):
        super().__init__(endog, exog, **kwargs)
        self.initialize()

    def initialize(self):
        """
        Initialize (possibly re-initialize) a Model instance.

        For example, if the the design matrix of a linear model changes then
        initialized can be used to recompute values using the modified design
        matrix.
        """

    # TODO: if the intent is to re-initialize the model with new data then this
    # method needs to take inputs...

    def loglike(self, params):
        """
        Log-likelihood of model.

        Parameters
        ----------
        params : ndarray
            The model parameters used to compute the log-likelihood.

        Notes
        -----
        Must be overridden by subclasses.
        """
        raise NotImplementedError

    def score(self, params):
        """
        Score vector of model.

        The gradient of logL with respect to each parameter.

        Parameters
        ----------
        params : ndarray
            The parameters to use when evaluating the Hessian.

        Returns
        -------
        ndarray
            The score vector evaluated at the parameters.
        """
        raise NotImplementedError

    def information(self, params):
        """
        Fisher information matrix of model.

        Returns -1 * Hessian of the log-likelihood evaluated at params.

        Parameters
        ----------
        params : ndarray
            The model parameters.
        """
        raise NotImplementedError

    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.
        """
        raise NotImplementedError

    def fit(
        self,
        start_params=None,
        method="newton",
        maxiter=100,
        full_output=True,
        disp=True,
        fargs=(),
        callback=None,
        retall=False,
        skip_hessian=False,
        **kwargs,
    ):
        """
        Fit method for likelihood based models

        Parameters
        ----------
        start_params : array_like, optional
            Initial guess of the solution for the loglikelihood maximization.
            The default is an array of zeros.
        method : str, optional
            The `method` determines which solver from `scipy.optimize`
            is used, and it can be chosen from among the following strings:

            - 'newton' for Newton-Raphson, 'nm' for Nelder-Mead
            - 'bfgs' for Broyden-Fletcher-Goldfarb-Shanno (BFGS)
            - 'lbfgs' for limited-memory BFGS with optional box constraints
            - 'powell' for modified Powell's method
            - 'cg' for conjugate gradient
            - 'ncg' for Newton-conjugate gradient
            - 'basinhopping' for global basin-hopping solver
            - 'minimize' for generic wrapper of scipy minimize (BFGS by default)

            The explicit arguments in `fit` are passed to the solver,
            with the exception of the basin-hopping solver. Each
            solver has several optional arguments that are not the same across
            solvers. See the notes section below (or scipy.optimize) for the
            available arguments and for the list of explicit arguments that the
            basin-hopping solver supports.
        maxiter : int, optional
            The maximum number of iterations to perform.
        full_output : bool, optional
            Set to True to have all available output in the Results object's
            mle_retvals attribute. The output is dependent on the solver.
            See LikelihoodModelResults notes section for more information.
        disp : bool, optional
            Set to True to print convergence messages.
        fargs : tuple, optional
            Extra arguments passed to the likelihood function, i.e.,
            loglike(x,*args)
        callback : callable callback(xk), optional
            Called after each iteration, as callback(xk), where xk is the
            current parameter vector.
        retall : bool, optional
            Set to True to return list of solutions at each iteration.
            Available in Results object's mle_retvals attribute.
        skip_hessian : bool, optional
            If False (default), then the negative inverse hessian is calculated
            after the optimization. If True, then the hessian will not be
            calculated. However, it will be available in methods that use the
            hessian in the optimization (currently only with `"newton"`).
        kwargs : keywords
            All kwargs are passed to the chosen solver with one exception. The
            following keyword controls what happens after the fit::

                warn_convergence : bool, optional
                    If True, checks the model for the converged flag. If the
                    converged flag is False, a ConvergenceWarning is issued.

        Notes
        -----
        The 'basinhopping' solver ignores `maxiter`, `retall`, `full_output`
        explicit arguments.

        Optional arguments for solvers (see returned Results.mle_settings)::

            'newton'
                tol : float
                    Relative error in params acceptable for convergence.
            'nm' -- Nelder Mead
                xtol : float
                    Relative error in params acceptable for convergence
                ftol : float
                    Relative error in loglike(params) acceptable for
                    convergence
                maxfun : int
                    Maximum number of function evaluations to make.
            'bfgs'
                gtol : float
                    Stop when norm of gradient is less than gtol.
                norm : float
                    Order of norm (np.inf is max, -np.inf is min)
                epsilon
                    If fprime is approximated, use this value for the step
                    size. Only relevant if LikelihoodModel.score is None.
            'lbfgs'
                m : int
                    This many terms are used for the Hessian approximation.
                factr : float
                    A stop condition that is a variant of relative error.
                pgtol : float
                    A stop condition that uses the projected gradient.
                epsilon
                    If fprime is approximated, use this value for the step
                    size. Only relevant if LikelihoodModel.score is None.
                maxfun : int
                    Maximum number of function evaluations to make.
                bounds : sequence
                    (min, max) pairs for each element in x,
                    defining the bounds on that parameter.
                    Use None for one of min or max when there is no bound
                    in that direction.
            'cg'
                gtol : float
                    Stop when norm of gradient is less than gtol.
                norm : float
                    Order of norm (np.inf is max, -np.inf is min)
                epsilon : float
                    If fprime is approximated, use this value for the step
                    size. Can be scalar or vector.  Only relevant if
                    Likelihoodmodel.score is None.
            'ncg'
                fhess_p : callable f'(x,*args)
                    Function which computes the Hessian of f times an arbitrary
                    vector, p.  Should only be supplied if
                    LikelihoodModel.hessian is None.
                avextol : float
                    Stop when the average relative error in the minimizer
                    falls below this amount.
                epsilon : float or ndarray
                    If fhess is approximated, use this value for the step size.
                    Only relevant if Likelihoodmodel.hessian is None.
            'powell'
                xtol : float
                    Line-search error tolerance
                ftol : float
                    Relative error in loglike(params) for acceptable for
                    convergence.
                maxfun : int
                    Maximum number of function evaluations to make.
                start_direc : ndarray
                    Initial direction set.
            'basinhopping'
                niter : int
                    The number of basin hopping iterations.
                niter_success : int
                    Stop the run if the global minimum candidate remains the
                    same for this number of iterations.
                T : float
                    The "temperature" parameter for the accept or reject
                    criterion. Higher "temperatures" mean that larger jumps
                    in function value will be accepted. For best results
                    `T` should be comparable to the separation (in function
                    value) between local minima.
                stepsize : float
                    Initial step size for use in the random displacement.
                interval : int
                    The interval for how often to update the `stepsize`.
                minimizer : dict
                    Extra keyword arguments to be passed to the minimizer
                    `scipy.optimize.minimize()`, for example 'method' - the
                    minimization method (e.g. 'L-BFGS-B'), or 'tol' - the
                    tolerance for termination. Other arguments are mapped from
                    explicit argument of `fit`:
                      - `args` <- `fargs`
                      - `jac` <- `score`
                      - `hess` <- `hess`
            'minimize'
                min_method : str, optional
                    Name of minimization method to use.
                    Any method specific arguments can be passed directly.
                    For a list of methods and their arguments, see
                    documentation of `scipy.optimize.minimize`.
                    If no method is specified, then BFGS is used.
        """
        Hinv = None  # JP error if full_output=0, Hinv not defined

        if start_params is None:
            if hasattr(self, "start_params"):
                start_params = self.start_params
            elif self.exog is not None:
                # fails for shape (K,)?
                start_params = [0.0] * self.exog.shape[1]
            else:
                raise ValueError(
                    "If exog is None, then start_params should be specified"
                )

        # TODO: separate args from nonarg taking score and hessian, ie.,
        # user-supplied and numerically evaluated estimate frprime does not take
        # args in most (any?) of the optimize function

        nobs = self.endog.shape[0]
        # f = lambda params, *args: -self.loglike(params, *args) / nobs

        def f(params, *args):
            return -self.loglike(params, *args) / nobs

        if method == "newton":
            # TODO: why are score and hess positive?
            def score(params, *args):
                return self.score(params, *args) / nobs

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

        else:

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

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

        warn_convergence = kwargs.pop("warn_convergence", True)

        # Remove covariance args before calling fir to allow strict checking
        if "cov_type" in kwargs:
            cov_kwds = kwargs.get("cov_kwds", {})
            kwds = {"cov_type": kwargs["cov_type"], "cov_kwds": cov_kwds}
            if cov_kwds:
                del kwargs["cov_kwds"]
            del kwargs["cov_type"]
        else:
            kwds = {}
        if "use_t" in kwargs:
            kwds["use_t"] = kwargs["use_t"]
            del kwargs["use_t"]

        optimizer = Optimizer()
        xopt, retvals, optim_settings = optimizer._fit(
            f,
            score,
            start_params,
            fargs,
            kwargs,
            hessian=hess,
            method=method,
            disp=disp,
            maxiter=maxiter,
            callback=callback,
            retall=retall,
            full_output=full_output,
        )
        # Restore cov_type, cov_kwds and use_t
        optim_settings.update(kwds)
        # NOTE: this is for fit_regularized and should be generalized
        cov_params_func = kwargs.setdefault("cov_params_func", None)
        if cov_params_func:
            Hinv = cov_params_func(self, xopt, retvals)
        elif method == "newton" and full_output:
            Hinv = np.linalg.inv(-retvals["Hessian"]) / nobs
        elif not skip_hessian:
            H = -1 * self.hessian(xopt)
            invertible = False
            if np.all(np.isfinite(H)):
                eigvals, eigvecs = np.linalg.eigh(H)
                if np.min(eigvals) > 0:
                    invertible = True

            if invertible:
                Hinv = eigvecs.dot(np.diag(1.0 / eigvals)).dot(eigvecs.T)
                Hinv = np.asfortranarray((Hinv + Hinv.T) / 2.0)
            else:
                warnings.warn(
                    "Inverting hessian failed, no bse or cov_params available",
                    HessianInversionWarning,
                    stacklevel=2,
                )
                Hinv = None

        # TODO: add Hessian approximation and change the above if needed
        mlefit = LikelihoodModelResults(self, xopt, Hinv, scale=1.0, **kwds)

        # TODO: hardcode scale?
        mlefit.mle_retvals = retvals
        if isinstance(retvals, dict):
            if warn_convergence and not retvals["converged"]:
                from statsmodels.tools.sm_exceptions import ConvergenceWarning

                warnings.warn(
                    "Maximum Likelihood optimization failed to "
                    "converge. Check mle_retvals",
                    ConvergenceWarning,
                    stacklevel=2,
                )

        mlefit.mle_settings = optim_settings
        return mlefit

    def _fit_zeros(
        self,
        keep_index=None,
        start_params=None,
        return_auxiliary=False,
        k_params=None,
        **fit_kwds,
    ):
        """experimental, fit the model subject to zero constraints

        Intended for internal use cases until we know what we need.
        API will need to change to handle models with two exog.
        This is not yet supported by all model subclasses.

        This is essentially a simplified version of `fit_constrained`, and
        does not need to use `offset`.

        The estimation creates a new model with transformed design matrix,
        exog, and converts the results back to the original parameterization.

        Some subclasses could use a more efficient calculation than using a
        new model.

        Parameters
        ----------
        keep_index : array_like (int or bool) or slice
            variables that should be dropped.
        start_params : None or array_like
            starting values for the optimization. `start_params` needs to be
            given in the original parameter space and are internally
            transformed.
        k_params : int or None
            If None, then we try to infer from start_params or model.
        **fit_kwds : keyword arguments
            fit_kwds are used in the optimization of the transformed model.

        Returns
        -------
        results : Results instance
        """
        # we need to append index of extra params to keep_index as in
        # NegativeBinomial
        if hasattr(self, "k_extra") and self.k_extra > 0:
            # we cannot change the original, TODO: should we add keep_index_params?
            keep_index = np.array(keep_index, copy=True)
            k = self.exog.shape[1]
            extra_index = np.arange(k, k + self.k_extra)
            keep_index_p = np.concatenate((keep_index, extra_index))
        else:
            keep_index_p = keep_index

        # not all models support start_params, drop if None, hide them in fit_kwds
        if start_params is not None:
            fit_kwds["start_params"] = start_params[keep_index_p]
            k_params = len(start_params)
            # ignore k_params in this case, or verify consisteny?

        # build auxiliary model and fit
        init_kwds = self._get_init_kwds()
        mod_constr = self.__class__(self.endog, self.exog[:, keep_index], **init_kwds)
        res_constr = mod_constr.fit(**fit_kwds)
        # switch name, only need keep_index for params below
        keep_index = keep_index_p

        if k_params is None:
            k_params = self.exog.shape[1]
            k_params += getattr(self, "k_extra", 0)

        params_full = np.zeros(k_params)
        params_full[keep_index] = res_constr.params

        # create dummy results Instance, TODO: wire up properly
        # TODO: this could be moved into separate private method if needed
        # discrete L1 fit_regularized doens't reestimate AFAICS
        # RLM does not have method, disp nor warn_convergence keywords
        # OLS, WLS swallows extra kwds with **kwargs, but does not have method='nm'
        try:
            # Note: addding full_output=False causes exceptions
            res = self.fit(
                maxiter=0,
                disp=0,
                method="nm",
                skip_hessian=True,
                warn_convergence=False,
                start_params=params_full,
            )
            # we get a wrapper back
        except (TypeError, ValueError):
            res = self.fit()

        # Warning: make sure we are not just changing the wrapper instead of
        # results #2400
        # TODO: do we need to change res._results.scale in some models?
        if hasattr(res_constr.model, "scale"):
            # Note: res.model is self
            # GLM problem, see #2399,
            # TODO: remove from model if not needed anymore
            res.model.scale = res._results.scale = res_constr.model.scale

        if hasattr(res_constr, "mle_retvals"):
            res._results.mle_retvals = res_constr.mle_retvals
            # not available for not scipy optimization, e.g. glm irls
            # TODO: what retvals should be required?
            # res.mle_retvals['fcall'] = res_constr.mle_retvals.get('fcall', np.nan)
            # res.mle_retvals['iterations'] = res_constr.mle_retvals.get(
            #                                                 'iterations', np.nan)
            # res.mle_retvals['converged'] = res_constr.mle_retvals['converged']
        # overwrite all mle_settings
        if hasattr(res_constr, "mle_settings"):
            res._results.mle_settings = res_constr.mle_settings

        res._results.params = params_full
        if (
            not hasattr(res._results, "normalized_cov_params")
            or res._results.normalized_cov_params is None
        ):
            res._results.normalized_cov_params = np.zeros((k_params, k_params))
        else:
            res._results.normalized_cov_params[...] = 0

        # fancy indexing requires integer array
        keep_index = np.array(keep_index)
        res._results.normalized_cov_params[keep_index[:, None], keep_index] = (
            res_constr.normalized_cov_params
        )
        k_constr = res_constr.df_resid - res._results.df_resid
        if hasattr(res_constr, "cov_params_default"):
            res._results.cov_params_default = np.zeros((k_params, k_params))
            res._results.cov_params_default[keep_index[:, None], keep_index] = (
                res_constr.cov_params_default
            )
        if hasattr(res_constr, "cov_type"):
            res._results.cov_type = res_constr.cov_type
            res._results.cov_kwds = res_constr.cov_kwds

        res._results.keep_index = keep_index
        res._results.df_resid = res_constr.df_resid
        res._results.df_model = res_constr.df_model

        res._results.k_constr = k_constr
        res._results.results_constrained = res_constr

        # special temporary workaround for RLM
        # need to be able to override robust covariances
        if hasattr(res.model, "M"):
            del res._results._cache["resid"]
            del res._results._cache["fittedvalues"]
            del res._results._cache["sresid"]
            cov = res._results._cache["bcov_scaled"]
            # inplace adjustment
            cov[...] = 0
            cov[keep_index[:, None], keep_index] = res_constr.bcov_scaled
            res._results.cov_params_default = cov

        return res

    def _fit_collinear(self, atol=1e-14, rtol=1e-13, **kwds):
        """experimental, fit of the model without collinear variables

        This currently uses QR to drop variables based on the given
        sequence.
        Options will be added in future, when the supporting functions
        to identify collinear variables become available.
        """

        # ------ copied from PR #2380 remove when merged
        x = self.exog
        tol = atol + rtol * x.var(0)
        r = np.linalg.qr(x, mode="r")
        mask = np.abs(r.diagonal()) < np.sqrt(tol)
        # TODO add to results instance
        # idx_collinear = np.where(mask)[0]
        idx_keep = np.where(~mask)[0]
        return self._fit_zeros(keep_index=idx_keep, **kwds)


# TODO: the below is unfinished
class GenericLikelihoodModel(LikelihoodModel):
    """
    Allows the fitting of any likelihood function via maximum likelihood.

    A subclass needs to specify at least the log-likelihood
    If the log-likelihood is specified for each observation, then results that
    require the Jacobian will be available. (The other case is not tested yet.)

    Notes
    -----
    Optimization methods that require only a likelihood function are 'nm' and
    'powell'

    Optimization methods that require a likelihood function and a
    score/gradient are 'bfgs', 'cg', and 'ncg'. A function to compute the
    Hessian is optional for 'ncg'.

    Optimization method that require a likelihood function, a score/gradient,
    and a Hessian is 'newton'

    If they are not overwritten by a subclass, then numerical gradient,
