Python statsmodels.api.GLM Examples

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

#    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 

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.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 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 doesn't count degree of freedom for scale with gamma
        cls.res2 = res2 

def test_formula_missing_exposure():
    # see 2083
    import statsmodels.formula.api as smf
    import pandas as pd

    d = {'Foo': [1, 2, 10, 149], 'Bar': [1, 2, 3, np.nan],
         'constant': [1] * 4, 'exposure' : np.random.uniform(size=4),
         'x': [1, 3, 2, 1.5]}
    df = pd.DataFrame(d)

    family = sm.families.Gaussian(link=sm.families.links.log())

    mod = smf.glm("Foo ~ Bar", data=df, exposure=df.exposure,
                  family=family)
    assert_(type(mod.exposure) is np.ndarray, msg='Exposure is not ndarray')

    exposure = pd.Series(np.random.uniform(size=5))
    df.loc[3, 'Bar'] = 4   # nan not relevant for Valueerror for shape mismatch
    assert_raises(ValueError, smf.glm, "Foo ~ Bar", data=df,
                  exposure=exposure, family=family)
    assert_raises(ValueError, GLM, df.Foo, df[['constant', 'Bar']],
                  exposure=exposure, family=family) 

def setup_class(cls):

        '''
        Tests Poisson family with canonical log link.
        '''
        super(TestWtdGlmPoissonHC0, cls).setup_class()

        start_params = np.array([1.82794424e-04, -4.76785037e-02,
                                 -9.48249717e-02, -2.92293226e-04,
                                 2.63728909e+00, -2.05934384e+01])

        fit_kwds = dict(cov_type='HC0')
        cls.res1 = GLM(cls.endog, cls.exog,
                        freq_weights=cls.weight,
                        family=sm.families.Poisson()).fit(**fit_kwds)
        fit_kwds = dict(cov_type='HC0', start_params=start_params)
        cls.res2 = GLM(cls.endog_big, cls.exog_big,
                        family=sm.families.Poisson()).fit(**fit_kwds) 

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.exog = add_constant(cls.data.exog, prepend=False)
        params = sm.OLS(cls.data.endog, cls.data.exog).fit().params
        cls.res1 = GLM(cls.data.endog, cls.data.exog,
                        family=sm.families.Gaussian()).fit(start_params=params)
        from .results.results_glm import Longley
        cls.res2 = Longley() 

def setup_class(cls):

        '''
        Tests Poisson family with canonical log link.
        '''
        super(TestWtdGlmPoissonClu, cls).setup_class()

        start_params = np.array([1.82794424e-04, -4.76785037e-02,
                                 -9.48249717e-02, -2.92293226e-04,
                                 2.63728909e+00, -2.05934384e+01])

        gid = np.arange(1, len(cls.endog) + 1) // 2
        fit_kwds = dict(cov_type='cluster', cov_kwds={'groups': gid, 'use_correction':False})

        import warnings
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            cls.res1 = GLM(cls.endog, cls.exog,
                            freq_weights=cls.weight,
                            family=sm.families.Poisson()).fit(**fit_kwds)
            gidr = np.repeat(gid, cls.weight)
            fit_kwds = dict(cov_type='cluster', cov_kwds={'groups': gidr, 'use_correction':False})
            cls.res2 = GLM(cls.endog_big, cls.exog_big,
                            family=sm.families.Poisson()).fit(start_params=start_params,
                                                              **fit_kwds) 

def test_score_test_OLS():
    # nicer example than Longley
    from statsmodels.regression.linear_model import OLS
    np.random.seed(5)
    nobs = 100
    sige = 0.5
    x = np.random.uniform(0, 1, size=(nobs, 5))
    x[:, 0] = 1
    beta = 1. / np.arange(1., x.shape[1] + 1)
    y = x.dot(beta) + sige * np.random.randn(nobs)

    res_ols = OLS(y, x).fit()
    res_olsc = OLS(y, x[:, :-2]).fit()
    co = res_ols.compare_lm_test(res_olsc, demean=False)

    res_glm = GLM(y, x[:, :-2], family=sm.families.Gaussian()).fit()
    co2 = res_glm.model.score_test(res_glm.params, exog_extra=x[:, -2:])
    # difference in df_resid versus nobs in scale see #1786
    assert_allclose(co[0] * 97 / 100., co2[0], rtol=1e-13) 

