"""
Test functions for models.GLM
"""

from statsmodels.compat.scipy import SP_LT_17

import os
import re
import warnings

import numpy as np
from numpy.testing import (
    assert_,
    assert_allclose,
    assert_almost_equal,
    assert_array_less,
    assert_equal,
)
import pandas as pd
from pandas.testing import assert_series_equal
import pytest
from scipy import stats

import statsmodels.api as sm
from statsmodels.datasets import cpunish, longley
from statsmodels.discrete import discrete_model as discrete
from statsmodels.formula._manager import FormulaManager
from statsmodels.genmod.generalized_linear_model import GLM, SET_USE_BIC_LLF
from statsmodels.tools.numdiff import (
    approx_fprime,
    approx_fprime_cs,
    approx_hess,
    approx_hess_cs,
)
from statsmodels.tools.sm_exceptions import (
    DomainWarning,
    PerfectSeparationWarning,
    ValueWarning,
)
from statsmodels.tools.tools import add_constant

# Test Precisions
DECIMAL_4 = 4
DECIMAL_3 = 3
DECIMAL_2 = 2
DECIMAL_1 = 1
DECIMAL_0 = 0

pdf_output = False

if pdf_output:
    from matplotlib.backends.backend_pdf import PdfPages

    pdf = PdfPages("test_glm.pdf")
else:
    pdf = None


def close_or_save(pdf, fig):
    if pdf_output:
        pdf.savefig(fig)


def teardown_module():
    if pdf_output:
        pdf.close()


@pytest.fixture(scope="module")
def iris():
    cur_dir = os.path.dirname(os.path.abspath(__file__))
    return np.genfromtxt(
        os.path.join(cur_dir, "results", "iris.csv"),
        delimiter=",",
        skip_header=1,
    )


