"""
Generalized linear models currently supports estimation using the one-parameter
exponential families

References
----------
Gill, Jeff. 2000. Generalized Linear Models: A Unified Approach.
    SAGE QASS Series.

Green, PJ. 1984.  "Iteratively reweighted least squares for maximum
    likelihood estimation, and some robust and resistant alternatives."
    Journal of the Royal Statistical Society, Series B, 46, 149-192.

Hardin, J.W. and Hilbe, J.M. 2007.  "Generalized Linear Models and
    Extensions."  2nd ed.  Stata Press, College Station, TX.

McCullagh, P. and Nelder, J.A.  1989.  "Generalized Linear Models." 2nd ed.
    Chapman & Hall, Boca Rotan.
"""

from statsmodels.compat.pandas import Appender

import warnings

import numpy as np
from numpy.linalg import LinAlgError

from statsmodels.base import _prediction_inference as pred
import statsmodels.base._parameter_inference as pinfer
from statsmodels.base._prediction_inference import PredictionResultsMean
import statsmodels.base.model as base
import statsmodels.base.wrapper as wrap
from statsmodels.graphics._regressionplots_doc import (
    _plot_added_variable_doc,
    _plot_ceres_residuals_doc,
    _plot_partial_residuals_doc,
)
import statsmodels.regression._tools as reg_tools
import statsmodels.regression.linear_model as lm
from statsmodels.tools.data import _as_array_with_name
from statsmodels.tools.decorators import (
    cache_readonly,
    cached_data,
    cached_value,
)
from statsmodels.tools.docstring import Docstring
from statsmodels.tools.sm_exceptions import (
    ConvergenceWarning,
    DomainWarning,
    HessianInversionWarning,
    PerfectSeparationWarning,
)
from statsmodels.tools.validation import float_like

# need import in module instead of lazily to copy `__doc__`
from . import families

__all__ = ["GLM", "PredictionResultsMean"]


def _check_convergence(criterion, iteration, atol, rtol):
    return np.allclose(
        criterion[iteration], criterion[iteration + 1], atol=atol, rtol=rtol
    )


# Remove after 0.13 when bic changes to bic llf
class _ModuleVariable:
    _value = None

    @property
    def use_bic_llf(self):
        return self._value

    def set_use_bic_llf(self, val):
        if val not in (True, False, None):
            raise ValueError("Must be True, False or None")
        self._value = bool(val) if val is not None else val


_use_bic_helper = _ModuleVariable()
SET_USE_BIC_LLF = _use_bic_helper.set_use_bic_llf