def setup_class(cls):
        fweights = [1, 1, 1, 2, 2, 2, 3, 3, 3, 1, 1, 1, 2, 2, 2, 3, 3]
        # faking aweights by using normalized freq_weights
        fweights = np.array(fweights)
        wsum = fweights.sum()
        nobs = len(cpunish_data.endog)
        aweights = fweights / wsum * nobs

        cls.res1 = GLM(cpunish_data.endog, cpunish_data.exog,
                    family=sm.families.Poisson(), var_weights=aweights).fit()
        # compare with discrete, start close to save time
        modd = discrete.Poisson(cpunish_data.endog, cpunish_data.exog)

        # Need to copy to avoid inplace adjustment
        from copy import copy
        cls.res2 = copy(res_stata.results_poisson_aweight_nonrobust)
        cls.res2.resids = cls.res2.resids.copy()

        # Need to adjust resids for pearson and deviance to add weights
        cls.res2.resids[:, 3:5] *= np.sqrt(aweights[:, np.newaxis])


# prob_weights fail with HC, not properly implemented yet 

def test_framing_example_formula():

    cur_dir = os.path.dirname(os.path.abspath(__file__))
    data = pd.read_csv(os.path.join(cur_dir, 'results', "framing.csv"))

    probit = sm.families.links.probit
    outcome_model = sm.GLM.from_formula("cong_mesg ~ emo + treat + age + educ + gender + income",
                                        data, family=sm.families.Binomial(link=probit()))

    mediator_model = sm.OLS.from_formula("emo ~ treat + age + educ + gender + income", data)

    med = Mediation(outcome_model, mediator_model, "treat", "emo",
                    outcome_fit_kwargs={'atol': 1e-11})

    np.random.seed(4231)
    med_rslt = med.fit(method='boot', n_rep=100)
    diff = np.asarray(med_rslt.summary() - framing_boot_4231)
    assert_allclose(diff, 0, atol=1e-6)

    np.random.seed(4231)
    med_rslt = med.fit(method='parametric', n_rep=100)
    diff = np.asarray(med_rslt.summary() - framing_para_4231)
    assert_allclose(diff, 0, atol=1e-6) 

def test_framing_example_moderator_formula():

    cur_dir = os.path.dirname(os.path.abspath(__file__))
    data = pd.read_csv(os.path.join(cur_dir, 'results', "framing.csv"))

    probit = sm.families.links.probit
    outcome_model = sm.GLM.from_formula("cong_mesg ~ emo + treat*age + emo*age + educ + gender + income",
                                        data, family=sm.families.Binomial(link=probit()))

    mediator_model = sm.OLS.from_formula("emo ~ treat*age + educ + gender + income", data)

    moderators = {"age" : 20}
    med = Mediation(outcome_model, mediator_model, "treat", "emo",
                    moderators=moderators)

    np.random.seed(4231)
    med_rslt = med.fit(method='parametric', n_rep=100)
    diff = np.asarray(med_rslt.summary() - framing_moderated_4231)
    assert_allclose(diff, 0, atol=1e-6) 

def setup_class(cls):

        # generate artificial data
        np.random.seed(98765678)
        nobs = 200
        rvs = np.random.randn(nobs,6)
        data_exog = rvs
        data_exog = sm.add_constant(data_exog, prepend=False)
        xbeta = 0.1 + 0.1*rvs.sum(1)
        data_endog = np.random.poisson(np.exp(xbeta))

        #estimate discretemod.Poisson as benchmark
        cls.res_discrete = Poisson(data_endog, data_exog).fit(disp=0)

        mod_glm = sm.GLM(data_endog, data_exog, family=sm.families.Poisson())
        cls.res_glm = mod_glm.fit()

        #estimate generic MLE
        cls.mod = PoissonGMLE(data_endog, data_exog)
        cls.res = cls.mod.fit(start_params=0.9 * cls.res_discrete.params,
                                method='bfgs', disp=0) 