class CheckModelResultsMixin:
    """
    res2 should be either the results from RModelWrap
    or the results as defined in model_results_data
    """

    decimal_params = DECIMAL_4

    def test_params(self):
        assert_almost_equal(self.res1.params, self.res2.params, self.decimal_params)

    decimal_bse = DECIMAL_4

    def test_standard_errors(self):
        assert_allclose(
            self.res1.bse,
            self.res2.bse,
            atol=10 ** (-self.decimal_bse),
            rtol=1e-5,
        )

    decimal_resids = DECIMAL_4

    def test_residuals(self):
        # fix incorrect numbers in resid_working results
        # residuals for Poisson are also tested in test_glm_weights.py
        import copy

        # new numpy would have copy method
        resid2 = copy.copy(self.res2.resids)
        resid2[:, 2] *= self.res1.family.link.deriv(self.res1.mu) ** 2

        atol = float(10 ** (-self.decimal_resids))
        resid_a = self.res1.resid_anscombe_unscaled
        resids = np.column_stack(
            (
                self.res1.resid_pearson,
                self.res1.resid_deviance,
                self.res1.resid_working,
                resid_a,
                self.res1.resid_response,
            )
        )
        assert_allclose(resids, resid2, rtol=1e-6, atol=atol)

    decimal_aic_R = DECIMAL_4

    def test_aic_R(self):
        # R includes the estimation of the scale as a lost dof
        # Does not with Gamma though
        if self.res1.scale != 1:
            dof = 2
        else:
            dof = 0
        if isinstance(self.res1.model.family, (sm.families.NegativeBinomial)):
            llf = self.res1.model.family.loglike(
                self.res1.model.endog,
                self.res1.mu,
                self.res1.model.var_weights,
                self.res1.model.freq_weights,
                scale=1,
            )
            aic = -2 * llf + 2 * (self.res1.df_model + 1)
        else:
            aic = self.res1.aic
        assert_almost_equal(aic + dof, self.res2.aic_R, self.decimal_aic_R)

    decimal_aic_Stata = DECIMAL_4

    def test_aic_Stata(self):
        # Stata uses the below llf for aic definition for these families
        if isinstance(
            self.res1.model.family,
            (
                sm.families.Gamma,
                sm.families.InverseGaussian,
                sm.families.NegativeBinomial,
            ),
        ):
            llf = self.res1.model.family.loglike(
                self.res1.model.endog,
                self.res1.mu,
                self.res1.model.var_weights,
                self.res1.model.freq_weights,
                scale=1,
            )
            aic = (-2 * llf + 2 * (self.res1.df_model + 1)) / self.res1.nobs
        else:
            aic = self.res1.aic / self.res1.nobs
        assert_almost_equal(aic, self.res2.aic_Stata, self.decimal_aic_Stata)

    decimal_deviance = DECIMAL_4

    def test_deviance(self):
        assert_almost_equal(
            self.res1.deviance, self.res2.deviance, self.decimal_deviance
        )

    decimal_scale = DECIMAL_4

    def test_scale(self):
        assert_almost_equal(self.res1.scale, self.res2.scale, self.decimal_scale)

    decimal_loglike = DECIMAL_4

    def test_loglike(self):
        # Stata uses the below llf for these families
        # We differ with R for them
        if isinstance(
            self.res1.model.family,
            (
                sm.families.Gamma,
                sm.families.InverseGaussian,
                sm.families.NegativeBinomial,
            ),
        ):
            llf = self.res1.model.family.loglike(
                self.res1.model.endog,
                self.res1.mu,
                self.res1.model.var_weights,
                self.res1.model.freq_weights,
                scale=1,
            )
        else:
            llf = self.res1.llf
        assert_almost_equal(llf, self.res2.llf, self.decimal_loglike)

    decimal_null_deviance = DECIMAL_4

    def test_null_deviance(self):
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", DomainWarning)

            assert_almost_equal(
                self.res1.null_deviance,
                self.res2.null_deviance,
                self.decimal_null_deviance,
            )

    decimal_bic = DECIMAL_4

    def test_bic(self):
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            assert_almost_equal(self.res1.bic, self.res2.bic_Stata, self.decimal_bic)

    def test_degrees(self):
        assert_equal(self.res1.model.df_resid, self.res2.df_resid)

    decimal_fittedvalues = DECIMAL_4

    def test_fittedvalues(self):
        assert_almost_equal(
            self.res1.fittedvalues,
            self.res2.fittedvalues,
            self.decimal_fittedvalues,
        )

    def test_tpvalues(self):
        # test comparing tvalues and pvalues with normal implementation
        # make sure they use normal distribution (inherited in results class)
        params = self.res1.params
        tvalues = params / self.res1.bse
        pvalues = stats.norm.sf(np.abs(tvalues)) * 2
        half_width = stats.norm.isf(0.025) * self.res1.bse
        conf_int = np.column_stack((params - half_width, params + half_width))
        if isinstance(tvalues, pd.Series):
            assert_series_equal(self.res1.tvalues, tvalues)
        else:
            assert_almost_equal(self.res1.tvalues, tvalues)
        assert_almost_equal(self.res1.pvalues, pvalues)
        assert_almost_equal(self.res1.conf_int(), conf_int)

    def test_pearson_chi2(self):
        if hasattr(self.res2, "pearson_chi2"):
            assert_allclose(
                self.res1.pearson_chi2,
                self.res2.pearson_chi2,
                atol=1e-6,
                rtol=1e-6,
            )

    def test_prsquared(self):
        if hasattr(self.res2, "prsquared"):
            assert_allclose(
                self.res1.pseudo_rsquared(kind="mcf"),
                self.res2.prsquared,
                rtol=0.05,
            )

        if hasattr(self.res2, "prsquared_cox_snell"):
            assert_allclose(
                float(self.res1.pseudo_rsquared(kind="cs")),
                self.res2.prsquared_cox_snell,
                rtol=0.05,
            )

    @pytest.mark.smoke
    def test_summary(self):
        self.res1.summary()

    @pytest.mark.smoke
    def test_summary2(self):
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", DomainWarning)
            self.res1.summary2()

    def test_get_distribution(self):
        res1 = self.res1
        if not hasattr(res1.model.family, "get_distribution"):
            # only Tweedie has not get_distribution
            pytest.skip("get_distribution not available")

        if isinstance(res1.model.family, sm.families.NegativeBinomial):
            res_scale = 1  # QMLE scale can differ from 1
        else:
            res_scale = res1.scale

        distr = res1.model.family.get_distribution(res1.fittedvalues, res_scale)
        var_endog = res1.model.family.variance(res1.fittedvalues) * res_scale
        m, v = distr.stats()
        assert_allclose(res1.fittedvalues, m, rtol=1e-13)
        assert_allclose(var_endog, v, rtol=1e-13)
        # check model method
        distr2 = res1.model.get_distribution(res1.params, res_scale)
        for k in distr2.kwds:
            assert_allclose(distr.kwds[k], distr2.kwds[k], rtol=1e-13)

        # compare var with predict
        var_ = res1.predict(which="var_unscaled")
        assert_allclose(var_ * res_scale, var_endog, rtol=1e-13)

        # check get_distribution of results instance
        if getattr(self, "has_edispersion", False):
            with pytest.warns(UserWarning, match="using scale=1"):
                distr3 = res1.get_distribution()
        else:
            distr3 = res1.get_distribution()
        for k in distr2.kwds:
            assert_allclose(distr3.kwds[k], distr2.kwds[k], rtol=1e-13)


