##' @importFrom MASS negative.binomial
##' @importFrom MASS theta.ml
## ==> user should use getME(object, "glmer.nb.theta")
getNBdisp <- function(object) environment(object@resp$family$aic)[[".Theta"]]
## Hidden, originally used at least once, well tested (!), but
if(FALSE) # not needed anymore currently
getNBdisp.fam <- function(familyString)
as.numeric(sub(".*([-+]?\\d+(\\.\\d*)?([Ee][-+]?\\d+)?){1}.*",
"\\1", familyString))
## Package "constants" {only on depending the glmResp definition in ./AllClass.R}:
glmResp.f.nms <- names(glmResp$fields())
glmNB.to.change <- setdiff(glmResp.f.nms, c("Ptr", "family"))
##' setNBdisp(object,theta) :=
##' NB-object with changed [DISP]ersion parameter 'theta' (and all that entails)
setNBdisp <- function(object,theta) {
## assign(".Theta",theta,envir=environment(object@resp$family$aic))
rr <- object@resp
newresp <- do.call(glmResp$new,
c(lapply(setNames(nm=glmNB.to.change), rr$field),
list(family = MASS::negative.binomial(theta=theta))))
newresp$setOffset(rr$offset)
newresp$updateMu(rr$eta - rr$offset)
object@resp <- newresp
object
}
refitNB <- function(object, theta, control = NULL) {
refit(setNBdisp(object, theta), control = control)
}
##' @title Optimize over the Negative Binomial Parameter Theta
##' @param object a "glmerMod" object, updated from poisson to negative.binomial()
##' @param interval and
##' @param tol are both passed to \code{\link{optimize}()}.
##' @param verbose ## show our own progress
##' @param control passed to \code{\link{refit}()}
##' @return the last fit, an object like 'object'
optTheta <- function(object,
interval = c(-5,5),
tol = .Machine$double.eps^0.25,
verbose = FALSE,
control = NULL)
{
lastfit <- object
it <- 0L
NBfun <- function(t) {
## Kluge to retain last value and evaluation count {good enough for ..}
f_refitNB <- factory(refitNB,types=c("message","warning"))
lastfit <<- f_refitNB(lastfit, theta = exp(t),
control = control)
dev <- -2*logLik(lastfit)
it <<- it+1L
if (verbose)
cat(sprintf("%2d: th=%#15.10g, dev=%#14.8f, beta[1]=%#14.8f\n",
it, exp(t), dev, lastfit@beta[1]))
dev
}
optval <- optimize(NBfun, interval=interval, tol=tol)
stopifnot(all.equal(optval$minimum, log(getNBdisp(lastfit))))
## FIXME: return eval count info somewhere else? MM: new slot there, why not?
attr(lastfit,"nevals") <- it
## fix up the 'th' expression, replacing it by the real number,
## so effects:::mer.to.glm() can eval() it:
lastfit@call$family[["theta"]] <- exp(optval$minimum)
## output warnings/messages from last fit
for (m in c("warning","message")) {
if (!is.null(x <- attr(lastfit,paste0("factory-",m)))) {
for (i in x) {
get(m,pos="package:base")(i)
}
}
}
lastfit
}
## use MASS machinery to estimate theta from residuals
est_theta <- function(object, limit = 20,
eps = .Machine$double.eps^0.25, trace = 0)
{
Y <- model.response(model.frame(object))
## may have NA values if na.exclude was used ...
mu <- na.omit(fitted(object))
theta.ml(Y, mu, weights = object@resp$weights,
limit = limit, eps = eps, trace = trace)
}
## -------> ../man/glmer.Rd
##' glmer() for Negative Binomial
##' @param ... formula, data, etc: the arguments for
##' \code{\link{glmer}(..)} (apart from \code{family}!).
glmer.nb <- function(..., 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))
{
## ?? E() is a placeholder; [-1] removes it
## left with a 'headless' list of elements in ...
dotE <- as.list(substitute(E(...))[-1])
## nE <- names(dotE <- as.list(substitute(E(...))[-1]))
## i <- match(c("formula",""), nE)
## i.frml <- i[!is.na(i)][[1]] # the index of "formula" in '...'
## dots <- list(...)
mc <- match.call()
if (is.null(th <- initCtrl$theta)) {
mc[[1]] <- quote(lme4::glmer)
mc$family <- quote(stats::poisson)
mc$verbose <- (verbose>=2)
mc$nb.control <- NULL
## ** FIXME: specifically add check.conv.singular="ignore"?
## suppress other warnings unless explicitly specified?
g0 <- suppressMessages(
eval(mc, parent.frame(1L))
)
th <- est_theta(g0, limit = initCtrl$limit,
eps = initCtrl$eps, trace = initCtrl$trace)
## using format() on purpose, influenced by options(digits = *) :
if(verbose) cat("th := est_theta(glmer(..)) =", format(th))
}
mc$initCtrl <- NULL ## clear to prevent infinite recursion
## in initCtrl$theta reference above ...
mc$family <- bquote(MASS::negative.binomial(theta=.(th)))
## ** see FIXME above
g1 <- suppressMessages(
eval(mc, parent.frame(1L))
)
if(verbose) cat(" --> dev.= -2*logLik(.) =", format(-2*logLik(g1)),"\n")
## fix the 'data' part (only now!)
if("data" %in% names(g1@call)) {
if (!is.null(dotE[["data"]])) {
g1@call[["data"]] <- dotE[["data"]]
}
} else
warning("no 'data = *' in glmer.nb() call ... Not much is guaranteed")
other.args <- c("verbose","control")
for (a in other.args) {
if (a %in% names(g1@call)) {
g1@call[[a]] <- dotE[[a]]
}
}
## FIXME: optTheta should also work by modifying mc directly,
## then re-evaluating, not via refit() ...
control <- eval.parent(g1@call$control)
res <- optTheta(g1, interval=interval, tol=tol, verbose=verbose,
control = c(control, nb.control))
res
}
## do we want to facilitate profiling on theta??
## save evaluations used in optimize() fit?
## ('memoise'?)
## Again, I think that a reference class object would be a better approach.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.