HSI
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
******************************************************************