class CheckComparisonMixin:

    def test_compare_discrete(self):
        res1 = self.res1
        resd = self.resd

        assert_allclose(res1.llf, resd.llf, rtol=1e-10)
        score_obs1 = res1.model.score_obs(res1.params * 0.98)
        score_obsd = resd.model.score_obs(resd.params * 0.98)
        assert_allclose(score_obs1, score_obsd, rtol=1e-10)

        # score
        score1 = res1.model.score(res1.params * 0.98)
        assert_allclose(score1, score_obs1.sum(0), atol=1e-20)
        score0 = res1.model.score(res1.params)
        assert_allclose(score0, np.zeros(score_obs1.shape[1]), atol=5e-7)

        hessian1 = res1.model.hessian(res1.params * 0.98, observed=False)
        hessiand = resd.model.hessian(resd.params * 0.98)
        assert_allclose(hessian1, hessiand, rtol=1e-10)

        hessian1 = res1.model.hessian(res1.params * 0.98, observed=True)
        hessiand = resd.model.hessian(resd.params * 0.98)
        assert_allclose(hessian1, hessiand, rtol=1e-9)

    def test_score_test(self):
        res1 = self.res1
        # fake example, should be zero, k_constraint should be 0
        st, pv, df = res1.model.score_test(res1.params, k_constraints=1)
        assert_allclose(st, 0, atol=1e-20)
        assert_allclose(pv, 1, atol=1e-10)
        assert_equal(df, 1)

        st, pv, df = res1.model.score_test(res1.params, k_constraints=0)
        assert_allclose(st, 0, atol=1e-20)
        assert_(np.isnan(pv), msg=repr(pv))
        assert_equal(df, 0)

        # TODO: no verified numbers largely SMOKE test
        exog_extra = res1.model.exog[:, 1] ** 2
        st, pv, df = res1.model.score_test(res1.params, exog_extra=exog_extra)
        assert_array_less(0.1, st)
        assert_array_less(0.1, pv)
        assert_equal(df, 1)

    def test_get_prediction(self):
        pred1 = self.res1.get_prediction()  # GLM
        predd = self.resd.get_prediction()  # discrete class
        assert_allclose(predd.predicted, pred1.predicted_mean, rtol=1e-11)
        assert_allclose(predd.se, pred1.se_mean, rtol=1e-6)
        assert_allclose(
            predd.summary_frame().values,
            pred1.summary_frame().values,
            rtol=1e-6,
        )

        pred1 = self.res1.get_prediction(which="mean")  # GLM
        predd = self.resd.get_prediction()  # discrete class
        assert_allclose(predd.predicted, pred1.predicted, rtol=1e-11)
        assert_allclose(predd.se, pred1.se, rtol=1e-6)
        assert_allclose(
            predd.summary_frame().values,
            pred1.summary_frame().values,
            rtol=1e-6,
        )


