## this file contains internal functions taken from lme4 v1.1-12
## which are altered as little as possible and preceded by 'lme4_'
lme4_safeDeparse <- function(x, collapse = " ") {
paste(deparse(x, 500L), collapse = collapse)
}
lme4_barnames <- function(bars) {
vapply(bars, function(x) lme4_safeDeparse(x[[3]]), "")
}
lme4_wmsg <- function(n, cmp.val, allow.n, msg1 = "", msg2 = "", msg3 = "") {
if (allow.n) {
unident <- n < cmp.val
cmp <- "<"
rstr <- ""
} else {
unident <- n <= cmp.val
cmp <- "<="
rstr <- " and the residual variance (or scale parameter)"
}
wstr <- sprintf("%s (=%d) %s %s (=%d)%s; the random-effects parameters%s are probably unidentifiable",
msg1, n, cmp, msg2, cmp.val, msg3, rstr)
list(unident = unident, wstr = wstr)
}
lme4_RHSForm <- function(form, as.form = FALSE) {
rhsf <- form[[length(form)]]
if (as.form) return(stats::reformulate(deparse(rhsf)))
return(rhsf)
}
`lme4_RHSForm<-` <- function(formula, value) {
formula[[length(formula)]] <- value
formula
}
lme4_getFixedFormula <- function(form) {
lme4_RHSForm(form) <- lme4::nobars(lme4_RHSForm(form))
form
}
lme4_reOnly <- function(f, response = FALSE) {
if (response && length(f) == 3) {
response <- f[[2]]
} else {
response <- NULL
}
stats::reformulate(paste0("(", vapply(lme4::findbars(f), lme4_safeDeparse,
""), ")"), response = response)
}
lme4_reFormHack <- function(re.form, ReForm, REForm, REform) {
warnDeprec <- function(name) {
warning(gettextf("'%s' is deprecated; use '%s' instead", name, "re.form"),
call. = FALSE, domain = NA)
}
if (!missing(ReForm)) {
warnDeprec("ReForm")
return(ReForm)
}
if (!missing(REForm)) {
warnDeprec("REForm")
return(REForm)
}
if (!missing(REform)) {
warnDeprec("REform")
return(REform)
}
re.form
}
lme4_isGLMM <- function(x, ...) {
as.logical(x@devcomp$dims[["GLMM"]])
}
lme4_isNLMM <- function(x, ...) {
as.logical(x@devcomp$dims[["NLMM"]])
}
lme4_isLMM <- function(x, ...) {
!lme4_isGLMM(x) && !lme4_isNLMM(x)
}
lme4_checkArgs <- function(type, ...) {
l... <- list(...)
if (isTRUE(l...[["sparseX"]])) {
warning("sparseX = TRUE has no effect at present", call. = FALSE)
}
if (length(l... <- list(...))) {
if (!is.null(l...[["family"]])) {
warning("calling lmer with family() is deprecated: please use glmer() instead",
call. = FALSE)
type <- "glmer"
}
if (!is.null(l...[["method"]])) {
msg <- paste("Argument", sQuote("method"), "is deprecated.")
if (type == "lmer") {
msg <- paste(msg, "Use the REML argument to specify ML or REML estimation.")
} else if (type == "glmer") {
msg <- paste(msg, "Use the nAGQ argument to specify Laplace (nAGQ=1) or adaptive",
"Gauss-Hermite quadrature (nAGQ>1). PQL is no longer available.")
}
warning(msg, call. = FALSE)
l... <- l...[names(l...) != "method"]
}
if (length(l...)) {
warning("extra argument(s) ", paste(sQuote(names(l...)),
collapse = ", "), " disregarded", call. = FALSE)
}
}
}
lme4_checkCtrlLevels <- function(cstr, val, smallOK = FALSE) {
bvals <- c("stop", "warning", "ignore")
if (smallOK) bvals <- outer(bvals, c("", "Small"), paste0)
if (!is.null(val) && !val %in% bvals) {
stop("invalid control level ", sQuote(val), " in ", cstr,
": valid options are {", paste(sapply(bvals, sQuote),
collapse = ","), "}")
}
invisible(NULL)
}
lme4_missDataFun <- function(d) {
ff <- sys.frames()
ex <- substitute(d)
ii <- rev(seq_along(ff))
for (i in ii) {
ex <- eval(substitute(substitute(x, env = sys.frames()[[n]]),
env = list(x = ex, n = i)))
}
return(is.symbol(ex) && !exists(deparse(ex)))
}
lme4_checkFormulaData <- function(formula, data, checkLHS = TRUE,
debug = FALSE) {
nonexist.data <- lme4_missDataFun(data)
wd <- tryCatch(eval(data), error = identity)
if (wrong.data <- inherits(wd, "simpleError")) {
wrong.data.msg <- wd$message
}
bad.data <- nonexist.data || wrong.data
if (bad.data || debug) {
varex <- function(v, env) exists(v, envir = env, inherits = FALSE)
allvars <- all.vars(stats::as.formula(formula))
allvarex <- function(env, vvec = allvars) all(vapply(vvec, varex, NA, env))
}
if (bad.data) {
if (allvarex(environment(formula))) {
stop("bad 'data', but variables found in environment of formula: ",
"try specifying 'formula' as a formula rather ",
"than a string in the original model", call. = FALSE)
} else {
if (nonexist.data) {
stop("'data' not found, and some variables missing from formula environment",
call. = FALSE)
}
else {
stop("bad 'data': ", wrong.data.msg)
}
}
} else {
if (is.null(data)) {
if (!is.null(ee <- environment(formula))) {
denv <- ee
} else {
denv <- parent.frame(2L)
}
} else {
denv <- list2env(data)
}
}
if (debug) {
cat("Debugging parent frames in checkFormulaData:\n")
glEnv <- 1L
while (!identical(parent.frame(glEnv), .GlobalEnv)) {
glEnv <- glEnv + 1L
}
for (i in 1:glEnv) {
OK <- allvarex(parent.frame(i))
cat("vars exist in parent frame ", i)
if (i == glEnv) cat(" (global)")
cat(" ", OK, "\n")
}
cat("vars exist in env of formula ", allvarex(denv), "\n")
}
stopifnot(!checkLHS || length(as.formula(formula, env = denv)) == 3)
return(denv)
}
lme4_doCheck <- function(x) {
is.character(x) && !any(x == "ignore")
}
lme4_checkNlevels <- function(flist, n, ctrl, allow.n = FALSE) {
stopifnot(is.list(ctrl), is.numeric(n))
nlevelVec <- unlist(lapply(flist, function(x) nlevels(droplevels(x))))
cstr <- "check.nlev.gtr.1"
lme4_checkCtrlLevels(cstr, cc <- ctrl[[cstr]])
if (lme4_doCheck(cc) && any(nlevelVec < 2)) {
wstr <- "grouping factors must have > 1 sampled level"
switch(cc, warning = warning(wstr, call. = FALSE), stop = stop(wstr,
call. = FALSE), stop(gettextf("unknown check level for '%s'",
cstr), domain = NA))
} else {
wstr <- character()
}
cstr <- "check.nobs.vs.nlev"
lme4_checkCtrlLevels(cstr, cc <- ctrl[[cstr]])
if (lme4_doCheck(cc) && any(if (allow.n) nlevelVec > n else nlevelVec >= n)) {
wst2 <- gettextf("number of levels of each grouping factor must be %s number of observations",
if (allow.n) "<=" else "<")
switch(cc, warning = warning(wst2, call. = FALSE), stop = stop(wst2,
call. = FALSE))
} else {
wst2 <- character()
}
cstr <- "check.nlev.gtreq.5"
lme4_checkCtrlLevels(cstr, cc <- ctrl[[cstr]])
if (lme4_doCheck(cc) && any(nlevelVec < 5)) {
wst3 <- "grouping factors with < 5 sampled levels may give unreliable estimates"
switch(cc, warning = warning(wst3, call. = FALSE), stop = stop(wst3,
call. = FALSE), stop(gettextf("unknown check level for '%s'",
cstr), domain = NA))
} else {
wst3 <- character()
}
c(wstr, wst2, wst3)
}
lme4_checkZdims <- function(Ztlist, n, ctrl, allow.n = FALSE) {
stopifnot(is.list(Ztlist), is.numeric(n))
cstr <- "check.nobs.vs.nRE"
lme4_checkCtrlLevels(cstr, cc <- ctrl[[cstr]])
term.names <- names(Ztlist)
rows <- vapply(Ztlist, nrow, 1L)
cols <- vapply(Ztlist, ncol, 1L)
stopifnot(all(cols == n))
if (lme4_doCheck(cc)) {
unique(unlist(lapply(seq_along(Ztlist), function(i) {
ww <- lme4_wmsg(cols[i], rows[i], allow.n, "number of observations",
"number of random effects", sprintf(" for term (%s)", term.names[i]))
if (ww$unident) {
switch(cc, warning = warning(ww$wstr, call. = FALSE),
stop = stop(ww$wstr, call. = FALSE), stop(gettextf("unknown check level for '%s'",
cstr), domain = NA))
ww$wstr
} else {
character()
}
})))
} else {
character()
}
}
#' @importFrom Matrix rankMatrix t
lme4_checkZrank <- function(Zt, n, ctrl, nonSmall = 1e+06, allow.n = FALSE) {
stopifnot(is.list(ctrl), is.numeric(n), is.numeric(nonSmall))
cstr <- "check.nobs.vs.rankZ"
if (lme4_doCheck(cc <- ctrl[[cstr]])) {
lme4_checkCtrlLevels(cstr, cc, smallOK = TRUE)
d <- dim(Zt)
doTr <- d[1L] < d[2L]
if (!(grepl("Small", cc) && prod(d) > nonSmall)) {
rankZ <- Matrix::rankMatrix(if (doTr) Matrix::t(Zt) else Zt,
method = "qr", sval = numeric(min(d)))
ww <- lme4_wmsg(n, rankZ, allow.n, "number of observations", "rank(Z)")
if (ww$unident) {
switch(cc, warningSmall = , warning = warning(ww$wstr,
call. = FALSE), stopSmall = , stop = stop(ww$wstr,
call. = FALSE), stop(gettextf("unknown check level for '%s'",
cstr), domain = NA))
ww$wstr
} else {
character()
}
} else {
character()
}
} else {
character()
}
}
lme4_chkRank.drop.cols <- function(X, kind, tol = 1e-07, method = "qr.R") {
stopifnot(is.matrix(X))
kinds <- eval(formals(lme4::lmerControl)[["check.rankX"]])
if (!kind %in% kinds) stop(gettextf("undefined option for 'kind': %s", kind))
if (kind == "ignore") return(X)
p <- ncol(X)
if (kind == "stop.deficient") {
if ((rX <- Matrix::rankMatrix(X, tol = tol, method = method)) < p) {
stop(gettextf(sub("\n +", "\n", "the fixed-effects model matrix is column rank deficient (rank(X) = %d < %d = p);\n the fixed effects will be jointly unidentifiable"), rX, p), call. = FALSE)
}
} else {
qr.X <- qr(X, tol = tol, LAPACK = FALSE)
rnkX <- qr.X$rank
if (rnkX == p) return(X)
msg <- sprintf(ngettext(p - rnkX, "fixed-effect model matrix is rank deficient so dropping %d column / coefficient",
"fixed-effect model matrix is rank deficient so dropping %d columns / coefficients"),
p - rnkX)
if (kind != "silent.drop.cols") {
if (kind == "warn+drop.cols") {
warning(msg, domain = NA)
} else {
message(msg, domain = NA)
}
}
contr <- attr(X, "contrasts")
asgn <- attr(X, "assign")
keep <- qr.X$pivot[seq_len(rnkX)]
dropped.names <- colnames(X[, -keep, drop = FALSE])
X <- X[, keep, drop = FALSE]
if (Matrix::rankMatrix(X, tol = tol, method = method) < ncol(X)) {
stop(gettextf("Dropping columns failed to produce full column rank design matrix"),
call. = FALSE)
}
if (!is.null(contr)) attr(X, "contrasts") <- contr
if (!is.null(asgn)) attr(X, "assign") <- asgn[keep]
attr(X, "msgRankdrop") <- msg
attr(X, "col.dropped") <- stats::setNames(qr.X$pivot[(rnkX + 1L):p],
dropped.names)
}
X
}
lme4_checkScaleX <- function (X, kind = "warning", tol = 1000) {
kinds <- eval(formals(lme4::lmerControl)[["check.scaleX"]])
if (!kind %in% kinds) stop(gettextf("unknown check-scale option: %s", kind))
if (is.null(kind) || kind == "ignore") return(X)
cont.cols <- apply(X, 2, function(z) !all(z %in% c(0, 1)))
col.sd <- apply(X[, cont.cols, drop = FALSE], 2L, sd)
sdcomp <- outer(col.sd, col.sd, "/")
logcomp <- abs(log(sdcomp[lower.tri(sdcomp)]))
logsd <- abs(log(col.sd))
if (any(c(logcomp, logsd) > log(tol))) {
wmsg <- "Some predictor variables are on very different scales:"
if (kind %in% c("warning", "stop")) {
wmsg <- paste(wmsg, "consider rescaling")
switch(kind, warning = warning(wmsg, call. = FALSE),
stop = stop(wmsg, call. = FALSE))
} else {
wmsg <- paste(wmsg, "auto-rescaled (results NOT adjusted)")
X[, cont.cols] <- sweep(X[, cont.cols, drop = FALSE], 2, col.sd, "/")
attr(X, "scaled:scale") <- stats::setNames(col.sd, colnames(X)[cont.cols])
if (kind == "warn+rescale") warning(wmsg, call. = FALSE)
}
} else {
wmsg <- character()
}
structure(X, msgScaleX = wmsg)
}
lme4_checkHess <- function(H, tol, hesstype = "") {
res <- list(code = numeric(0), messages = list())
evd <- tryCatch(eigen(H, symmetric = TRUE, only.values = TRUE)$values,
error = function(e) e)
if (inherits(evd, "error")) {
res$code <- -6L
res$messages <- gettextf("Problem with Hessian check (infinite or missing values?)")
} else {
negative <- sum(evd < -tol)
if (negative) {
res$code <- -3L
res$messages <- gettextf(paste("Model failed to converge:",
"degenerate", hesstype, "Hessian with %d negative eigenvalues"),
negative)
} else {
zero <- sum(abs(evd) < tol)
if (zero || inherits(tryCatch(chol(H), error = function(e) e), "error")) {
res$code <- -4L
res$messages <- paste(hesstype, "Hessian is numerically singular: parameters are not uniquely determined")
} else {
res$cond.H <- max(evd) / min(evd)
if (max(evd) * tol > 1) {
res$code <- c(res$code, 2L)
res$messages <- c(res$messages, paste("Model is nearly unidentifiable: ",
"very large eigenvalue", "\n - Rescale variables?", sep = ""))
}
if ((min(evd)/max(evd)) < tol) {
res$code <- c(res$code, 3L)
if (!5L %in% res$code) {
res$messages <- c(res$messages, paste("Model is nearly unidentifiable: ",
"large eigenvalue ratio", "\n - Rescale variables?", sep = ""))
}
}
}
}
}
if (length(res$code) == 0) res$code <- 0
res
}
lme4_checkConv <- function(derivs, coefs, ctrl, lbound, debug = FALSE) {
if (is.null(derivs)) return(NULL)
if (anyNA(derivs$gradient)) {
return(list(code = -5L, messages = gettextf("Gradient contains NAs")))
}
ntheta <- length(lbound)
res <- list()
ccl <- ctrl[[cstr <- "check.conv.grad"]]
lme4_checkCtrlLevels(cstr, cc <- ccl[["action"]])
wstr <- NULL
if (lme4_doCheck(cc)) {
scgrad <- tryCatch(with(derivs, solve(chol(Hessian),
gradient)), error = function(e) e)
if (inherits(scgrad, "error")) {
wstr <- "unable to evaluate scaled gradient"
res$code <- -1L
} else {
mingrad <- pmin(abs(scgrad), abs(derivs$gradient))
maxmingrad <- max(mingrad)
if (maxmingrad > ccl$tol) {
w <- which.max(maxmingrad)
res$code <- -1L
wstr <- gettextf("Model failed to converge with max|grad| = %g (tol = %g, component %d)",
maxmingrad, ccl$tol, w)
}
}
if (!is.null(wstr)) {
res$messages <- wstr
switch(cc, warning = warning(wstr), stop = stop(wstr),
stop(gettextf("unknown check level for '%s'", cstr), domain = NA))
}
if (!is.null(ccl$relTol) && (max.rel.grad <- max(abs(derivs$gradient/coefs))) > ccl$relTol) {
res$code <- -2L
wstr <- gettextf("Model failed to converge with max|relative grad| = %g (tol = %g)",
max.rel.grad, ccl$relTol)
res$messages <- wstr
switch(cc, warning = warning(wstr), stop = stop(wstr),
stop(gettextf("unknown check level for '%s'", cstr), domain = NA))
}
}
ccl <- ctrl[[cstr <- "check.conv.singular"]]
lme4_checkCtrlLevels(cstr, cc <- ccl[["action"]])
if (lme4_doCheck(cc)) {
bcoefs <- seq(ntheta)[lbound == 0]
if (any(coefs[bcoefs] < ccl$tol)) {
wstr <- "singular fit"
res$messages <- c(res$messages, wstr)
switch(cc, warning = warning(wstr), stop = stop(wstr),
stop(gettextf("unknown check level for '%s'", cstr), domain = NA))
}
}
ccl <- ctrl[[cstr <- "check.conv.hess"]]
lme4_checkCtrlLevels(cstr, cc <- ccl[["action"]])
if (lme4_doCheck(cc)) {
if (length(coefs) > ntheta) {
H.beta <- derivs$Hessian[-seq(ntheta), -seq(ntheta)]
resHess <- lme4_checkHess(H.beta, ccl$tol, "fixed-effect")
if (any(resHess$code != 0)) {
res$code <- resHess$code
res$messages <- c(res$messages, resHess$messages)
wstr <- paste(resHess$messages, collapse = ";")
switch(cc, warning = warning(wstr), stop = stop(wstr),
stop(gettextf("unknown check level for '%s'", cstr), domain = NA))
}
}
resHess <- lme4_checkHess(derivs$Hessian, ccl$tol)
if (any(resHess$code != 0)) {
res$code <- resHess$code
res$messages <- c(res$messages, resHess$messages)
wstr <- paste(resHess$messages, collapse = ";")
switch(cc, warning = warning(wstr), stop = stop(wstr),
stop(gettextf("unknown check level for '%s'", cstr), domain = NA))
}
}
if (debug && length(res$messages) > 0) print(res$messages)
res
}
lme4_updateStart <- function(start, theta) {
if (is.null(start)) return(NULL)
if (is.numeric(start)) {
start <- theta
} else if (!is.null(start$theta)) {
start$theta <- theta
}
start
}
lme4_est_theta <- function(object, limit = 20, eps = .Machine$double.eps^0.25,
trace = 0) {
Y <- stats::model.response(model.frame(object))
mu <- fitted(object)
MASS::theta.ml(Y, mu, weights = object@resp$weights, limit = limit,
eps = eps, trace = trace)
}
# Not sure if this one will work
lme4_setNBdisp <- function(object, theta) {
rr <- object@resp
newresp <- do.call(lme4::glmResp$new, c(lapply(
stats::setNames(nm = c("mu", "offset", "sqrtXwt", "sqrtrwt", "weights",
"wtres", "y", "eta", "n")), rr$field),
list(family = MASS::negative.binomial(theta = theta))))
newresp$setOffset(rr$offset)
newresp$updateMu(rr$eta - rr$offset)
object@resp <- newresp
object
}
lme4_refitNB <- function(object, theta, control = NULL) {
lme4::refit(lme4_setNBdisp(object, theta), control = control)
}
lme4_getNBdisp <- function(object) {
environment(object@resp$family$aic)[[".Theta"]]
}
lme4_optTheta <- function(object, interval = c(-5, 5),
tol = .Machine$double.eps^0.25, verbose = FALSE,
control = NULL) {
lastfit <- object
it <- 0L
NBfun <- function(t) {
dev <- -2 * logLik(lastfit <<- lme4_refitNB(lastfit, theta = exp(t),
control = control))
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 <- stats::optimize(NBfun, interval = interval, tol = tol)
stopifnot(all.equal(optval$minimum, log(lme4_getNBdisp(lastfit))))
attr(lastfit, "nevals") <- it
lastfit@call$family[["theta"]] <- exp(optval$minimum)
lastfit
}
lme4_anova.merMod <- function(object, ..., refit = TRUE, model.names = NULL) {
mCall <- match.call(expand.dots = TRUE)
dots <- list(...)
.sapply <- function(L, FUN, ...) unlist(lapply(L, FUN, ...))
modp <- (as.logical(vapply(dots, is, NA, "merMod")) | as.logical(vapply(dots,
is, NA, "lm")))
if (any(modp)) {
mods <- c(list(object), dots[modp])
nobs.vec <- vapply(mods, nobs, 1L)
if (var(nobs.vec) > 0) {
stop("models were not all fitted to the same size of dataset")
}
if (is.null(mNms <- model.names)) {
mNms <- vapply(as.list(mCall)[c(FALSE, TRUE, modp)], lme4_safeDeparse, "")
}
if (any(duplicated(mNms))) {
warning("failed to find unique model names, assigning generic names")
mNms <- paste0("MODEL", seq_along(mNms))
}
if (length(mNms) != length(mods)) {
stop("model names vector and model list have different lengths")
}
names(mods) <- sub("@env$", "", mNms)
models.reml <- vapply(mods,
function(x) is(x, "merMod") && lme4::isREML(x), NA)
models.GHQ <- vapply(mods, function(x) is(x, "glmerMod") &&
lme4::getME(x, "devcomp")$dims["nAGQ"] > 1, NA)
if (any(models.GHQ) && any(vapply(mods, function(x) is(x, "glm"), NA))) {
stop("GLMMs with nAGQ>1 have log-likelihoods incommensurate with glm() objects")
}
if (refit) {
if (any(models.reml)) {
message("refitting model(s) with ML (instead of REML)")
}
mods[models.reml] <- lapply(mods[models.reml], lme4::refitML)
} else {
if (any(models.reml) && any(!models.reml)) {
warning("some models fit with REML = TRUE, some not")
}
}
llks <- lapply(mods, logLik)
ii <- order(Df <- vapply(llks, attr, FUN.VALUE = numeric(1), "df"))
mods <- mods[ii]
llks <- llks[ii]
Df <- Df[ii]
calls <- lapply(mods, getCall)
data <- lapply(calls, `[[`, "data")
if (!all(vapply(data, identical, NA, data[[1]]))) {
stop("all models must be fit to the same data object")
}
header <- paste("Data:", lme4_abbrDeparse(data[[1]]))
subset <- lapply(calls, `[[`, "subset")
if (!all(vapply(subset, identical, NA, subset[[1]]))) {
stop("all models must use the same subset")
}
if (!is.null(subset[[1]])) {
header <- c(header, paste("Subset:", lme4_abbrDeparse(subset[[1]])))
}
llk <- unlist(llks)
chisq <- 2 * pmax(0, c(NA, diff(llk)))
dfChisq <- c(NA, diff(Df))
val <- data.frame(Df = Df, AIC = .sapply(llks, AIC),
BIC = .sapply(llks, BIC), logLik = llk, deviance = -2 * llk,
Chisq = chisq, `Chi Df` = dfChisq, `Pr(>Chisq)` = pchisq(chisq,
dfChisq, lower.tail = FALSE), row.names = names(mods),
check.names = FALSE)
class(val) <- c("anova", class(val))
forms <- lapply(lapply(calls, `[[`, "formula"), deparse)
structure(val, heading = c(header, "Models:", paste(rep(names(mods),
times = lengths(forms)), unlist(forms), sep = ": ")))
} else {
dc <- lme4::getME(object, "devcomp")
X <- lme4::getME(object, "X")
stopifnot(length(asgn <- attr(X, "assign")) == dc$dims[["p"]])
ss <- as.vector(object@pp$RX() %*% object@beta)^2
names(ss) <- colnames(X)
terms <- terms(object)
nmeffects <- attr(terms, "term.labels")[unique(asgn)]
if ("(Intercept)" %in% names(ss)) {
nmeffects <- c("(Intercept)", nmeffects)
}
ss <- unlist(lapply(split(ss, asgn), sum))
stopifnot(length(ss) == length(nmeffects))
df <- lengths(split(asgn, asgn))
ms <- ss/df
f <- ms/(sigma(object)^2)
table <- data.frame(df, ss, ms, f)
dimnames(table) <- list(nmeffects, c("Df", "Sum Sq",
"Mean Sq", "F value"))
if ("(Intercept)" %in% nmeffects) {
table <- table[-match("(Intercept)", nmeffects), ]
}
structure(table, heading = "Analysis of Variance Table",
class = c("anova", "data.frame"))
}
}
lme4_abbrDeparse <- function(x, width = 60) {
r <- deparse(x, width)
if (length(r) > 1) {
paste(r[1], "...")
} else {
r
}
}
lme4_simfunList <- function() {
simfunList <- list()
simfunList$gaussian <- function(object, nsim, ftd = fitted(object),
wts = weights(object)) {
if (any(wts != 1)) warning("ignoring prior weights")
rnorm(nsim * length(ftd), ftd, sd = sigma(object))
}
simfunList$binomial <- function(object, nsim, ftd = fitted(object),
wts = weights(object)) {
n <- length(ftd)
ntot <- n * nsim
if (any(wts%%1 != 0)) {
stop("cannot simulate from non-integer prior.weights")
}
if (!is.null(m <- model.frame(object))) {
y <- model.response(m)
if (is.factor(y)) {
yy <- factor(1 + rbinom(ntot, size = 1, prob = ftd), labels = levels(y))
split(yy, rep(seq_len(nsim), each = n))
} else if (is.matrix(y) && ncol(y) == 2) {
yy <- vector("list", nsim)
for (i in seq_len(nsim)) {
Y <- rbinom(n, size = wts, prob = ftd)
YY <- cbind(Y, wts - Y)
colnames(YY) <- colnames(y)
yy[[i]] <- YY
}
yy
} else rbinom(ntot, size = wts, prob = ftd)/wts
} else rbinom(ntot, size = wts, prob = ftd)/wts
}
simfunList$poisson <- function(object, nsim, ftd = fitted(object),
wts = weights(object)) {
wts <- weights(object)
if (any(wts != 1)) warning("ignoring prior weights")
rpois(nsim * length(ftd), ftd)
}
simfunList$Gamma <- function(object, nsim, ftd = fitted(object),
wts = weights(object)) {
if (any(wts != 1)) message("using weights to scale shape parameter")
shape <- sigma(object) * wts
rgamma(nsim * length(ftd), shape = shape, rate = shape/ftd)
}
simfunList$negative.binomial <- function(object, nsim, ftd = fitted(object),
wts = weights(object)) {
if (any(wts != 1)) warning("ignoring prior weights")
theta <- lme4_getNBdisp(object)
rnbinom(nsim * length(ftd), mu = ftd, size = theta)
}
return(simfunList)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.