Repository: statsmodels/statsmodels (Python)

Issue:
ENH: family: add option to turn of valid link check, or warn. 
see #1509

add option: `invalid_link : string in ['raise' 'warn', 'ignore']`   to turn exception into warning or ignore

needed to allow the GLM equivalent for Probit in #1870
#1509 also includes other options to improve how links are specified, that are not addressed here.

Also, binomial should allow probit and CDFLink, which are also not included here.



Hints:
@kshedden  This will also change which links can be used in GEE by turning off the exception.

test failure in GEE Multinomial   (I only ran test_glm.py locally)

```
======================================================================

ERROR: statsmodels.genmod.tests.test_gee.TestGEE.test_nominal
----------------------------------------------------------------------
Traceback (most recent call last):
File "/home/travis/miniconda/envs/statsmodels-test/lib/python2.7/site-packages/nose/case.py", line 197, in runTest
self.test(*self.arg)
File "/home/travis/miniconda/envs/statsmodels-test/lib/python2.7/site-packages/statsmodels-0.6.0-py2.7-linux-x86_64.egg/statsmodels/genmod/tests/test_gee.py", line 522, in test_nominal
family = Multinomial(3)
File "/home/travis/miniconda/envs/statsmodels-test/lib/python2.7/site-packages/statsmodels-0.6.0-py2.7-linux-x86_64.egg/statsmodels/genmod/generalized_estimating_equations.py", line 1946, in __init__
self.link = MultinomialLogit(self.ncut)
File "/home/travis/miniconda/envs/statsmodels-test/lib/python2.7/site-packages/statsmodels-0.6.0-py2.7-linux-x86_64.egg/statsmodels/genmod/families/family.py", line 59, in _setlink
invalid_link = self.invalid_link
AttributeError: 'Multinomial' object has no attribute 'invalid_link'
```

[![Coverage Status](https://coveralls.io/builds/1067981/badge)](https://coveralls.io/builds/1067981)

Coverage decreased (-0.02%) when pulling **723bdb8a51065257a101f01754bb124bfbc031c6 on josef-pkt:REF_genmod_link_check** into **2e806fc4203d528a17307f4fb08ab8bfe7ee43dd on statsmodels:master**.

about naming: `invalid_link` doesn't indicate that it is a flag, #2212 has `safe_links` as a list of links which sounds similar

aside:
CDFLink is valid for Binomial family because CDFLink is a subclass of Logit, so isinstance with Logit is true.
trying to rebase this
too many merge conflict, all `__init__`
`__init__` was refactored in main since this PR was opened.

AFAICS, it will be faster starting a new PR from scratch instead of editing all merge conflicts

only planned change is to use `check_link` as keyword and flag.

Patch:
diff --git a/statsmodels/discrete/discrete_model.py b/statsmodels/discrete/discrete_model.py
--- a/statsmodels/discrete/discrete_model.py
+++ b/statsmodels/discrete/discrete_model.py
@@ -2594,7 +2594,7 @@ class Probit(BinaryModel):
     @cache_readonly
     def link(self):
         from statsmodels.genmod.families import links
-        link = links.probit()
+        link = links.Probit()
         return link
 
     def cdf(self, X):
diff --git a/statsmodels/genmod/families/family.py b/statsmodels/genmod/families/family.py
--- a/statsmodels/genmod/families/family.py
+++ b/statsmodels/genmod/families/family.py
@@ -14,6 +14,9 @@
 from scipy import special, stats
 
 from statsmodels.compat.scipy import SP_LT_17
+from statsmodels.tools.sm_exceptions import (
+    ValueWarning,
+    )
 from . import links as L, varfuncs as V
 
 FLOAT_EPS = np.finfo(float).eps
@@ -59,13 +62,14 @@ def _setlink(self, link):
         #  <statsmodels.family.links.Power object at 0x9a4236c>]
         # for Poisson...
         self._link = link
-        if not isinstance(link, L.Link):
-            raise TypeError("The input should be a valid Link object.")
-        if hasattr(self, "links"):
-            validlink = max([isinstance(link, _) for _ in self.links])
-            if not validlink:
-                errmsg = "Invalid link for family, should be in %s. (got %s)"
-                raise ValueError(errmsg % (repr(self.links), link))
+        if self._check_link:
+            if not isinstance(link, L.Link):
+                raise TypeError("The input should be a valid Link object.")
+            if hasattr(self, "links"):
+                validlink = max([isinstance(link, _) for _ in self.links])
+                if not validlink:
+                    msg = "Invalid link for family, should be in %s. (got %s)"
+                    raise ValueError(msg % (repr(self.links), link))
 
     def _getlink(self):
         """
@@ -76,15 +80,16 @@ def _getlink(self):
     # link property for each family is a pointer to link instance
     link = property(_getlink, _setlink, doc="Link function for family")
 
-    def __init__(self, link, variance):
+    def __init__(self, link, variance, check_link=True):
+        self._check_link = check_link
         if inspect.isclass(link):
             warnmssg = (
                 "Calling Family(..) with a link class is not allowed. Use an "
                 "instance of a link class instead."
             )
             raise TypeError(warnmssg)
-        else:
-            self.link = link
+
+        self.link = link
         self.variance = variance
 
     def starting_mu(self, y):
@@ -397,10 +402,14 @@ class Poisson(Family):
     valid = [0, np.inf]
     safe_links = [L.Log, ]
 
-    def __init__(self, link=None):
+    def __init__(self, link=None, check_link=True):
         if link is None:
             link = L.Log()
-        super(Poisson, self).__init__(link=link, variance=Poisson.variance)
+        super(Poisson, self).__init__(
+            link=link,
+            variance=Poisson.variance,
+            check_link=check_link
+            )
 
     def _resid_dev(self, endog, mu):
         r"""
@@ -545,10 +554,14 @@ class Gaussian(Family):
     variance = V.constant
     safe_links = links
 
-    def __init__(self, link=None):
+    def __init__(self, link=None, check_link=True):
         if link is None:
             link = L.Identity()
-        super(Gaussian, self).__init__(link=link, variance=Gaussian.variance)
+        super(Gaussian, self).__init__(
+            link=link,
+            variance=Gaussian.variance,
+            check_link=check_link
+            )
 
     def _resid_dev(self, endog, mu):
         r"""
@@ -710,10 +723,14 @@ class Gamma(Family):
     variance = V.mu_squared
     safe_links = [L.Log, ]
 
-    def __init__(self, link=None):
+    def __init__(self, link=None, check_link=True):
         if link is None:
-            link = L.inverse_power()
-        super(Gamma, self).__init__(link=link, variance=Gamma.variance)
+            link = L.InversePower()
+        super(Gamma, self).__init__(
+            link=link,
+            variance=Gamma.variance,
+            check_link=check_link
+            )
 
     def _resid_dev(self, endog, mu):
         r"""
@@ -884,7 +901,7 @@ class Binomial(Family):
     # Other safe links, e.g. cloglog and probit are subclasses
     safe_links = [L.Logit, L.CDFLink]
 
-    def __init__(self, link=None):  # , n=1.):
+    def __init__(self, link=None, check_link=True):  # , n=1.):
         if link is None:
             link = L.Logit()
         # TODO: it *should* work for a constant n>1 actually, if freq_weights
@@ -892,8 +909,11 @@ def __init__(self, link=None):  # , n=1.):
         self.n = 1
         # overwritten by initialize if needed but always used to initialize
         # variance since endog is assumed/forced to be (0,1)
-        super(Binomial, self).__init__(link=link,
-                                       variance=V.Binomial(n=self.n))
+        super(Binomial, self).__init__(
+            link=link,
+            variance=V.Binomial(n=self.n),
+            check_link=check_link
+            )
 
     def starting_mu(self, y):
         r"""
@@ -1114,7 +1134,7 @@ class InverseGaussian(Family):
     link : a link instance, optional
         The default link for the inverse Gaussian family is the
         inverse squared link.
-        Available links are inverse_squared, inverse, log, and identity.
+        Available links are InverseSquared, Inverse, Log, and Identity.
         See statsmodels.genmod.families.links for more information.
 
     Attributes
@@ -1140,11 +1160,14 @@ class InverseGaussian(Family):
     variance = V.mu_cubed
     safe_links = [L.InverseSquared, L.Log, ]
 
-    def __init__(self, link=None):
+    def __init__(self, link=None, check_link=True):
         if link is None:
-            link = L.inverse_squared()
+            link = L.InverseSquared()
         super(InverseGaussian, self).__init__(
-            link=link, variance=InverseGaussian.variance)
+            link=link,
+            variance=InverseGaussian.variance,
+            check_link=check_link
+            )
 
     def _resid_dev(self, endog, mu):
         r"""
@@ -1310,15 +1333,19 @@ class NegativeBinomial(Family):
     variance = V.nbinom
     safe_links = [L.Log, ]
 
-    def __init__(self, link=None, alpha=1.):
+    def __init__(self, link=None, alpha=1., check_link=True):
         self.alpha = 1. * alpha  # make it at least float
         if alpha is self.__init__.__defaults__[1]:  # `is` is intentional
             warnings.warn("Negative binomial dispersion parameter alpha not "
-                          f"set. Using default value alpha={alpha}.")
+                          f"set. Using default value alpha={alpha}.",
+                          ValueWarning)
         if link is None:
             link = L.Log()
         super(NegativeBinomial, self).__init__(
-            link=link, variance=V.NegativeBinomial(alpha=self.alpha))
+            link=link,
+            variance=V.NegativeBinomial(alpha=self.alpha),
+            check_link=check_link
+            )
 
     def _resid_dev(self, endog, mu):
         r"""
@@ -1518,7 +1545,7 @@ class Tweedie(Family):
     variance = V.Power(power=1.5)
     safe_links = [L.Log, L.Power]
 
-    def __init__(self, link=None, var_power=1., eql=False):
+    def __init__(self, link=None, var_power=1., eql=False, check_link=True):
         self.var_power = var_power
         self.eql = eql
         if eql and (var_power < 1 or var_power > 2):
@@ -1527,7 +1554,10 @@ def __init__(self, link=None, var_power=1., eql=False):
         if link is None:
             link = L.Log()
         super(Tweedie, self).__init__(
-            link=link, variance=V.Power(power=var_power * 1.))
+            link=link,
+            variance=V.Power(power=var_power * 1.),
+            check_link=check_link
+            )
 
     def _resid_dev(self, endog, mu):
         r"""