def setup_class(cls):
        fweights = [1, 1, 1, 2, 2, 2, 3, 3, 3, 1, 1, 1, 2, 2, 2, 3, 3]
        # faking aweights by using normalized freq_weights
        fweights = np.array(fweights)
        wsum = fweights.sum()
        nobs = len(cpunish_data.endog)
        aweights = fweights / wsum * nobs
        cls.corr_fact = np.sqrt((wsum - 1.) / wsum)

        cls.res1 = GLM(cpunish_data.endog, cpunish_data.exog,
                        family=sm.families.Poisson(), freq_weights=fweights
                        ).fit(cov_type='HC0') #, cov_kwds={'use_correction':False})
        # compare with discrete, start close to save time
        #modd = discrete.Poisson(cpunish_data.endog, cpunish_data.exog)
        cls.res2 = res_stata.results_poisson_fweight_hc1


# var_weights (aweights fail with HC, not properly implemented yet 

def setup_class(cls):
        fweights = [1, 1, 1, 2, 2, 2, 3, 3, 3, 1, 1, 1, 2, 2, 2, 3, 3]
        # faking aweights by using normalized freq_weights
        fweights = np.array(fweights)
        wsum = fweights.sum()
        nobs = len(cpunish_data.endog)
        aweights = fweights / wsum * nobs

        # This is really close when corr_fact = (wsum - 1.) / wsum, but to
        # avoid having loosen precision of the assert_allclose, I'm doing this
        # manually. Its *possible* lowering the IRLS convergence criterion
        # in stata and here will make this less sketchy.
        cls.corr_fact = np.sqrt((wsum - 1.) / wsum) * 0.98518473599905609
        cls.res1 = GLM(cpunish_data.endog, cpunish_data.exog,
                        family=sm.families.Poisson(), var_weights=aweights
                        ).fit(cov_type='HC0') #, cov_kwds={'use_correction':False})
        # compare with discrete, start close to save time
        # modd = discrete.Poisson(cpunish_data.endog, cpunish_data.exog)
        cls.res2 = res_stata.results_poisson_aweight_hc1 

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

#    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 

def setup_class(cls):
        np.random.seed(4321)
        n = 10000
        p = 5
        exog = np.empty((n, p))
        exog[:, 0] = 1
        exog[:, 1] = np.random.randint(low=-5, high=5, size=n)
        x = np.repeat(np.array([1, 2, 3, 4]), n / 4)
        exog[:, 2:] = get_dummies(x)
        beta = np.array([-1, 0.1, -0.05, .2, 0.35])
        lin_pred = (exog * beta).sum(axis=1)
        family = sm.families.Binomial
        link = sm.families.links.logit
        endog = gen_endog(lin_pred, family, link, binom_version=0)
        wt = np.random.randint(1, 5, n)
        mod1 = sm.GLM(endog, exog, family=family(link=link), freq_weights=wt)
        cls.res1 = mod1.fit()

        exog_dup = np.repeat(exog, wt, axis=0)
        endog_dup = np.repeat(endog, wt)
        mod2 = sm.GLM(endog_dup, exog_dup, family=family(link=link))
        cls.res2 = mod2.fit() 

def test_warnings_raised():
    if sys.version_info < (3, 4):
        raise SkipTest
    weights = [1, 1, 1, 2, 2, 2, 3, 3, 3, 1, 1, 1, 2, 2, 2, 3, 3]
    # faking aweights by using normalized freq_weights
    weights = np.array(weights)

    gid = np.arange(1, 17 + 1) // 2

    cov_kwds = {'groups': gid, 'use_correction': False}
    with warnings.catch_warnings(record=True) as w:
        res1 = GLM(cpunish_data.endog, cpunish_data.exog,
                   family=sm.families.Poisson(), freq_weights=weights
                   ).fit(cov_type='cluster', cov_kwds=cov_kwds)
        res1.summary()
        assert len(w) >= 1

    with warnings.catch_warnings(record=True) as w:
        res1 = GLM(cpunish_data.endog, cpunish_data.exog,
                   family=sm.families.Poisson(), var_weights=weights
                   ).fit(cov_type='cluster', cov_kwds=cov_kwds)
        res1.summary()
        assert len(w) >= 1 

