glmer.nb: Fitting Negative Binomial GLMMs

Description Usage Arguments Value Note See Also Examples

View source: R/nbinom.R

Description

Fits a generalized linear mixed-effects model (GLMM) for the negative binomial family, building on glmer, and initializing via theta.ml from MASS.

Usage

1
2
3
4
glmer.nb(..., interval = log(th) + c(-3, 3),
         tol = 5e-5, verbose = FALSE, nb.control = NULL,
         initCtrl = list(limit = 20, eps = 2*tol, trace = verbose,
                         theta = NULL))

Arguments

...

arguments as for glmer(.) such as formula, data, control, etc, but not family!

interval

interval in which to start the optimization. The default is symmetric on log scale around the initially estimated theta.

tol

tolerance for the optimization via optimize.

verbose

logical indicating how much progress information should be printed during the optimization. Use verbose = 2 (or larger) to enable verbose=TRUE in the glmer() calls.

nb.control

optional list, like glmerControl(), used in refit(*, control = control.nb) during the optimization.

initCtrl

(experimental, do not rely on this:) a list with named components as in the default, passed to theta.ml (package MASS) for the initial value of the negative binomial parameter theta. May also include a theta component, in which case the initial estimation step is skipped

Value

An object of class glmerMod, for which many methods are available (e.g. methods(class="glmerMod")), see glmer.

Note

For historical reasons, the shape parameter of the negative binomial and the random effects parameters in our (G)LMM models are both called theta (θ), but are unrelated here.

The negative binomial θ can be extracted from a fit g <- glmer.nb() by getME(g, "glmer.nb.theta").

Parts of glmer.nb() are still experimental and methods are still missing or suboptimal. In particular, there is no inference available for the dispersion parameter θ, yet.

To fit a negative binomial model with known overdispersion parameter (e.g. as part of a model comparison exercise, use glmer with the negative.binomial family from the MASS package, e.g. glmer(...,family=MASS::negative.binomial(theta=1.75)).

See Also

glmer; from package MASS, negative.binomial (which we re-export currently) and theta.ml, the latter for initialization of optimization.

The ‘Details’ of pnbinom for the definition of the negative binomial distribution.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
set.seed(101)
dd <- expand.grid(f1 = factor(1:3),
                  f2 = LETTERS[1:2], g=1:9, rep=1:15,
          KEEP.OUT.ATTRS=FALSE)
summary(mu <- 5*(-4 + with(dd, as.integer(f1) + 4*as.numeric(f2))))
dd$y <- rnbinom(nrow(dd), mu = mu, size = 0.5)
str(dd)
require("MASS")## and use its glm.nb() - as indeed we have zero random effect:
## Not run: 
m.glm <- glm.nb(y ~ f1*f2, data=dd, trace=TRUE)
summary(m.glm)
m.nb <- glmer.nb(y ~ f1*f2 + (1|g), data=dd, verbose=TRUE)
m.nb
## The neg.binomial theta parameter:
getME(m.nb, "glmer.nb.theta")
LL <- logLik(m.nb)
## mixed model has 1 additional parameter (RE variance)
stopifnot(attr(LL,"df")==attr(logLik(m.glm),"df")+1)
plot(m.nb, resid(.) ~ g)# works, as long as data 'dd' is found

## End(Not run)

Example output

Loading required package: Matrix
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      5      10      20      20      30      35 
'data.frame':	810 obs. of  5 variables:
 $ f1 : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...
 $ f2 : Factor w/ 2 levels "A","B": 1 1 1 2 2 2 1 1 1 2 ...
 $ g  : int  1 1 1 1 1 1 2 2 2 2 ...
 $ rep: int  1 1 1 1 1 1 1 1 1 1 ...
 $ y  : num  3 16 31 6 51 14 19 31 0 15 ...
Loading required package: MASS
Theta(1) = 0.479770, 2(Ls - Lm) = 941.735000

Call:
glm.nb(formula = y ~ f1 * f2, data = dd, trace = TRUE, init.theta = 0.4797700505, 
    link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0340  -1.2213  -0.5140   0.1935   2.6476  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   1.6501     0.1299  12.707  < 2e-16 ***
f12           0.7672     0.1816   4.225 2.38e-05 ***
f13           1.0115     0.1812   5.583 2.36e-08 ***
f2B           1.5124     0.1806   8.375  < 2e-16 ***
f12:f2B      -0.6151     0.2538  -2.423   0.0154 *  
f13:f2B      -0.6104     0.2534  -2.409   0.0160 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

    Null deviance: 1074.76  on 809  degrees of freedom
Residual deviance:  941.74  on 804  degrees of freedom
AIC: 5983.3

Number of Fisher Scoring iterations: 1


              Theta:  0.4798 
          Std. Err.:  0.0243 

 2 x log-likelihood:  -5969.2690 
theta.ml: iter 0 'theta = 0.470115'
theta.ml: iter1 theta =0.482275
theta.ml: iter2 theta =0.482681
theta.ml: iter3 theta =0.482681
th := est_theta(glmer(..)) = 0.4826813 --> dev.= -2*logLik(.) = 5969.284 
 1: th=   0.2377340479, dev= 6140.60808029, beta[1]=    1.65008211
 2: th=   0.9800077384, dev= 6192.26657489, beta[1]=    1.64682447
 3: th=  0.09906383105, dev= 6721.48253778, beta[1]=    1.65007943
 4: th=   0.4547081059, dev= 5970.38221149, beta[1]=    1.65008438
 5: th=   0.4605939030, dev= 5969.91396596, beta[1]=    1.65008438
 6: th=   0.4826727026, dev= 5969.28348056, beta[1]=    1.65008438
 7: th=   0.4798351618, dev= 5969.26927968, beta[1]=    1.65008438
 8: th=   0.4797599933, dev= 5969.26927267, beta[1]=    1.65008438
 9: th=   0.4797700781, dev= 5969.26927250, beta[1]=    1.65008438
10: th=   0.4797780796, dev= 5969.26927261, beta[1]=    1.65008438
11: th=   0.4797700781, dev= 5969.26927250, beta[1]=    1.65008438
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Negative Binomial(0.4798)  ( log )
Formula: y ~ f1 * f2 + (1 | g)
   Data: dd
      AIC       BIC    logLik  deviance  df.resid 
 5985.269  6022.846 -2984.635  5969.269       802 
Random effects:
 Groups Name        Std.Dev.
 g      (Intercept) 4.29e-06
Number of obs: 810, groups:  g, 9
Fixed Effects:
(Intercept)          f12          f13          f2B      f12:f2B      f13:f2B  
     1.6501       0.7672       1.0115       1.5124      -0.6151      -0.6104  
[1] 0.4797701

lme4 documentation built on Nov. 10, 2018, 5:04 p.m.