diff --git a/statsmodels/genmod/families/links.py b/statsmodels/genmod/families/links.py
--- a/statsmodels/genmod/families/links.py
+++ b/statsmodels/genmod/families/links.py
@@ -15,6 +15,7 @@ def _link_deprecation_warning(old, new):
         f"link alias will be removed after the 0.15.0 release.",
         FutureWarning
     )
+    # raise
 
 
 class Link:
@@ -268,10 +269,10 @@ class Power(Link):
     Notes
     -----
     Aliases of Power:
-    inverse = Power(power=-1)
-    sqrt = Power(power=.5)
-    inverse_squared = Power(power=-2.)
-    identity = Power(power=1.)
+    Inverse = Power(power=-1)
+    Sqrt = Power(power=.5)
+    InverseSquared = Power(power=-2.)
+    Identity = Power(power=1.)
     """
 
     def __init__(self, power=1.):
diff --git a/statsmodels/genmod/generalized_estimating_equations.py b/statsmodels/genmod/generalized_estimating_equations.py
--- a/statsmodels/genmod/generalized_estimating_equations.py
+++ b/statsmodels/genmod/generalized_estimating_equations.py
@@ -2985,7 +2985,7 @@ class _Multinomial(families.Family):
     variance = varfuncs.binary
     safe_links = [_MultinomialLogit, ]
 
-    def __init__(self, nlevels):
+    def __init__(self, nlevels, check_link=True):
         """
         Parameters
         ----------
