Synthetic negative binomial data and model in R

From: Bayesian Models for Astrophysical Data, Cambridge Univ. Press

you are kindly asked to include the complete citation if you used this material in a publication

Code 6.9 Synthetic negative binomial data and model in R

====================================================

library(MASS)

set.seed(141)
nobs <- 2500

x1 <- rbinom(nobs,size=1, prob=0.6)
x2 <- runif(nobs)
xb <- 1 + 2.0*x1 - 1.5*x2

a <- 3.3

theta <- 0.303 # 1/a

exb <- exp(xb)

nby <- rnegbin(n=nobs, mu=exb, theta=theta)
negbml <- data.frame(nby, x1, x2)

nb2 <- glm.nb(nby ~ x1 + x2, data=negbml)

summary(nb2)

====================================================

Output on screen:

Call:

glm.nb(formula = nby ~ x1 + x2, data = negbml, init.theta = 0.295684269,

link = log)

Deviance Residuals:

Min 1Q Median 3Q Max

-1.58697 -1.12176 -0.78956 0.06661 3.05031

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) 0.98613 0.08845 11.15 <2e-16 ***

x1 2.03995 0.08081 25.24 <2e-16 ***

x2 -1.61330 0.13544 -11.91 <2e-16 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.2957) family taken to be 1)

Null deviance: 3030.9 on 2499 degrees of freedom

Residual deviance: 2355.5 on 2497 degrees of freedom

AIC: 11553

Number of Fisher Scoring iterations: 1

Theta: 0.2957

Std. Err.: 0.0108

2 x log-likelihood: -11544.6130

HSI