inst/tinytest/test-binomial.R

In brglm2: Bias Reduction in Generalized Linear Models

library("brglm")
data("lizards", package = "brglm2")
links <- lapply(c("logit", "probit", "cloglog", "cauchit"), make.link)

tol <- 1e-08
for (l in seq_along(links)) {
    expect_warning(
        lizardsBRlegacy <- brglm(cbind(grahami, opalinus) ~ height + diameter +
                                     light + time, family = binomial(links[[l]]),
                                 data = lizards, method = "brglm.fit",
                                 br.epsilon = 1e-10, br.maxit = 1000)
    )

    expect_warning(
        lizardsBR <- glm(cbind(grahami, opalinus) ~ height + diameter +
                             light + time, family = binomial(links[[l]]), data=lizards,
                         method = "brglmFit", epsilon = 1e-10, maxit = 1000)
    )
    ## glm with brglm.fit method and brglm_0 return the same coefficients for the lizards when link is links[[l]]$name)
    expect_equal(coef(lizardsBR), coef(lizardsBRlegacy), tolerance = tol)
}


## Performance comparisons BR versus ML
## link1 <- "cauchit"
## system.time(replicate(100, {lizardsBR <- glm(cbind(grahami, opalinus) ~ height + diameter + light + time, family = binomial(link1), data=lizards, method = "brglmFit", epsilon = 1e-10, maxit = 1000, type = "ML")}))

## system.time(replicate(100, {lizardsBR <- glm(cbind(grahami, opalinus) ~ height + diameter + light + time, family = binomial(link1), data=lizards, method = "glm.fit", epsilon = 1e-10, maxit = 1000)}))

## link1 <- "cauchit"
## system.time(replicate(100, {lizardsBR <- glm(cbind(grahami, opalinus) ~ height + diameter + light + time, family = binomial(link1), data=lizards, method = "brglmFit", epsilon = 1e-10, maxit = 1000)}))

## system.time(replicate(100, {lizardsBR <- brglm(cbind(grahami, opalinus) ~ height + diameter + light + time, family = binomial(link1), data=lizards, br.epsilon = 1e-10, br.maxit = 1000)}))