class TestGlmGaussian(CheckModelResultsMixin):
    @classmethod
    def setup_class(cls):
        """
        Test Gaussian family with canonical identity link
        """
        # Test Precisions
        cls.decimal_resids = DECIMAL_3
        cls.decimal_params = DECIMAL_2
        cls.decimal_bic = DECIMAL_0
        cls.decimal_bse = DECIMAL_3

        from statsmodels.datasets.longley import load

        cls.data = load()
        cls.data.endog = np.require(cls.data.endog, requirements="W")
        cls.data.exog = np.require(cls.data.exog, requirements="W")
        cls.data.exog = add_constant(cls.data.exog, prepend=False)
        cls.res1 = GLM(
            cls.data.endog, cls.data.exog, family=sm.families.Gaussian()
        ).fit()
        from .results.results_glm import Longley

        cls.res2 = Longley()

    def test_compare_OLS(self):
        res1 = self.res1
        # OLS does not define score_obs
        from statsmodels.regression.linear_model import OLS

        resd = OLS(self.data.endog, self.data.exog).fit(use_t=False)
        self.resd = resd  # attach to access from the outside

        assert_allclose(res1.llf, resd.llf, rtol=1e-10)
        score_obs1 = res1.model.score_obs(res1.params, scale=None)
        score_obsd = resd.resid[:, None] / resd.scale * resd.model.exog
        # low precision because of badly scaled exog
        assert_allclose(score_obs1, score_obsd, rtol=1e-8)

        score_obs1 = res1.model.score_obs(res1.params, scale=1)
        score_obsd = resd.resid[:, None] * resd.model.exog
        assert_allclose(score_obs1, score_obsd, rtol=1e-8)

        hess_obs1 = res1.model.hessian(res1.params, scale=None)
        hess_obsd = -1.0 / resd.scale * resd.model.exog.T.dot(resd.model.exog)
        # low precision because of badly scaled exog
        assert_allclose(hess_obs1, hess_obsd, rtol=1e-8)

        pred1 = res1.get_prediction()  # GLM
        predd = resd.get_prediction()  # discrete class
        assert_allclose(predd.predicted, pred1.predicted_mean, rtol=1e-11)
        assert_allclose(predd.se, pred1.se_mean, rtol=1e-6)
        assert_allclose(
            predd.summary_frame().values[:, :4],
            pred1.summary_frame().values,
            rtol=1e-6,
        )

        pred1 = self.res1.get_prediction(which="mean")  # GLM
        predd = self.resd.get_prediction()  # discrete class
        assert_allclose(predd.predicted, pred1.predicted, rtol=1e-11)
        assert_allclose(predd.se, pred1.se, rtol=1e-6)
        assert_allclose(
            predd.summary_frame().values[:, :4],
            pred1.summary_frame().values,
            rtol=1e-6,
        )


# FIXME: enable or delete
#    def setup_method(self):
#        if skipR:
#            raise SkipTest, "Rpy not installed."
#        Gauss = r.gaussian
#        self.res2 = RModel(self.data.endog, self.data.exog, r.glm, family=Gauss)
#        self.res2.resids = np.array(self.res2.resid)[:,None]*np.ones((1,5))
#        self.res2.null_deviance = 185008826 # taken from R. Rpy bug?


class TestGlmGaussianGradient(TestGlmGaussian):
    @classmethod
    def setup_class(cls):
        """
        Test Gaussian family with canonical identity link
        """
        # Test Precisions
        cls.decimal_resids = DECIMAL_3
        cls.decimal_params = DECIMAL_2
        cls.decimal_bic = DECIMAL_0
        cls.decimal_bse = DECIMAL_2

        from statsmodels.datasets.longley import load

        cls.data = load()
        cls.data.endog = np.require(cls.data.endog, requirements="W")
        cls.data.exog = np.require(cls.data.exog, requirements="W")
        cls.data.exog = add_constant(cls.data.exog, prepend=False)
        cls.res1 = GLM(
            cls.data.endog, cls.data.exog, family=sm.families.Gaussian()
        ).fit(method="newton")
        from .results.results_glm import Longley

        cls.res2 = Longley()


class TestGaussianLog(CheckModelResultsMixin):
    @classmethod
    def setup_class(cls):
        # Test Precision
        cls.decimal_aic_R = DECIMAL_0
        cls.decimal_aic_Stata = DECIMAL_2
        cls.decimal_loglike = DECIMAL_0
        cls.decimal_null_deviance = DECIMAL_1

        nobs = 100
        x = np.arange(nobs)
        np.random.seed(54321)
        #        y = 1.0 - .02*x - .001*x**2 + 0.001 * np.random.randn(nobs)
        cls.X = np.c_[np.ones((nobs, 1)), x, x**2]
        cls.lny = np.exp(
            -(-1.0 + 0.02 * x + 0.0001 * x**2)
        ) + 0.001 * np.random.randn(nobs)

        GaussLog_Model = GLM(
            cls.lny,
            cls.X,
            family=sm.families.Gaussian(sm.families.links.Log()),
        )
        cls.res1 = GaussLog_Model.fit()
        from .results.results_glm import GaussianLog

        cls.res2 = GaussianLog()


