R/fromLme4.R

In kollerma/robustlmm: Robust Linear Mixed Effects Models

## This file contains functions that have been copied
## from lme4 to enable a release of this package onto CRAN.
## Copied from Revision 1785
## Updated from commit 9c7edde3f081af68f84b1e9b61be1f5cc52cfaa0
## The authors of lme4 are:
## Douglas Bates <bates@stat.wisc.edu>,
## Martin Maechler <maechler@R-project.org> and
## Ben Bolker <bolker@mcmaster.ca>
## Steven Walker

## follows the file lme4/R/lmer.R

## minimal changes: merMod -> rlmerMod
## and replaced with a stop message if not defined

## use getS3method to copy methods from lme4
##' @importFrom utils getS3method

##' @importFrom stats coef
coefMer <- getS3method("coef", "merMod")
##' @export
coef.rlmerMod <- function(object, ...) {
    val <- coefMer(object, ...)
    class(val) <- "coef.rlmerMod"
    val
}

## deviance is in accessors.R

## FIXME what about drop1

##' @importFrom stats extractAIC
##' @export
extractAIC.rlmerMod <- function(fit, scale = 0, k = 2, ...)
   stop("AIC is not defined for rlmerMod objects")

##' @importFrom stats family
##' @export
family.rlmerMod <- function(object, ...) gaussian()

##' @importFrom stats fitted
##' @export
fitted.rlmerMod <- getS3method("fitted", "merMod")

## Extract the fixed-effects estimates.
##
## Extract the estimates of the fixed-effects parameters from a fitted model.
## @name fixef
## @title Extract fixed-effects estimates
## @aliases fixef fixed.effects fixef.rlmerMod
## @docType methods
## @param object any fitted model object from which fixed effects estimates can
## be extracted.
## @param \dots optional additional arguments. Currently none are used in any
## methods.
## @return a named, numeric vector of fixed-effects estimates.
## @keywords models
## @examples
## ## doFit = FALSE to speed up example
## fixef(rlmer(Reaction ~ Days + (Days|Subject), sleepstudy, doFit=FALSE))
##' @importFrom nlme fixef
##' @export
fixef.rlmerMod <- function(object, ...)
    structure(object@beta, names = dimnames(.X(object))[[2]])

##' @importFrom lme4 nobars
getFixedFormula <- function(form) {
    form[[3]] <- if (is.null(nb <- nobars(form[[3]]))) 1 else nb
    form
}

##' @importFrom stats formula
##' @export
formula.rlmerMod <- getS3method("formula", "merMod")

##' @importFrom lme4 isREML
##' @export
isREML.rlmerMod <- function(x, ...) .isREML(x, ...)

## needed for predict():
##' @importFrom lme4 isGLMM
##' @export
isGLMM.rlmerMod <- function(x, ...) FALSE

##' @importFrom lme4 isLMM
##' @export
isLMM.rlmerMod <- function(x, ...) TRUE

##' @importFrom lme4 isNLMM
##' @export
isNLMM.rlmerMod <- function(x, ...) FALSE

##' @importFrom stats logLik
##' @export
logLik.rlmerMod <- function(object, REML = NULL, ...)
   stop("log-likelihood is not defined for rlmerMod objects")

##' @importFrom stats model.frame
##' @export
model.frame.rlmerMod <- getS3method("model.frame", "merMod")

##' @importFrom stats model.matrix
##' @export
model.matrix.rlmerMod <- function(object, ...) .X(object)

## we have our own nobs.rlmerMod method

## ranef function is in accessors.R

## no refit methods

##' @importFrom lme4 REMLcrit
## don't add a function for now.

## residuals is in accessors.R

## sigma in in accessors.R

## no simulate method

##' @importFrom stats terms
##' @export
terms.rlmerMod <- getS3method("terms", "merMod")

## update is in helpers.R

## ...

.prt.resids <- function(resids, digits, title="Scaled residuals:", ...) {
    cat(title,"\n")
    rq <- setNames(zapsmall(quantile(resids), digits + 1L),
                   c("Min", "1Q", "Median", "3Q", "Max"))
    print(rq, digits=digits, ...)
    cat("\n")
}

.prt.call <- function(call, long=TRUE) {
    if (!is.null(cc <- call$formula))
	cat("Formula:", deparse(cc),"\n")
    if (!is.null(cc <- call$data))
	cat("   Data:", deparse(cc), "\n")
    if (!is.null(cc <- call$weights))
        cat("Weights:", deparse(cc), "\n")
    if (!is.null(cc <- call$offset))
        cat(" Offset:", deparse(cc), "\n")
    if (long && length(cc <- call$control) &&
	!identical((dc <- deparse(cc)), "lmerControl()"))
	## && !identical(eval(cc), lmerControl()))
	cat("Control:", dc, "\n")
    if (!is.null(cc <- call$subset))
	cat(" Subset:", deparse(asOneSidedFormula(cc)[[2]]),"\n")
}

