
R version 3.6.2 Patched (2020-02-12 r77795) -- "Dark and Stormy Night"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ## Example in which the fit for the null deviance fails to converge:
> # https://stat.ethz.ch/pipermail/r-help/2012-May/313161.html
> Y <- c(rep(0,35),1,2,0,6,8,16,43)
> beta <- 42:1
> cst <- lchoose(42, beta)
> tau <- (beta^2)/2
> fit <- glm(formula = Y ~ offset(cst) + beta + tau, family = poisson)
Warning messages:
1: glm.fit: algorithm did not converge 
2: In glm(formula = Y ~ offset(cst) + beta + tau, family = poisson) :
  fitting to calculate the null deviance did not converge -- increase 'maxit'?
> 
> ## Ensure make.link() consistency:
> linkNames <- c("logit", "probit", "cauchit", "cloglog",
+                "identity",
+                "log",  "sqrt",  "1/mu^2", "inverse")
> links <- lapply(setNames(,linkNames), make.link)
> fns <- c("linkfun", "linkinv", "mu.eta", "valideta")
> stopifnot(exprs = {
+     is.matrix(nms <- sapply(links, names)) # matching number & type
+     is.character(nms)
+     nms[,1] == nms ## all columns are the same
+     identical(setNames(,linkNames), vapply(links, `[[`, "", "name"))
+     fns %in% nms[,1]
+ })
> links <- lapply(links, `[`, fns) # functions only
> stopifnot(unlist(lapply(links, function(L) vapply(L, is.function, NA))))
> ## all functions having consistent arguments :
> lf <- lapply(links, function(L) lapply(L, formals))
> stopifnot(exprs = { ## all functions have 1 argument
+     unlist(lapply(lf, lengths), recursive=FALSE) == 1L
+     is.matrix(argNms <- sapply(lf, function(L) vapply(L, names, "")))
+     argNms[,1] == argNms ## all columns are the same
+ })
> noquote(t(argNms))
         linkfun linkinv mu.eta valideta
logit    mu      eta     eta    eta     
probit   mu      eta     eta    eta     
cauchit  mu      eta     eta    eta     
cloglog  mu      eta     eta    eta     
identity mu      eta     eta    eta     
log      mu      eta     eta    eta     
sqrt     mu      eta     eta    eta     
1/mu^2   mu      eta     eta    eta     
inverse  mu      eta     eta    eta     
> 
> ## Calling all functions
> ## 1. valideta
> stopifnot(vv <- vapply(links, function(L) L$valideta((1:3)/4), NA))
> ## 2. all others
> other <- fns != "valideta"
> str(linkO <- lapply(links, function(L) L[other]))
List of 9
 $ logit   :List of 3
  ..$ linkfun:function (mu)  
  ..$ linkinv:function (eta)  
  ..$ mu.eta :function (eta)  
 $ probit  :List of 3
  ..$ linkfun:function (mu)  
  ..$ linkinv:function (eta)  
  ..$ mu.eta :function (eta)  
 $ cauchit :List of 3
  ..$ linkfun:function (mu)  
  ..$ linkinv:function (eta)  
  ..$ mu.eta :function (eta)  
 $ cloglog :List of 3
  ..$ linkfun:function (mu)  
  ..$ linkinv:function (eta)  
  ..$ mu.eta :function (eta)  
 $ identity:List of 3
  ..$ linkfun:function (mu)  
  ..$ linkinv:function (eta)  
  ..$ mu.eta :function (eta)  
 $ log     :List of 3
  ..$ linkfun:function (mu)  
  ..$ linkinv:function (eta)  
  ..$ mu.eta :function (eta)  
 $ sqrt    :List of 3
  ..$ linkfun:function (mu)  
  ..$ linkinv:function (eta)  
  ..$ mu.eta :function (eta)  
 $ 1/mu^2  :List of 3
  ..$ linkfun:function (mu)  
  ..$ linkinv:function (eta)  
  ..$ mu.eta :function (eta)  
 $ inverse :List of 3
  ..$ linkfun:function (mu)  
  ..$ linkinv:function (eta)  
  ..$ mu.eta :function (eta)  
> v <- sapply(linkO, function(L) sapply(L, function(F) F((0:4)/4)),
+             simplify = "array")
> stopifnot(exprs = {
+     is.numeric(v)
+     identical(dim(v), c(5L, sum(other), length(links)))
+     identical(dimnames(v)[[2]], fns[other])
+     ## check that all functions are monotone (incr. _or_ decr.) <==>
+     ## signs of differences are constant <==> var(*) == 0
+     apply(v, 2:3, function(f) var(sign(diff(f))) == 0)
+ })
> 
> ## Could further check  [for 'okLinks' of given families]:
> ##	<family>(          "<linkname>")  ==
> ##      <family>(make.link("<linkname>"))
> 
> 
> ## <family>$aic() vs logLik() vs AIC() -- for Gamma:
> # From example(glm) :
> clotting <- data.frame(
+     u    = c( 5, 10,15,20,30,40,60,80,100),
+     lot1 = c(118,58,42,35,27,25,21,19,18),
+     lot2 = c(69, 35,26,21,18,16,13,12,12))
> summary(fm1 <- glm(lot1 ~ log(u), data = clotting, family = Gamma))

Call:
glm(formula = lot1 ~ log(u), family = Gamma, data = clotting)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.04008  -0.03756  -0.02637   0.02905   0.08641  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.0165544  0.0009275  -17.85 4.28e-07 ***
log(u)       0.0153431  0.0004150   36.98 2.75e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Gamma family taken to be 0.002446059)

    Null deviance: 3.51283  on 8  degrees of freedom
Residual deviance: 0.01673  on 7  degrees of freedom
AIC: 37.99

Number of Fisher Scoring iterations: 3

> summary(fm2 <- glm(lot2 ~ log(u), data = clotting, family = Gamma))

Call:
glm(formula = lot2 ~ log(u), family = Gamma, data = clotting)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.05574  -0.02925   0.01030   0.01714   0.06371  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.0239085  0.0013265  -18.02 4.00e-07 ***
log(u)       0.0235992  0.0005768   40.91 1.36e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Gamma family taken to be 0.001813354)

    Null deviance: 3.118557  on 8  degrees of freedom
Residual deviance: 0.012672  on 7  degrees of freedom
AIC: 27.032

Number of Fisher Scoring iterations: 3

> 
> hasDisp <- 1 # have dispersion (here, but not e.g., for binomial, poisson)
> for(fm in list(fm1, fm2)) {
+     print(ll <- logLik(fm))
+     p <- attr(ll, "df")
+     A0 <- AIC(fm)
+     A1 <- -2*c(ll) + 2*p
+     aic.v <- fm$family$aic(y  = fm$y, mu = fitted(fm),
+                            wt = weights(fm), dev= deviance(fm))
+     stopifnot(p == (p. <- length(coef(fm))) + hasDisp,
+               all.equal(-2*c(ll) + 2*hasDisp, aic.v)) # <fam>$aic() = -2 * loglik + 2s
+     A2 <- aic.v + 2 * p.
+     stopifnot(exprs = {
+         all.equal(A0, A1)
+         all.equal(A1, A2)
+         all.equal(A1, fm$aic)
+     })
+ }
'log Lik.' -15.99496 (df=3)
'log Lik.' -10.51608 (df=3)
> 
> 
> 
> proc.time()
   user  system elapsed 
  0.250   0.051   0.277