def test_wtd_patsy_missing():
    from statsmodels.datasets.cpunish import load
    import pandas as pd
    data = load()
    data.exog[0, 0] = np.nan
    data.endog[[2, 4, 6, 8]] = np.nan
    data.pandas = pd.DataFrame(data.exog, columns=data.exog_name)
    data.pandas['EXECUTIONS'] = data.endog
    weights = np.arange(1, len(data.endog)+1)
    formula = """EXECUTIONS ~ INCOME + PERPOVERTY + PERBLACK + VC100k96 +
                 SOUTH + DEGREE"""
    mod_misisng = GLM.from_formula(formula, data=data.pandas,
                                   freq_weights=weights)
    assert_equal(mod_misisng.freq_weights.shape[0],
                 mod_misisng.endog.shape[0])
    assert_equal(mod_misisng.freq_weights.shape[0],
                 mod_misisng.exog.shape[0])
    assert_equal(mod_misisng.freq_weights.shape[0], 12)
    keep_weights = np.array([2,  4,  6,  8, 10, 11, 12, 13, 14, 15, 16, 17])
    assert_equal(mod_misisng.freq_weights, keep_weights) 

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_incompatible_input():
    weights = [1, 1, 1, 2, 2, 2, 3, 3, 3, 1, 1, 1, 2, 2, 2, 3, 3]
    exog = cpunish_data.exog
    endog = cpunish_data.endog
    family = sm.families.Poisson()
    # Too short
    assert_raises(ValueError, GLM, endog, exog, family=family,
                  freq_weights=weights[:-1])
    assert_raises(ValueError, GLM, endog, exog, family=family,
                  var_weights=weights[:-1])
    # Too long
    assert_raises(ValueError, GLM, endog, exog, family=family,
                  freq_weights=weights + [3])
    assert_raises(ValueError, GLM, endog, exog, family=family,
                  var_weights=weights + [3])

    # Too many dimensions
    assert_raises(ValueError, GLM, endog, exog, family=family,
                  freq_weights=[weights, weights])
    assert_raises(ValueError, GLM, endog, exog, family=family,
                  var_weights=[weights, weights]) 

def setup_class(cls):
        '''
        Tests Tweedie family with Power(1) link and var_power=2.
        '''
        from statsmodels.datasets.cpunish import load_pandas
        cls.data = load_pandas()
        cls.endog = cls.data.endog
        cls.exog = cls.data.exog[['INCOME', 'SOUTH']]
        np.random.seed(1234)
        cls.weight = np.random.randint(5, 100, len(cls.endog))
        cls.endog_big = np.repeat(cls.endog.values, cls.weight)
        cls.exog_big = np.repeat(cls.exog.values, cls.weight, axis=0)
        link = sm.families.links.Power
        family_link = sm.families.Tweedie(link=link, var_power=2)
        cls.res1 = GLM(cls.endog, cls.exog,
                        freq_weights=cls.weight,
                        family=family_link).fit()
        cls.res2 = GLM(cls.endog_big, cls.exog_big,
                        family=family_link).fit() 

def setup_class(cls):
        '''
        Tests Poisson family with canonical log link.

        Test results were obtained by R.
        '''
        from .results.results_glm import Cpunish
        from statsmodels.datasets.cpunish import load
        cls.data = load()
        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 

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")
            cls.res1 = GLM(res2.endog, res2.exog,
                           family=sm.families.Gamma(
                                link=sm.families.links.identity())
                            ).fit()
        cls.res2 = res2

#    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 