.prt.VC <- function(varcor, digits, comp, formatter=format, ...) {
    cat("Random effects:\n")
    fVC <- if(missing(comp))
	formatVC(varcor, digits=digits, formatter=formatter)
    else
	formatVC(varcor, digits=digits, formatter=formatter, comp=comp)
    print(fVC, quote = FALSE, digits = digits, ...)
}

.prt.grps <- function(ngrps, nobs) {
    cat(sprintf("Number of obs: %d, groups: ", nobs))
    cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))
    cat("\n")
}

## .printRlmerMod is in helpers.R

## print.rlmerMod is in helpers.R

##' @exportMethod show
setMethod("show", "rlmerMod", function(object) print.rlmerMod(object))

## print.summary.rlmerMod is in helpers.R

## can import tnames if required

## getME is in accessors.R

globalVariables("forceSymmetric", add=TRUE)

##' @importFrom stats vcov
##' @export
vcov.rlmerMod <- getS3method("vcov", "merMod")

##' @importFrom stats vcov
##' @export
vcov.summary.rlmerMod <- function(object, correlation = TRUE, ...) {
    if(is.null(object$vcov)) stop("logic error in summary of rlmerMod object")
    object$vcov
}

##' @importFrom nlme VarCorr
##' @method VarCorr rlmerMod
##' @export
VarCorrMer <- getS3method("VarCorr", "merMod")
VarCorr.rlmerMod <- function(x, ...)# <- 3 args from nlme
{
    val <- VarCorrMer(x, ...)
    class(val) <- "VarCorr.rlmerMod"
    val
}

##' @export
VarCorr.summary.rlmerMod <- function(x, ...) x$varcor

##' @export
print.VarCorr.rlmerMod <- getS3method("print", "VarCorr.merMod")

##' @export
as.data.frame.VarCorr.rlmerMod <-
    getS3method("as.data.frame", "VarCorr.merMod")

## __NOT YET EXPORTED__
## "format()" the 'VarCorr' matrix of the random effects -- for
## print()ing and show()ing
##
## @title Format the 'VarCorr' Matrix of Random Effects
## @param varc a \code{\link{VarCorr}} (-like) matrix with attributes.
## @param digits the number of significant digits.
## @param comp character vector of length one or two indicating which
## columns out of "Variance" and "Std.Dev." should be shown in the
## formatted output.
## @param formatter the \code{\link{function}} to be used for
## formatting the standard deviations and or variances (but
## \emph{not} the correlations which (currently) are always formatted
## as "0.nnn"
## @param ... optional arguments for \code{formatter(*)} in addition
## to the first (numeric vector) and \code{digits}.
## @return a character matrix of formatted VarCorr entries from \code{varc}.
formatVC <- function(varcor, digits = max(3, getOption("digits") - 2),
                     comp = "Std.Dev.", formatter = format,
                     useScale = attr(varcor, "useSc"),
                     ...)
{
  c.nms <- c("Groups", "Name", "Variance", "Std.Dev.")
  avail.c <- c.nms[-(1:2)]
  if(anyNA(mcc <- pmatch(comp, avail.c)))
    stop("Illegal 'comp': ", comp[is.na(mcc)])
  nc <- length(colnms <- c(c.nms[1:2], (use.c <- avail.c[mcc])))
  if(length(use.c) == 0)
    stop("Must *either* show variances or standard deviations")
  reStdDev <- c(lapply(varcor, attr, "stddev"),
                if(useScale) list(Residual = unname(attr(varcor, "sc"))))
  reLens <- lengths(reStdDev)
  nr <- sum(reLens)
  reMat <- array('', c(nr, nc), list(rep.int('', nr), colnms))
  reMat[1+cumsum(reLens)-reLens, "Groups"] <- names(reLens)
  reMat[,"Name"] <- c(unlist(lapply(varcor, colnames)), if(useScale) "")
  if(any("Variance" == use.c))
    reMat[,"Variance"] <- formatter(unlist(reStdDev)^2, digits = digits, ...)
  if(any("Std.Dev." == use.c))
    reMat[,"Std.Dev."] <- formatter(unlist(reStdDev),   digits = digits, ...)
  if (any(reLens > 1)) {
    maxlen <- max(reLens)
    recorr <- lapply(varcor, attr, "correlation")
    corr <-
      do.call(rbind,
              lapply(recorr,
                     function(x) {
                       x <- as.matrix(x)
                       dig <- max(2, digits - 2) # use 'digits' !
                       ## not using formatter() for correlations
                       cc <- format(round(x, dig), nsmall = dig)
                       cc[!lower.tri(cc)] <- ""
                       nr <- nrow(cc)
                       if (nr >= maxlen) return(cc)
                       cbind(cc, matrix("", nr, maxlen-nr))
                     }))[, -maxlen, drop = FALSE]
    if (nrow(corr) < nrow(reMat))
      corr <- rbind(corr, matrix("", nrow(reMat) - nrow(corr), ncol(corr)))
    colnames(corr) <- c("Corr", rep.int("", max(0L, ncol(corr)-1L)))
    cbind(reMat, corr)
  } else reMat
}