# FIXME: enable or delete
#    def setup(cls):
#        if skipR:
#            raise SkipTest, "Rpy not installed"
#        GaussLogLink = r.gaussian(link = "log")
#        GaussLog_Res_R = RModel(cls.lny, cls.X, r.glm, family=GaussLogLink)
#        cls.res2 = GaussLog_Res_R


class TestGaussianInverse(CheckModelResultsMixin):
    @classmethod
    def setup_class(cls):
        # Test Precisions
        cls.decimal_bic = DECIMAL_1
        cls.decimal_aic_R = DECIMAL_1
        cls.decimal_aic_Stata = DECIMAL_3
        cls.decimal_loglike = DECIMAL_1
        cls.decimal_resids = DECIMAL_3

        nobs = 100
        x = np.arange(nobs)
        np.random.seed(54321)
        cls.y = 1.0 + 2.0 * x + x**2 + 0.1 * np.random.randn(nobs)
        cls.X = np.c_[np.ones((nobs, 1)), x, x**2]
        cls.y_inv = (1.0 + 0.02 * x + 0.001 * x**2) ** -1 + 0.001 * np.random.randn(
            nobs
        )
        InverseLink_Model = GLM(
            cls.y_inv,
            cls.X,
            family=sm.families.Gaussian(sm.families.links.InversePower()),
        )
        InverseLink_Res = InverseLink_Model.fit()
        cls.res1 = InverseLink_Res
        from .results.results_glm import GaussianInverse

        cls.res2 = GaussianInverse()


# FIXME: enable or delete
#    def setup(cls):
#        if skipR:
#            raise SkipTest, "Rpy not installed."
#        InverseLink = r.gaussian(link = "inverse")
#        InverseLink_Res_R = RModel(cls.y_inv, cls.X, r.glm, family=InverseLink)
#        cls.res2 = InverseLink_Res_R


class TestGlmBinomial(CheckModelResultsMixin):
    @classmethod
    def setup_class(cls):
        """
        Test Binomial family with canonical logit link using star98 dataset.
        """
        cls.decimal_resids = DECIMAL_1
        cls.decimal_bic = DECIMAL_2

        from statsmodels.datasets.star98 import load

        from .results.results_glm import Star98

        data = load()
        data.endog = np.require(data.endog, requirements="W")
        data.exog = np.require(data.exog, requirements="W")
        data.exog = add_constant(data.exog, prepend=False)
        cls.res1 = GLM(data.endog, data.exog, family=sm.families.Binomial()).fit()
        # NOTE: if you want to replicate with RModel
        # res2 = RModel(data.endog[:,0]/trials, data.exog, r.glm,
        #        family=r.binomial, weights=trials)

        cls.res2 = Star98()

    def test_endog_dtype(self):
        from statsmodels.datasets.star98 import load

        data = load()
        data.exog = add_constant(data.exog, prepend=False)
        endog = data.endog.astype(int)
        res2 = GLM(endog, data.exog, family=sm.families.Binomial()).fit()
        assert_allclose(res2.params, self.res1.params)
        endog = data.endog.astype(np.double)
        res3 = GLM(endog, data.exog, family=sm.families.Binomial()).fit()
        assert_allclose(res3.params, self.res1.params)

    def test_invalid_endog(self, reset_randomstate):
        # GH2733 inspired check
        endog = np.random.randint(0, 100, size=(1000, 3))
        exog = np.random.standard_normal((1000, 2))
        with pytest.raises(ValueError, match="endog has more than 2 columns"):
            GLM(endog, exog, family=sm.families.Binomial())

    def test_invalid_endog_formula(self, reset_randomstate):
        # GH2733
        n = 200
        exog = np.random.normal(size=(n, 2))
        endog = np.random.randint(0, 3, size=n).astype(str)
        # formula interface
        data = pd.DataFrame({"y": endog, "x1": exog[:, 0], "x2": exog[:, 1]})
        with pytest.raises(ValueError, match="array with multiple columns"):
            sm.GLM.from_formula("y ~ x1 + x2", data, family=sm.families.Binomial())

    def test_get_distribution_binom_count(self):
        # test for binomial counts with n_trials > 1
        res1 = self.res1
        res_scale = 1  # QMLE scale can differ from 1

        mu_prob = res1.fittedvalues
        n = res1.model.n_trials
        distr = res1.model.family.get_distribution(mu_prob, res_scale, n_trials=n)
        var_endog = res1.model.family.variance(mu_prob) * res_scale
        m, v = distr.stats()
        assert_allclose(mu_prob * n, m, rtol=1e-13)
        assert_allclose(var_endog * n, v, rtol=1e-13)

        # check model method
        distr2 = res1.model.get_distribution(res1.params, res_scale, n_trials=n)
        for k in distr2.kwds:
            assert_allclose(distr.kwds[k], distr2.kwds[k], rtol=1e-13)


