top of page

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

(c) 2017,  Joseph M. Hilbe, Rafael S. de Souza and Emille E. O. Ishida

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

Code 8.19 Random Intercept negative binomial data in R

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

library(MASS)

â€‹

N <- 2000                                                         # 10 groups, each with 200 observations
NGroups <- 10

â€‹

x1 <- runif(N)
x2 <- runif(N)

Groups <- rep(1:10, each = 200)

a <- rnorm(NGroups, mean = 0, sd = 0.5)
eta <- 1 + 0.2 * x1 - 0.75 * x2 + a[Groups]

mu <- exp(eta)
y <- rnegbin(mu, theta=2)

â€‹

nbri <- data.frame(
y = y,
x1 = x1,
x2 = x2,
Groups = Groups,
RE = a[Groups]
)

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

â€‹

â€‹

â€‹

â€‹

Code 8.20 Random intercept negative binomial mixed model in R

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

library(gamlss.mx)

nbrani <- gamlssNP(y ~ x1 + x2,
data = nbri,
random = ~ 1 | Groups,
family = NBI,
mixture = "gq", k = 20)

â€‹

summary(nbrani)

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

Anchor 1

Output on screen:

â€‹

Family: "NBI Mixture with NO"

Call: gamlssNP(formula = y ~ x1 + x2, random = ~1 | Groups, family = NBI,

data = nbri, mixture = "gq", k = 20)

Fitting method: RS()

------------------------------------------------------------------

Mu Coefficients:

â€‹

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

(Intercept)       1.09115                0.06212       17.565           <2e-16 ***

x1                    0.14289                0.08189        1.745               0.081 .

x2                   -0.92644                0.08209     -11.286             <2e-16 ***

z                      0.42583                 0.02226      19.134             <2e-16 ***

---

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

â€‹

------------------------------------------------------------------

Sigma Coefficients:

â€‹

Estimate         Std. Er        t value           Pr(>|t|)

(Intercept)      -0.69134            0.07339     -9.421          <2e-16 ***

---

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

â€‹

------------------------------------------------------------------

No. of observations in the fit: 8000

Degrees of Freedom for the fit: 5

Residual Deg. of Freedom: 1995

at cycle: 1

Global Deviance: 7036.691

AIC: 7046.691

SBC: 7074.696

******************************************************************

bottom of page