## summary is in helpers.R

## plots are in plots.R

##' @importFrom stats weights
##' @export
weights.rlmerMod <- function(object, ...) {
  object@resp$weights
}

## no optimizer options in rlmer

#######################################################
## predict method                                    ##
#######################################################

reFormHack <- function(re.form,ReForm,REForm,REform) {
    if (!missing(ReForm)) {
        message(shQuote("re.form")," is now preferred to ",shQuote("ReForm"))
        return(ReForm)
    }
    if (!missing(REForm)) {
        message(shQuote("re.form")," is now preferred to ",shQuote("REForm"))
        return(REForm)
    }
    if (!missing(REform)) {
        message(shQuote("re.form")," is now preferred to ",shQuote("REform"))
        return(REform)
    }
    re.form
}

##' @importFrom stats predict
##' @importFrom lme4 mkReTrms findbars
##' @export
## FIXME: This is slightly modified from lme4 version 1.0-6.
##        Newer versions have an improved version.
predict.rlmerMod <- function(object, newdata=NULL,
                             re.form=NULL,
                             ReForm,
                             REForm,
                             REform,
                             terms=NULL, type=c("link","response"),
                             allow.new.levels=FALSE, na.action=na.pass, ...) {
    ## FIXME: appropriate names for result vector?
    ## FIXME: make sure behaviour is entirely well-defined for NA in grouping factors

    REform <- reFormHack(re.form,ReForm,REForm,REform)

    if (length(list(...)>0)) warning("unused arguments ignored")
    if (isLMM(object) && !missing(type)) warning("type argument ignored for linear mixed models")
    fit.na.action <- attr(object@frame,"na.action")
    type <- match.arg(type)
    if (!is.null(terms)) stop("terms functionality for predict not yet implemented")
    ## FIXME/WARNING: how do we/can we do this in an eval-safe way???
    form_orig <- formula(object)
    if (is.null(newdata) && is.null(REform)) {
        ## raw predict() call, just return fitted values (inverse-link if appropriate)
        if (isLMM(object) || isNLMM(object)) {
            pred <- fitted(object)
        } else { ## inverse-link
            pred <-  switch(type,response=object@resp$mu, ## fitted(object),
                            link=object@resp$eta)  ## fixme: getME() ?
        }
        if (!is.null(fit.na.action)) {
            pred <- napredict(fit.na.action,pred)
        }
        return(pred)

      } else { ## newdata and/or REform specified
        X_orig <- getME(object, "X")
        if (is.null(newdata)) {
            X <- X_orig
        } else {
            ## evaluate new fixed effect
            RHS <- formula(object,fixed.only=TRUE)[-2]
            Terms <- terms(object,fixed.only=TRUE)
            X <- model.matrix(RHS, mfnew <- model.frame(delete.response(Terms),
                                                        newdata, na.action=na.action),
                              contrasts.arg=attr(X_orig,"contrasts"))
        }
        pred <- drop(X %*% fixef(object))
        ## modified from predict.glm ...
        offset <- rep(0, nrow(X))
        tt <- terms(object)
        ## FIXME:: need to unname()  ?
        if (!is.null(off.num <- attr(tt, "offset"))) {
            for (i in off.num) offset <- offset + eval(attr(tt,"variables")[[i + 1]], newdata)
        }
        ## FIXME: is this redundant??
        if (!is.null(frOffset <- attr(object@frame,"offset")))
            offset <- offset + eval(frOffset, newdata)
        pred <- pred+offset
        if (is.null(REform)) {
            REform <- form_orig[-2]
        }
        ## FIXME: ??? can't apply is.na() to a 'language' object?
        ##  what's the appropriate test?
        if (is.language(REform)) {
            na.action.name <- deparse(match.call()$na.action) ## ugh
            if (!is.null(newdata) && na.action.name %in% c("na.exclude","na.omit")) {
                ## strip NAs from data for random-effects matrix construction
                if (length(nadrop <- attr(mfnew,"na.action"))>0) {
                    newdata <- newdata[-nadrop,]
                }
            }
            re <- ranef(object)
            ## ok? -- newdata used even though it was just tested for null
            if(is.null(newdata)) rfd <- object@frame else rfd <- newdata # get data for REform
            ReTrms <- mkReTrms(findbars(REform[[2]]), rfd)
            if (!allow.new.levels && any(sapply(ReTrms$flist,function(x) any(is.na(x)))))
                stop("NAs are not allowed in prediction data for grouping variables unless allow.new.levels is TRUE")
            unames <- unique(sort(names(ReTrms$cnms)))  ## FIXME: same as names(ReTrms$flist) ?
            ## convert numeric grouping variables to factors as necessary
            for (i in all.vars(REform[[2]])) {
                newdata[[i]] <- factor(newdata[[i]])
            }
            Rfacs <- setNames(lapply(unames,function(x) factor(eval(parse(text=x),envir=newdata))),
                              unames)
            new_levels <- lapply(Rfacs,function(x) levels(droplevels(factor(x))))
            ## FIXME: should this be unique(as.character(x)) instead?
            ##   (i.e., what is the proper way to protect against non-factors?)
            levelfun <- function(x,n) {
                ## find and deal with new levels
                if (any(!new_levels[[n]] %in% rownames(x))) {
                    if (!allow.new.levels) stop("new levels detected in newdata")
                    ## create an all-zero data frame corresponding to the new set of levels ...
                    newx <- as.data.frame(matrix(0,nrow=length(new_levels[[n]]),ncol=ncol(x),
                                                 dimnames=list(new_levels[[n]],names(x))))
                    ## then paste in the matching RE values from the original fit/set of levels
                    newx[rownames(x),] <- x
                    x <- newx
                }
                ## find and deal with missing old levels
                if (any(!rownames(x) %in% new_levels[[n]])) {
                    x <- x[rownames(x) %in% new_levels[[n]],,drop=FALSE]
                }
                x
            }
            ## fill in/delete levels as appropriate
            re_x <- mapply(levelfun,re,names(re),SIMPLIFY=FALSE)
            re_new <- list()
            if (any(!names(ReTrms$cnms) %in% names(re)))
                stop("grouping factors specified in REform that were not present in original model")
            ## pick out random effects values that correspond to
            ##  random effects incorporated in REform ...
            for (i in seq_along(ReTrms$cnms)) {
                rname <- names(ReTrms$cnms)[i]
                if (any(!ReTrms$cnms[[rname]] %in% names(re[[rname]])))
                    stop("random effects specified in REform that were not present in original model")
                re_new[[i]] <- re_x[[rname]][,ReTrms$cnms[[rname]]]
            }
            re_newvec <- unlist(lapply(re_new,t))  ## must TRANSPOSE RE matrices before unlisting
            if(!is.null(newdata)) pred <- pred + drop(as.matrix(re_newvec %*% ReTrms$Zt))
        } ## predictions with REform!=NA
        if (isGLMM(object) && type=="response") {
            pred <- object@resp$family$linkinv(pred)
        }
        ## fill in NAs as appropriate
        if (is.null(newdata) && !is.null(fit.na.action)) {
            pred <- napredict(fit.na.action,pred)
        } else {
            pred <- napredict(na.action,pred)
        }
        return(pred)
    }
}