# FIXME: enable/xfail/skip or delete
# TODO:
# Non-Canonical Links for the Binomial family require the algorithm to be
# slightly changed
# class TestGlmBinomialLog(CheckModelResultsMixin):
#    pass

# class TestGlmBinomialLogit(CheckModelResultsMixin):
#    pass

# class TestGlmBinomialProbit(CheckModelResultsMixin):
#    pass

# class TestGlmBinomialCloglog(CheckModelResultsMixin):
#    pass

# class TestGlmBinomialPower(CheckModelResultsMixin):
#    pass

# class TestGlmBinomialLoglog(CheckModelResultsMixin):
#    pass

# class TestGlmBinomialLogc(CheckModelResultsMixin):
# TODO: need include logc link
#    pass


class TestGlmBernoulli(CheckModelResultsMixin, CheckComparisonMixin):
    @classmethod
    def setup_class(cls):
        from .results.results_glm import Lbw

        cls.res2 = Lbw()
        cls.res1 = GLM(
            cls.res2.endog, cls.res2.exog, family=sm.families.Binomial()
        ).fit()

        modd = discrete.Logit(cls.res2.endog, cls.res2.exog)
        cls.resd = modd.fit(start_params=cls.res1.params * 0.9, disp=False)

    def test_score_r(self):
        res1 = self.res1
        res2 = self.res2
        st, pv, df = res1.model.score_test(
            res1.params, exog_extra=res1.model.exog[:, 1] ** 2
        )
        st_res = 0.2837680293459376  # (-0.5326988167303712)**2
        assert_allclose(st, st_res, rtol=1e-4)

        st, pv, df = res1.model.score_test(
            res1.params, exog_extra=res1.model.exog[:, 0] ** 2
        )
        st_res = 0.6713492821514992  # (-0.8193590679009413)**2
        assert_allclose(st, st_res, rtol=1e-4)

        select = list(range(9))
        select.pop(7)

        res1b = GLM(
            res2.endog,
            res2.exog.iloc[:, select],
            family=sm.families.Binomial(),
        ).fit()
        tres = res1b.model.score_test(res1b.params, exog_extra=res1.model.exog[:, -2])
        tres = np.asarray(tres[:2]).ravel()
        tres_r = (2.7864148487452, 0.0950667)
        assert_allclose(tres, tres_r, rtol=1e-4)

        cmd_r = """\
        data = read.csv("...statsmodels\\statsmodels\\genmod\\tests\\results\\stata_lbw_glm.csv")

        data["race_black"] = data["race"] == "black"
        data["race_other"] = data["race"] == "other"
        mod = glm(low ~ age + lwt + race_black + race_other + smoke + ptl + ht + ui, family=binomial, data=data)
        options(digits=16)
        anova(mod, test="Rao")

        library(statmod)
        s = glm.scoretest(mod, data["age"]**2)
        s**2
        s = glm.scoretest(mod, data["lwt"]**2)
        s**2
        """
        assert isinstance(cmd_r, str)


# class TestGlmBernoulliIdentity(CheckModelResultsMixin):
#    pass

# class TestGlmBernoulliLog(CheckModelResultsMixin):
#    pass

# class TestGlmBernoulliProbit(CheckModelResultsMixin):
#    pass

# class TestGlmBernoulliCloglog(CheckModelResultsMixin):
#    pass

# class TestGlmBernoulliPower(CheckModelResultsMixin):
#    pass

# class TestGlmBernoulliLoglog(CheckModelResultsMixin):
#    pass

# class test_glm_bernoulli_logc(CheckModelResultsMixin):
#    pass


