R Under development (unstable) (2022-03-19 r81942) -- "Unsuffered Consequences"
Copyright (C) 2022 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)
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)
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.535 0.082 0.593