The R Project SVN
R
(root)/branches/R-ucrt/src/library/stats/tests/glm.Rout.save – Rev 81322
Rev 81048 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
R Under development (unstable) (2021-10-13 r81047) -- "Unsuffered Consequences"Copyright (C) 2021 The R Foundation for Statistical ComputingPlatform: 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 converge2: 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 validetalogit mu eta eta etaprobit mu eta eta etacauchit mu eta eta etacloglog mu eta eta etaidentity mu eta eta etalog mu eta eta etasqrt mu eta eta eta1/mu^2 mu eta eta etainverse 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.08641Coefficients: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 freedomResidual deviance: 0.01673 on 7 degrees of freedomAIC: 37.99Number 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.06371Coefficients: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 freedomResidual deviance: 0.012672 on 7 degrees of freedomAIC: 27.032Number 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 elapsed0.287 0.044 0.312