class TestGlmGamma(CheckModelResultsMixin):

    @classmethod
    def setup_class(cls):
        """
        Tests Gamma family with canonical inverse link (power -1)
        """
        # Test Precisions
        cls.decimal_aic_R = -1  # TODO: off by about 1, we are right with Stata
        cls.decimal_resids = DECIMAL_2

        from statsmodels.datasets.scotland import load

        from .results.results_glm import Scotvote

        data = load()
        data.exog = add_constant(data.exog, prepend=False)
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            res1 = GLM(data.endog, data.exog, family=sm.families.Gamma()).fit()
        cls.res1 = res1
        #        res2 = RModel(data.endog, data.exog, r.glm, family=r.Gamma)
        res2 = Scotvote()
        res2.aic_R += 2  # R does not count degree of freedom for scale with gamma
        cls.res2 = res2


class TestGlmGammaLog(CheckModelResultsMixin):
    @classmethod
    def setup_class(cls):
        # Test Precisions
        cls.decimal_resids = DECIMAL_3
        cls.decimal_aic_R = DECIMAL_0
        cls.decimal_fittedvalues = DECIMAL_3

        from .results.results_glm import CancerLog

        res2 = CancerLog()
        cls.res1 = GLM(
            res2.endog,
            res2.exog,
            family=sm.families.Gamma(link=sm.families.links.Log()),
        ).fit()
        cls.res2 = res2


# FIXME: enable or delete
#    def setup(cls):
#        if skipR:
#            raise SkipTest, "Rpy not installed."
#        cls.res2 = RModel(cls.data.endog, cls.data.exog, r.glm,
#            family=r.Gamma(link="log"))
#        cls.res2.null_deviance = 27.92207137420696 # From R (bug in rpy)
#        cls.res2.bic = -154.1582089453923 # from Stata


class TestGlmGammaIdentity(CheckModelResultsMixin):
    @classmethod
    def setup_class(cls):
        # Test Precisions
        cls.decimal_resids = -100  # TODO Very off from Stata?
        cls.decimal_params = DECIMAL_2
        cls.decimal_aic_R = DECIMAL_0
        cls.decimal_loglike = DECIMAL_1

        from .results.results_glm import CancerIdentity

        res2 = CancerIdentity()
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            fam = sm.families.Gamma(link=sm.families.links.Identity())
            cls.res1 = GLM(res2.endog, res2.exog, family=fam).fit()
        cls.res2 = res2


# FIXME: enable or delete
#    def setup(cls):
#        if skipR:
#            raise SkipTest, "Rpy not installed."
#        cls.res2 = RModel(cls.data.endog, cls.data.exog, r.glm,
#            family=r.Gamma(link="identity"))
#        cls.res2.null_deviance = 27.92207137420696 # from R, Rpy bug


class TestGlmPoisson(CheckModelResultsMixin, CheckComparisonMixin):
    @classmethod
    def setup_class(cls):
        """
        Tests Poisson family with canonical log link.

        Test results were obtained by R.
        """
        from .results.results_glm import Cpunish

        cls.data = cpunish.load()
        cls.data.endog = np.require(cls.data.endog, requirements="W")
        cls.data.exog = np.require(cls.data.exog, requirements="W")
        cls.data.exog[:, 3] = np.log(cls.data.exog[:, 3])
        cls.data.exog = add_constant(cls.data.exog, prepend=False)
        cls.res1 = GLM(
            cls.data.endog, cls.data.exog, family=sm.families.Poisson()
        ).fit()
        cls.res2 = Cpunish()
        # compare with discrete, start close to save time
        modd = discrete.Poisson(cls.data.endog, cls.data.exog)
        cls.resd = modd.fit(start_params=cls.res1.params * 0.9, disp=False)


# class TestGlmPoissonIdentity(CheckModelResultsMixin):
#    pass

# class TestGlmPoissonPower(CheckModelResultsMixin):
#    pass