@@ -2993,6 +2993,7 @@ def __init__(self, nlevels):
             The number of distinct categories for the multinomial
             distribution.
         """
+        self._check_link = check_link
         self.initialize(nlevels)
 
     def initialize(self, nlevels):


Test patch:
diff --git a/statsmodels/discrete/tests/test_predict.py b/statsmodels/discrete/tests/test_predict.py
--- a/statsmodels/discrete/tests/test_predict.py
+++ b/statsmodels/discrete/tests/test_predict.py
@@ -385,7 +385,7 @@ def test_distr(case):
     if cls_model in models_influ:
         with warnings.catch_warnings():
             warnings.simplefilter("ignore", category=UserWarning)
-            # ZI models warn about missint hat_matrix_diag
+            # ZI models warn about missing hat_matrix_diag
             influ = res.get_influence()
             influ.summary_frame()
 
@@ -405,7 +405,7 @@ def test_distr(case):
 
         try:
             with warnings.catch_warnings():
-                # ZI models warn about missint hat_matrix_diag
+                # ZI models warn about missing hat_matrix_diag
                 warnings.simplefilter("ignore", category=UserWarning)
                 influ.plot_influence()
         except ImportError:
diff --git a/statsmodels/discrete/tests/test_truncated_model.py b/statsmodels/discrete/tests/test_truncated_model.py
--- a/statsmodels/discrete/tests/test_truncated_model.py
+++ b/statsmodels/discrete/tests/test_truncated_model.py
@@ -1,10 +1,15 @@
 
+import warnings
+
 import numpy as np
 from numpy.testing import assert_allclose, assert_equal
 
 from statsmodels import datasets
 from statsmodels.tools.tools import add_constant
 from statsmodels.tools.testing import Holder
+from statsmodels.tools.sm_exceptions import (
+    ConvergenceWarning,
+    )
 
 from statsmodels.distributions.discrete import (
     truncatedpoisson,
@@ -45,9 +50,12 @@ def test_bic(self):
 
     def test_fit_regularized(self):
         model = self.res1.model
-
         alpha = np.ones(len(self.res1.params))
-        res_reg = model.fit_regularized(alpha=alpha*0.01, disp=0)
+
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", category=ConvergenceWarning)
+            # This does not catch all Convergence warnings, why?
+            res_reg = model.fit_regularized(alpha=alpha*0.01, disp=0)
 
         assert_allclose(res_reg.params, self.res1.params,
                         rtol=1e-3, atol=5e-3)
diff --git a/statsmodels/genmod/families/tests/test_family.py b/statsmodels/genmod/families/tests/test_family.py
--- a/statsmodels/genmod/families/tests/test_family.py
+++ b/statsmodels/genmod/families/tests/test_family.py
@@ -11,6 +11,9 @@
 from scipy import integrate
 
 from statsmodels.compat.scipy import SP_LT_17
+from statsmodels.tools.sm_exceptions import (
+    ValueWarning,
+    )
 import statsmodels.genmod.families as F
 from statsmodels.genmod.families.family import Tweedie
 import statsmodels.genmod.families.links as L
@@ -56,6 +59,8 @@ def test_invalid_family_link(family, links):
                    "Using default value alpha=1.0.")
             warnings.filterwarnings("ignore", message=msg,
                                     category=UserWarning)
+            warnings.filterwarnings("ignore",
+                                    category=FutureWarning)
             for link in invalid_links:
                 family(link())
 
@@ -66,11 +71,26 @@ def test_family_link(family, links):
         msg = ("Negative binomial dispersion parameter alpha not set. "
                "Using default value alpha=1.0.")
         warnings.filterwarnings("ignore", message=msg,
-                                category=UserWarning)
+                                category=ValueWarning)
+        warnings.filterwarnings("ignore",
+                                category=FutureWarning)
         for link in links:
             assert family(link())
 
 
+@pytest.mark.parametrize("family, links", link_cases)
+def test_family_link_check(family, links):
+    # check that we can turn of all link checks
+    class Hugo():
+        pass
+    with warnings.catch_warnings():
+        msg = ("Negative binomial dispersion parameter alpha not set. "
+               "Using default value alpha=1.0.")
+        warnings.filterwarnings("ignore", message=msg,
+                                category=ValueWarning)
+        assert family(Hugo(), check_link=False)
+
+
 @pytest.mark.skipif(SP_LT_17, reason="Scipy too old, function not available")
 @pytest.mark.parametrize("power", (1.1, 1.5, 1.9))
 def test_tweedie_loglike_obs(power):
diff --git a/statsmodels/genmod/families/tests/test_link.py b/statsmodels/genmod/families/tests/test_link.py
--- a/statsmodels/genmod/families/tests/test_link.py
+++ b/statsmodels/genmod/families/tests/test_link.py
@@ -12,14 +12,14 @@
 # Family instances
 links = families.links
 logit = links.Logit()
-inverse_power = links.inverse_power()
-sqrt = links.sqrt()
-inverse_squared = links.inverse_squared()
-identity = links.identity()
-log = links.log()
-logc = links.logc()
-probit = links.probit()
-cauchy = links.cauchy()
+inverse_power = links.InversePower()
+sqrt = links.Sqrt()
+inverse_squared = links.InverseSquared()
+identity = links.Identity()
+log = links.Log()
+logc = links.LogC()
+probit = links.Probit()
+cauchy = links.Cauchy()
 cloglog = links.CLogLog()
 loglog = links.LogLog()
 negbinom = links.NegativeBinomial()
@@ -72,7 +72,7 @@ def test_deriv():
     for link in Links:
         for k in range(10):
             p = np.random.uniform(0, 1)
-            if isinstance(link, links.cauchy):
+            if isinstance(link, links.Cauchy):
                 p = np.clip(p, 0.03, 0.97)
             d = link.deriv(p)
             da = nd.approx_fprime(np.r_[p], link)
@@ -166,15 +166,15 @@ class CasesCDFLink():
     # just as namespace to hold cases for test_cdflink
 
     link_pairs = [
-        (links.CDFLink(dbn=stats.gumbel_l), links.cloglog()),
-        (links.CDFLink(dbn=stats.gumbel_r), links.loglog()),
-        (links.CDFLink(dbn=stats.norm), links.probit()),
-        (links.CDFLink(dbn=stats.logistic), links.logit()),
-        (links.CDFLink(dbn=stats.t(1)), links.cauchy()),
+        (links.CDFLink(dbn=stats.gumbel_l), links.CLogLog()),
+        (links.CDFLink(dbn=stats.gumbel_r), links.LogLog()),
+        (links.CDFLink(dbn=stats.norm), links.Probit()),
+        (links.CDFLink(dbn=stats.logistic), links.Logit()),
+        (links.CDFLink(dbn=stats.t(1)), links.Cauchy()),
         # approximation of t by normal is not good enough for rtol, atol
         # (links.CDFLink(dbn=stats.t(1000000)), links.probit()),
 
-        (MyCLogLog(), links.cloglog()),  # not a cdflink, but compares
+        (MyCLogLog(), links.CLogLog()),  # not a cdflink, but compares
         ]
 
     methods = ['__call__', 'deriv', 'inverse', 'inverse_deriv', 'deriv2',
diff --git a/statsmodels/genmod/tests/test_gee_glm.py b/statsmodels/genmod/tests/test_gee_glm.py
--- a/statsmodels/genmod/tests/test_gee_glm.py
+++ b/statsmodels/genmod/tests/test_gee_glm.py
@@ -110,7 +110,7 @@ class TestCompareGamma(CheckGEEGLM):
     def setup_class(cls):
         # adjusted for Gamma, not in test_gee.py
         vs = Independence()
-        family = families.Gamma(link=links.log())
+        family = families.Gamma(link=links.Log())
         np.random.seed(987126)
         #Y = np.random.normal(size=100)**2
         Y = np.exp(0.1 + np.random.normal(size=100))   # log-normal
diff --git a/statsmodels/genmod/tests/test_glm.py b/statsmodels/genmod/tests/test_glm.py
--- a/statsmodels/genmod/tests/test_glm.py
+++ b/statsmodels/genmod/tests/test_glm.py
@@ -427,7 +427,7 @@ def setup_class(cls):
                         0.001 * np.random.randn(nobs)
 
         GaussLog_Model = GLM(cls.lny, cls.X,
-                             family=sm.families.Gaussian(sm.families.links.log()))
+                             family=sm.families.Gaussian(sm.families.links.Log()))
         cls.res1 = GaussLog_Model.fit()
         from .results.results_glm import GaussianLog
         cls.res2 = GaussianLog()
@@ -457,7 +457,7 @@ def setup_class(cls):
         cls.X = np.c_[np.ones((nobs,1)),x,x**2]
         cls.y_inv = (1. + .02*x + .001*x**2)**-1 + .001 * np.random.randn(nobs)
         InverseLink_Model = GLM(cls.y_inv, cls.X,
-                family=sm.families.Gaussian(sm.families.links.inverse_power()))
+                family=sm.families.Gaussian(sm.families.links.InversePower()))
         InverseLink_Res = InverseLink_Model.fit()
         cls.res1 = InverseLink_Res
         from .results.results_glm import GaussianInverse
@@ -683,7 +683,7 @@ def setup_class(cls):
         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()
+            family=sm.families.Gamma(link=sm.families.links.Log())).fit()
         cls.res2 = res2
 
 # FIXME: enable or delete
@@ -709,7 +709,7 @@ def setup_class(cls):
         res2 = CancerIdentity()
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
-            fam = sm.families.Gamma(link=sm.families.links.identity())
+            fam = sm.families.Gamma(link=sm.families.links.Identity())
             cls.res1 = GLM(res2.endog, res2.exog, family=fam).fit()
         cls.res2 = res2
 
@@ -793,7 +793,7 @@ def setup_class(cls):
         res2 = InvGaussLog()
         cls.res1 = GLM(res2.endog, res2.exog,
             family=sm.families.InverseGaussian(
-                link=sm.families.links.log())).fit()
+                link=sm.families.links.Log())).fit()
         cls.res2 = res2
 
 # FIXME: enable or delete
@@ -820,7 +820,7 @@ def setup_class(cls):
             warnings.simplefilter("ignore")
             cls.res1 = GLM(data.endog, data.exog,
                             family=sm.families.InverseGaussian(
-                                link=sm.families.links.identity())).fit()
+                                link=sm.families.links.Identity())).fit()
         from .results.results_glm import InvGaussIdentity
         cls.res2 = InvGaussIdentity()
 
@@ -1098,7 +1098,7 @@ def test_formula_missing_exposure():
          'x': [1, 3, 2, 1.5]}
     df = pd.DataFrame(d)
 
-    family = sm.families.Gaussian(link=sm.families.links.log())
+    family = sm.families.Gaussian(link=sm.families.links.Log())
 
     mod = smf.glm("Foo ~ Bar", data=df, exposure=df.exposure,
                   family=family)
@@ -1235,12 +1235,12 @@ def test_gradient_irls():
 
     fam = sm.families
     lnk = sm.families.links
-    families = [(fam.Binomial, [lnk.logit, lnk.probit, lnk.cloglog, lnk.log, lnk.cauchy]),
-                (fam.Poisson, [lnk.log, lnk.identity, lnk.sqrt]),
-                (fam.Gamma, [lnk.log, lnk.identity, lnk.inverse_power]),
-                (fam.Gaussian, [lnk.identity, lnk.log, lnk.inverse_power]),
-                (fam.InverseGaussian, [lnk.log, lnk.identity, lnk.inverse_power, lnk.inverse_squared]),
-                (fam.NegativeBinomial, [lnk.log, lnk.inverse_power, lnk.inverse_squared, lnk.identity])]
+    families = [(fam.Binomial, [lnk.Logit, lnk.Probit, lnk.CLogLog, lnk.Log, lnk.Cauchy]),
+                (fam.Poisson, [lnk.Log, lnk.Identity, lnk.Sqrt]),
+                (fam.Gamma, [lnk.Log, lnk.Identity, lnk.InversePower]),
+                (fam.Gaussian, [lnk.Identity, lnk.Log, lnk.InversePower]),
+                (fam.InverseGaussian, [lnk.Log, lnk.Identity, lnk.InversePower, lnk.InverseSquared]),
+                (fam.NegativeBinomial, [lnk.Log, lnk.InversePower, lnk.InverseSquared, lnk.Identity])]
 
     n = 100
     p = 3
@@ -1255,36 +1255,36 @@ def test_gradient_irls():
                 if family_class != fam.Binomial and binom_version == 1:
                     continue
 
-                if (family_class, link) == (fam.Poisson, lnk.identity):
+                if (family_class, link) == (fam.Poisson, lnk.Identity):
                     lin_pred = 20 + exog.sum(1)
-                elif (family_class, link) == (fam.Binomial, lnk.log):
+                elif (family_class, link) == (fam.Binomial, lnk.Log):
                     lin_pred = -1 + exog.sum(1) / 8
-                elif (family_class, link) == (fam.Poisson, lnk.sqrt):
+                elif (family_class, link) == (fam.Poisson, lnk.Sqrt):
                     lin_pred = 2 + exog.sum(1)
-                elif (family_class, link) == (fam.InverseGaussian, lnk.log):
+                elif (family_class, link) == (fam.InverseGaussian, lnk.Log):
                     #skip_zero = True
                     lin_pred = -1 + exog.sum(1)
-                elif (family_class, link) == (fam.InverseGaussian, lnk.identity):
+                elif (family_class, link) == (fam.InverseGaussian, lnk.Identity):
                     lin_pred = 20 + 5*exog.sum(1)
                     lin_pred = np.clip(lin_pred, 1e-4, np.inf)
-                elif (family_class, link) == (fam.InverseGaussian, lnk.inverse_squared):
+                elif (family_class, link) == (fam.InverseGaussian, lnk.InverseSquared):
                     lin_pred = 0.5 + exog.sum(1) / 5
                     continue # skip due to non-convergence
-                elif (family_class, link) == (fam.InverseGaussian, lnk.inverse_power):
+                elif (family_class, link) == (fam.InverseGaussian, lnk.InversePower):
                     lin_pred = 1 + exog.sum(1) / 5
-                elif (family_class, link) == (fam.NegativeBinomial, lnk.identity):
+                elif (family_class, link) == (fam.NegativeBinomial, lnk.Identity):
                     lin_pred = 20 + 5*exog.sum(1)
                     lin_pred = np.clip(lin_pred, 1e-4, np.inf)
-                elif (family_class, link) == (fam.NegativeBinomial, lnk.inverse_squared):
+                elif (family_class, link) == (fam.NegativeBinomial, lnk.InverseSquared):
                     lin_pred = 0.1 + np.random.uniform(size=exog.shape[0])
                     continue # skip due to non-convergence
-                elif (family_class, link) == (fam.NegativeBinomial, lnk.inverse_power):
+                elif (family_class, link) == (fam.NegativeBinomial, lnk.InversePower):
                     lin_pred = 1 + exog.sum(1) / 5
 
-                elif (family_class, link) == (fam.Gaussian, lnk.inverse_power):
+                elif (family_class, link) == (fam.Gaussian, lnk.InversePower):
                     # adding skip because of convergence failure
                     skip_one = True
-                # the following fails with identity link, because endog < 0
+                # the following fails with Identity link, because endog < 0
                 # elif family_class == fam.Gamma:
                 #     lin_pred = 0.5 * exog.sum(1) + np.random.uniform(size=exog.shape[0])
                 else:
@@ -1297,9 +1297,9 @@ def test_gradient_irls():
                     mod_irls = sm.GLM(endog, exog, family=family_class(link=link()))
                 rslt_irls = mod_irls.fit(method="IRLS")
 
-                if not (family_class, link) in [(fam.Poisson, lnk.sqrt),
-                                                (fam.Gamma, lnk.inverse_power),
-                                                (fam.InverseGaussian, lnk.identity)
+                if not (family_class, link) in [(fam.Poisson, lnk.Sqrt),
+                                                (fam.Gamma, lnk.InversePower),
+                                                (fam.InverseGaussian, lnk.Identity)
                                                 ]:
                     check_score_hessian(rslt_irls)
 
@@ -1348,16 +1348,16 @@ def test_gradient_irls_eim():
 
     fam = sm.families
     lnk = sm.families.links
-    families = [(fam.Binomial, [lnk.logit, lnk.probit, lnk.cloglog, lnk.log,
-                                lnk.cauchy]),
-                (fam.Poisson, [lnk.log, lnk.identity, lnk.sqrt]),
-                (fam.Gamma, [lnk.log, lnk.identity, lnk.inverse_power]),
-                (fam.Gaussian, [lnk.identity, lnk.log, lnk.inverse_power]),
-                (fam.InverseGaussian, [lnk.log, lnk.identity,
-                                       lnk.inverse_power,
-                                       lnk.inverse_squared]),
-                (fam.NegativeBinomial, [lnk.log, lnk.inverse_power,
-                                        lnk.inverse_squared, lnk.identity])]
+    families = [(fam.Binomial, [lnk.Logit, lnk.Probit, lnk.CLogLog, lnk.Log,
+                                lnk.Cauchy]),
+                (fam.Poisson, [lnk.Log, lnk.Identity, lnk.Sqrt]),
+                (fam.Gamma, [lnk.Log, lnk.Identity, lnk.InversePower]),
+                (fam.Gaussian, [lnk.Identity, lnk.Log, lnk.InversePower]),
+                (fam.InverseGaussian, [lnk.Log, lnk.Identity,
+                                       lnk.InversePower,
+                                       lnk.InverseSquared]),
+                (fam.NegativeBinomial, [lnk.Log, lnk.InversePower,
+                                        lnk.InverseSquared, lnk.Identity])]
 
     n = 100
     p = 3
@@ -1372,39 +1372,39 @@ def test_gradient_irls_eim():
                 if family_class != fam.Binomial and binom_version == 1:
                     continue
 
-                if (family_class, link) == (fam.Poisson, lnk.identity):
+                if (family_class, link) == (fam.Poisson, lnk.Identity):
                     lin_pred = 20 + exog.sum(1)
-                elif (family_class, link) == (fam.Binomial, lnk.log):
+                elif (family_class, link) == (fam.Binomial, lnk.Log):
                     lin_pred = -1 + exog.sum(1) / 8
-                elif (family_class, link) == (fam.Poisson, lnk.sqrt):
+                elif (family_class, link) == (fam.Poisson, lnk.Sqrt):
                     lin_pred = 2 + exog.sum(1)
-                elif (family_class, link) == (fam.InverseGaussian, lnk.log):
+                elif (family_class, link) == (fam.InverseGaussian, lnk.Log):
                     # skip_zero = True
                     lin_pred = -1 + exog.sum(1)
                 elif (family_class, link) == (fam.InverseGaussian,
-                                              lnk.identity):
+                                              lnk.Identity):
                     lin_pred = 20 + 5*exog.sum(1)
                     lin_pred = np.clip(lin_pred, 1e-4, np.inf)
                 elif (family_class, link) == (fam.InverseGaussian,
-                                              lnk.inverse_squared):
+                                              lnk.InverseSquared):
                     lin_pred = 0.5 + exog.sum(1) / 5
                     continue  # skip due to non-convergence
                 elif (family_class, link) == (fam.InverseGaussian,
-                                              lnk.inverse_power):
+                                              lnk.InversePower):
                     lin_pred = 1 + exog.sum(1) / 5
                 elif (family_class, link) == (fam.NegativeBinomial,
-                                              lnk.identity):
+                                              lnk.Identity):
                     lin_pred = 20 + 5*exog.sum(1)
                     lin_pred = np.clip(lin_pred, 1e-4, np.inf)
                 elif (family_class, link) == (fam.NegativeBinomial,
-                                              lnk.inverse_squared):
+                                              lnk.InverseSquared):
                     lin_pred = 0.1 + np.random.uniform(size=exog.shape[0])
                     continue  # skip due to non-convergence
                 elif (family_class, link) == (fam.NegativeBinomial,
-                                              lnk.inverse_power):
+                                              lnk.InversePower):
                     lin_pred = 1 + exog.sum(1) / 5
 
-                elif (family_class, link) == (fam.Gaussian, lnk.inverse_power):
+                elif (family_class, link) == (fam.Gaussian, lnk.InversePower):
                     # adding skip because of convergence failure
                     skip_one = True
                 else:
@@ -1701,7 +1701,7 @@ def setup_class(cls):
         with warnings.catch_warnings():
             warnings.simplefilter("ignore", category=DomainWarning)
             family_link = sm.families.NegativeBinomial(
-                link=sm.families.links.nbinom(alpha=alpha),
+                link=sm.families.links.NegativeBinomial(alpha=alpha),
                 alpha=alpha)
             cls.res1 = GLM(cls.endog, cls.exog,
                            freq_weights=cls.weight,
@@ -1718,7 +1718,7 @@ def setup_class(cls):
         Tests Gamma family with log link.
         '''
         super(TestWtdGlmGamma, cls).setup_class()
-        family_link = sm.families.Gamma(sm.families.links.log())
+        family_link = sm.families.Gamma(sm.families.links.Log())
         cls.res1 = GLM(cls.endog, cls.exog,
                        freq_weights=cls.weight,
                        family=family_link).fit()
@@ -1733,7 +1733,7 @@ def setup_class(cls):
         Tests Gaussian family with log link.
         '''
         super(TestWtdGlmGaussian, cls).setup_class()
-        family_link = sm.families.Gaussian(sm.families.links.log())
+        family_link = sm.families.Gaussian(sm.families.links.Log())
         cls.res1 = GLM(cls.endog, cls.exog,
                        freq_weights=cls.weight,
                        family=family_link).fit()
@@ -1748,7 +1748,7 @@ def setup_class(cls):
         Tests InverseGaussian family with log link.
         '''
         super(TestWtdGlmInverseGaussian, cls).setup_class()
-        family_link = sm.families.InverseGaussian(sm.families.links.log())
+        family_link = sm.families.InverseGaussian(sm.families.links.Log())
         cls.res1 = GLM(cls.endog, cls.exog,
                        freq_weights=cls.weight,
                        family=family_link).fit()
@@ -1763,7 +1763,7 @@ def setup_class(cls):
         Tests Gamma family with log link.
         '''
         super(TestWtdGlmGammaNewton, cls).setup_class()
-        family_link = sm.families.Gamma(sm.families.links.log())
+        family_link = sm.families.Gamma(sm.families.links.Log())
         cls.res1 = GLM(cls.endog, cls.exog,
                        freq_weights=cls.weight,
                        family=family_link
@@ -1773,7 +1773,7 @@ def setup_class(cls):
                        ).fit(method='newton')
 
     def test_init_kwargs(self):
-        family_link = sm.families.Gamma(sm.families.links.log())
+        family_link = sm.families.Gamma(sm.families.links.Log())
 
         with pytest.warns(ValueWarning, match="unknown kwargs"):
             GLM(self.endog, self.exog, family=family_link,
@@ -1788,7 +1788,7 @@ def setup_class(cls):
         Tests Gamma family with log link.
         '''
         super(TestWtdGlmGammaScale_X2, cls).setup_class()
-        family_link = sm.families.Gamma(sm.families.links.log())
+        family_link = sm.families.Gamma(sm.families.links.Log())
         cls.res1 = GLM(cls.endog, cls.exog,
                        freq_weights=cls.weight,
                        family=family_link,
@@ -1805,7 +1805,7 @@ def setup_class(cls):
         Tests Gamma family with log link.
         '''
         super(TestWtdGlmGammaScale_dev, cls).setup_class()
-        family_link = sm.families.Gamma(sm.families.links.log())
+        family_link = sm.families.Gamma(sm.families.links.Log())
         cls.res1 = GLM(cls.endog, cls.exog,
                        freq_weights=cls.weight,
                        family=family_link,
@@ -1837,7 +1837,7 @@ def setup_class(cls):
         Tests Tweedie family with log link and var_power=1.
         '''
         super(TestWtdTweedieLog, cls).setup_class()
-        family_link = sm.families.Tweedie(link=sm.families.links.log(),
+        family_link = sm.families.Tweedie(link=sm.families.links.Log(),
                                           var_power=1)
         cls.res1 = GLM(cls.endog, cls.exog,
                         freq_weights=cls.weight,
@@ -2000,7 +2000,7 @@ def setup_class(cls):
         cls.data = cpunish.load_pandas()
         cls.exog = cls.data.exog[['INCOME', 'SOUTH']]
         cls.endog = cls.data.endog
-        family_link = sm.families.Tweedie(link=sm.families.links.log(),
+        family_link = sm.families.Tweedie(link=sm.families.links.Log(),
                                           var_power=1.)
         cls.res1 = sm.GLM(endog=cls.data.endog,
                           exog=cls.data.exog[['INCOME', 'SOUTH']],
@@ -2015,7 +2015,7 @@ def setup_class(cls):
 
         from .results.results_glm import FairTweedieLog15
         data = load_pandas()
-        family_link = sm.families.Tweedie(link=sm.families.links.log(),
+        family_link = sm.families.Tweedie(link=sm.families.links.Log(),
                                           var_power=1.5)
         cls.res1 = sm.GLM(endog=data.endog,
                           exog=data.exog[['rate_marriage', 'age',
@@ -2048,11 +2048,11 @@ def setup_class(cls):
         cls.data = cpunish.load_pandas()
         cls.exog = cls.data.exog[['INCOME', 'SOUTH']]
         cls.endog = cls.data.endog
-        family1 = sm.families.Gaussian(link=sm.families.links.log())
+        family1 = sm.families.Gaussian(link=sm.families.links.Log())
         cls.res1 = sm.GLM(endog=cls.data.endog,
                           exog=cls.data.exog[['INCOME', 'SOUTH']],
                           family=family1).fit()
-        family2 = sm.families.Tweedie(link=sm.families.links.log(),
+        family2 = sm.families.Tweedie(link=sm.families.links.Log(),
                                       var_power=0)
         cls.res2 = sm.GLM(endog=cls.data.endog,
                           exog=cls.data.exog[['INCOME', 'SOUTH']],
@@ -2065,11 +2065,11 @@ def setup_class(cls):
         cls.data = cpunish.load_pandas()
         cls.exog = cls.data.exog[['INCOME', 'SOUTH']]
         cls.endog = cls.data.endog
-        family1 = sm.families.Poisson(link=sm.families.links.log())
+        family1 = sm.families.Poisson(link=sm.families.links.Log())
         cls.res1 = sm.GLM(endog=cls.data.endog,
                           exog=cls.data.exog[['INCOME', 'SOUTH']],
                           family=family1).fit()
-        family2 = sm.families.Tweedie(link=sm.families.links.log(),
+        family2 = sm.families.Tweedie(link=sm.families.links.Log(),
                                       var_power=1)
         cls.res2 = sm.GLM(endog=cls.data.endog,
                           exog=cls.data.exog[['INCOME', 'SOUTH']],
@@ -2082,11 +2082,11 @@ def setup_class(cls):
         cls.data = cpunish.load_pandas()
         cls.exog = cls.data.exog[['INCOME', 'SOUTH']]
         cls.endog = cls.data.endog
-        family1 = sm.families.Gamma(link=sm.families.links.log())
+        family1 = sm.families.Gamma(link=sm.families.links.Log())
         cls.res1 = sm.GLM(endog=cls.data.endog,
                           exog=cls.data.exog[['INCOME', 'SOUTH']],
                           family=family1).fit()
-        family2 = sm.families.Tweedie(link=sm.families.links.log(),
+        family2 = sm.families.Tweedie(link=sm.families.links.Log(),
                                       var_power=2)
         cls.res2 = sm.GLM(endog=cls.data.endog,
                           exog=cls.data.exog[['INCOME', 'SOUTH']],
@@ -2099,11 +2099,11 @@ def setup_class(cls):
         cls.data = cpunish.load_pandas()
         cls.exog = cls.data.exog[['INCOME', 'SOUTH']]
         cls.endog = cls.data.endog
-        family1 = sm.families.InverseGaussian(link=sm.families.links.log())
+        family1 = sm.families.InverseGaussian(link=sm.families.links.Log())
         cls.res1 = sm.GLM(endog=cls.data.endog,
                           exog=cls.data.exog[['INCOME', 'SOUTH']],
                           family=family1).fit()
-        family2 = sm.families.Tweedie(link=sm.families.links.log(),
+        family2 = sm.families.Tweedie(link=sm.families.links.Log(),
                                       var_power=3)
         cls.res2 = sm.GLM(endog=cls.data.endog,
                           exog=cls.data.exog[['INCOME', 'SOUTH']],
@@ -2278,12 +2278,12 @@ def testTweediePowerEstimate():
          1.56354040e-09,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00]
     model1 = sm.GLM(y, data.exog[['INCOME', 'SOUTH']],
-                    family=sm.families.Tweedie(link=sm.families.links.log(),
+                    family=sm.families.Tweedie(link=sm.families.links.Log(),
                                                var_power=1.5))
     res1 = model1.fit()
     model2 = sm.GLM((y - res1.mu) ** 2,
                     np.column_stack((np.ones(len(res1.mu)), np.log(res1.mu))),
-                    family=sm.families.Gamma(sm.families.links.log()))
+                    family=sm.families.Gamma(sm.families.links.Log()))
     res2 = model2.fit()
     # Sample may be too small for this...
     # assert_allclose(res1.scale, np.exp(res2.params[0]), rtol=0.25)
diff --git a/statsmodels/genmod/tests/test_glm_weights.py b/statsmodels/genmod/tests/test_glm_weights.py
--- a/statsmodels/genmod/tests/test_glm_weights.py
+++ b/statsmodels/genmod/tests/test_glm_weights.py
@@ -305,7 +305,7 @@ def setup_class(cls):
                 data=data,
                 family=sm.families.Tweedie(
                     var_power=1.55,
-                    link=sm.families.links.log()
+                    link=sm.families.links.Log()
                     ),
                 var_weights=aweights
         )
@@ -326,7 +326,7 @@ def setup_class(cls):
         aweights[::5] = 5
         aweights[::13] = 3
         model = sm.GLM(endog, exog,
-                       family=sm.families.Gamma(link=sm.families.links.log()),
+                       family=sm.families.Gamma(link=sm.families.links.Log()),
                        var_weights=aweights)
         cls.res1 = model.fit(rtol=1e-25, atol=0)
         cls.res2 = res_r.results_gamma_aweights_nonrobust
@@ -355,7 +355,7 @@ def setup_class(cls):
         model = smf.glm(
                 'EXECUTIONS ~ INCOME + SOUTH - 1',
                 data=data,
-                family=sm.families.Gaussian(link=sm.families.links.log()),
+                family=sm.families.Gaussian(link=sm.families.links.Log()),
                 var_weights=aweights
         )
         cls.res1 = model.fit(rtol=1e-25, atol=0)
@@ -432,16 +432,16 @@ def test_wtd_gradient_irls():
 
     fam = sm.families
     lnk = sm.families.links
-    families = [(fam.Binomial, [lnk.logit, lnk.probit, lnk.cloglog, lnk.log,
-                                lnk.cauchy]),
-                (fam.Poisson, [lnk.log, lnk.identity, lnk.sqrt]),
-                (fam.Gamma, [lnk.log, lnk.identity, lnk.inverse_power]),
-                (fam.Gaussian, [lnk.identity, lnk.log, lnk.inverse_power]),
-                (fam.InverseGaussian, [lnk.log, lnk.identity,
-                                       lnk.inverse_power,
-                                       lnk.inverse_squared]),
-                (fam.NegativeBinomial, [lnk.log, lnk.inverse_power,
-                                        lnk.inverse_squared, lnk.identity])]
+    families = [(fam.Binomial, [lnk.Logit, lnk.Probit, lnk.CLogLog, lnk.Log,
+                                lnk.Cauchy]),
+                (fam.Poisson, [lnk.Log, lnk.Identity, lnk.Sqrt]),
+                (fam.Gamma, [lnk.Log, lnk.Identity, lnk.InversePower]),
+                (fam.Gaussian, [lnk.Identity, lnk.Log, lnk.InversePower]),
+                (fam.InverseGaussian, [lnk.Log, lnk.Identity,
+                                       lnk.InversePower,
+                                       lnk.InverseSquared]),
+                (fam.NegativeBinomial, [lnk.Log, lnk.InversePower,
+                                        lnk.InverseSquared, lnk.Identity])]
 
     n = 100
     p = 3
@@ -456,67 +456,67 @@ def test_wtd_gradient_irls():
 
                 if family_class != fam.Binomial and binom_version == 1:
                     continue
-                elif family_class == fam.Binomial and link == lnk.cloglog:
+                elif family_class == fam.Binomial and link == lnk.CLogLog:
                     # Cannot get gradient to converage with var_weights here
                     continue
-                elif family_class == fam.Binomial and link == lnk.log:
+                elif family_class == fam.Binomial and link == lnk.Log:
                     # Cannot get gradient to converage with var_weights here
                     continue
-                elif (family_class, link) == (fam.Poisson, lnk.identity):
+                elif (family_class, link) == (fam.Poisson, lnk.Identity):
                     lin_pred = 20 + exog.sum(1)
-                elif (family_class, link) == (fam.Binomial, lnk.log):
+                elif (family_class, link) == (fam.Binomial, lnk.Log):
                     lin_pred = -1 + exog.sum(1) / 8
-                elif (family_class, link) == (fam.Poisson, lnk.sqrt):
+                elif (family_class, link) == (fam.Poisson, lnk.Sqrt):
                     lin_pred = -2 + exog.sum(1)
-                elif (family_class, link) == (fam.Gamma, lnk.log):
+                elif (family_class, link) == (fam.Gamma, lnk.Log):
                     # Cannot get gradient to converge with var_weights here
                     continue
-                elif (family_class, link) == (fam.Gamma, lnk.identity):
+                elif (family_class, link) == (fam.Gamma, lnk.Identity):
                     # Cannot get gradient to converage with var_weights here
                     continue
-                elif (family_class, link) == (fam.Gamma, lnk.inverse_power):
+                elif (family_class, link) == (fam.Gamma, lnk.InversePower):
                     # Cannot get gradient to converage with var_weights here
                     continue
-                elif (family_class, link) == (fam.Gaussian, lnk.log):
+                elif (family_class, link) == (fam.Gaussian, lnk.Log):
                     # Cannot get gradient to converage with var_weights here
                     continue
-                elif (family_class, link) == (fam.Gaussian, lnk.inverse_power):
+                elif (family_class, link) == (fam.Gaussian, lnk.InversePower):
                     # Cannot get gradient to converage with var_weights here
                     continue
-                elif (family_class, link) == (fam.InverseGaussian, lnk.log):
+                elif (family_class, link) == (fam.InverseGaussian, lnk.Log):
                     # Cannot get gradient to converage with var_weights here
                     lin_pred = -1 + exog.sum(1)
                     continue
                 elif (family_class, link) == (fam.InverseGaussian,
-                                              lnk.identity):
+                                              lnk.Identity):
                     # Cannot get gradient to converage with var_weights here
                     lin_pred = 20 + 5*exog.sum(1)
                     lin_pred = np.clip(lin_pred, 1e-4, np.inf)
                     continue
                 elif (family_class, link) == (fam.InverseGaussian,
-                                              lnk.inverse_squared):
+                                              lnk.InverseSquared):
                     lin_pred = 0.5 + exog.sum(1) / 5
                     continue  # skip due to non-convergence
                 elif (family_class, link) == (fam.InverseGaussian,
-                                              lnk.inverse_power):
+                                              lnk.InversePower):
                     lin_pred = 1 + exog.sum(1) / 5
                     method = 'newton'
                 elif (family_class, link) == (fam.NegativeBinomial,
-                                              lnk.identity):
+                                              lnk.Identity):
                     lin_pred = 20 + 5*exog.sum(1)
                     lin_pred = np.clip(lin_pred, 1e-3, np.inf)
                     method = 'newton'
                 elif (family_class, link) == (fam.NegativeBinomial,
-                                              lnk.inverse_squared):
+                                              lnk.InverseSquared):
                     lin_pred = 0.1 + np.random.uniform(size=exog.shape[0])
                     continue  # skip due to non-convergence
                 elif (family_class, link) == (fam.NegativeBinomial,
-                                              lnk.inverse_power):
+                                              lnk.InversePower):
                     # Cannot get gradient to converage with var_weights here
                     lin_pred = 1 + exog.sum(1) / 5
                     continue
 
-                elif (family_class, link) == (fam.Gaussian, lnk.inverse_power):
+                elif (family_class, link) == (fam.Gaussian, lnk.InversePower):
                     # adding skip because of convergence failure
                     skip_one = True
                 else:
@@ -601,7 +601,7 @@ def setup_class(cls):
         beta = np.array([-1, 0.1, -0.05, .2, 0.35])
         lin_pred = (exog * beta).sum(axis=1)
         family = sm.families.Poisson
-        link = sm.families.links.log
+        link = sm.families.links.Log
         endog = gen_endog(lin_pred, family, link)
         mod1 = sm.GLM(endog, exog, family=family(link=link()))
         cls.res1 = mod1.fit()
@@ -632,7 +632,7 @@ def setup_class(cls):
         beta = np.array([-1, 0.1, -0.05, .2, 0.35])
         lin_pred = (exog * beta).sum(axis=1)
         family = sm.families.Poisson
-        link = sm.families.links.log
+        link = sm.families.links.Log
         endog = gen_endog(lin_pred, family, link)
         mod1 = sm.GLM(endog, exog, family=family(link=link()))
         cls.res1 = mod1.fit()
@@ -663,7 +663,7 @@ def setup_class(cls):
         beta = np.array([7, 0.1, -0.05, .2, 0.35])
         lin_pred = (exog * beta).sum(axis=1)
         family = sm.families.Tweedie
-        link = sm.families.links.log
+        link = sm.families.links.Log
         endog = gen_endog(lin_pred, family, link)
         mod1 = sm.GLM(endog, exog, family=family(link=link(), var_power=1.5))
         cls.res1 = mod1.fit(rtol=1e-20, atol=0, tol_criterion='params')
@@ -695,7 +695,7 @@ def setup_class(cls):
         beta = np.array([7, 0.1, -0.05, .2, 0.35])
         lin_pred = (exog * beta).sum(axis=1)
         family = sm.families.Tweedie
-        link = sm.families.links.log
+        link = sm.families.links.Log
         endog = gen_endog(lin_pred, family, link)
         mod1 = sm.GLM(endog, exog, family=family(link=link(), var_power=1.5))
         cls.res1 = mod1.fit(rtol=1e-10, atol=0, tol_criterion='params',
@@ -728,7 +728,7 @@ def setup_class(cls):
         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
+        link = sm.families.links.Logit
         endog = gen_endog(lin_pred, family, link, binom_version=0)
         mod1 = sm.GLM(endog, exog, family=family(link=link()))
         cls.res1 = mod1.fit(rtol=1e-10, atol=0, tol_criterion='params',
@@ -761,7 +761,7 @@ def setup_class(cls):
         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
+        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)
@@ -855,7 +855,7 @@ def setup_class(cls):
         model = smf.glm(
                 'EXECUTIONS ~ INCOME + SOUTH - 1',
                 data=data,
-                family=sm.families.Gaussian(link=sm.families.links.identity()),
+                family=sm.families.Gaussian(link=sm.families.links.Identity()),
                 var_weights=aweights
         )
         wlsmodel = smf.wls(