The R Project SVN R

Rev

Rev 59360 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 59360 Rev 74908
Line 3... Line 3...
3 
Y <- c(rep(0,35),1,2,0,6,8,16,43)
 3 
Y <- c(rep(0,35),1,2,0,6,8,16,43)
4 
beta <- 42:1
 4 
beta <- 42:1
5 
cst <- lchoose(42, beta)
 5 
cst <- lchoose(42, beta)
6 
tau <- (beta^2)/2
 6 
tau <- (beta^2)/2
7 
fit <- glm(formula = Y ~ offset(cst) + beta + tau, family = poisson)
 7 
fit <- glm(formula = Y ~ offset(cst) + beta + tau, family = poisson)
- 
 
 8 
 
- 
 
 9 
## Ensure make.link() consistency:
- 
 
 10 
linkNames <- c("logit", "probit", "cauchit", "cloglog",
- 
 
 11 
               "identity",
- 
 
 12 
               "log",  "sqrt",  "1/mu^2", "inverse")
- 
 
 13 
links <- lapply(setNames(,linkNames), make.link)
- 
 
 14 
fns <- c("linkfun", "linkinv", "mu.eta", "valideta")
- 
 
 15 
stopifnot(exprs = {
- 
 
 16 
    is.matrix(nms <- sapply(links, names)) # matching number & type
- 
 
 17 
    is.character(nms)
- 
 
 18 
    nms[,1] == nms ## all columns are the same
- 
 
 19 
    identical(setNames(,linkNames), vapply(links, `[[`, "", "name"))
- 
 
 20 
    fns %in% nms[,1]
- 
 
 21 
})
- 
 
 22 
links <- lapply(links, `[`, fns) # functions only
- 
 
 23 
stopifnot(unlist(lapply(links, function(L) vapply(L, is.function, NA))))
- 
 
 24 
## all functions having consistent arguments :
- 
 
 25 
lf <- lapply(links, function(L) lapply(L, formals))
- 
 
 26 
stopifnot(exprs = { ## all functions have 1 argument
- 
 
 27 
    unlist(lapply(lf, lengths), recursive=FALSE) == 1L
- 
 
 28 
    is.matrix(argNms <- sapply(lf, function(L) vapply(L, names, "")))
- 
 
 29 
    argNms[,1] == argNms ## all columns are the same
- 
 
 30 
})
- 
 
 31 
noquote(t(argNms))
- 
 
 32 
 
- 
 
 33 
## Calling all functions
- 
 
 34 
## 1. valideta
- 
 
 35 
stopifnot(vv <- vapply(links, function(L) L$valideta((1:3)/4), NA))
- 
 
 36 
## 2. all others
- 
 
 37 
other <- fns != "valideta"
- 
 
 38 
str(linkO <- lapply(links, function(L) L[other]))
- 
 
 39 
v <- sapply(linkO, function(L) sapply(L, function(F) F((0:4)/4)),
- 
 
 40 
            simplify = "array")
- 
 
 41 
stopifnot(exprs = {
- 
 
 42 
    is.numeric(v)
- 
 
 43 
    identical(dim(v), c(5L, sum(other), length(links)))
- 
 
 44 
    identical(dimnames(v)[[2]], fns[other])
- 
 
 45 
    ## check that all functions are monotone (incr. _or_ decr.) <==>
- 
 
 46 
    ## signs of differences are constant <==> var(*) == 0
- 
 
 47 
    apply(v, 2:3, function(f) var(sign(diff(f))) == 0)
- 
 
 48 
})
- 
 
 49 
 
- 
 
 50 
## Could further check  [for 'okLinks' of given families]:
- 
 
 51 
##	<family>(          "<linkname>")  ==
- 
 
 52 
##      <family>(make.link("<linkname>"))
- 
 
 53 
 
- 
 
 54 
 
- 
 
 55 
## <family>$aic() vs logLik() vs AIC() -- for Gamma:
- 
 
 56 
# From example(glm) :
- 
 
 57 
clotting <- data.frame(
- 
 
 58 
    u    = c( 5, 10,15,20,30,40,60,80,100),
- 
 
 59 
    lot1 = c(118,58,42,35,27,25,21,19,18),
- 
 
 60 
    lot2 = c(69, 35,26,21,18,16,13,12,12))
- 
 
 61 
summary(fm1 <- glm(lot1 ~ log(u), data = clotting, family = Gamma))
- 
 
 62 
summary(fm2 <- glm(lot2 ~ log(u), data = clotting, family = Gamma))
- 
 
 63 
 
- 
 
 64 
hasDisp <- 1 # have dispersion (here, but not e.g., for binomial, poisson)
- 
 
 65 
for(fm in list(fm1, fm2)) {
- 
 
 66 
    print(ll <- logLik(fm))
- 
 
 67 
    p <- attr(ll, "df")
- 
 
 68 
    A0 <- AIC(fm)
- 
 
 69 
    A1 <- -2*c(ll) + 2*p
- 
 
 70 
    aic.v <- fm$family$aic(y  = fm$y, mu = fitted(fm),
- 
 
 71 
                           wt = weights(fm), dev= deviance(fm))
- 
 
 72 
    stopifnot(p == (p. <- length(coef(fm))) + hasDisp,
- 
 
 73 
              all.equal(-2*c(ll) + 2*hasDisp, aic.v)) # <fam>$aic() = -2 * loglik + 2s
- 
 
 74 
    A2 <- aic.v + 2 * p.
- 
 
 75 
    stopifnot(exprs = {
- 
 
 76 
        all.equal(A0, A1)
- 
 
 77 
        all.equal(A1, A2)
- 
 
 78 
        all.equal(A1, fm$aic)
- 
 
 79 
    })
- 
 
 80 
}
- 
 
 81 
 
- 
 
 82