negbin1.R
In negbin1: Negative Binomial 1 Regression

### R code from vignette source 'negbin1.Rnw'

###################################################
### code chunk number 1: preliminaries
###################################################
library("negbin1")
library("VGAM")
library("gamlss")
library("MASS")
library("lmtest")
library("Formula")
library("numDeriv")
library("memisc")
options(prompt = "R> ", continue = "+   ")


###################################################
### code chunk number 2: dgp
###################################################
### data-generating process
dgp <- function(n = 10000, coef = c(0.2, 1, 0.3, 0.4)) {
  d <- data.frame(
    x1 = runif(n, -1, 1),
    x2 = runif(n, -1, 1)
    )
  d$mu <- exp(coef[1] + coef[2] * d$x1 + coef[3] * d$x2)
  d$y <- rnbinom(n, mu = d$mu, size = d$mu / exp(coef[4]))
  return(d)
}

## simulate data
set.seed(2017-05-15)
d <- dgp()


###################################################
### code chunk number 3: predictions
###################################################
mod <- negbin1(y ~ x1 + x2, data = d)
newdata <- data.frame(x1 = c(0, 0.4, 0.1, 2), x2 = c(0.1, 0, -0.1, 0.1))
predict(mod, newdata, type = "response")
predict(mod, newdata, type = "quantile", at = 0.95)
predict(mod, newdata, type = "probability", at = 0)  


###################################################
### code chunk number 4: models
###################################################
st <- rbind(system.time(negbin1(y ~ x1 + x2, data = d, grad = FALSE)),
system.time(negbin1(y ~ x1 + x2, data = d, grad = TRUE)),
system.time(vglm(y ~ x1 + x2, negbinomial(parallel = TRUE, zero = NULL),
                       data = d, trace = FALSE)),
system.time(gamlss(y ~ x1 + x2, family = NBI, data = d, trace = FALSE))
)
st <- st[, -c(4,5)]
rownames(st) <- c("negbin1(grad = FALSE)", "negbin1(grad = TRUE)", "vglm", "gamlss")
st


###################################################
### code chunk number 5: data
###################################################
data("publications", package = "negbin1")


###################################################
### code chunk number 6: poisson
###################################################
poi <- glm(articles ~ gender + married + kids + prestige + mentor,
  family = "poisson", data = publications)
summary(poi)


###################################################
### code chunk number 7: sum
###################################################
var(publications$articles) / mean(publications$articles)


###################################################
### code chunk number 8: nb1
###################################################
nb1 <- negbin1(articles ~ gender + married + kids + prestige + mentor,
   data = publications, grad = TRUE)
summary(nb1)


###################################################
### code chunk number 9: se
###################################################
cbind("Poisson" = sqrt(diag(vcov(poi))), "NB1" = sqrt(diag(vcov(nb1)))[-7])


###################################################
### code chunk number 10: nb2
###################################################
nb2 <- vglm(articles ~ gender + married + kids + prestige + mentor, negbinomial(parallel = TRUE, zero = NULL),
            data = publications, trace = TRUE)    
summary(nb2)


###################################################
### code chunk number 11: nb3
###################################################
nb3 <- gamlss(articles ~ gender + married + kids + prestige + mentor,
              family = NBII(sigma.link = "identity"), data = publications)
summary(nb3)


###################################################
### code chunk number 12: mtable (eval = FALSE)
###################################################
## library("memisc")
## mtable("Poisson" = poi, "NB1(negbin1)" = nb1, 
##        summary.stats = c("Log-likelihood", "AIC", "BIC", "N"))


###################################################
### code chunk number 13: mtable-latex
###################################################
options(factor.style = "($f):  ($l)")
toLatex(mtable("Poisson" = poi, "NB1" = nb1, summary.stats = c("Log-likelihood", "AIC", "BIC", "N")))