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)

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

 

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 link function: log

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 link function: log

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

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