R Under development (unstable) (2018-10-14 r75439) -- "Unsuffered Consequences" Copyright (C) 2018 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]: > ## ( "") == > ## (make.link("")) > > > ## $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)) # $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.261 0.040 0.278