--- branches/trunk-lme4/R/lmer.R 2005/04/02 14:41:46 689 +++ pkg/R/lmer.R 2006/01/14 23:41:55 1165 @@ -1,226 +1,539 @@ -## Methods for lmer and for the objects that it produces +# Methods for lmer and for the objects that it produces ## Some utilities -contr.SAS <- function(n, contrasts = TRUE) -## Eliminate this function after R-2.1.0 is released - contr.treatment(n, - if (is.numeric(n) && length(n) == 1) n else length(n), - contrasts) - -Lind <- function(i,j) { - if (i < j) stop(paste("Index i=", i,"must be >= index j=", j)) - ((i - 1) * i)/2 + j +## Return the pairs of expressions separated by vertical bars +findbars <- function(term) +{ + if (is.name(term) || is.numeric(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]])) +} + +## Return the formula omitting the pairs of expressions +## that are separated by vertical bars +nobars <- function(term) +{ + 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 +} + +## Substitute the '+' function for the '|' function +subbars <- function(term) +{ + if (is.name(term) || is.numeric(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('+') + term[[2]] <- subbars(term[[2]]) + term[[3]] <- subbars(term[[3]]) + term +} + +## Expand an expression with colons to the sum of the lhs +## and the current expression +colExpand <- function(term) +{ + if (is.name(term) || is.numeric(term)) return(term) + if (length(term) == 2) { + term[[2]] <- colExpand(term[[2]]) + return(term) + } + stopifnot(length(term) == 3) + if (is.call(term) && term[[1]] == as.name(':')) { + return(substitute(A+B, list(A = term, B = colExpand(term[[2]])))) + } + term[[2]] <- colExpand(term[[2]]) + term[[3]] <- colExpand(term[[3]]) + term } -Dhalf <- function(from) { - D <- from@D - nf <- length(D) - Gp <- from@Gp - res <- array(0, rep(Gp[nf+1],2)) - for (i in 1:nf) { - DD <- D[[i]] - dd <- dim(DD) - for (k in 1:dd[3]) { - mm <- array(DD[ , , k], dd[1:2]) - base <- Gp[i] + (k - 1)*dd[1] - res[cbind(c(base + row(mm)), c(base + col(mm)))] <- c(mm) - } +abbrvNms <- function(gnm, cnms) +{ + 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 = '.')) } - res + ans } -lmerControl <- # Control parameters for lmer - function(maxIter = 50, - msMaxIter = 50, - tolerance = sqrt((.Machine$double.eps)), - niterEM = 20, - msTol = sqrt(.Machine$double.eps), +## Control parameters for lmer +lmerControl <- + function(maxIter = 200, # used in ../src/lmer.c only + tolerance = sqrt(.Machine$double.eps),# ditto + msMaxIter = 200, + ## msTol = sqrt(.Machine$double.eps), + ## FIXME: should be able to pass tolerances to nlminb() msVerbose = getOption("verbose"), - PQLmaxIt = 20, + niterEM = 15, EMverbose = getOption("verbose"), + PQLmaxIt = 30,# FIXME: unused; PQL currently uses 'maxIter' instead analyticGradient = TRUE, - analyticHessian=FALSE) + analyticHessian = FALSE # unused _FIXME_ + ) { - list(maxIter = maxIter, - msMaxIter = msMaxIter, - tolerance = tolerance, - niterEM = niterEM, - msTol = msTol, - msVerbose = msVerbose, - PQLmaxIt = PQLmaxIt, - EMverbose=EMverbose, - analyticHessian=analyticHessian, - analyticGradient=analyticGradient) + list(maxIter = as.integer(maxIter), + tolerance = as.double(tolerance), + msMaxIter = as.integer(msMaxIter), + ## msTol = as.double(msTol), + msVerbose = as.integer(msVerbose),# "integer" on purpose + niterEM = as.integer(niterEM), + EMverbose = as.logical(EMverbose), + PQLmaxIt = as.integer(PQLmaxIt), + analyticGradient = as.logical(analyticGradient), + analyticHessian = as.logical(analyticHessian)) } -setMethod("lmer", signature(formula = "formula", family = "missing"), +setMethod("coef", signature(object = "lmer"), + function(object, ...) + { + fef <- data.frame(rbind(object@fixed), check.names = FALSE) + ref <- as(ranef(object), "list") + names(ref) <- names(object@flist) + val <- lapply(ref, function(x) fef[rep(1, nrow(x)),]) + for (i in seq(a = val)) { + refi <- ref[[i]] + row.names(val[[i]]) <- row.names(refi) + if (!all(names(refi) %in% names(fef))) + stop("unable to align random and fixed effects") + val[[i]][ , names(refi)] <- val[[i]][ , names(refi)] + refi + } + new("lmer.coef", val) + }) + +## setMethod("plot", signature(x = "lmer.coef"), +## function(x, y, ...) +## { +## varying <- unique(do.call("c", +## lapply(x, function(el) +## names(el)[sapply(el, +## function(col) +## any(col != col[1]))]))) +## gf <- do.call("rbind", lapply(x, "[", j = varying)) +## gf$.grp <- factor(rep(names(x), sapply(x, nrow))) +## switch(min(length(varying), 3), +## qqmath(eval(substitute(~ x | .grp, +## list(x = as.name(varying[1])))), gf, ...), +## xyplot(eval(substitute(y ~ x | .grp, +## list(y = as.name(varying[1]), +## x = as.name(varying[2])))), gf, ...), +## splom(~ gf | .grp, ...)) +## }) + +## setMethod("plot", signature(x = "lmer.ranef"), +## function(x, y, ...) +## { +## lapply(x, function(x) { +## cn <- lapply(colnames(x), as.name) +## switch(min(ncol(x), 3), +## qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...), +## xyplot(eval(substitute(y ~ x, +## list(y = cn[[1]], +## x = cn[[2]]))), x, ...), +## splom(~ x, ...)) +## }) +## }) + +setMethod("with", signature(data = "lmer"), + function(data, expr, ...) { + dat <- eval(data@call$data) + if (!is.null(na.act <- attr(data@frame, "na.action"))) + dat <- dat[-na.act, ] + lst <- c(list(. = data), data@flist, data@frame, dat) + eval(substitute(expr), lst[unique(names(lst))]) + }) + +setMethod("terms", signature(x = "lmer"), + function(x, ...) x@terms) + +rWishart <- function(n, df, invScal) + .Call("Matrix_rWishart", n, df, invScal, PACKAGE = "Matrix") + + +setMethod("lmer", signature(formula = "formula"), function(formula, data, family, method = c("REML", "ML", "PQL", "Laplace", "AGQ"), - control = list(), + control = list(), start, subset, weights, na.action, offset, - model = TRUE, x = FALSE, y = FALSE, ...) + model = TRUE, x = FALSE, y = FALSE , ...) { - # match and check parameters - REML <- match.arg(method) == "REML" - controlvals <- do.call("lmerControl", control) - controlvals$REML <- REML + ## match and check parameters if (length(formula) < 3) stop("formula must be a two-sided formula") + cv <- do.call("lmerControl", control) - mf <- match.call() # create the model frame as frm - m <- match(c("data", "subset", "weights", "na.action", "offset"), - names(mf), 0) - mf <- mf[c(1, m)] - mf[[1]] <- as.name("model.frame") - frame.form <- subbars(formula) - environment(frame.form) <- environment(formula) + ## Must evaluate the model frame first and then fit the glm using + ## that frame. Otherwise missing values in the grouping factors + ## cause inconsistent numbers of observations. + mf <- match.call() + m <- match(c("family", "data", "subset", "weights", + "na.action", "offset"), names(mf), 0) + mf <- fe <- mf[c(1, m)] + frame.form <- subbars(formula) # substitute +' for |' + fixed.form <- nobars(formula) # remove any terms with `|' + if (inherits(fixed.form, "name")) # RHS is empty - use a constant + fixed.form <- substitute(foo ~ 1, list(foo = fixed.form)) + environment(fixed.form) <- environment(frame.form) <- environment(formula) + + ## evaluate a model frame for fixed and random effects mf$formula <- frame.form + mf$family <- NULL mf$drop.unused.levels <- TRUE + mf[[1]] <- as.name("model.frame") frm <- eval(mf, parent.frame()) - - ## grouping factors and model matrices for random effects + + ## fit a glm model to the fixed formula + fe$formula <- fixed.form + fe$subset <- NULL # subset has already been created in call to data.frame + fe$data <- frm + fe$x <- fe$model <- fe$y <- TRUE + fe[[1]] <- as.name("glm") + glm.fit <- eval(fe, parent.frame()) + x <- glm.fit$x + y <- as.double(glm.fit$y) + family <- glm.fit$family + + ## check for a linear mixed model + lmm <- family$family == "gaussian" && family$link == "identity" + if (lmm) { # linear mixed model + method <- match.arg(method) + if (method %in% c("PQL", "Laplace", "AGQ")) { + warning(paste('Argument method = "', method, + '" is not meaningful for a linear mixed model.\n', + 'Using method = "REML".\n', sep = '')) + method <- "REML" + } + } else { # generalized linear mixed model + if (missing(method)) method <- "PQL" + else { + method <- match.arg(method) + if (method == "ML") method <- "PQL" + if (method == "REML") + warning('Argument method = "REML" is not meaningful ', + 'for a generalized linear mixed model.', + '\nUsing method = "PQL".\n') + } + } + ## create factor list for the random effects bars <- findbars(formula[[3]]) - random <- - lapply(bars, - function(x) list(model.matrix(eval(substitute(~term, - list(term=x[[2]]))), - frm), - eval(substitute(as.factor(fac)[,drop = TRUE], - list(fac = x[[3]])), frm))) - names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]]))) - + names(bars) <- unlist(lapply(bars, function(x) deparse(x[[3]]))) + fl <- lapply(bars, + function(x) + eval(substitute(as.factor(fac)[,drop = TRUE], + list(fac = x[[3]])), frm)) ## order factor list by decreasing number of levels - nlev <- sapply(random, function(x) length(levels(x[[2]]))) + nlev <- sapply(fl, function(x) length(levels(x))) if (any(diff(nlev) > 0)) { - random <- random[rev(order(nlev))] + ord <- rev(order(nlev)) + bars <- bars[ord] + fl <- fl[ord] + } + ## create list of transposed model matrices for random effects + Ztl <- lapply(bars, function(x) + t(model.matrix(eval(substitute(~ expr, + list(expr = x[[2]]))), + frm))) + ## Create the mixed-effects representation (mer) object + mer <- .Call("mer_create", fl, + .Call("Zt_create", fl, Ztl, PACKAGE = "Matrix"), + x, y, method, sapply(Ztl, nrow), + c(lapply(Ztl, rownames), list(.fixed = colnames(x))), + !(family$family %in% c("binomial", "poisson")), + match.call(), family, + PACKAGE = "Matrix") + if (lmm) { + .Call("mer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix") + LMEoptimize(mer) <- cv + return(mer) } - fixed.form <- nobars(formula) - if (!inherits(fixed.form, "formula")) fixed.form <- ~ 1 # default formula - Xmat <- model.matrix(fixed.form, frm) - mmats <- c(lapply(random, "[[", 1), - .fixed = list(cbind(Xmat, .response = model.response(frm)))) - ## FIXME: Use Xfrm and Xmat to get the terms and assign - ## slots, pass them to lmer_create, then destroy them - obj <- .Call("lmer_create", lapply(random, "[[", 2), - mmats, PACKAGE = "Matrix") - slot(obj, "terms") <- attr(model.frame(fixed.form, data), "terms") - slot(obj, "assign") <- attr(Xmat, "assign") - slot(obj, "call") <- match.call() - slot(obj, "REML") <- REML - rm(Xmat) - .Call("lmer_initial", obj, PACKAGE="Matrix") - .Call("lmer_ECMEsteps", obj, - controlvals$niterEM, - controlvals$REML, - controlvals$EMverbose, - PACKAGE = "Matrix") - LMEoptimize(obj) <- controlvals - slot(obj, "residuals") <- - unname(model.response(frm) - - (slot(obj, "fitted") <- - .Call("lmer_fitted", obj, mmats, TRUE, PACKAGE = "Matrix"))) - obj + + ## The rest of the function applies to generalized linear mixed models + gVerb <- getOption("verbose") + eta <- glm.fit$linear.predictors + wts <- glm.fit$prior.weights + wtssqr <- wts * wts + offset <- glm.fit$offset + if (is.null(offset)) offset <- numeric(length(eta)) + linkinv <- quote(family$linkinv(eta)) + mu.eta <- quote(family$mu.eta(eta)) + mu <- family$linkinv(eta) + variance <- quote(family$variance(mu)) + dev.resids <- quote(family$dev.resids(y, mu, wtssqr)) + LMEopt <- get("LMEoptimize<-") + doLMEopt <- quote(LMEopt(x = mer, value = cv)) + + GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix") + .Call("glmer_PQL", GSpt, PACKAGE = "Matrix") # obtain PQL estimates + fixInd <- seq(ncol(x)) + ## pars[fixInd] == beta, pars[-fixInd] == theta + PQLpars <- c(fixef(mer), + .Call("mer_coef", mer, 2, PACKAGE = "Matrix")) + .Call("glmer_devLaplace", PQLpars, GSpt, PACKAGE = "Matrix") + ## indicator of constrained parameters + const <- c(rep(FALSE, length(fixInd)), + unlist(lapply(mer@nc[seq(along = fl)], + function(k) 1:((k*(k+1))/2) <= k) + )) + devLaplace <- function(pars) + .Call("glmer_devLaplace", pars, GSpt, PACKAGE = "Matrix") + + optimRes <- + nlminb(PQLpars, devLaplace, + lower = ifelse(const, 5e-10, -Inf), + control = list(trace = getOption("verbose"), + iter.max = cv$msMaxIter)) + .Call("glmer_finalize", GSpt, PACKAGE = "Matrix") + return(mer) + deviance <- devAGQ(PQLpars, 1) + +### FIXME: For nf == 1 change this to an AGQ evaluation. Needs +### AGQ for nc > 1 first. + fxd <- PQLpars[fixInd] + loglik <- logLik(mer) + + if (method %in% c("Laplace", "AGQ")) { + nAGQ <- 1 + if (method == "AGQ") { # determine nAGQ at PQL estimates + dev11 <- devAGQ(PQLpars, 11) + ## FIXME: Should this be an absolute or a relative tolerance? + devTol <- sqrt(.Machine$double.eps) * abs(dev11) + for (nAGQ in c(9, 7, 5, 3, 1)) + if (abs(dev11 - devAGQ(PQLpars, nAGQ - 2)) > devTol) break + nAGQ <- nAGQ + 2 + if (gVerb) + cat(paste("Using", nAGQ, "quadrature points per column\n")) + } + obj <- function(pars) + .Call("glmer_devAGQ", pars, GSpt, nAGQ, PACKAGE = "Matrix") + optimRes <- + nlminb(PQLpars, obj, + lower = ifelse(const, 5e-10, -Inf), + control = list(trace = getOption("verbose"), + iter.max = cv$msMaxIter)) + optpars <- optimRes$par + if (optimRes$convergence != 0) + warning("nlminb failed to converge") + deviance <- optimRes$objective + if (gVerb) + cat(paste("convergence message", optimRes$message, "\n")) + fxd[] <- optpars[fixInd] ## preserve the names + .Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix") + } + + .Call("glmer_finalize", GSpt, PACKAGE = "Matrix") + loglik[] <- -deviance/2 + new("lmer", mer, + frame = if (model) frm else data.frame(), + terms = glm.fit$terms, + assign = attr(glm.fit$x, "assign"), + call = match.call(), family = family, + logLik = loglik, fixed = fxd) }) -setReplaceMethod("LMEoptimize", signature(x="lmer", value="list"), + +## Extract the permutation +setAs("mer", "pMatrix", function(from) + .Call("mer_pMatrix", from, PACKAGE = "Matrix")) + +## Extract the L matrix +setAs("mer", "dtCMatrix", function(from) + .Call("mer_dtCMatrix", from, PACKAGE = "Matrix")) + +## Extract the fixed effects +setMethod("fixef", signature(object = "mer"), + function(object, ...) + .Call("mer_fixef", object, PACKAGE = "Matrix")) + +## Extract the random effects +setMethod("ranef", signature(object = "mer"), + function(object, ...) + .Call("mer_ranef", object, PACKAGE = "Matrix") + ) + +## Optimization for mer objects +setReplaceMethod("LMEoptimize", signature(x="mer", value="list"), function(x, value) { if (value$msMaxIter < 1) return(x) - st <- ccoef(x) # starting values nc <- x@nc - nc <- nc[1:(length(nc) - 2)] constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k)) - fn <- function(pars) { - ccoef(x) <- pars - deviance(x, REML = value$REML) - } + fn <- function(pars) + deviance(.Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix")) gr <- if (value$analyticGradient) function(pars) { - if (!identical(TRUE,all.equal(pars, ccoef(x)))) ccoef(x) <- pars - grad <- gradient(x, REML = value$REML, unconst = TRUE) - grad[constr] <- -grad[constr]/pars[constr] - grad - } else NULL - optimRes <- optim(st, fn, gr, - method = "L-BFGS-B", - lower = ifelse(constr, 1e-10, -Inf), - control = list(maxit = value$msMaxIter, - trace = as.integer(value$msVerbose))) + if (!isTRUE(all.equal(pars, + .Call("mer_coef", x, + 2, PACKAGE = "Matrix")))) + .Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix") + .Call("mer_gradient", x, 2, PACKAGE = "Matrix") + } + else NULL + optimRes <- nlminb(.Call("mer_coef", x, 2, PACKAGE = "Matrix"), + fn, gr, + lower = ifelse(constr, 5e-10, -Inf), + control = list(iter.max = value$msMaxIter, + trace = as.integer(value$msVerbose))) + .Call("mer_coefGets", x, optimRes$par, 2, PACKAGE = "Matrix") if (optimRes$convergence != 0) { - warning(paste("optim returned message",optimRes$message,"\n")) + warning(paste("nlminb returned message", + optimRes$message,"\n")) } - ccoef(x) <- optimRes$par return(x) }) -setMethod("ranef", signature(object = "lmer"), - function(object, accumulate = FALSE, ...) { - val <- new("lmer.ranef", - lapply(.Call("lmer_ranef", object, PACKAGE = "Matrix"), - data.frame, check.names = FALSE), - varFac = object@bVar, - stdErr = .Call("lmer_sigma", object, - object@REML, PACKAGE = "Matrix")) - if (!accumulate || length(val@varFac) == 1) return(val) - ## check for nested factors - L <- object@L - if (any(sapply(seq(a = val), function(i) length(L[[Lind(i,i)]]@i)))) - error("Require nested grouping factors to accumulate random effects") - val - }) - -setMethod("fixef", signature(object = "lmer"), +setMethod("deviance", signature(object = "mer"), function(object, ...) { - val <- .Call("lmer_fixef", object, PACKAGE = "Matrix") - val[-length(val)] + .Call("mer_factor", object, PACKAGE = "Matrix") + object@deviance[[ifelse(object@method == "REML", "REML", "ML")]] }) -setMethod("VarCorr", signature(x = "lmer"), - function(x, REML = TRUE, useScale = TRUE, ...) { - val <- .Call("lmer_variances", x, PACKAGE = "Matrix") - for (i in seq(along = val)) { - dimnames(val[[i]]) = list(x@cnames[[i]], x@cnames[[i]]) - val[[i]] = as(as(val[[i]], "pdmatrix"), "corrmatrix") - } - new("VarCorr", - scale = .Call("lmer_sigma", x, REML, PACKAGE = "Matrix"), - reSumry = val, - useScale = useScale) - }) +setMethod("mcmcsamp", signature(object = "mer"), + function(object, n = 1, verbose = FALSE, saveb = FALSE, + trans = TRUE, ...) + { + ans <- t(.Call("mer_MCMCsamp", object, saveb, n, + trans, PACKAGE = "Matrix")) + attr(ans, "mcpar") <- as.integer(c(1, n, 1)) + class(ans) <- "mcmc" + glmer <- FALSE + gnms <- names(object@flist) + cnms <- object@cnames + ff <- fixef(object) + colnms <- c(names(ff), if (glmer) character(0) else "sigma^2", + unlist(lapply(seq(along = gnms), + function(i) + abbrvNms(gnms[i],cnms[[i]])))) + if (trans) { + ## parameter type: 0 => fixed effect, 1 => variance, + ## 2 => covariance + ptyp <- c(integer(length(ff)), if (glmer) integer(0) else 1:1, + unlist(lapply(seq(along = gnms), + function(i) + { + k <- length(cnms[[i]]) + rep(1:2, c(k, (k*(k-1))/2)) + }))) + colnms[ptyp == 1] <- + paste("log(", colnms[ptyp == 1], ")", sep = "") + colnms[ptyp == 2] <- + paste("atanh(", colnms[ptyp == 2], ")", sep = "") + } + colnames(ans) <- colnms + ans + }) -setMethod("gradient", signature(x = "lmer"), - function(x, REML, unconst, ...) - .Call("lmer_gradient", x, REML, unconst, PACKAGE = "Matrix")) +setMethod("simulate", signature(object = "mer"), + function(object, nsim = 1, seed = NULL, ...) + { + if(!exists(".Random.seed", envir = .GlobalEnv)) + runif(1) # initialize the RNG if necessary + if(is.null(seed)) + RNGstate <- .Random.seed + else { + R.seed <- .Random.seed + set.seed(seed) + RNGstate <- structure(seed, kind = as.list(RNGkind())) + on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv)) + } -setMethod("summary", signature(object = "lmer"), - function(object, ...) - new("summary.lmer", object, useScale = TRUE, showCorrelation = TRUE)) + family <- object@family + if (family$family != "gaussian" || + family$link != "identity") + stop("simulation of generalized linear mixed models not yet implemented") + ## similate the linear predictors + lpred <- .Call("mer_simulate", object, nsim, PACKAGE = "Matrix") + sc <- 1 + if (object@useScale) + sc <- .Call("mer_sigma", object, object@method == "REML", + PACKAGE = "Matrix") + ## add fixed-effects contribution and per-observation noise term + lpred <- as.data.frame(lpred + drop(object@X %*% fixef(object)) + + rnorm(prod(dim(lpred)), sd = sc)) + ## save the seed + attr(lpred, "seed") <- RNGstate + lpred + }) -setMethod("show", signature(object = "lmer"), - function(object) - show(new("summary.lmer", object, useScale = TRUE, - showCorrelation = FALSE)) - ) -setMethod("show", "summary.lmer", +setMethod("show", "mer", function(object) { - fcoef <- fixef(object) + vcShow <- function(varc, useScale) + { + digits <- max(3, getOption("digits") - 2) + sc <- attr(varc, "sc") + recorr <- lapply(varc, function(el) el@factors$correlation) + reStdDev <- c(lapply(recorr, slot, "sd"), list(Residual = sc)) + reLens <- unlist(c(lapply(reStdDev, length))) + reMat <- array('', c(sum(reLens), 4), + list(rep('', sum(reLens)), + c("Groups", "Name", "Variance", "Std.Dev."))) + reMat[1+cumsum(reLens)-reLens, 1] <- names(reLens) + reMat[,2] <- c(unlist(lapply(reStdDev, names)), "") + 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("", maxlen - 1)) + reMat <- cbind(reMat, rbind(corr, rep("", ncol(corr)))) + } + if (!useScale) reMat <- reMat[-nrow(reMat),] + print(reMat, quote = FALSE) + } + + fcoef <- .Call("mer_fixef", object, PACKAGE = "Matrix") useScale <- object@useScale - corF <- as(as(vcov(object, useScale = useScale), "pdmatrix"), - "corrmatrix") + corF <- vcov(object)@factors$correlation DF <- getFixDF(object) - coefs <- cbind(fcoef, corF@stdDev, DF) - nc <- object@nc + coefs <- cbind(fcoef, corF@sd, DF) dimnames(coefs) <- list(names(fcoef), c("Estimate", "Std. Error", "DF")) - digits <- max(3, getOption("digits") - 2) - REML <- length(object@REML) > 0 && object@REML[1] - llik <- logLik(object) + digits <- max(3, getOption("digits") - 2) + REML <- object@method == "REML" + llik <- logLik(object, REML) dev <- object@deviance - + devc <- object@devComp + rdig <- 5 - cat("Linear mixed-effects model fit by ") - cat(ifelse(object@REML, "REML\n", "maximum likelihood\n") ) + if (glz <- !(object@method %in% c("REML", "ML"))) { + cat(paste("Generalized linear mixed model fit using", + object@method, "\n")) + } else { + cat("Linear mixed-effects model fit by ") + cat(if(REML) "REML\n" else "maximum likelihood\n") + } if (!is.null(object@call$formula)) { cat("Formula:", deparse(object@call$formula),"\n") } @@ -231,20 +544,29 @@ cat(" Subset:", deparse(asOneSidedFormula(object@call$subset)[[2]]),"\n") } - print(data.frame(AIC = AIC(llik), BIC = BIC(llik), - logLik = c(llik), - MLdeviance = dev["ML"], - REMLdeviance = dev["REML"], - row.names = "")) + if (glz) { + cat(" Family: ", object@family$family, "(", + object@family$link, " link)\n", sep = "") + print(data.frame(AIC = AIC(llik), BIC = BIC(llik), + logLik = c(llik), + deviance = -2*llik, + row.names = "")) + } else { + print(data.frame(AIC = AIC(llik), BIC = BIC(llik), + logLik = c(llik), + MLdeviance = dev["ML"], + REMLdeviance = dev["REML"], + row.names = "")) + } cat("Random effects:\n") - show(VarCorr(object)) + vcShow(VarCorr(object), useScale) ngrps <- lapply(object@flist, function(x) length(levels(x))) - cat(sprintf("# of obs: %d, groups: ", object@nc[length(object@nc)])) + cat(sprintf("# of obs: %d, groups: ", devc[1])) cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; ")) cat("\n") if (!useScale) cat("\nEstimated scale (compare to 1) ", - .Call("lmer_sigma", object, object@REML, PACKAGE = "Matrix"), + .Call("mer_sigma", object, FALSE, PACKAGE = "Matrix"), "\n") if (nrow(coefs) > 0) { if (useScale) { @@ -263,181 +585,82 @@ } cat("\nFixed effects:\n") printCoefmat(coefs, tst.ind = 4, zap.ind = 3) - if (length(object@showCorrelation) > 0 && object@showCorrelation[1]) { - rn <- rownames(coefs) - dimnames(corF) <- list( - abbreviate(rn, minlen=11), - abbreviate(rn, minlen=6)) - if (!is.null(corF)) { - p <- NCOL(corF) - if (p > 1) { - cat("\nCorrelation of Fixed Effects:\n") - corF <- format(round(corF, 3), nsmall = 3) - corF[!lower.tri(corF)] <- "" - print(corF[-1, -p, drop=FALSE], quote = FALSE) - } + rn <- rownames(coefs) + if (!is.null(corF)) { + p <- ncol(corF) + if (p > 1) { + cat("\nCorrelation of Fixed Effects:\n") + corF <- matrix(format(round(corF@x, 3), nsmall = 3), + nc = p) + dimnames(corF) <- list( + abbreviate(rn, minlen=11), + abbreviate(rn, minlen=6)) + corF[!lower.tri(corF)] <- "" + print(corF[-1, -p, drop=FALSE], quote = FALSE) } } } invisible(object) }) -setMethod("lmer", signature(formula = "formula"), - function(formula, family, data, - method = c("REML", "ML", "PQL", "Laplace", "AGQ"), - control = list(), - subset, weights, na.action, offset, - model = TRUE, x = FALSE, y = FALSE, ...) - { - gVerb <- getOption("verbose") - # match and check parameters - controlvals <- do.call("lmerControl", control) - controlvals$REML <- FALSE - if (length(formula) < 3) stop("formula must be a two-sided formula") - - ## initial glm fit - mf <- match.call() - m <- match(c("family", "data", "subset", "weights", - "na.action", "offset"), - names(mf), 0) - mf <- mf[c(1, m)] - mf[[1]] <- as.name("glm") - fixed.form <- nobars(formula) - if (!inherits(fixed.form, "formula")) fixed.form <- ~ 1 # default formula - environment(fixed.form) <- environment(formula) - mf$formula <- fixed.form - mf$x <- mf$model <- mf$y <- TRUE - glm.fit <- eval(mf, parent.frame()) - family <- glm.fit$family - ## Note: offset is on the linear scale - offset <- glm.fit$offset - if (is.null(offset)) offset <- 0 - weights <- sqrt(abs(glm.fit$prior.weights)) - ## initial 'fitted' values on linear scale - etaold <- eta <- glm.fit$linear.predictors - - ## evaluation of model frame - mf$x <- mf$model <- mf$y <- mf$family <- NULL - mf$drop.unused.levels <- TRUE - this.form <- subbars(formula) - environment(this.form) <- environment(formula) - mf$formula <- this.form - mf[[1]] <- as.name("model.frame") - frm <- eval(mf, parent.frame()) - - ## grouping factors and model matrices for random effects - bars <- findbars(formula[[3]]) - random <- - lapply(bars, - function(x) list(model.matrix(eval(substitute(~term, - list(term=x[[2]]))), - frm), - eval(substitute(as.factor(fac)[,drop = TRUE], - list(fac = x[[3]])), frm))) - names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]]))) - - ## order factor list by decreasing number of levels - nlev <- sapply(random, function(x) length(levels(x[[2]]))) - if (any(diff(nlev) > 0)) { - random <- random[rev(order(nlev))] - } - mmats <- c(lapply(random, "[[", 1), - .fixed = list(cbind(glm.fit$x, .response = glm.fit$y))) - ## FIXME: Use Xfrm and Xmat to get the terms and assign - ## slots, pass these to lmer_create, then destroy Xfrm, Xmat, etc. - obj <- .Call("lmer_create", lapply(random, "[[", 2), mmats, PACKAGE = "Matrix") - obj@terms <- attr(glm.fit$model, "terms") - obj@assign <- attr(glm.fit$x, "assign") - obj@call <- match.call() - obj@REML <- FALSE - rm(glm.fit) - .Call("lmer_initial", obj, PACKAGE="Matrix") - mmats.unadjusted <- mmats - mmats[[1]][1,1] <- mmats[[1]][1,1] - conv <- FALSE - firstIter <- TRUE - msMaxIter.orig <- controlvals$msMaxIter - responseIndex <- ncol(mmats$.fixed) - - for (iter in seq(length = controlvals$PQLmaxIt)) - { - mu <- family$linkinv(eta) - dmu.deta <- family$mu.eta(eta) - ## weights (note: weights is already square-rooted) - w <- weights * dmu.deta / sqrt(family$variance(mu)) - ## adjusted response (should be comparable to X \beta, not including offset - z <- eta - offset + (mmats.unadjusted$.fixed[, responseIndex] - mu) / dmu.deta - .Call("nlme_weight_matrix_list", - mmats.unadjusted, w, z, mmats, PACKAGE="Matrix") - .Call("lmer_update_mm", obj, mmats, PACKAGE="Matrix") - if (firstIter) { - .Call("lmer_initial", obj, PACKAGE="Matrix") - if (gVerb) cat(" PQL iterations convergence criterion\n") - } - .Call("lmer_ECMEsteps", obj, - controlvals$niterEM, - FALSE, - controlvals$EMverbose, - PACKAGE = "Matrix") - LMEoptimize(obj) <- controlvals - eta[] <- offset + ## FIXME: should the offset be here ? - .Call("lmer_fitted", obj, - mmats.unadjusted, TRUE, PACKAGE = "Matrix") - crit <- max(abs(eta - etaold)) / (0.1 + max(abs(eta))) - if (gVerb) cat(sprintf("%03d: %#11g\n", as.integer(iter), crit)) - ## use this to determine convergence - if (crit < controlvals$tolerance) { - conv <- TRUE - break - } - etaold[] <- eta - - ## Changing number of iterations on second and - ## subsequent iterations. - if (firstIter) - { - controlvals$niterEM <- 2 - controlvals$msMaxIter <- 10 - firstIter <- FALSE - } - } - if (!conv) warning("IRLS iterations for glmm did not converge") - obj - }) +setMethod("vcov", signature(object = "mer"), + function(object, REML = object@method == "REML", + useScale = object@useScale,...) { + sc <- if (object@useScale) { + .Call("mer_sigma", object, REML, PACKAGE = "Matrix") + } else { 1 } + rr <- as(sc^2 * tcrossprod(solve(object@RXX)), "dpoMatrix") + rr@factors$correlation <- as(rr, "correlation") + rr + }) ## calculates degrees of freedom for fixed effects Wald tests ## This is a placeholder. The answers are generally wrong. It will ## be very tricky to decide what a 'right' answer should be with ## crossed random effects. -setMethod("getFixDF", signature(object="lmer"), - function(object, ...) - { - nc <- object@nc[-seq(along = object@Omega)] - p <- nc[1] - 1 - n <- nc[2] - rep(n-p, p) - }) +setMethod("getFixDF", signature(object="mer"), + function(object, ...) { + devc <- object@devComp + rep(as.integer(devc[1]- devc[2]), devc[2]) + }) -setMethod("logLik", signature(object="lmer"), - function(object, REML = object@REML, ...) { +setMethod("logLik", signature(object="mer"), + function(object, REML = object@method == "REML", ...) { val <- -deviance(object, REML = REML)/2 - nc <- object@nc[-seq(a = object@Omega)] - attr(val, "nall") <- attr(val, "nobs") <- nc[2] - attr(val, "df") <- nc[1] + length(ccoef(object)) - attr(val, "REML") <- REML + devc <- as.integer(object@devComp[1:2]) + attr(val, "nall") <- attr(val, "nobs") <- devc[1] + attr(val, "df") <- abs(devc[2]) + + length(.Call("mer_coef", object, 0, PACKAGE = "Matrix")) + attr(val, "REML") <- REML class(val) <- "logLik" val }) -setMethod("anova", signature(object = "lmer"), +setMethod("VarCorr", signature(x = "mer"), + function(x, REML = x@method == "REML", useScale = x@useScale, ...) + { + sc <- 1 + if (useScale) + sc <- .Call("mer_sigma", x, REML, PACKAGE = "Matrix") + sc2 <- sc * sc + ans <- lapply(x@Omega, function(el) { + el <- as(sc2 * solve(el), "dpoMatrix") + el@factors$correlation <- as(el, "correlation") + el + }) + attr(ans, "sc") <- sc + ans + }) + +setMethod("anova", signature(object = "mer"), function(object, ...) { mCall <- match.call(expand.dots = TRUE) dots <- list(...) modp <- logical(0) if (length(dots)) - modp <- sapply(dots, inherits, "lmer") | sapply(dots, inherits, "lm") + modp <- sapply(dots, inherits, "mer") | sapply(dots, inherits, "lm") if (any(modp)) { # multiple models - form table opts <- dots[!modp] mods <- c(list(object), dots[modp]) @@ -467,19 +690,18 @@ check.names = FALSE) class(val) <- c("anova", class(val)) attr(val, "heading") <- - c(header, "", "Models:", + c(header, "Models:", paste(names(mods), unlist(lapply(lapply(calls, "[[", "formula"), deparse)), - sep = ": "),"") + sep = ": ")) return(val) } else { foo <- object foo@status["factored"] <- FALSE - .Call("lmer_factor", foo, PACKAGE="Matrix") + .Call("mer_factor", foo, PACKAGE="Matrix") dfr <- getFixDF(foo) - rcol <- ncol(foo@RXX) - ss <- foo@RXX[ , rcol]^2 - ssr <- ss[[rcol]] + ss <- foo@rXy^2 + ssr <- exp(foo@devComp["logryy2"]) ss <- ss[seq(along = dfr)] names(ss) <- object@cnames[[".fixed"]][seq(along = dfr)] asgn <- foo@assign @@ -489,14 +711,16 @@ 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])) + #dfr <- unlist(lapply(split(dfr, asgn), function(x) x[1])) ms <- ss/df - f <- ms/(ssr/dfr) - P <- pf(f, df, dfr, lower.tail = FALSE) - table <- data.frame(df, ss, ms, dfr, f, P) + #f <- ms/(ssr/dfr) + #P <- pf(f, df, dfr, lower.tail = FALSE) + #table <- data.frame(df, ss, ms, dfr, f, P) + table <- data.frame(df, ss, ms) dimnames(table) <- list(nmeffects, - c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)")) +# c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)")) + c("Df", "Sum Sq", "Mean Sq")) if ("(Intercept)" %in% nmeffects) table <- table[-1,] attr(table, "heading") <- "Analysis of Variance Table" class(table) <- c("anova", "data.frame") @@ -504,208 +728,53 @@ } }) -setMethod("update", signature(object = "lmer"), - function(object, formula., ..., evaluate = TRUE) - { - call <- object@call - if (is.null(call)) - stop("need an object with call component") - extras <- match.call(expand.dots = FALSE)$... - if (!missing(formula.)) - call$formula <- update.formula(formula(object), formula.) - if (length(extras) > 0) { - existing <- !is.na(match(names(extras), names(call))) - for (a in names(extras)[existing]) call[[a]] <- extras[[a]] - if (any(!existing)) { - call <- c(as.list(call), extras[!existing]) - call <- as.call(call) - } - } - if (evaluate) - eval(call, parent.frame()) - else call - }) - - -setMethod("confint", signature(object = "lmer"), - function (object, parm, level = 0.95, ...) - { - cf <- fixef(object) - pnames <- names(cf) - if (missing(parm)) - parm <- seq(along = pnames) - else if (is.character(parm)) - parm <- match(parm, pnames, nomatch = 0) - a <- (1 - level)/2 - a <- c(a, 1 - a) - pct <- paste(round(100 * a, 1), "%") - ci <- array(NA, dim = c(length(parm), 2), - dimnames = list(pnames[parm], pct)) - ses <- sqrt(diag(vcov(object)))[parm] - ci[] <- cf[parm] + ses * t(outer(a, getFixDF(object)[parm], qt)) - ci - }) - -setMethod("param", signature(object = "lmer"), - function(object, unconst = FALSE, ...) { - .Call("lmer_coef", object, unconst, PACKAGE = "Matrix") - }) - -setMethod("deviance", "lmer", - function(object, REML = NULL, ...) { - .Call("lmer_factor", object, PACKAGE = "Matrix") - if (is.null(REML)) - REML <- if (length(oR <- object@REML)) oR else FALSE - object@deviance[[ifelse(REML, "REML", "ML")]] - }) - -setMethod("chol", signature(x = "lmer"), - function(x, pivot = FALSE, LINPACK = pivot) { - x@status["factored"] <- FALSE # force a decomposition - .Call("lmer_factor", x, PACKAGE = "Matrix") - }) - -setMethod("solve", signature(a = "lmer", b = "missing"), - function(a, b, ...) - .Call("lmer_invert", a, PACKAGE = "Matrix") +setMethod("confint", signature(object = "mer"), + function(object, parm, level = 0.95, ...) + .NotYetImplemented() ) -setMethod("formula", "lmer", function(x, ...) x@call$formula) - -setMethod("vcov", signature(object = "lmer"), - function(object, REML = object@REML, useScale = TRUE,...) { - sc <- .Call("lmer_sigma", object, REML, PACKAGE = "Matrix") - rr <- object@RXX - nms <- object@cnames[[".fixed"]] - dimnames(rr) <- list(nms, nms) - nr <- nrow(rr) - rr <- rr[-nr, -nr, drop = FALSE] - rr <- rr %*% t(rr) - if (useScale) { - rr = sc^2 * rr - } - rr - }) - -## Extract the L matrix -setAs("lmer", "dtTMatrix", - function(from) - { - ## force a refactorization if the factors have been inverted - if (from@status["inverted"]) from@status["factored"] <- FALSE - .Call("lmer_factor", from, PACKAGE = "Matrix") - L <- lapply(from@L, as, "dgTMatrix") - nf <- length(from@D) - Gp <- from@Gp - nL <- Gp[nf + 1] - Li <- integer(0) - Lj <- integer(0) - Lx <- double(0) - for (i in 1:nf) { - for (j in 1:i) { - Lij <- L[[Lind(i, j)]] - Li <- c(Li, Lij@i + Gp[i]) - Lj <- c(Lj, Lij@j + Gp[j]) - Lx <- c(Lx, Lij@x) - } - } - new("dtTMatrix", Dim = as.integer(c(nL, nL)), i = Li, j = Lj, x = Lx, - uplo = "L", diag = "U") - }) - -## Extract the ZZX matrix -setAs("lmer", "dsTMatrix", - function(from) - { - .Call("lmer_inflate", from, PACKAGE = "Matrix") - ZZpO <- lapply(from@ZZpO, as, "dgTMatrix") - ZZ <- lapply(from@ZtZ, as, "dgTMatrix") - nf <- length(ZZpO) - Gp <- from@Gp - nZ <- Gp[nf + 1] - Zi <- integer(0) - Zj <- integer(0) - Zx <- double(0) - for (i in 1:nf) { - ZZpOi <- ZZpO[[i]] - Zi <- c(Zi, ZZpOi@i + Gp[i]) - Zj <- c(Zj, ZZpOi@j + Gp[i]) - Zx <- c(Zx, ZZpOi@x) - if (i > 1) { - for (j in 1:(i-1)) { - ZZij <- ZZ[[Lind(i, j)]] - ## off-diagonal blocks are transposed - Zi <- c(Zi, ZZij@j + Gp[j]) - Zj <- c(Zj, ZZij@i + Gp[i]) - Zx <- c(Zx, ZZij@x) - } - } - } - new("dsTMatrix", Dim = as.integer(c(nZ, nZ)), i = Zi, j = Zj, x = Zx, - uplo = "U") - }) +setMethod("fitted", signature(object = "mer"), + function(object, ...) + .Call("mer_fitted", object, PACKAGE = "Matrix") + ) -setMethod("fitted", signature(object = "lmer"), - function(object, ...) object@fitted) +setMethod("formula", signature(x = "mer"), + function(x, ...) + x@callformula + ) -setMethod("residuals", signature(object = "lmer"), - function(object, ...) object@residuals) +setMethod("residuals", signature(object = "mer"), + function(object, ...) + .NotYetImplemented() + ) -setMethod("resid", signature(object = "lmer"), - function(object, ...) do.call("residuals", c(list(object), list(...)))) +setMethod("resid", signature(object = "mer"), + function(object, ...) + .NotYetImplemented() + ) -setMethod("coef", signature(object = "lmer"), +setMethod("summary", signature(object = "mer"), function(object, ...) - { - fef <- data.frame(rbind(fixef(object)), check.names = FALSE) - ref <- as(ranef(object), "list") - names(ref) <- names(object@flist) - val <- lapply(ref, function(x) fef[rep(1, nrow(x)),]) - for (i in seq(a = val)) { - refi <- ref[[i]] - row.names(val[[i]]) <- row.names(refi) - if (!all(names(refi) %in% names(fef))) - stop("unable to align random and fixed effects") - val[[i]][ , names(refi)] <- val[[i]][ , names(refi)] + refi - } - new("lmer.coef", val) - }) + .NotYetImplemented() + ) -setMethod("plot", signature(x = "lmer.coef"), - function(x, y, ...) - { - varying <- unique(do.call("c", - lapply(x, function(el) - names(el)[sapply(el, - function(col) - any(col != col[1]))]))) - gf <- do.call("rbind", lapply(x, "[", j = varying)) - gf.grp <- factor(rep(names(x), sapply(x, nrow))) - switch(min(length(varying), 3), - qqmath(eval(substitute(~ x | .grp, - list(x = as.name(varying[1])))), gf, ...), - xyplot(eval(substitute(y ~ x | .grp, - list(y = as.name(varying[1]), - x = as.name(varying[2])))), gf, ...), - splom(~ gf | .grp, ...)) - }) +setMethod("update", signature(object = "mer"), + function(object, ...) + .NotYetImplemented() + ) -setMethod("plot", signature(x = "lmer.ranef"), - function(x, y, ...) - { - lapply(x, function(x) { - cn <- lapply(colnames(x), as.name) - switch(min(ncol(x), 3), - qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...), - xyplot(eval(substitute(y ~ x, list(y = cn[[1]], x = cn[[2]]))), - x, ...), - splom(~ x, ...)) - }) - }) -setMethod("with", signature(data = "lmer"), - function(data, expr, ...) - eval(substitute(substitute(expr, list(. = quote(data))), - append(data@flist, eval(data@call\$data)), - enclos = parent.frame()) - }) +simss <- function(fm0, fma, nsim) +{ + ysim <- simulate(fm0, nsim) + cv <- list(analyticGradient = FALSE, msMaxIter = 200:200, + msVerbose = 0:0) + sapply(ysim, function(yy) { + .Call("mer_update_y", fm0, yy, PACKAGE = "Matrix") + LMEoptimize(fm0) <- cv + .Call("mer_update_y", fma, yy, PACKAGE = "Matrix") + LMEoptimize(fma) <- cv + exp(c(H0 = fm0@devComp[["logryy2"]], + Ha = fma@devComp[["logryy2"]])) + }) +}