glmer.nb.RdFits a generalized linear mixed-effects model (GLMM) for the negative
binomial family, building on glmer, and initializing via
theta.ml from MASS.
arguments as for glmer(.) such as formula,
data, control, etc, but not family!
interval in which to start the optimization. The default is symmetric on log scale around the initially estimated theta.
tolerance for the optimization via optimize.
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.
optional list, like the output of glmerControl(),
used in refit(*, control = control.nb) during the
optimization (control, if included in ...,
will be used in the initial-stage glmer(...,family=poisson)
fit, and passed on to the later optimization stages as well)
(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
An object of class glmerMod, for which many
methods are available (e.g. methods(class="glmerMod")), see
glmer.
For historical reasons, the shape parameter of the negative
binomial and the random effects parameters in our (G)LMM models are
both called theta (\(\theta\)), but are unrelated here.
The negative binomial \(\theta\) 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 \(\theta\), 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)).
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.
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))))
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 5 10 20 20 30 35
dd$y <- rnbinom(nrow(dd), mu = mu, size = 0.5)
str(dd)
#> '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 ...
require("MASS")## and use its glm.nb() - as indeed we have zero random effect:
#> Loading required package: MASS
if (FALSE) { # \dontrun{
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
} # }