ENH/REF GLM simplify computation if is_canonical

[josef-pkt]

triggered by looking at the computation of IRLS weights in #4267

We currently don't use any special casing and simplified computation in GLM when a canonical link function is used.

This would create some speed improvements, but the more important point is that it will most likely increase numerical stability and precision.

Example Pregibon 1981 on regression diagnostics points out that the weights in IRLS are zero for the corner case p=0 and p=1 in Binomial. However, we divide by the variance which divides by zero in this case, most likely 0/0. I guess we get currently around this with clipping away from the boundaries.

Other maybe more important cases are gamma and negative binomial where canonical link functions don't work so well, and the non-canonical log link works better.
I haven't checked, but my guess is that some of the messy instabilities in corner cases will likely disappear.

(a bit related: thequackdaddy added an EIM option for the hessian in fit to avoid the complicated OIM which is generic for non-canonical cases. However, we don't skip the unnecessary part for OIM in the canonical link case, where OIM==EIM, IIRC)

[josef-pkt]

for gamma: link, link_deriv and varfunc are all power functions
#3316 which refers to discussion and comments in #2208 (about the DomainWarning that we added for the canonical link in gamma)

[josef-pkt]

implementation for irls weights

• family has the link as argument and attribute and could check for is_canonical
• family computes the weights in generic method inherited by specific families
• this could delegate to a family specific _weights_canonical and delegate to it if is_canonical

for gradient, nonlinear optimizers

• loglike/loglikeobs is defined in families and it doesn't look like there is a straightforward simplification, argument is mu which already incorporates the link, so computation is independent of link.
• score_factor has a product self.family.link.deriv(mu) * self.family.variance(mu) which will might in the canonical case. score uses score_obs which uses score_factor, so score_factor is the only relevant part (no code duplication here)
• hessian_factor, there is a shortcut for observed=False, which AFAIR should apply for the canonical link case. the computation for eim_factor
  eim_factor = 1 / (self.family.link.deriv(mu)**2 * self.family.variance(mu))
  is essentially the same as the irls weights and will simplify in canonical case.
  oim has some messier terms which however won't be relevant for the canonical case

the remaining parts, including deviance, are mainly for post-estimation where speed or numerical robustness is less relevant.

so it looks like we might want or need to new special methods in the individual families

weights = eim_factor = 1 / (self.family.link.deriv(mu)**2 * self.family.variance(mu))
and likely also
??? = self.family.variance.deriv(mu) * self.family.link.deriv(mu)

Caveat:
is_canonical assumes most likely that the lin_pred is linear, however derivatives are mostly with respect to mu
We need to avoid too much special casing or make the condition for special casing narrow enough, so that it also does the right thing if lin_pred is actually nonlinear, e.g. nonlinear GLM as subclass of GLM
#2179 (comment)

[josef-pkt]

adding label good as first pr

I think the irls weights part is easy to implement, but requires a reference with the list of weight functions for the canonical cases, maybe SPSS which has a good list of formulas in the technical documentation. (I don't have it in my current pile of papers to quickly check.)