class GLM(base.LikelihoodModel):
    __doc__ = """
    Generalized Linear Models

    GLM inherits from statsmodels.base.model.LikelihoodModel

    Parameters
    ----------
    endog : array_like
        1d array of endogenous response variable.  This array can be 1d or 2d.
        Binomial family models accept a 2d array with two columns. If
        supplied, each observation is expected to be [success, failure].
    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 (models specified using a formula
        include an intercept by default). See `statsmodels.tools.add_constant`.
    family : family class instance
        The default is Gaussian.  To specify the binomial distribution
        family = sm.family.Binomial()
        Each family can take a link instance as an argument.  See
        statsmodels.family.family for more information.
    offset : array_like or None
        An offset to be included in the model.  If provided, must be
        an array whose length is the number of rows in exog.
    exposure : array_like or None
        Log(exposure) will be added to the linear prediction in the model.
        Exposure is only valid if the log link is used. If provided, it must be
        an array with the same length as endog.
    freq_weights : array_like
        1d array of frequency weights. The default is None. If None is selected
        or a blank value, then the algorithm will replace with an array of 1's
        with length equal to the endog.
        WARNING: Using weights is not verified yet for all possible options
        and results, see Notes.
    var_weights : array_like
        1d array of variance (analytic) weights. The default is None. If None
        is selected or a blank value, then the algorithm will replace with an
        array of 1's with length equal to the endog.
        WARNING: Using weights is not verified yet for all possible options
        and results, see Notes.
    {extra_params}

    Attributes
    ----------
    df_model : float
        Model degrees of freedom is equal to p - 1, where p is the number
        of regressors.  Note that the intercept is not reported as a
        degree of freedom.
    df_resid : float
        Residual degrees of freedom is equal to the number of observation n
        minus the number of regressors p.
    endog : ndarray
        See Notes.  Note that `endog` is a reference to the data so that if
        data is already an array and it is changed, then `endog` changes
        as well.
    exposure : array_like
        Include ln(exposure) in model with coefficient constrained to 1. Can
        only be used if the link is the logarithm function.
    exog : ndarray
        See Notes.  Note that `exog` is a reference to the data so that if
        data is already an array and it is changed, then `exog` changes
        as well.
    freq_weights : ndarray
        See Notes. Note that `freq_weights` is a reference to the data so that
        if data is already an array and it is changed, then `freq_weights`
        changes as well.
    var_weights : ndarray
        See Notes. Note that `var_weights` is a reference to the data so that
        if data is already an array and it is changed, then `var_weights`
        changes as well.
    iteration : int
        The number of iterations that fit has run.  Initialized at 0.
    family : family class instance
        The distribution family of the model. Can be any family in
        statsmodels.families.  Default is Gaussian.
    mu : ndarray
        The mean response of the transformed variable.  `mu` is the value of
        the inverse of the link function at lin_pred, where lin_pred is the
        linear predicted value of the WLS fit of the transformed variable.
        `mu` is only available after fit is called.  See
        statsmodels.families.family.fitted of the distribution family for more
        information.
    n_trials : ndarray
        See Notes. Note that `n_trials` is a reference to the data so that if
        data is already an array and it is changed, then `n_trials` changes
        as well. `n_trials` is the number of binomial trials and only available
        with that distribution. See statsmodels.families.Binomial for more
        information.
    normalized_cov_params : ndarray
        The p x p normalized covariance of the design / exogenous data.
        This is approximately equal to (X.T X)^(-1)
    offset : array_like
        Include offset in model with coefficient constrained to 1.
    scale : float
        The estimate of the scale / dispersion of the model fit.  Only
        available after fit is called.  See GLM.fit and GLM.estimate_scale
        for more information.
    scaletype : str
        The scaling used for fitting the model.  This is only available after
        fit is called.  The default is None.  See GLM.fit for more information.
    weights : ndarray
        The value of the weights after the last iteration of fit.  Only
        available after fit is called.  See statsmodels.families.family for
        the specific distribution weighting functions.

    Examples
    --------
    >>> import statsmodels.api as sm
    >>> data = sm.datasets.scotland.load()
    >>> data.exog = sm.add_constant(data.exog)

    Instantiate a gamma family model with the default link function.

    >>> gamma_model = sm.GLM(data.endog, data.exog,
    ...                      family=sm.families.Gamma())

    >>> gamma_results = gamma_model.fit()
    >>> gamma_results.params
    array([-0.01776527,  0.00004962,  0.00203442, -0.00007181,  0.00011185,
           -0.00000015, -0.00051868, -0.00000243])
    >>> gamma_results.scale
    0.0035842831734919055
    >>> gamma_results.deviance
    0.087388516416999198
    >>> gamma_results.pearson_chi2
    0.086022796163805704
    >>> gamma_results.llf
    -83.017202161073527

    See Also
    --------
    statsmodels.genmod.families.family.Family
    :ref:`families`
    :ref:`links`

    Notes
    -----
    Note: PerfectSeparationError exception has been converted to a
    PerfectSeparationWarning and perfect separation or perfect prediction will
    not raise an exception by default. (changed in version 0.14)

    Only the following combinations make sense for family and link:

     ============= ===== === ===== ====== ======= === ==== ====== ====== ====
     Family        ident log logit probit cloglog pow opow nbinom loglog logc
     ============= ===== === ===== ====== ======= === ==== ====== ====== ====
     Gaussian      x     x   x     x      x       x   x     x      x
     inv Gaussian  x     x                        x
     binomial      x     x   x     x      x       x   x           x      x
     Poisson       x     x                        x
     neg binomial  x     x                        x        x
     gamma         x     x                        x
     Tweedie       x     x                        x
     ============= ===== === ===== ====== ======= === ==== ====== ====== ====

    Not all of these link functions are currently available.

    Endog and exog are references so that if the data they refer to are already
    arrays and these arrays are changed, endog and exog will change.

    statsmodels supports two separate definitions of weights: frequency weights
    and variance weights.

    Frequency weights produce the same results as repeating observations by the
    frequencies (if those are integers). Frequency weights will keep the number
    of observations consistent, but the degrees of freedom will change to
    reflect the new weights.

    Variance weights (referred to in other packages as analytic weights) are
    used when ``endog`` represents an an average or mean. This relies on the
    assumption that that the inverse variance scales proportionally to the
    weight--an observation that is deemed more credible should have less
    variance and therefore have more weight. For the ``Poisson`` family--which
    assumes that occurrences scale proportionally with time--a natural practice
    would be to use the amount of time as the variance weight and set ``endog``
    to be a rate (occurrences per period of time). Similarly, using a
    compound Poisson family, namely ``Tweedie``, makes a similar assumption
    about the rate (or frequency) of occurrences having variance proportional to
    time.

    Both frequency and variance weights are verified for all basic results with
    nonrobust or heteroscedasticity robust ``cov_type``. Other robust
    covariance types have not yet been verified, and at least the small sample
    correction is currently not based on the correct total frequency count.

    Currently, all residuals are not weighted by frequency, although they may
    incorporate ``n_trials`` for ``Binomial`` and ``var_weights``

    +---------------+----------------------------------+
    | Residual Type | Applicable weights               |
    +===============+==================================+
    | Anscombe      | ``var_weights``                  |
    +---------------+----------------------------------+
    | Deviance      | ``var_weights``                  |
    +---------------+----------------------------------+
    | Pearson       | ``var_weights`` and ``n_trials`` |
    +---------------+----------------------------------+
    | Reponse       | ``n_trials``                     |
    +---------------+----------------------------------+
    | Working       | ``n_trials``                     |
    +---------------+----------------------------------+

    WARNING: Loglikelihood and deviance are not valid in models where
    scale is equal to 1 (i.e., ``Binomial``, ``NegativeBinomial``, and
    ``Poisson``). If variance weights are specified, then results such as
    ``loglike`` and ``deviance`` are based on a quasi-likelihood
    interpretation. The loglikelihood is not correctly specified in this case,
    and statistics based on it, such AIC or likelihood ratio tests, are not
    appropriate.
    """.format(
        extra_params=base._missing_param_doc
    )
    # Maximum number of endogenous variables when using a formula
    _formula_max_endog = 2

    def __init__(
        self,
        endog,
        exog,
        family=None,
        offset=None,
        exposure=None,
        freq_weights=None,
        var_weights=None,
        missing="none",
        **kwargs,
    ):

        if type(self) is GLM:
            self._check_kwargs(kwargs, ["n_trials"])

        if (family is not None) and not isinstance(
            family.link, tuple(family.safe_links)
        ):

            warnings.warn(
                (
                    f"The {type(family.link).__name__} link function "
                    "does not respect the domain of the "
                    f"{type(family).__name__} family."
                ),
                DomainWarning,
                stacklevel=2,
            )

        self._exposure_name = None
        self._offset_name = None
        self._freq_weights_name = None
        self._var_weights_name = None

        if exposure is not None:
            exposure_array, self._exposure_name = _as_array_with_name(
                exposure, "exposure"
            )
            exposure = np.log(exposure_array)
        if offset is not None:  # this should probably be done upstream
            offset, self._offset_name = _as_array_with_name(offset, "offset")

        if freq_weights is not None:
            freq_weights, self._freq_weights_name = _as_array_with_name(
                freq_weights, "freq_weights"
            )
        if var_weights is not None:
            var_weights, self._var_weights_name = _as_array_with_name(
                var_weights, "var_weights"
            )

        self.freq_weights = freq_weights
        self.var_weights = var_weights

        super().__init__(
            endog,
            exog,
            missing=missing,
            offset=offset,
            exposure=exposure,
            freq_weights=freq_weights,
            var_weights=var_weights,
            **kwargs,
        )
        self._check_inputs(
            family,
            self.offset,
            self.exposure,
            self.endog,
            self.freq_weights,
            self.var_weights,
        )
        if offset is None:
            del self.offset
        if exposure is None:
            del self.exposure

        self.nobs = self.endog.shape[0]

        # things to remove_data
        self._data_attr.extend(
            [
                "weights",
                "mu",
                "freq_weights",
                "var_weights",
                "iweights",
                "_offset_exposure",
                "n_trials",
            ]
        )
        # register kwds for __init__, offset and exposure are added by super
        self._init_keys.append("family")

        self._setup_binomial()
        # internal usage for recreating a model
        if "n_trials" in kwargs:
            self.n_trials = kwargs["n_trials"]

        # Construct a combined offset/exposure term.  Note that
        # exposure has already been logged if present.
        offset_exposure = 0.0
        if hasattr(self, "offset"):
            offset_exposure = self.offset
        if hasattr(self, "exposure"):
            offset_exposure = offset_exposure + self.exposure
        self._offset_exposure = offset_exposure

        self.scaletype = None

    def initialize(self):
        """
        Initialize a generalized linear model.
        """
        self.df_model = np.linalg.matrix_rank(self.exog) - 1

        if (self.freq_weights is not None) and (
            self.freq_weights.shape[0] == self.endog.shape[0]
        ):
            self.wnobs = self.freq_weights.sum()
            self.df_resid = self.wnobs - self.df_model - 1
        else:
            self.wnobs = self.exog.shape[0]
            self.df_resid = self.exog.shape[0] - self.df_model - 1

    def _check_inputs(self, family, offset, exposure, endog, freq_weights, var_weights):

        # Default family is Gaussian
        if family is None:
            family = families.Gaussian()
        self.family = family

        if exposure is not None:
            if not isinstance(self.family.link, families.links.Log):
                raise ValueError(
                    "exposure can only be used with the log link function"
                )
            elif exposure.shape[0] != endog.shape[0]:
                raise ValueError("exposure is not the same length as endog")

        if offset is not None:
            if offset.shape[0] != endog.shape[0]:
                raise ValueError("offset is not the same length as endog")

        if freq_weights is not None:
            if freq_weights.shape[0] != endog.shape[0]:
                raise ValueError("freq weights not the same length as endog")
            if len(freq_weights.shape) > 1:
                raise ValueError("freq weights has too many dimensions")

        # internal flag to store whether freq_weights were not None
        self._has_freq_weights = self.freq_weights is not None
        if self.freq_weights is None:
            self.freq_weights = np.ones(endog.shape[0])
            # TODO: check do we want to keep None as sentinel for freq_weights

        if np.shape(self.freq_weights) == () and self.freq_weights > 1:
            self.freq_weights = self.freq_weights * np.ones(endog.shape[0])

        if var_weights is not None:
            if var_weights.shape[0] != endog.shape[0]:
                raise ValueError("var weights not the same length as endog")
            if len(var_weights.shape) > 1:
                raise ValueError("var weights has too many dimensions")

        # internal flag to store whether var_weights were not None
        self._has_var_weights = var_weights is not None
        if var_weights is None:
            self.var_weights = np.ones(endog.shape[0])
            # TODO: check do we want to keep None as sentinel for var_weights
        self.iweights = np.asarray(self.freq_weights * self.var_weights)

    def _get_init_kwds(self):
        # this is a temporary fixup because exposure has been transformed
        # see #1609, copied from discrete_model.CountModel
        kwds = super()._get_init_kwds()
        if "exposure" in kwds and kwds["exposure"] is not None:
            kwds["exposure"] = np.exp(kwds["exposure"])
        return kwds

    def loglike_mu(self, mu, scale=1.0):
        """
        Evaluate the log-likelihood for a generalized linear model.
        """
        scale = float_like(scale, "scale")
        return self.family.loglike(
            self.endog, mu, self.var_weights, self.freq_weights, scale
        )

    def loglike(self, params, scale=None):
        """
        Evaluate the log-likelihood for a generalized linear model.
        """
        scale = float_like(scale, "scale", optional=True)
        lin_pred = np.dot(self.exog, params) + self._offset_exposure
        expval = self.family.link.inverse(lin_pred)
        if scale is None:
            scale = self.estimate_scale(expval)
        llf = self.family.loglike(
            self.endog, expval, self.var_weights, self.freq_weights, scale
        )
        return llf

    def score_obs(self, params, scale=None):
        """score first derivative of the loglikelihood for each observation.

        Parameters
        ----------
        params : ndarray
            Parameter at which score is evaluated.
        scale : None or float
            If scale is None, then the default scale will be calculated.
            Default scale is defined by `self.scaletype` and set in fit.
            If scale is not None, then it is used as a fixed scale.

        Returns
        -------
        score_obs : ndarray, 2d
            The first derivative of the loglikelihood function evaluated at
            params for each observation.
        """
        scale = float_like(scale, "scale", optional=True)
        score_factor = self.score_factor(params, scale=scale)
        return score_factor[:, None] * self.exog

    def score(self, params, scale=None):
        """score, first derivative of the loglikelihood function

        Parameters
        ----------
        params : ndarray
            Parameter at which score is evaluated.
        scale : None or float
            If scale is None, then the default scale will be calculated.
            Default scale is defined by `self.scaletype` and set in fit.
            If scale is not None, then it is used as a fixed scale.

        Returns
        -------
        score : ndarray_1d
            The first derivative of the loglikelihood function calculated as
            the sum of `score_obs`
        """
        scale = float_like(scale, "scale", optional=True)
        score_factor = self.score_factor(params, scale=scale)
        return np.dot(score_factor, self.exog)

    def score_factor(self, params, scale=None):
        """weights for score for each observation

        This can be considered as score residuals.

        Parameters
        ----------
        params : ndarray
            parameter at which score is evaluated
        scale : None or float
            If scale is None, then the default scale will be calculated.
            Default scale is defined by `self.scaletype` and set in fit.
            If scale is not None, then it is used as a fixed scale.

        Returns
        -------
        score_factor : ndarray_1d
            A 1d weight vector used in the calculation of the score_obs.
            The score_obs are obtained by `score_factor[:, None] * exog`
        """
        scale = float_like(scale, "scale", optional=True)
        mu = self.predict(params)
        if scale is None:
            scale = self.estimate_scale(mu)

        score_factor = (self.endog - mu) / self.family.link.deriv(mu)
        score_factor /= self.family.variance(mu)
        score_factor *= self.iweights * self.n_trials

        if not scale == 1:
            score_factor /= scale

        return score_factor

    def hessian_factor(self, params, scale=None, observed=True):
        """Weights for calculating Hessian

        Parameters
        ----------
        params : ndarray
            parameter at which Hessian is evaluated
        scale : None or float
            If scale is None, then the default scale will be calculated.
            Default scale is defined by `self.scaletype` and set in fit.
            If scale is not None, then it is used as a fixed scale.
        observed : bool
            If True, then the observed Hessian is returned. If false then the
            expected information matrix is returned.

        Returns
        -------
        hessian_factor : ndarray, 1d
            A 1d weight vector used in the calculation of the Hessian.
            The hessian is obtained by `(exog.T * hessian_factor).dot(exog)`
        """

        # calculating eim_factor
        mu = self.predict(params)
        if scale is None:
            scale = self.estimate_scale(mu)

        eim_factor = 1 / (self.family.link.deriv(mu) ** 2 * self.family.variance(mu))
        eim_factor *= self.iweights * self.n_trials

        if not observed:
            if not scale == 1:
                eim_factor /= scale
            return eim_factor

        # calculating oim_factor, eim_factor is with scale=1

        score_factor = self.score_factor(params, scale=1.0)
        if eim_factor.ndim > 1 or score_factor.ndim > 1:
            raise RuntimeError("something wrong")

        tmp = self.family.variance(mu) * self.family.link.deriv2(mu)
        tmp += self.family.variance.deriv(mu) * self.family.link.deriv(mu)

        tmp = score_factor * tmp
        # correct for duplicatee iweights in oim_factor and score_factor
        tmp /= self.iweights * self.n_trials
        oim_factor = eim_factor * (1 + tmp)

        if tmp.ndim > 1:
            raise RuntimeError("something wrong")

        if not scale == 1:
            oim_factor /= scale

        return oim_factor

    def hessian(self, params, scale=None, observed=None):
        """Hessian, second derivative of loglikelihood function

        Parameters
        ----------
        params : ndarray
            parameter at which Hessian is evaluated
        scale : None or float
            If scale is None, then the default scale will be calculated.
            Default scale is defined by `self.scaletype` and set in fit.
            If scale is not None, then it is used as a fixed scale.
        observed : bool
            If True, then the observed Hessian is returned (default).
            If False, then the expected information matrix is returned.

        Returns
        -------
        hessian : ndarray
            Hessian, i.e. observed information, or expected information matrix.
        """
        if observed is None:
            if getattr(self, "_optim_hessian", None) == "eim":
                observed = False
            else:
                observed = True
        scale = float_like(scale, "scale", optional=True)
        tmp = getattr(self, "_tmp_like_exog", np.empty_like(self.exog, dtype=float))

        factor = self.hessian_factor(params, scale=scale, observed=observed)
        np.multiply(self.exog.T, factor, out=tmp.T)
        return -tmp.T.dot(self.exog)

    def information(self, params, scale=None):
        """
        Fisher information matrix.
        """
        scale = float_like(scale, "scale", optional=True)
        return self.hessian(params, scale=scale, observed=False)

    def _derivative_exog(
        self,
        params,
        exog=None,
        transform="dydx",
        dummy_idx=None,
        count_idx=None,
        offset=None,
        exposure=None,
    ):
        """
        Derivative of mean, expected endog with respect to the parameters
        """
        if exog is None:
            exog = self.exog
        if (offset is not None) or (exposure is not None):
            raise NotImplementedError("offset and exposure not supported")

        lin_pred = self.predict(
            params, exog, which="linear", offset=offset, exposure=exposure
        )

        k_extra = getattr(self, "k_extra", 0)
        params_exog = params if k_extra == 0 else params[:-k_extra]

        margeff = self.family.link.inverse_deriv(lin_pred)[:, None] * params_exog
        if "ex" in transform:
            margeff *= exog
        if "ey" in transform:
            mean = self.family.link.inverse(lin_pred)
            margeff /= mean[:, None]

        return self._derivative_exog_helper(
            margeff, params, exog, dummy_idx, count_idx, transform
        )

    def _derivative_exog_helper(
        self, margeff, params, exog, dummy_idx, count_idx, transform
    ):
        """
        Helper for _derivative_exog to wrap results appropriately
        """
        from statsmodels.discrete.discrete_margins import (
            _get_count_effects,
            _get_dummy_effects,
        )

        if count_idx is not None:
            margeff = _get_count_effects(
                margeff, exog, count_idx, transform, self, params
            )
        if dummy_idx is not None:
            margeff = _get_dummy_effects(
                margeff, exog, dummy_idx, transform, self, params
            )

        return margeff

    def _derivative_predict(
        self, params, exog=None, transform="dydx", offset=None, exposure=None
    ):
        """
        Derivative of the expected endog with respect to the parameters.

        Parameters
        ----------
        params : ndarray
            parameter at which score is evaluated
        exog : ndarray or None
            Explanatory variables at which derivative are computed.
            If None, then the estimation exog is used.
        offset, exposure : None
            Not yet implemented.

        Returns
        -------
        The value of the derivative of the expected endog with respect
        to the parameter vector.
        """
        # core part is same as derivative_mean_params
        # additionally handles exog and transform
        if exog is None:
            exog = self.exog
        if (
            (offset is not None)
            or (exposure is not None)
            or (getattr(self, "offset", None) is not None)
        ):
            raise NotImplementedError("offset and exposure not supported")

        lin_pred = self.predict(params, exog=exog, which="linear")
        idl = self.family.link.inverse_deriv(lin_pred)
        dmat = exog * idl[:, None]
        if "ey" in transform:
            mean = self.family.link.inverse(lin_pred)
            dmat /= mean[:, None]

        return dmat

    def _deriv_mean_dparams(self, params):
        """
        Derivative of the expected endog with respect to the parameters.

        Parameters
        ----------
        params : ndarray
            parameter at which score is evaluated

        Returns
        -------
        The value of the derivative of the expected endog with respect
        to the parameter vector.
        """
        lin_pred = self.predict(params, which="linear")
        idl = self.family.link.inverse_deriv(lin_pred)
        dmat = self.exog * idl[:, None]
        return dmat

    def _deriv_score_obs_dendog(self, params, scale=None):
        """derivative of score_obs w.r.t. endog

        Parameters
        ----------
        params : ndarray
            parameter at which score is evaluated
        scale : None or float
            If scale is None, then the default scale will be calculated.
            Default scale is defined by `self.scaletype` and set in fit.
            If scale is not None, then it is used as a fixed scale.

        Returns
        -------
        derivative : ndarray_2d
            The derivative of the score_obs with respect to endog. This
            can is given by `score_factor0[:, None] * exog` where
            `score_factor0` is the score_factor without the residual.
        """
        scale = float_like(scale, "scale", optional=True)
        mu = self.predict(params)
        if scale is None:
            scale = self.estimate_scale(mu)

        score_factor = 1 / self.family.link.deriv(mu)
        score_factor /= self.family.variance(mu)
        score_factor *= self.iweights * self.n_trials

        if not scale == 1:
            score_factor /= scale

        return score_factor[:, None] * self.exog

    def score_test(
        self, params_constrained, k_constraints=None, exog_extra=None, observed=True
    ):
        """score test for restrictions or for omitted variables

        The covariance matrix for the score is based on the Hessian, i.e.
        observed information matrix or optionally on the expected information
        matrix..

        Parameters
        ----------
        params_constrained : array_like
            estimated parameter of the restricted model. This can be the
            parameter estimate for the current when testing for omitted
            variables.
        k_constraints : int or None
            Number of constraints that were used in the estimation of params
            restricted relative to the number of exog in the model.
            This must be provided if no exog_extra are given. If exog_extra is
            not None, then k_constraints is assumed to be zero if it is None.
        exog_extra : None or array_like
            Explanatory variables that are jointly tested for inclusion in the
            model, i.e. omitted variables.
        observed : bool
            If True, then the observed Hessian is used in calculating the
            covariance matrix of the score. If false then the expected
            information matrix is used.

        Returns
        -------
        chi2_stat : float
            chisquare statistic for the score test
        p-value : float
            P-value of the score test based on the chisquare distribution.
        df : int
            Degrees of freedom used in the p-value calculation. This is equal
            to the number of constraints.

        Notes
        -----