class TestGlmInvgauss(CheckModelResultsMixin):
    @classmethod
    def setup_class(cls):
        """
        Tests the Inverse Gaussian family in GLM.

        Notes
        -----
        Used the rndivgx.ado file provided by Hardin and Hilbe to
        generate the data.  Results are read from model_results, which
        were obtained by running R_ig.s
        """
        # Test Precisions
        cls.decimal_aic_R = DECIMAL_0
        cls.decimal_loglike = DECIMAL_0

        from .results.results_glm import InvGauss

        res2 = InvGauss()
        res1 = GLM(res2.endog, res2.exog, family=sm.families.InverseGaussian()).fit()
        cls.res1 = res1
        cls.res2 = res2

    def test_get_distribution(self):
        res1 = self.res1
        distr = res1.model.family.get_distribution(res1.fittedvalues, res1.scale)
        var_endog = res1.model.family.variance(res1.fittedvalues) * res1.scale
        m, v = distr.stats()
        assert_allclose(res1.fittedvalues, m, rtol=1e-13)
        assert_allclose(var_endog, v, rtol=1e-13)


class TestGlmInvgaussLog(CheckModelResultsMixin):
    @classmethod
    def setup_class(cls):
        # Test Precisions
        cls.decimal_aic_R = -10  # Big difference vs R.
        cls.decimal_resids = DECIMAL_3

        from .results.results_glm import InvGaussLog

        res2 = InvGaussLog()
        cls.res1 = GLM(
            res2.endog,
            res2.exog,
            family=sm.families.InverseGaussian(link=sm.families.links.Log()),
        ).fit()
        cls.res2 = res2


# FIXME: enable or delete
#    def setup(cls):
#        if skipR:
#            raise SkipTest, "Rpy not installed."
#        cls.res2 = RModel(cls.data.endog, cls.data.exog, r.glm,
#            family=r.inverse_gaussian(link="log"))
#        cls.res2.null_deviance = 335.1539777981053 # from R, Rpy bug
#        cls.res2.llf = -12162.72308 # from Stata, R's has big rounding diff


class TestGlmInvgaussIdentity(CheckModelResultsMixin):
    @classmethod
    def setup_class(cls):
        # Test Precisions
        cls.decimal_aic_R = -10  # TODO: Big difference vs R
        cls.decimal_fittedvalues = DECIMAL_3
        cls.decimal_params = DECIMAL_3

        from .results.results_glm import Medpar1

        data = Medpar1()
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            cls.res1 = GLM(
                data.endog,
                data.exog,
                family=sm.families.InverseGaussian(link=sm.families.links.Identity()),
            ).fit()
        from .results.results_glm import InvGaussIdentity

        cls.res2 = InvGaussIdentity()


# FIXME: enable or delete
#    def setup(cls):
#        if skipR:
#            raise SkipTest, "Rpy not installed."
#        cls.res2 = RModel(cls.data.endog, cls.data.exog, r.glm,
#            family=r.inverse_gaussian(link="identity"))
#        cls.res2.null_deviance = 335.1539777981053 # from R, Rpy bug
#        cls.res2.llf = -12163.25545    # from Stata, big diff with R


class TestGlmNegbinomial(CheckModelResultsMixin):
    @classmethod
    def setup_class(cls):
        """
        Test Negative Binomial family with log link
        """
        # Test Precision
        cls.decimal_resid = DECIMAL_1
        cls.decimal_params = DECIMAL_3
        cls.decimal_resids = -1  # 1 % mismatch at 0
        cls.decimal_fittedvalues = DECIMAL_1

        from statsmodels.datasets.committee import load

        cls.data = load()
        cls.data.endog = np.require(cls.data.endog, requirements="W")
        cls.data.exog = np.require(cls.data.exog, requirements="W")
        cls.data.exog[:, 2] = np.log(cls.data.exog[:, 2])
        interaction = cls.data.exog[:, 2] * cls.data.exog[:, 1]
        cls.data.exog = np.column_stack((cls.data.exog, interaction))
        cls.data.exog = add_constant(cls.data.exog, prepend=False)
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", category=DomainWarning)
            with pytest.warns(ValueWarning, match="Negative binomial dispersio"):
                fam = sm.families.NegativeBinomial()

        cls.res1 = GLM(cls.data.endog, cls.data.exog, family=fam).fit(scale="x2")
        from .results.results_glm import Committee

        res2 = Committee()
        res2.aic_R += 2  # They do not count a degree of freedom for the scale
        cls.res2 = res2
        cls.has_edispersion = True


# FIXME: enable or delete
#    def setup_method(self):
#        if skipR:
#            raise SkipTest, "Rpy not installed"