def setup_class(cls):

        '''
        Tests Negative Binomial family with canonical link
        g(p) = log(p/(p + 1/alpha))
        '''
        super(TestWtdGlmNegativeBinomial, cls).setup_class()
        alpha = 1.

        with warnings.catch_warnings():
            warnings.simplefilter("ignore", category=DomainWarning)
            family_link = sm.families.NegativeBinomial(
                link=sm.families.links.nbinom(alpha=alpha),
                alpha=alpha)
            cls.res1 = GLM(cls.endog, cls.exog,
                           freq_weights=cls.weight,
                           family=family_link).fit()
            cls.res2 = GLM(cls.endog_big, cls.exog_big,
                           family=family_link).fit() 

def setup_class(cls):

        '''
        Tests Binomial family with canonical logit link.
        '''
        super(TestWtdGlmBinomial, cls).setup_class()
        cls.endog = cls.endog / 100
        cls.endog_big = cls.endog_big / 100
        cls.res1 = GLM(cls.endog, cls.exog,
                        freq_weights=cls.weight,
                        family=sm.families.Binomial()).fit()
        cls.res2 = GLM(cls.endog_big, cls.exog_big,
                        family=sm.families.Binomial()).fit() 

def setup_class(cls):

        '''
        Tests Gamma family with log link.
        '''
        super(TestWtdGlmGamma, cls).setup_class()
        family_link = sm.families.Gamma(sm.families.links.log())
        cls.res1 = GLM(cls.endog, cls.exog,
                        freq_weights=cls.weight,
                        family=family_link).fit()
        cls.res2 = GLM(cls.endog_big, cls.exog_big,
                        family=family_link).fit() 

def test_glm_irls_method():
    nobs, k_vars = 50, 4
    np.random.seed(987126)
    x = np.random.randn(nobs, k_vars - 1)
    exog = add_constant(x, has_constant='add')
    y = exog.sum(1) + np.random.randn(nobs)

    mod = GLM(y, exog)
    res1 = mod.fit()
    res2 = mod.fit(wls_method='pinv', attach_wls=True)
    res3 = mod.fit(wls_method='qr', attach_wls=True)
    # fit_gradient does not attach mle_settings
    res_g1 = mod.fit(start_params=res1.params, method='bfgs')

    for r in [res1, res2, res3]:
        assert_equal(r.mle_settings['optimizer'], 'IRLS')
        assert_equal(r.method, 'IRLS')

    assert_equal(res1.mle_settings['wls_method'], 'lstsq')
    assert_equal(res2.mle_settings['wls_method'], 'pinv')
    assert_equal(res3.mle_settings['wls_method'], 'qr')

    assert_(hasattr(res2.results_wls.model, 'pinv_wexog'))
    assert_(hasattr(res3.results_wls.model, 'exog_Q'))

    # fit_gradient currently does not attach mle_settings
    assert_equal(res_g1.method, 'bfgs') 

def setup_class(cls):
        '''
        Tests InverseGuassian family with log link.
        '''
        super(TestWtdGlmInverseGaussian, cls).setup_class()
        family_link = sm.families.InverseGaussian(sm.families.links.log())
        cls.res1 = GLM(cls.endog, cls.exog,
                        freq_weights=cls.weight,
                        family=family_link).fit()
        cls.res2 = GLM(cls.endog_big, cls.exog_big,
                        family=family_link).fit() 

def setup_class(cls):
        '''
        Tests Poisson family with canonical log link.
        '''
        super(TestWtdGlmPoisson, cls).setup_class()
        cls.res1 = GLM(cls.endog, cls.exog,
                        freq_weights=cls.weight,
                        family=sm.families.Poisson()).fit()
        cls.res2 = GLM(cls.endog_big, cls.exog_big,
                        family=sm.families.Poisson()).fit() 

def test_perfect_pred():
    cur_dir = os.path.dirname(os.path.abspath(__file__))
    iris = np.genfromtxt(os.path.join(cur_dir, 'results', 'iris.csv'),
                         delimiter=",", skip_header=1)
    y = iris[:, -1]
    X = iris[:, :-1]
    X = X[y != 2]
    y = y[y != 2]
    X = add_constant(X, prepend=True)
    glm = GLM(y, X, family=sm.families.Binomial())
    with warnings.catch_warnings():
        warnings.simplefilter("ignore", category=RuntimeWarning)
        assert_raises(PerfectSeparationError, glm.fit)