# lmer, glmer and nlmer plus methods and utilities
### Utilities for parsing the mixed model formula
findbars <- function(term)
### Return the pairs of expressions that separated by vertical bars
{
if (is.name(term) || !is.language(term)) return(NULL)
if (term[[1]] == as.name("(")) return(findbars(term[[2]]))
if (!is.call(term)) stop("term must be of class call")
if (term[[1]] == as.name('|')) return(term)
if (length(term) == 2) return(findbars(term[[2]]))
c(findbars(term[[2]]), findbars(term[[3]]))
}
nobars <- function(term)
### Return the formula omitting the pairs of expressions that are
### separated by vertical bars
{
if (!('|' %in% all.names(term))) return(term)
if (is.call(term) && term[[1]] == as.name('|')) return(NULL)
if (length(term) == 2) {
nb <- nobars(term[[2]])
if (is.null(nb)) return(NULL)
term[[2]] <- nb
return(term)
}
nb2 <- nobars(term[[2]])
nb3 <- nobars(term[[3]])
if (is.null(nb2)) return(nb3)
if (is.null(nb3)) return(nb2)
term[[2]] <- nb2
term[[3]] <- nb3
term
}
subbars <- function(term)
### Substitute the '+' function for the '|' function
{
if (is.name(term) || !is.language(term)) return(term)
if (length(term) == 2) {
term[[2]] <- subbars(term[[2]])
return(term)
}
stopifnot(length(term) >= 3)
if (is.call(term) && term[[1]] == as.name('|'))
term[[1]] <- as.name('+')
for (j in 2:length(term)) term[[j]] <- subbars(term[[j]])
term
}
subnms <- function(term, nlist)
### Substitute any names from nlist in term with 1
{
if (!is.language(term)) return(term)
if (is.name(term)) {
if (any(unlist(lapply(nlist, get("=="), term)))) return(1)
return(term)
}
stopifnot(length(term) >= 2)
for (j in 2:length(term)) term[[j]] <- subnms(term[[j]], nlist)
term
}
slashTerms <- function(x)
### Return the list of '/'-separated terms in an expression that
### contains slashes
{
if (!("/" %in% all.names(x))) return(x)
if (x[[1]] != as.name("/"))
stop("unparseable formula for grouping factor")
list(slashTerms(x[[2]]), slashTerms(x[[3]]))
}
makeInteraction <- function(x)
### from a list of length 2 return recursive interaction terms
{
if (length(x) < 2) return(x)
trm1 <- makeInteraction(x[[1]])
trm11 <- if(is.list(trm1)) trm1[[1]] else trm1
list(substitute(foo:bar, list(foo=x[[2]], bar = trm11)), trm1)
}
expandSlash <- function(bb)
### expand any slashes in the grouping factors returned by findbars
{
if (!is.list(bb)) return(expandSlash(list(bb)))
## I really do mean lapply(unlist(... - unlist returns a
## flattened list in this case
unlist(lapply(bb, function(x) {
if (length(x) > 2 && is.list(trms <- slashTerms(x[[3]])))
return(lapply(unlist(makeInteraction(trms)),
function(trm) substitute(foo|bar,
list(foo = x[[2]],
bar = trm))))
x
}))
}
### Utilities used in lmer, glmer and nlmer
createCm <- function(A, s)
### Create the nonzero pattern for the sparse matrix Cm from A.
### ncol(A) is s * ncol(Cm). The s groups of ncol(Cm) consecutive
### columns in A are overlaid to produce Cm.
{
stopifnot(is(A, "dgCMatrix"))
s <- as.integer(s)[1]
if (s == 1L) return(A)
if ((nc <- ncol(A)) %% s)
stop(gettextf("ncol(A) = %d is not a multiple of s = %d",
nc, s))
ncC <- as.integer(nc / s)
TA <- as(A, "TsparseMatrix")
as(new("dgTMatrix", Dim = c(nrow(A), ncC),
i = TA@i, j = as.integer(TA@j %% ncC), x = TA@x),
"CsparseMatrix")
}
### FIXME: somehow the environment of the mf formula does not have
### .globalEnv in its parent list. example(Mmmec, package = "mlmRev")
### used to have a formula of ~ offset(log(expected)) + ... and the
### offset function was not found in eval(mf, parent.frame(2))
lmerFrames <- function(mc, formula, contrasts, vnms = character(0))
### Create the model frame, X, Y, wts, offset and terms
### mc - matched call of calling function
### formula - two-sided formula
### contrasts - contrasts argument
### vnms - names of variables to be included in the model frame
{
mf <- mc
m <- match(c("data", "subset", "weights", "na.action", "offset"),
names(mf), 0)
mf <- mf[c(1, m)]
## The model formula for evaluation of the model frame. It looks
## like a linear model formula but includes any random effects
## terms and any names of parameters used in a nonlinear mixed model.
frame.form <- subbars(formula) # substitute `+' for `|'
if (length(vnms) > 0) # add the variables names for nlmer
frame.form[[3]] <-
substitute(foo + bar,
list(foo = parse(text = paste(vnms, collapse = ' + '))[[1]],
bar = frame.form[[3]]))
## The model formula for the fixed-effects terms only.
fixed.form <- nobars(formula) # remove any terms with `|'
if (!inherits(fixed.form, "formula"))
## RHS is empty - use `y ~ 1'
fixed.form <- as.formula(substitute(foo ~ 1, list(foo = fixed.form)))
## attach the correct environment
environment(fixed.form) <- environment(frame.form) <- environment(formula)
## evaluate a model frame
mf$formula <- frame.form
mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
fe <- mf # save a copy of the call
mf <- eval(mf, parent.frame(2))
## evaluate the terms for the fixed-effects only (used in anova)
fe$formula <- fixed.form
fe <- eval(fe, parent.frame(2)) # allow model.frame to update them
## response vector
Y <- model.response(mf, "any")
## avoid problems with 1D arrays, but keep names
if(length(dim(Y)) == 1) {
nm <- rownames(Y)
dim(Y) <- NULL
if(!is.null(nm)) names(Y) <- nm
}
mt <- attr(fe, "terms")
## Extract X checking for a null model. This check shouldn't be
## needed because an empty formula is changed to ~ 1 but it can't hurt.
X <- if (!is.empty.model(mt))
model.matrix(mt, mf, contrasts) else matrix(,NROW(Y),0)
storage.mode(X) <- "double" # when ncol(X) == 0, X is logical
fixef <- numeric(ncol(X))
names(fixef) <- colnames(X)
dimnames(X) <- NULL
## Extract the weights and offset. For S4 classes we want the
## `not used' condition to be numeric(0) instead of NULL
wts <- model.weights(mf); if (is.null(wts)) wts <- numeric(0)
off <- model.offset(mf); if (is.null(off)) off <- numeric(0)
## check weights and offset
if (any(wts <= 0))
stop(gettextf("negative weights or weights of zero are not allowed"))
if(length(off) && length(off) != NROW(Y))
stop(gettextf("number of offsets is %d should equal %d (number of observations)",
length(off), NROW(Y)))
## remove the terms attribute from mf
attr(mf, "terms") <- mt
list(Y = Y, X = X, wts = as.double(wts), off = as.double(off), mf = mf, fixef = fixef)
}
##' Is f1 nested within f2?
##'
##' Does every level of f1 occur in conjunction with exactly one level
##' of f2? The function is based on converting a triplet sparse matrix
##' to a compressed column-oriented form in which the nesting can be
##' quickly evaluated.
##'
##' @param f1 factor 1
##' @param f2 factor 2
##' @return TRUE if factor 1 is nested within factor 2
isNested <- function(f1, f2)
{
f1 <- as.factor(f1)
f2 <- as.factor(f2)
stopifnot(length(f1) == length(f2))
sm <- as(new("ngTMatrix",
i = as.integer(f2) - 1L,
j = as.integer(f1) - 1L,
Dim = c(length(levels(f2)),
length(levels(f1)))),
"CsparseMatrix")
all(diff(sm@p) < 2)
}
isREML <- function(x, ...) UseMethod("isREML")
isLMM <- function(x, ...) UseMethod("isLMM")
isNLMM <- function(x, ...) UseMethod("isNLMM")
isGLMM <- function(x, ...) UseMethod("isGLMM")
##' @S3method isREML mer
isREML.mer <- function(x, ...) as.logical(x@dims["REML"])
##' @S3method isGLMM mer
isGLMM.mer <- function(x,...) {
length(x@muEta) > 0
## or: is(x@resp,"glmResp")
}
##' @S3method isNLMM mer
isNLMM.mer <- function(x,...) {
## or: is(x@resp,"nlsResp")
!isLMM.mer(x) & !isGLMM.mer(x)
}
##' @S3method isLMM mer
isLMM.mer <- function(x,...) as.logical(x@dims["LMM"])
## or: is(x@resp,"lmerResp") ?
##' dimsNames and devNames are in the package's namespace rather than
##' in the function lmerFactorList because the function sparseRasch
##' needs to access them.
dimsNames <- c("nt", "n", "p", "q", "s", "np", "LMM", "REML",
"fTyp", "lTyp", "vTyp", "nest", "useSc", "nAGQ",
"verb", "mxit", "mxfn", "cvg")
dimsDefault <- list(s = 1L, # identity mechanistic model
mxit= 300L, # maximum number of iterations
mxfn= 900L, # maximum number of function evaluations
verb= 0L, # no verbose output
np= 0L, # number of parameters in ST
LMM= 0L, # not a linear mixed model
REML= 0L, # glmer and nlmer don't use REML
fTyp= 2L, # default family is "gaussian"
lTyp= 5L, # default link is "identity"
vTyp= 1L, # default variance function is "constant"
useSc= 1L, # default is to use the scale parameter
nAGQ= 1L, # default is Laplace
cvg = 0L) # no optimization yet attempted
devNames <- c("ML", "REML", "ldL2", "ldRX2", "sigmaML",
"sigmaREML", "pwrss", "disc", "usqr", "wrss",
"dev", "llik", "NULLdev")
##' Create model matrices from r.e. terms.
##'
##' Create the list of model matrices from the random-effects terms in
##' the formula and the model frame.
##'
##' @param formula model formula
##' @param fr: list with '$mf': model frame; '$X': .. matrix
##' @param rmInt logical scalar - should the `(Intercept)` column
##' be removed before creating Zt
##' @param drop logical scalar indicating if elements with numeric
##' value 0 should be dropped from the sparse model matrices
##'
##' @return a list with components named \code{"trms"}, \code{"fl"}
##' and \code{"dims"}
lmerFactorList <- function(formula, fr, rmInt, drop)
{
mf <- fr$mf
## record dimensions and algorithm settings
## create factor list for the random effects
bars <- expandSlash(findbars(formula[[3]]))
if (!length(bars)) stop("No random effects terms specified in formula")
names(bars) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
fl <- lapply(bars,
function(x)
{
ff <- eval(substitute(as.factor(fac)[,drop = TRUE],
list(fac = x[[3]])), mf)
im <- as(ff, "sparseMatrix") # transpose of indicators
## Could well be that we should rather check earlier .. :
if(!isTRUE(validObject(im, test=TRUE)))
stop("invalid conditioning factor in random effect: ", format(x[[3]]))
mm <- model.matrix(eval(substitute(~ expr, # model matrix
list(expr = x[[2]]))),
mf)
if (rmInt) {
if (is.na(icol <- match("(Intercept)", colnames(mm)))) break
if (ncol(mm) < 2)
stop("lhs of a random-effects term cannot be an intercept only")
mm <- mm[ , -icol , drop = FALSE]
}
ans <- list(f = ff,
A = do.call(rBind,
lapply(seq_len(ncol(mm)), function(j) im)),
Zt = do.call(rBind,
lapply(seq_len(ncol(mm)),
function(j) {im@x <- mm[,j]; im})),
ST = matrix(0, ncol(mm), ncol(mm),
dimnames = list(colnames(mm), colnames(mm))))
if (drop) {
## This is only used for nlmer models.
## Need to do something more complicated for A
## here. Essentially you need to create a copy
## of im for each column of mm, im@x <- mm[,j],
## create the appropriate number of copies,
## prepend matrices of zeros, then rBind and drop0.
ans$A@x <- rep(0, length(ans$A@x))
ans$Zt <- drop0(ans$Zt)
}
ans
})
dd <-
VecFromNames(dimsNames, "integer",
c(list(n = nrow(mf), p = ncol(fr$X), nt = length(fl),
q = sum(sapply(fl, function(el) nrow(el$Zt)))),
dimsDefault))
## order terms by decreasing number of levels in the factor but don't
## change the order if this is already true
nlev <- sapply(fl, function(el) length(levels(el$f)))
## determine the number of random effects at this point
if (any(diff(nlev)) > 0) fl <- fl[rev(order(nlev))]
## separate the terms from the factor list
trms <- lapply(fl, "[", -1)
names(trms) <- NULL
fl <- lapply(fl, "[[", "f")
attr(fl, "assign") <- seq_along(fl)
## check for repeated factors
fnms <- names(fl)
if (length(fnms) > length(ufn <- unique(fnms))) {
## check that the lengths of the number of levels coincide
fl <- fl[match(ufn, fnms)]
attr(fl, "assign") <- match(fnms, ufn)
}
names(fl) <- ufn
## check for nesting of factors
dd["nest"] <- all(sapply(seq_along(fl)[-1],
function(i) isNested(fl[[i-1]], fl[[i]])))
list(trms = trms, fl = fl, dims = dd)
}
checkSTform <- function(ST, STnew)
### Check that the 'STnew' argument matches the form of ST.
{
stopifnot(is.list(STnew), length(STnew) == length(ST),
all.equal(names(ST), names(STnew)))
lapply(seq_along(STnew), function (i)
stopifnot(class(STnew[[i]]) == class(ST[[i]]),
all.equal(dim(STnew[[i]]), dim(ST[[i]]))))
all(unlist(lapply(STnew, function(m) all(diag(m) > 0))))
}
lmerControl <- function(msVerbose = getOption("verbose"),
maxIter = 300L, maxFN = 900L)
### Control parameters for lmer, glmer and nlmer
{
stopifnot(maxIter >= 0, maxFN >= 0)
list(
maxIter = as.integer(maxIter),
maxFN = as.integer(maxFN),
msVerbose = as.integer(msVerbose))# "integer" on purpose
}
##' Generate a named vector of the given mode.
##' NB: If \code{defaults} contains more than one entry of a given name,
##' the *last* one wins
VecFromNames <- function(nms, mode = "numeric", defaults = list())
{
ans <- vector(mode = mode, length = length(nms))
names(ans) <- nms
ans[] <- NA
if ((nd <- length(defaults <- as.list(defaults))) > 0) {
if (length(dnms <- names(defaults)) < nd)
stop("defaults must be a named list")
stopifnot(all(dnms %in% nms))
ans[dnms] <- as(unlist(defaults), mode)
}
ans
}
mkZt <- function(FL, start, s = 1L)
### Create the standard versions of flist, Zt, Gp, ST, A, Cm,
### Cx, and L. Update dd.
{
dd <- FL$dims
fl <- FL$fl
asgn <- attr(fl, "assign")
trms <- FL$trms
ST <- lapply(trms, `[[`, "ST")
Ztl <- lapply(trms, `[[`, "Zt")
Zt <- do.call(rBind, Ztl)
Zt@Dimnames <- vector("list", 2)
Gp <- c(0L, cumsum(vapply(Ztl, nrow, 1L, USE.NAMES=FALSE)))
.Call("mer_ST_initialize", ST, Gp, Zt)
A <- do.call(rBind, lapply(trms, `[[`, "A"))
rm(Ztl, FL) # because they could be large
nc <- sapply(ST, ncol) # of columns in els of ST
Cm <- createCm(A, s)
L <- .Call("mer_create_L", Cm)
if (s < 2) Cm <- new("dgCMatrix")
if (!is.null(start) && checkSTform(ST, start)) ST <- start
nvc <- sapply(nc, function (qi) (qi * (qi + 1))/2) # no. of var. comp.
### FIXME: Check number of variance components versus number of
### levels in the factor for each term. Warn or stop as appropriate
dd["np"] <- as.integer(sum(nvc)) # number of parameters in optimization
dev <- VecFromNames(devNames, "numeric")
fl <- do.call(data.frame, c(fl, check.names = FALSE))
attr(fl, "assign") <- asgn
list(Gp = Gp, ST = ST, A = A, Cm = Cm, L = L, Zt = Zt,
dd = dd, dev = dev, flist = fl)
}
famNms <- c("binomial", "gaussian", "Gamma", "inverse.gaussian",
"poisson")
linkNms <- c("logit", "probit", "cauchit", "cloglog", "identity",
"log", "sqrt", "1/mu^2", "inverse")
varNms <- c("constant", "mu(1-mu)", "mu", "mu^2", "mu^3")
famType <- function(family)
{
if (!(fTyp <- match(family$family, famNms, nomatch = 0)))
stop(gettextf("unknown GLM family: %s",
sQuote(family$family), domain = "R-lme4"))
if (!(lTyp <- match(family$link, linkNms, nomatch = 0)))
stop(gettextf("unknown link: %s",
sQuote(family$link), domain = "R-lme4"))
vNam <- switch(fTyp,
"mu(1-mu)", # binomial
"constant", # gaussian
"mu^2", # Gamma
"mu^3", # inverse.gaussian
"mu") # poisson
if (!(vTyp <- match(vNam, varNms, nomatch = 0)))
stop(gettextf("unknown GLM family: %s",
sQuote(family$family), domain = "R-lme4"))
c(fTyp = fTyp, lTyp = lTyp, vTyp = vTyp)
}
convergenceMessage <- function(cvg)
### Create the convergence message
{
msg <- switch(as.character(cvg),
"3" = "X-convergence (3)",
"4" = "relative convergence (4)",
"5" = "both X-convergence and relative convergence (5)",
"6" = "absolute function convergence (6)",
"7" = "singular convergence (7)",
"8" = "false convergence (8)",
"9" = "function evaluation limit reached without convergence (9)",
"10" = "iteration limit reached without convergence (9)",
"14" = "storage has been allocated (?) (14)",
"15" = "LIV too small (15)",
"16" = "LV too small (16)",
"63" = "fn cannot be computed at initial par (63)",
"65" = "gr cannot be computed at initial par (65)")
if (is.null(msg))
msg <- paste("See PORT documentation. Code (", cvg, ")", sep = "")
msg
}
#### Extractors specific to mixed-effects models
coef.mer <- function(object, ...)
{
if (length(list(...)))
warning(paste('arguments named "',
paste(names(list(...)), collapse = ", "),
'" ignored', sep = ''))
fef <- data.frame(rbind(fixef(object)), check.names = FALSE)
ref <- ranef(object)
## check for variables in RE but missing from FE, fill in zeros in FE accordingly
refnames <- unlist(lapply(ref,colnames))
nmiss <- length(missnames <- setdiff(refnames,names(fef)))
if (nmiss >0) {
fillvars <- setNames(data.frame(rbind(rep(0,nmiss))),missnames)
fef <- cbind(fillvars,fef)
}
val <- lapply(ref, function(x) fef[rep(1, nrow(x)),,drop = FALSE])
for (i in seq(a = val)) {
refi <- ref[[i]]
row.names(val[[i]]) <- row.names(refi)
nmsi <- colnames(refi)
if (!all(nmsi %in% names(fef)))
stop("unable to align random and fixed effects")
for (nm in nmsi) val[[i]][[nm]] <- val[[i]][[nm]] + refi[,nm]
}
class(val) <- "coef.mer"
val
}
setMethod("coef", signature(object = "cpglmm"), coef.mer)
setAs("cpglmm", "dtCMatrix", function(from)
### Extract the L matrix
as(from@L, "sparseMatrix"))
setMethod("fixef", signature(object = "cpglmm"),
function(object, ...)
### Extract the fixed effects
object@fixef)
##' Create a list of lists from multiple parallel lists
##' @param A a list
##' @param ... other, parallel lists
##' @return a list of lists
plist <- function(A, ...)
{
dots <- list(...)
stopifnot(is.list(A), all(sapply(dots, is.list)),
all(sapply(dots, length) == length(A)))
dots <- c(list(A), dots)
ans <- A
for (i in seq_along(A)) ans[[i]] <- lapply(dots, "[[", i)
ans
}
##' Extract the random effects.
##'
##' Extract the conditional modes, which for a linear mixed model are
##' also the conditional means, of the random effects, given the
##' observed responses. These also depend on the model parameters.
##'
##' @param object an object that inherits from the \code{\linkS4class{mer}} class
##' @param postVar logical scalar - should the posterior variance be returned
##' @param drop logical scalar - drop dimensions of single extent
##' @param whichel - vector of names of factors for which to return results
##' @return a named list of arrays or vectors, aligned to the factor list
setMethod("ranef", signature(object = "cpglmm"),
function(object, postVar = FALSE, drop = FALSE, whichel = names(wt), ...)
{
rr <- object@ranef
## nt is the number of terms, cn is the list of column names
nt <- length(cn <- lapply(object@ST, colnames))
lterm <- lapply(plist(reinds(object@Gp), cn),
function(el) {
cni <- el[[2]]
matrix(rr[ el[[1]] ], ncol = length(cni),
dimnames = list(NULL, cni))
})
wt <- whichterms(object)
ans <- lapply(plist(wt, object@flist),
function(el) {
ans <- do.call(cbind, lterm[ el[[1]] ])
rownames(ans) <- levels(el[[2]])
data.frame(ans, check.names = FALSE)
})
## Process whichel
stopifnot(is(whichel, "character"))
whchL <- names(wt) %in% whichel
ans <- ans[whchL]
if (postVar) {
pV <- .Call("mer_postVar", object, whchL)
for (i in seq_along(ans))
attr(ans[[i]], "postVar") <- pV[[i]]
}
if (drop)
ans <- lapply(ans, function(el)
{
if (ncol(el) > 1) return(el)
pv <- drop(attr(el, "postVar"))
el <- drop(as.matrix(el))
if (!is.null(pv))
attr(el, "postVar") <- pv
el
})
class(ans) <- "ranef.mer"
ans
})
print.ranef.mer <- function(x, ...) print(unclass(x), ...)
print.coef.mer <- function(x, ...) print(unclass(x), ...)
setGeneric("sigma", function(object, ...) standardGeneric("sigma"))
setMethod("sigma", signature(object = "cpglmm"),
function (object, ...) {
dd <- object@dims
if(!dd[["useSc"]]) return(1)
object@deviance[[if(dd[["REML"]]) "sigmaREML" else "sigmaML"]]
})
#### Methods for standard extractors for fitted models
setMethod("anova", signature(object = "cpglmm"),
function(object, ...)
{
mCall <- match.call(expand.dots = TRUE)
dots <- list(...)
modp <- if (length(dots))
sapply(dots, is, "cpglmm") | sapply(dots, is, "lm") else logical(0)
if (any(modp)) { # multiple models - form table
opts <- dots[!modp]
mods <- c(list(object), dots[modp])
names(mods) <- sapply(as.list(mCall)[c(FALSE, TRUE, modp)],
as.character)
mods <- mods[order(sapply(lapply(mods, logLik, REML = FALSE),
attr, "df"))]
calls <- lapply(mods, slot, "call")
data <- lapply(calls, "[[", "data")
if (any(data != data[[1]]))
stop("all models must be fit to the same data object")
header <- paste("Data:", data[[1]])
subset <- lapply(calls, "[[", "subset")
if (any(subset != subset[[1]]))
stop("all models must use the same subset")
if (!is.null(subset[[1]]))
header <-
c(header, paste("Subset", deparse(subset[[1]]), sep = ": "))
llks <- lapply(mods, logLik, REML = FALSE)
Df <- sapply(llks, attr, "df")
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,
"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))
attr(val, "heading") <-
c(header, "Models:",
paste(rep(names(mods), times = unlist(lapply(lapply(lapply(calls,
"[[", "formula"), deparse), length))),
unlist(lapply(lapply(calls, "[[", "formula"), deparse)),
sep = ": "))
return(val)
}
else { ## ------ single model ---------------------
if (length(object@muEta))
stop("single argument anova for GLMMs not yet implemented")
if (length(object@V))
stop("single argument anova for NLMMs not yet implemented")
p <- object@dims[["p"]]
ss <- (.Call("mer_update_projection", object)[[2]])^2
names(ss) <- names(object@fixef)
asgn <- attr(object@X, "assign")
terms <- terms(object)
nmeffects <- attr(terms, "term.labels")
if ("(Intercept)" %in% names(ss))
nmeffects <- c("(Intercept)", nmeffects)
ss <- unlist(lapply(split(ss, asgn), sum))
df <- unlist(lapply(split(asgn, asgn), length))
## dfr <- unlist(lapply(split(dfr, asgn), function(x) x[1]))
ms <- ss/df
f <- ms/(sigma(object)^2)
## P <- pf(f, df, dfr, lower.tail = FALSE)
## table <- data.frame(df, ss, ms, dfr, f, P)
table <- data.frame(df, ss, ms, f)
dimnames(table) <-
list(nmeffects,
## c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))
c("Df", "Sum Sq", "Mean Sq", "F value"))
if ("(Intercept)" %in% nmeffects)
table <- table[-match("(Intercept)", nmeffects), ]
attr(table, "heading") <- "Analysis of Variance Table"
class(table) <- c("anova", "data.frame")
table
}
})
setMethod("fitted", signature(object = "cpglmm"),
function(object, ...)
napredict(attr(object@frame, "na.action"), object@mu))
setMethod("residuals", signature(object = "cpglmm"),
function(object, ...)
napredict(attr(object@frame, "na.action"), object@resid))
setMethod("resid", signature(object = "cpglmm"),
function(object, ...)
napredict(attr(object@frame, "na.action"), object@resid))
### Show and print methods and utilities for them
formatVC <- function(varc, digits = max(3, getOption("digits") - 2))
### "format()" the 'VarCorr' matrix of the random effects -- for show()ing
{
sc <- unname(attr(varc, "sc"))
recorr <- lapply(varc, attr, "correlation")
reStdDev <- c(lapply(varc, attr, "stddev"), list(Residual = sc))
reLens <- unlist(c(lapply(reStdDev, length)))
nr <- sum(reLens)
reMat <- array('', c(nr, 4),
list(rep.int('', nr),
c("Groups", "Name", "Variance", "Std.Dev.")))
reMat[1+cumsum(reLens)-reLens, 1] <- names(reLens)
reMat[,2] <- c(unlist(lapply(varc, colnames)), "")
reMat[,3] <- format(unlist(reStdDev)^2, digits = digits)
reMat[,4] <- format(unlist(reStdDev), digits = digits)
if (any(reLens > 1)) {
maxlen <- max(reLens)
corr <-
do.call("rBind",
lapply(recorr,
function(x, maxlen) {
x <- as(x, "matrix")
cc <- format(round(x, 3), nsmall = 3)
cc[!lower.tri(cc)] <- ""
nr <- dim(cc)[1]
if (nr >= maxlen) return(cc)
cbind(cc, matrix("", nr, maxlen-nr))
}, maxlen))
colnames(corr) <- c("Corr", rep.int("", maxlen - 1))
cbind(reMat, rBind(corr, rep.int("", ncol(corr))))
} else reMat
}
BlockDiagonal <- function(lst)
{
stopifnot(is(lst, "list"))
lst <- lapply(lapply(lst, as, Class = "generalMatrix"),
as, Class = "TsparseMatrix")
isSquare <- function(x) nrow(x) == ncol(x)
stopifnot(all(sapply(lst, isSquare)),
all(sapply(lst, is, class2 = "dMatrix")))
if ((nl <- length(lst)) == 1) return(lst[[1]])
offsets <- c(0L, cumsum(sapply(lst, ncol)))
new("dgTMatrix", Dim = rep.int(offsets[nl + 1], 2),
i = unlist(lapply(1:nl, function(i) lst[[i]]@i + offsets[i])),
j = unlist(lapply(1:nl, function(i) lst[[i]]@j + offsets[i])),
x = unlist(lapply(lst, slot, "x")))
}
abbrvNms <- function(gnm, cnms)
### Abbreviate names of columns in grouping factors
### gnm - group name
### cnms - column names
{
ans <- paste(abbreviate(gnm), abbreviate(cnms), sep = '.')
if (length(cnms) > 1) {
anms <- lapply(cnms, abbreviate, minlength = 3)
nmmat <- outer(anms, anms, paste, sep = '.')
ans <- c(ans, paste(abbreviate(gnm, minlength = 3),
nmmat[upper.tri(nmmat)], sep = '.'))
}
ans
}
ST2Omega <- function(ST)
### Temporary function to convert the ST representation of the
### relative variance-covariance matrix returned by lmer into the
### Omega representation required by lmer
{
if (nrow(ST) == 1) return(as(1/ST^2, "dpoMatrix"))
dd <- diag(ST)
T <- as(ST, "dtrMatrix")
T@diag <- "U"
crossprod(solve(T)/dd)
}
## Utilities for the fitted mer object
slotsz <- function(obj)
rev(sort(sapply(slotNames(obj), function(s) object.size(slot(obj, s)))))
slotApply <- function(object, f, ..., simplify = FALSE) {
.localFun <- function(what, ...) f(slot(object, what), ...)
sapply(slotNames(object), .localFun, ..., simplify = simplify)
}
yfrm <- function(fm)
{
stopifnot(is(fm, "cpglmm"))
snr <- slotApply(fm, function(x)
{
if (is(x, "matrix") ||
is(x, "data.frame") ||
is(x, "numeric")) return (NROW(x))
0
}, simplify = TRUE)
snr <- snr[snr > 0 & !(names(snr) %in%
c("Gp", "dims", "deviance", "frame", "flist", "X"))]
fr <- cbind(fm@frame, fm@flist[1:NROW(fm@frame), !(names(fm@flist) %in%
names(fm@frame))])
n <- NROW(fr)
if (NROW(fm@X) == n)
fr <- cbind(fr, X = fm@X, Xbeta = fm@X %*% fm@fixef,
Zb = crossprod(fm@Zt, fm@ranef)@x)
do.call(cbind, c(list(fr), sapply(names(which(snr == NROW(fr))),
slot, object = fm, simplify = FALSE)))
}
##' Find terms associated with grouping factor names.
##' Determine the random-effects associated with particular grouping
##' factors.
##' @param fm a fitted model object of S4 class "mer"
##' @param fnm one or more grouping factor names, as a character vector
##' @return a list of indices of terms
##' @keywords models
##' @export
##' @examples
##' fm1 <- lmer(strength ~ (1|batch) + (1|sample), Pastes)
##' whichterms(fm1)
whichterms <- function(fm, fnm = names(fm@flist))
{
stopifnot(is(fm, "cpglmm"), is.character(fnm))
fl <- fm@flist
asgn <- attr(fl, "assign")
fnms <- names(fl)
stopifnot(all(fnm %in% fnms))
if (is.null(names(fnm))) names(fnm) <- fnm
lapply(fnm, function(nm) which(asgn == match(nm, fnms)))
}
##' Random-effects indices by term
##' Returns a list of indices into the ranef vector by random-effects
##' terms.
##' @param Gp the Gp slot from an mer object
##' @return a list of random-effects indices
##' @keywords models
reinds <- function(Gp)
{
lens <- diff(Gp)
lapply(seq_along(lens), function(i) Gp[i] + seq_len(lens[i]))
}
##' Random-effects indices associated with grouping factor names
##' Determine the random-effects indices with particular grouping
##' factors.
##' @param fm a fitted model object of S4 class "mer"
##' @param fnm one or more grouping factor names, as a character vector
##' @return a list of indices of terms
##' @keywords models
##' @export
##' @examples
##' fm1 <- lmer(strength ~ (1|batch) + (1|sample), Pastes)
##' whichreind(fm1)
whichreind <- function(fm, fnm = names(fm@flist))
lapply(whichterms(fm, fnm),
function (ind) unlist(reinds(fm@Gp)[ind]))