forceCopy <- function(x) deepcopy(x)

# these functions pick up where findbars leaves off, in terms of sugar
# @param REtrm an element of the result of findbars
# @param REtrms the result of findbars
# @return \code{reexpr} gives a one-sided formula with the linear
# model formula for the raw model matrix. \code{grpfact} gives an
# expression with the name of the grouping factor associated with
# the raw model matrix. \code{termnms} gives a character vector with
# the names of the random effects terms.
reexpr <- function(REtrm) substitute( ~ foo, list(foo = REtrm[[2]]))
grpfact <- function(REtrm) substitute(factor(fac), list(fac = REtrm[[3]]))
termnms <- function(REtrms) vapply(REtrms, deparse1, "")

# mmList(): list of model matrices
# ------    called from getME()
mmList <- function(object, ...) UseMethod("mmList")
mmList.rlmerMod <- function(object, ...) mmList(formula(object), model.frame(object))
mmList.formula <- function(object, frame, ...) {
    bars <- findbars(object)
    mm <- setNames(lapply(bars, function(b) model.matrix(eval(reexpr(b), frame), frame)),
                   termnms(bars))
    grp <- lapply(lapply(bars, grpfact), eval, frame)
    nl <- vapply(grp, nlevels, 1L)
    if (any(diff(nl) > 0))
        mm[order(nl, decreasing=TRUE)]
    else
        mm
}