| |
| 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 |
