R/lmList.R

In nlme: Linear and Nonlinear Mixed Effects Models

###                  Create a list of lm objects
###
### Copyright 2005-2022  The R Core team
### Copyright 1997-2003  Jose C. Pinheiro,
###                      Douglas M. Bates <bates@stat.wisc.edu>
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available at
#  http://www.r-project.org/Licenses/
#

lmList <-
  ## a list of lm objects from a formula or a groupedData object
  function(object, data, level, subset, na.action = na.fail, pool = TRUE, warn.lm = TRUE)
  UseMethod("lmList")

lmList.groupedData <-
  function(object, data, level, subset, na.action = na.fail, pool = TRUE, warn.lm = TRUE)
{
  ### object will provide the formula, the data, and the groups
  form <- formula(object)
  args <- as.list(match.call())[-1L]
  args[["object"]] <- eval(substitute(Y ~ RHS,
                                      list(Y  = form[[2]],
                                           RHS= form[[3]][[2]])))
  if (!missing(data)) {
    message("'data' argument not used, but taken from groupedData object")
    args[["data"]] <- substitute(object)
  } else {
    args <- c(args, list(data = substitute(object)))
  }
  do.call(lmList.formula, args)
}

lmList.formula <- function(object, data, level, subset, na.action = na.fail,
                           pool = TRUE, warn.lm = TRUE)
{
  Call <- match.call()
  if (!missing(subset)) {
    data <-
      data[eval(asOneSidedFormula(Call[["subset"]])[[2]], data),, drop = FALSE]
  }
  if (!inherits(data, "data.frame")) data <- as.data.frame(data)
  data <- na.action(data)
  if (is.null(grpForm <- getGroupsFormula(object))) {
    if (inherits(data, "groupedData")) {
      if (missing(level))
        level <- length(getGroupsFormula(data, asList = TRUE))
      else if (length(level) > 1) {
	stop("multiple levels not allowed")
      }
      groups <- getGroups(data, level = level)[drop = TRUE]
      grpForm <- getGroupsFormula(data)
      Call$object <-
        eval(parse(text=paste(deparse(Call$object),
                              deparse(grpForm[[2]]), sep = "|")))
    } else {
      stop ("'data' must be a \"groupedData\" object if 'groups' argument is missing")
    }
  } else {
    if (missing(level))
      level <- length(getGroupsFormula(object, asList = TRUE))
    else if (length(level) > 1) {
      stop("multiple levels not allowed")
    }
    groups <- getGroups(data, form = grpForm, level = level)[drop = TRUE]
    object <- eval(substitute(Y ~ X, list(Y = getResponseFormula (object)[[2]],
                                          X = getCovariateFormula(object)[[2]])))
  }
  val <- lapply(split(data, groups),
		function(dat)
		    tryCatch(lm(object, data = dat, na.action = na.action),
			     error = function(e) e))
  val <- warnErrList(val, warn = warn.lm)
  if (inherits(data, "groupedData")) {
    ## saving labels and units for plots
    attr(val, "units") <- attr(data, "units")
    attr(val, "labels") <- attr(data, "labels")
  }

  structure(val, class = "lmList",
            dims = list(N = nrow(data), M = length(val)),
            call = Call,
            groupsForm = grpForm,
            groups = ordered(groups, levels = names(val)),
            origOrder = match(unique(as.character(groups)), names(val)),
            level = level,
            pool = pool)
}

###*# Methods for standard generics

augPred.lmList <-
  function(object, primary = NULL, minimum = min(primary),
	   maximum = max(primary), length.out = 51, ...)
{
  data <- eval(attr(object, "call")[["data"]])
  if (!inherits(data, "data.frame")) {
      stop(gettextf("'data' in %s call must evaluate to a data frame",
                     sQuote(substitute(object))), domain = NA)
  }
  if(is.null(primary)) {
    if (!inherits(data, "groupedData")) {
        stop(gettextf("%s without \"primary\" can only be used with fits of \"groupedData\" objects",
                      sys.call()[[1]]), domain = NA)
    }
    primary <- getCovariate(data)
    pr.var <- getCovariateFormula(data)[[2L]]
  } else {
    pr.var <- asOneSidedFormula(primary)[[2L]]
    primary <- eval(pr.var, data)
  }
  prName <- c_deparse(pr.var)
  newprimary <- seq(from = minimum, to = maximum, length.out = length.out)
  groups <- getGroups(object)
  grName <- deparse(gr.v <- getGroupsFormula(object)[[2]])
  ugroups <- unique(groups)
  value <- data.frame(rep(newprimary, length(ugroups)),
		      rep(ugroups, rep(length(newprimary), length(ugroups))))
  names(value) <- c(prName, grName)
  ## recovering other variables in data that may be needed for predictions
  ## varying variables will be replaced by their means
  summData <- gsummary(data, groups = groups)
  if (any(toAdd <- is.na(match(names(summData), names(value))))) {
    summData <- summData[, toAdd, drop = FALSE]
  }
  value[, names(summData)] <- summData[value[, 2], ]
  pred <- c(predict(object, value, asList = FALSE))
  newvals <- cbind(value[, 1:2], pred)
  names(newvals)[3] <- respName <-
    deparse(resp.v <- getResponseFormula(object)[[2]])
  orig <- data.frame(primary, groups, getResponse(object))
  names(orig) <- names(newvals)
  value <- rbind(orig, newvals)
  attributes(value[, 2]) <- attributes(groups)
  value[, ".type"] <- ordered(c(rep("original", nrow(data)),
				rep("predicted", nrow(newvals))),
			      levels = c("predicted", "original"))
  labs <- list(x = prName, y = respName)
  unts <- list(x = "", y = "")
  if(inherits(data, "groupedData")) {
    labs[names(attr(data, "labels"))] <- attr(data, "labels")
    unts[names(attr(data, "units"))] <- attr(data, "units")
    attr(value, "units") <- attr(data, "units")
  }
  structure(value, class = c("augPred", class(value)),
	    labels = labs,
	    units = unts,
	    formula = eval(substitute(Y ~ X | G,
				      list(Y = resp.v, X = pr.var, G = gr.v))))
}

coef.lmList <-
  ## Extract the coefficients and form a  data.frame if possible
  function(object, augFrame = FALSE, data = NULL,
           which = NULL, FUN = mean, omitGroupingFactor = TRUE, ...)
{
  coefs <- lapply(object, coef)
  non.null <- !vapply(coefs, is.null, NA)
  ## size the data frame to cope with combined levels for factors
  ## and name the columns so can fill by name
  if (any(non.null)) {
    coefNames <- unique(unlist(lapply(coefs[non.null], names)))
    co <- matrix(NA_real_,
                 ncol=length(coefNames),
                 nrow=length(coefs),
                 dimnames=list(names(object), coefNames))
    for (i in which(non.null)) {
      co[i, names(coefs[[i]])] <- coefs[[i]]
    }
    coefs <- as.data.frame(co)
    effectNames <- names(coefs)
    if(augFrame) {
      if (is.null(data)) {
        data <- getData(object)
      }
      data <- as.data.frame(data)
      if (is.null(which)) {
        which <- 1:ncol(data)
      }
      data <- data[, which, drop = FALSE]
      ## eliminating columns with same names as effects
      data <- data[, is.na(match(names(data), effectNames)), drop = FALSE]
      data <- gsummary(data, FUN = FUN, groups = getGroups(object))
      if (omitGroupingFactor) {
        data <- data[, is.na(match(names(data),
                                   names(getGroupsFormula(object, asList = TRUE)))),
                     drop = FALSE]
      }
      if (length(data) > 0) {
        coefs <- cbind(coefs, data[row.names(coefs),,drop = FALSE])
      }
    }
    attr(coefs, "level") <- attr(object, "level")
    attr(coefs, "label") <- "Coefficients"
    attr(coefs, "effectNames") <- effectNames
    attr(coefs, "standardized") <- FALSE
    attr(coefs, "grpNames") <- deparse(getGroupsFormula(object)[[2]])
    class(coefs) <- c("coef.lmList", "ranef.lmList", class(coefs))
  }
  coefs
}

fitted.lmList <-
  function(object, subset = NULL, asList = FALSE, ...)
{
  if(!is.null(subset)) {
    if(is.character(subset)) {
      if (any(is.na(match(subset, names(object))))) {
        stop("nonexistent groups requested in 'subset'")
      }
    } else {
      if (is.integer(subset)) {
        if (any(is.na(match(subset, seq_along(object))))) {
          stop("nonexistent groups requested in 'subset'")
        }
      } else {
        stop("'subset' can only be character or integer")
      }
    }
    oclass <- class(object)
    oatt <- attr(object, "call")
    object <- object[subset]
    attr(object, "call") <- oatt
    class(object) <- oclass
  }
  val <- lapply(object, fitted)
  if(!asList) {				#convert to array
    ngrps <- table(getGroups(object))[names(object)]
    if(any(aux <- vapply(object, is.null, NA))) {
      for(i in names(ngrps[aux])) {
	val[[i]] <- rep(NA, ngrps[i])
      }
    }
    val <- val[attr(object, "origOrder")] # putting in original order
    namVal <- names(val)
    val <- unlist(val)
    names(val) <- rep(namVal, ngrps)
  }
  lab <- "Fitted values"
  if (!is.null(aux <- attr(object, "units")$y)) {
    lab <- paste(lab, aux)
  }
  attr(val, "label") <- lab
  val
}

fixef.lmList <- function(object, ...)
{
  coeff <- coef(object)
  if(is.matrix(coeff) || is.data.frame(coeff)) {
    colMeans(coeff, na.rm = TRUE)
  } # else NULL
}

formula.lmList <-
  function(x, ...) eval(attr(x, "call")[["object"]])

getData.lmList <-
  function(object)
{
  mCall <- attr(object, "call")
  data <- eval(mCall$data)
  if (is.null(data)) return(data)
  naAct <- eval(mCall$na.action)
  if (!is.null(naAct)) {
    data <- naAct(data)
  }
  subset <- mCall$subset
  if (!is.null(subset)) {
    subset <- eval(asOneSidedFormula(subset)[[2]], data)
    data <- data[subset, ]
  }
  return(data)
}

getGroups.lmList <-  function(object, form, level, data, sep)
  attr(object, "groups")

getGroupsFormula.lmList <- function(object, asList = FALSE, sep) {
  val <- attr(object, "groupsForm")
  getGroupsFormula(eval(substitute(~ 1 | GR, list(GR = val[[2]]))),
		   asList = asList)
}

getResponse.lmList <- function(object, form) fitted(object) + resid(object)

intervals.lmList <-
  function(object, level = 0.95, pool = attr(object, "pool"), ...)
{
  smry <- summary(object, pool = pool)
  coeff <- coef(smry)
  out <- coeff[ , 1:3 , ]
  dn <- dimnames(out)
  dimnames(out) <-
    if(is.null(dn))
      list(NULL, c("lower", "est.", "upper"))
    else {
      dn[[2]] <- c("lower", "est.", "upper")
      dn
    }
  mult <- sqrt(qf(level, 1, smry$df.residual))
  out[ , "est.", ] <- coeff[ , "Estimate",  ]
  out[ , "lower", ] <- out[ , "est.", ] - mult * coeff[ , "Std. Error", ]
  out[ , "upper", ] <- out[ , "est.", ] + mult * coeff[ , "Std. Error", ]
  attr(out, "groupsName") <- deparse(attr(object, "groupsForm")[[2]])
  class(out) <- "intervals.lmList"
  out
}

logLik.lmList <-
  function(object, REML = FALSE, pool = attr(object, "pool"), ...)
{
  if(any(vapply(object, is.null, NA))) {
    stop("log-likelihood not available with NULL fits")
  }
  if(pool) {
    aux <- rowSums(sapply(object, function(el) {
      res <- resid(el)
      p <- el$rank
      n <- length(res)
      if (is.null(w <- el$weights)) w <- rep(1, n)
      else {
        excl <- w == 0
        if (any(excl)) {
          res <- res[!excl]
          n <- length(res)
          w <- w[!excl]
        }
      }
      c(n, sum(w * res^2), p, sum(log(w)),
        sum(log(abs(diag(el$qr$qr)[1:p]))))
    }))
    N <- aux[1] - REML * aux[3]
    val <- (aux[4] - N * (log(2 * pi) + 1 - log(N) + log(aux[2])))/2 -
      REML * aux[5]
    attr(val, "nall") <- aux[1]
    attr(val, "nobs") <- N
    attr(val, "df") <- aux[3] + 1
  } else {
    aux <- lapply(object, logLik, REML)
    val <- sum(unlist(aux))
    attr(val, "nobs") <- sum(sapply(aux, function(x) attr(x, "nobs")))
    attr(val, "df") <- sum(sapply(aux, function(x) attr(x, "df")))
  }
  class(val) <- "logLik"
  val
}

pairs.lmList <-
  function(x, form = ~ coef(.), label, id = NULL, idLabels = NULL,
	   grid = FALSE, ...)
{
  object <- x
  ## scatter plot matrix plots, generally based on coef or random.effects
  if (!inherits(form, "formula")) {
    stop("'form' must be a formula")
  }
  if (length(form) != 2) {
    stop("'form' must be a one-sided formula")
  }
  ## constructing data
  allV <- all.vars(asOneFormula(form, id, idLabels))
  allV <- allV[is.na(match(allV,c("T","F","TRUE","FALSE")))]
  if (length(allV) > 0) {
    data <- getData(object)
    if (is.null(data)) {		# try to construct data
      alist <- lapply(as.list(allV), as.name)
      names(alist) <- allV
      alist <- c(as.list(quote(data.frame)), alist)
      mode(alist) <- "call"
      data <- eval(alist, sys.parent(1))
    } else {
      if (any(naV <- is.na(match(allV, names(data))))) {
        stop(sprintf(ngettext(sum(naV),
                              "%s not found in data",
                              "%s not found in data"),
                     allV[naV]), domain = NA)
     }
    }
  } else data <- NULL

  ## argument list
  dots <- list(...)
  args <- if (length(dots) > 0) dots else list()

  ## covariate - must be present as a data.frame
  covF <- getCovariateFormula(form)
  .x <- eval(covF[[2]], list(. = object)) # only function of "."
  if (!inherits(.x, "data.frame")) {
    stop("covariate must be a data frame")
  }
  if (!is.null(effNams <- attr(.x, "effectNames"))) {
    .x <- .x[, effNams, drop = FALSE]
  }
  ## eliminating constant effects
  isFixed <- vapply(.x, function(el) length(unique(el)) == 1, NA)
  .x <- .x[, !isFixed, drop = FALSE]
  nc <- ncol(.x)
  if (nc == 1) {
    stop("cannot do pairs of just one variable")
  }
  if (!missing(label)) {
    names(.x) <- labels
  }
  if (nc == 2) {
    ## will use xyplot
    argForm <- .y ~ .x
    argData <- .x
    names(argData) <- c(".x", ".y")
    if (is.null(args$xlab)) {
      args$xlab <- names(.x)[1]
    }
    if (is.null(args$ylab)) {
      args$ylab <- names(.x)[2]
    }
  } else {				# splom
    argForm <- ~ .x
    argData <- list(.x = .x)
  }

  auxData <- list()
  ## groups - need not be present
  grpsF <- getGroupsFormula(form)
  if (!is.null(grpsF)) {
    gr <- splitFormula(grpsF, sep = "*")
    for(i in seq_along(gr)) {
      auxData[[deparse(gr[[i]][[2]])]] <- eval(gr[[i]][[2]], data)
    }
    argForm <- eval(substitute(if(length(argForm) == 2) ~ .x | R else .y ~ .x | R,
			       list(R = grpsF[[2L]])))
  }

  ## id and idLabels - need not be present
  if (!is.null(id)) {			# identify points in plot
    N <- attr(object, "dims")$N
    id <-
      switch(mode(id),
	     numeric = {
	       if ((id <= 0) || (id >= 1)) {
		 stop("'id' must be between 0 and 1")
	       }
	       aux <- as.matrix(na.omit(ranef(object)))
	       auxV <- t(chol(var(aux)))
	       as.logical(colSums((solve(auxV, t(aux)))^2) >
                          qchisq(1 - id, dim(aux)[2]))
	     },
	     call = eval(asOneSidedFormula(id)[[2]], data),
	     stop("'id' can only be a formula or numeric")
	     )
    if (length(id) == N) {
      ## id as a formula evaluated in data
      auxData[[".id"]] <- id
    }

    if (is.null(idLabels)) {
      idLabels <- row.names(.x)
    } else {
      if (mode(idLabels) == "call") {
	idLabels <-
	  as.character(eval(asOneSidedFormula(idLabels)[[2]], data))
      } else if (is.vector(idLabels)) {
	if (length(idLabels <- unlist(idLabels)) != N) {
	  stop("'idLabels' of incorrect length")
	}
	idLabels <- as.character(idLabels)
      } else {
	stop("'idLabels' can only be a formula or a vector")
      }
    }
    if (length(idLabels) == N) {
      ## idLabels as a formula evaluated in data
      auxData[[".Lid"]] <- idLabels
    }
  }

  if (length(auxData)) {		# need collapsing
    auxData <- gsummary(data.frame(auxData),
			groups = getGroups(object))
    auxData <- auxData[row.names(.x), , drop = FALSE]
    if (!is.null(auxData[[".g"]])) {
      argData[[".g"]] <- auxData[[".g"]]
    }

    if (!is.null(auxData[[".id"]])) {
      id <- auxData[[".id"]]
    }

    if (!is.null(auxData[[".Lid"]])) {
      idLabels <- auxData[[".Lid"]]
    }
    wchDat <- is.na(match(names(auxData), c(".id", ".idLabels")))
    if (any(wchDat)) {
      argData <- cbind(argData, auxData[, wchDat, drop = FALSE])

    }
  }
  if (is.null(id)) assign("id", FALSE)
  else assign("id", as.logical(as.character(id)) )
  assign("idLabels", as.character(idLabels))
  assign("grid", grid)

  ## adding to args list
  args <- c(list(argForm, data = argData), args)
##   if (is.null(args$strip)) {
##     args$strip <- function(...) strip.default(..., style = 1)
##   }
  if (is.null(args$cex)) args$cex <- par("cex")
  if (is.null(args$adj)) args$adj <- par("adj")

  ## defining the type of plot
  if (length(argForm) == 3) {		# xyplot
    plotFun <- "xyplot"
    args <- c(args,
	      panel = list(function(x, y, subscripts, ...)
		  {
                    x <- as.numeric(x)
                    y <- as.numeric(y)
                    dots <- list(...)
		    if (grid) panel.grid()
		    panel.xyplot(x, y)
                    if (any(id) & any(ids <- id[subscripts])) {
                        ltext(x[ids], y[ids], idLabels[subscripts][ids],
                             cex = dots$cex, adj = dots$adj)
		    }
		  }))
  } else {				# splom
    plotFun <- "splom"
    args <- c(args,
	      panel = list(function(x, y, subscripts, ...)
		  {
                    x <- as.numeric(x)
                    y <- as.numeric(y)
                    dots <- list(...)
		    if (grid) panel.grid()
		    panel.xyplot(x, y)
                    if (any(id) & any(ids <- id[subscripts])) {
                        ltext(x[ids], y[ids], idLabels[subscripts][ids],
                              cex = dots$cex, adj = dots$adj)
		    }
		  }))
  }
  do.call(plotFun, args)
} ## {pairs.lmList}

plot.intervals.lmList <-
    function(x,
             xlab = "", ylab = attr(x, "groupsName"),
             strip = function(...) strip.default(..., style = 1),
             ...)
{
  dims <- dim(x)
  dn <- dimnames(x)
  ## changed definition of what to ordered to preserve order of parameters
  tt <- data.frame(group = ordered(rep(dn[[1]], dims[2] * dims[3]),
                   levels = dn[[1]]),
		   intervals = as.vector(x),
		   what = ordered(rep(dn[[3]],
                   rep(dims[1] * dims[2], dims[3])), levels = dn[[3]]))
  dotplot(group ~ intervals | what,
	  data = tt,
	  scales = list(x="free"),
	  strip = strip,
	  xlab = xlab, ylab = ylab,
	  panel = function(x, y, pch = dot.symbol$pch,
	      col = dot.symbol$col, cex = dot.symbol$cex,
	      font = dot.symbol$font, ...)
	  {
            x <- as.numeric(x)
            y <- as.numeric(y)
	    ok <- !is.na(x) & !is.na(y)
	    yy <- y[ok]
	    xx <- x[ok]
	    dot.symbol <- trellis.par.get("dot.symbol")
	    dot.line <- trellis.par.get("dot.line")
	    panel.abline(h = yy, lwd = dot.line$lwd, lty = dot.line$lty, col =
                         dot.line$col)
	    lpoints(xx, yy, pch = "|", col = col, cex = cex, font = font, ...)
	    lower <- tapply(xx, yy, min)
	    upper <- tapply(xx, yy, max)
	    nams <- as.numeric(names(lower))
	    lsegments(lower, nams, upper, nams, col = 1, lty = 1,
                      lwd = if (dot.line$lwd) dot.line$lwd else 2)
	  }, ...)
} ## {plot.intervals.lmList}

plot.ranef.lmList <- function(x, form = NULL, grid = TRUE, control, ...)
{
    plot.ranef.lme(x, form=form, grid=grid, control=control, ...)
}

plot.lmList <-
  function(x, form = resid(., type = "pool") ~ fitted(.), abline,
	   id = NULL, idLabels = NULL,  grid, ...)
  ## Diagnostic plots based on residuals and/or fitted values
{
  if(!inherits(form, "formula")) stop("'form' must be a formula")

  ## constructing data
  allV <- all.vars(asOneFormula(form, id, idLabels))
  allV <- allV[is.na(match(allV,c("T","F","TRUE","FALSE")))]
  if (length(allV) > 0) {
    data <- getData(x)
    if (is.null(data)) {		# try to construct data
      alist <- lapply(as.list(allV), as.name)
      names(alist) <- allV
      alist <- c(as.list(quote(data.frame)), alist)
      mode(alist) <- "call"
      data <- eval(alist, sys.parent(1))
    } else {
      if (any(naV <- is.na(match(allV, names(data))))) {
        stop(sprintf(ngettext(sum(naV),
                              "%s not found in data",
                              "%s not found in data"),
                      allV[naV]), domain = NA)
     }
    }
  } else data <- NULL

  if (inherits(data, "groupedData")) {	# save labels and units, if present
    ff <- formula(data)
    rF <- deparse(getResponseFormula(ff)[[2]])
    cF <- c_deparse(getCovariateFormula(ff)[[2]])
    lbs <- attr(data, "labels")
    unts <- attr(data, "units")
    if (!is.null(lbs$x)) cL <- paste(lbs$x, unts$x) else cF <- NULL
    if (!is.null(lbs$y)) rL <- paste(lbs$y, unts$y) else rF <- NULL
  } else {
    rF <- cF <- NULL
  }

  ## argument list
  dots <- list(...)
  args <- if(length(dots) > 0) dots else list()
  ## appending object to data
  data <- as.list(c(as.list(data), . = list(x)))

  ## covariate - must always be present
  covF <- getCovariateFormula(form)
  .x <- eval(covF[[2]], data)
  if (!is.numeric(.x)) {
    stop("covariate must be numeric")
  }
  argForm <- ~ .x
  argData <- as.data.frame(.x)
  if (is.null(xlab <- attr(.x, "label"))) {
    xlab <- c_deparse(covF[[2]])
    if (!is.null(cF) && (xlab == cF)) xlab <- cL
    else if (!is.null(rF) && (xlab == rF)) xlab <- rL
  }
  if (is.null(args$xlab)) args$xlab <- xlab

  ## response - need not be present
  respF <- getResponseFormula(form)
  if (!is.null(respF)) {
    .y <- eval(respF[[2]], data)
    if (is.null(ylab <- attr(.y, "label"))) {
      ylab <- deparse(respF[[2]])
      if      (!is.null(cF) && (ylab == cF)) ylab <- cL
      else if (!is.null(rF) && (ylab == rF)) ylab <- rL
    }
    argForm <- .y ~ .x
    argData[, ".y"] <- .y
    if (is.null(args$ylab)) args$ylab <- ylab
  }

  ## groups - need not be present
  grpsF <- getGroupsFormula(form)
  if (!is.null(grpsF)) {
    gr <- splitFormula(grpsF, sep = "*")
    for(i in seq_along(gr)) {
      argData[[deparse(gr[[i]][[2]])]] <- eval(gr[[i]][[2]], data)
    }
    argForm <- eval(substitute(if(length(argForm) == 2) ~ .x | R else .y ~ .x | R,
			       list(R = grpsF[[2L]])))
  }
  ## adding to args list
  args <- c(list(argForm, data = argData), args)
##   if (is.null(args$strip)) {
##     args$strip <- function(...) strip.default(..., style = 1)
##   }
  if (is.null(args$cex)) args$cex <- par("cex")
  if (is.null(args$adj)) args$adj <- par("adj")

  if (!is.null(id)) {			# identify points in plot
    id <-
      switch(mode(id),
	     numeric = {
	       if ((id <= 0) || (id >= 1)) {
		 stop("'id' must be between 0 and 1")
	       }
	       as.logical(abs(resid(x, type = "pooled")) > -qnorm(id / 2))
	     },
	     call = eval(asOneSidedFormula(id)[[2]], data),
	     stop("'id' can only be a formula or numeric")
	     )
    if (is.null(idLabels)) {
      idLabels <- getGroups(x)
      if (length(idLabels) == 0) idLabels <- 1:x$dims$N
      idLabels <- as.character(idLabels)
    } else {
      if (mode(idLabels) == "call") {
	idLabels <-
	  as.character(eval(asOneSidedFormula(idLabels)[[2]], data))
      } else if (is.vector(idLabels)) {
	if (length(idLabels <- unlist(idLabels)) != length(id)) {
	  stop("'idLabels' of incorrect length")
	}
	idLabels <- as.character(idLabels)
      } else {
	stop("'idLabels' can only be a formula or a vector")
      }
    }
  }

  ## defining abline, if needed
  if(missing(abline)) abline <- if(missing(form)) c(0, 0) # else NULL

  ## defining the type of plot
  if (length(argForm) == 3) {
    if (is.numeric(.y)) {		# xyplot
      plotFun <- "xyplot"
      args$panel <- function(x, y, subscripts, ...)
		    {
                      x <- as.numeric(x)
                      y <- as.numeric(y)
                      dots <- list(...)
		      if (grid) panel.grid()
		      panel.xyplot(x, y)
                      if (any(ids <- id[subscripts])) {
                          ltext(x[ids], y[ids], idLabels[subscripts][ids],
                                cex = dots$cex, adj = dots$adj)
                      }
		      if (!is.null(abline)) {
                          if (length(abline) == 2)
                               panel.abline(a = abline, ...)
                          else panel.abline(h = abline, ...)
		      }
		    }
    } else {				# assume factor or character
      plotFun <- "bwplot"
      args$panel <- function(x, y, ...)
		    {
		      if (grid) panel.grid()
		      panel.bwplot(x, y)
		      if (!is.null(abline)) {
			panel.abline(v = abline[1], ...)
		      }
		    }
    }
  } else { ## '~ x'
    plotFun <- "histogram"
    args$panel <- function(x, y, ...)
		  {
		    if (grid) panel.grid()
		    panel.histogram(x, y)
		    if (!is.null(abline)) {
		      panel.abline(v = abline[1], ...)
		    }
		  }
  }
  ## defining grid (seen from panel()s defined here):
  if (missing(grid)) grid <- (plotFun == "xyplot") # T/F
  do.call(plotFun, args)
} ## {plot.lmList}

predict.lmList <-
  function(object, newdata, subset = NULL, pool = attr(object, "pool"),
	   asList = FALSE, se.fit = FALSE, ...)
{
  if(missing(newdata)) {
    if (!se.fit) return(fitted(object, subset, asList))
    myData <- getData(object)
    grps <- getGroups(object)
    myData <- split(myData, grps)
    newdata <- NULL
    sameData <- FALSE
  } else {
    newdata <- as.data.frame(newdata)
    sameData <- TRUE
    ## checking if same newdata for all groups
    formGrps <- getGroupsFormula(object)
    if(all(match(all.vars(formGrps), names(newdata), 0))) {
      ## newdata contains groups definition
      grps <- getGroups(newdata, getGroupsFormula(object, asList = TRUE),
			level = attr(object, "level"))
      grps <- grps[drop = TRUE]
      subset <- as.character(unique(grps))
      if(any(is.na(match(subset, names(object))))) {
	stop("nonexistent group in 'newdata'")
      }
      myData <- split(newdata, grps)
      newdata <- NULL
      sameData <- FALSE
    }
  }
  if(!is.null(subset)) {
    if(any(is.na(match(subset, names(object)))))
      stop("nonexistent group requested in 'subset'")
    oclass <- class(object)
    ## fix for PR#13788
    oatt <- attributes(object)[c("call", "groupsForm", "pool")]
    object <- object[subset]
    attributes(object) <- c(attributes(object), oatt)
    class(object) <- oclass
    if(is.null(newdata))
      myData <- myData[subset]
  }
  nmGrps <- names(object)
  noNull <- !vapply(object, is.null, NA)
  val <- vector("list", length(nmGrps))
  names(val) <- nmGrps
  if(!sameData) {
    if(!se.fit) {
      for(i in nmGrps[noNull]) {
        val[[i]] <- predict(object[[i]], myData[[i]])
      }
    } else {
      if(pool) {
	poolSD <- pooledSD(object)
      }
      for(i in nmGrps[noNull]) {
	aux <- predict(object[[i]], myData[[i]], se.fit = TRUE)
	if(pool) {
	  val[[i]] <- data.frame(fit = aux$fit,
				 se.fit = aux$se.fit*poolSD/aux$residual.scale)
	} else {
	  val[[i]] <- data.frame(fit = aux$fit, se.fit = aux$se.fit)
	}
      }
    }
  } else {
    if(pool) {
      poolSD <- pooledSD(object)
      val[noNull] <-
	lapply(object[noNull],
	       function(el, newdata, se.fit, poolSD) {
		 aux <- predict(el, newdata, se.fit = se.fit)
		 if(se.fit) {
		   data.frame(fit = aux$fit,
			      se.fit = aux$se.fit*poolSD/aux$residual.scale)
		 } else {
		   aux
		 }
	       }, newdata = newdata, se.fit = se.fit, poolSD = poolSD)
    } else {
      val[noNull] <-
	lapply(object[noNull],
	       function(el, newdata, se.fit) {
		 aux <- predict(el, newdata, se.fit = se.fit)
		 if(se.fit) {
		   data.frame(fit = aux$fit, se.fit = aux$se.fit)
		 } else {
		   aux
		 }
	       }, newdata = newdata, se.fit = se.fit)
    }
  }
  if(!asList) {				#convert to array
    if(is.null(newdata)) {
      ngrps <- table(grps)[names(object)]
    } else {
      ngrps <- rep(dim(newdata)[1], length(object))
      names(ngrps) <- names(object)
    }
    if(any(aux <- vapply(object, is.null, NA))) {
      for(i in names(ngrps[aux])) {
	aux1 <- rep(NA, ngrps[i])
	if(se.fit) {
	  val[[i]] <- data.frame(fit = aux1, se.fit = aux1)
	} else {
	  val[[i]] <- aux1
	}
      }
    }
    if(se.fit) {
      val <- do.call("rbind", val)
      val[, as.character(getGroupsFormula(object)[[2]])] <-
	rep(names(ngrps), ngrps)
      val <- val[, c(3,1,2)]
      row.names(val) <- 1:nrow(val)
    } else {
      val <- unlist(val)
      names(val) <- rep(names(ngrps), ngrps)
      attr(val, "label") <- "Predicted values"
      if (!is.null(aux <- attr(object, "units")$y)) {
        attr(val, "label") <- paste(attr(val, "label"), aux)
      }
    }
  }
  val
}

print.intervals.lmList <-
  function(x, ...)
{
  ox <- x
  x <- unclass(x)
  attr(x, "groupsName") <- NULL
  print(x, ...)
  invisible(ox)
}

print.lmList <- function(x, pool = attr(x, "pool"), ...)
{
  mCall <- attr(x, "call")
  cat("Call:\n")
  form <- formula(x)
  cat("  Model:", deparse(getResponseFormula(form)[[2]]),
      "~", c_deparse(getCovariateFormula(form)[[2]]), "|",
      deparse(getGroupsFormula(x)[[2]]), "\n")
  if (!is.null(mCall$level)) {
    cat(" Level:", mCall$level, "\n")
  }
  cat("   Data:", deparse(mCall$data),"\n\n")
  cat("Coefficients:\n")
  print(coef(x), ...)
  if(pool) {
    cat("\n")
    poolSD <- pooledSD(x)
    dfRes <- attr(poolSD, "df")
    RSE <- c(poolSD)
    cat("Degrees of freedom: ", length(unlist(lapply(x, fitted))),
	" total; ", dfRes, " residual\n", sep = "")
    cat("Residual standard error:", format(RSE))
    cat("\n")
  }
  invisible(x)
}

print.summary.lmList <- function(x, ...)
{
  cat("Call:\n")
  form <- formula(x)
  cat("  Model:", deparse(getResponseFormula(form)[[2]]),
      "~", c_deparse(getCovariateFormula(form)[[2]]), "|",
      deparse(attr(x, "groupsForm")[[2]]), "\n")
  if (!is.null(x$call$level)) {
    cat(" Level:", x$call$level, "\n")
  }
  cat("   Data:", deparse(x$call$data),"\n\n")
  cat("Coefficients:\n")
  for(i in dimnames(coef(x))[[3]]) {
    cat("  ",i,"\n")
    print(coef(x)[,,i], ...)
  }
  if(x$pool) {
    cat("\n")
    cat("Residual standard error:", format(x$RSE), "on",
	x$df.residual, "degrees of freedom\n")
  }
  cat("\n")
  invisible(x)
}

qqnorm.lmList <-
  function(y, form = ~ resid(., type = "pooled"), abline = NULL,
           id = NULL, idLabels = NULL, grid = FALSE, resType = "pool", ...)
  ## normal probability plots for residuals and random effects
{
  object <- y
  if (!inherits(form, "formula")) {
    stop("'form' must be a formula")
  }
  ## constructing data
  allV <- all.vars(asOneFormula(form, id, idLabels))
  allV <- allV[is.na(match(allV,c("T","F","TRUE","FALSE")))]
  if (length(allV) > 0) {
    data <- getData(object)
    if (is.null(data)) {		# try to construct data
      alist <- lapply(as.list(allV), as.name)
      names(alist) <- allV
      alist <- c(as.list(quote(data.frame)), alist)
      mode(alist) <- "call"
      data <- eval(alist, sys.parent(1))
    } else if (any(naV <- is.na(match(allV, names(data))))) {
        stop(sprintf(ngettext(sum(naV),
                              "%s not found in data",
                              "%s not found in data"),
                     allV[naV]), domain = NA)
    }
  } else data <- NULL
  ## argument list
  dots <- list(...)
  args <- if (length(dots) > 0) dots else list()
  ## appending object to data
  data <- as.list(c(as.list(data), . = list(object)))

  ## covariate - must always be present
  covF <- getCovariateFormula(form)
  .x <- eval(covF[[2]], data)
  labs <- attr(.x, "label")
  type <-
    if (inherits(.x, "ranef.lmList"))
      "reff" # random effects
    else if (!is.null(labs) && ((labs == "Standardized residuals") ||
                                (substr(labs, 1, 9) == "Residuals")))
      "res" # residuals
    else
      stop("only residuals and random effects allowed")
  if (is.null(args$xlab)) args$xlab <- labs
  if (is.null(args$ylab)) args$ylab <- "Quantiles of standard normal"
  if(type == "res") { # case 1: -------- residuals ------------------------
    fData <- qqnorm(.x, plot.it = FALSE)
    data[[".y"]] <- fData$x
    data[[".x"]] <- fData$y
    dform <- ".y ~ .x"
    if (!is.null(grp <- getGroupsFormula(form))) {
      dform <- paste(dform, deparse(grp[[2]]), sep = "|")
    }
    if (!is.null(id)) {			# identify points in plot
      id <-
        switch(mode(id),
               numeric = {
                 if ((id <= 0) || (id >= 1)) {
                   stop("'id' must be between 0 and 1")
                 }
                 as.logical(abs(resid(object, type=resType))
                            > -qnorm(id / 2))
               },
               call = eval(asOneSidedFormula(id)[[2]], data),
               stop("'id' can only be a formula or numeric")
               )
      if (is.null(idLabels)) {
        idLabels <- getGroups(object)
        if (length(idLabels) == 0) idLabels <- seq_len(object$dims$N)
        idLabels <- as.character(idLabels)
      } else {
        if (mode(idLabels) == "call") {
          idLabels <-
            as.character(eval(asOneSidedFormula(idLabels)[[2]], data))
        } else if (is.vector(idLabels)) {
          if (length(idLabels <- unlist(idLabels)) != length(id)) {
            stop("'idLabels' of incorrect length")
          }
          idLabels <- as.character(idLabels)
        } else {
          stop("'idLabels' can only be a formula or a vector")
        }
      }
    }
  } else { # case 2: ------- random.effects -------------------------------
    level <- attr(.x, "level")
    std <- attr(.x, "standardized")
    if (!is.null(effNams <- attr(.x, "effectNames"))) {
      .x <- .x[, effNams, drop = FALSE]
    }
    nc <- ncol(.x)
    nr <- nrow(.x)
    fData <- lapply(as.data.frame(.x), qqnorm, plot.it = FALSE)
    fData <- data.frame(.x = unlist(lapply(fData, function(x) x[["y"]])),
			.y = unlist(lapply(fData, function(x) x[["x"]])),
			.g = ordered(rep(names(fData),rep(nr, nc)),
				     levels = names(fData)))
    dform <- ".y ~ .x | .g"
    auxData <-
	if (!is.null(grp <- getGroupsFormula(form))) {
	    dform <- paste(dform, deparse(grp[[2]]), sep = "*")
	    data[is.na(match(names(data), "."))]
	} else
	    list()
    ## id and idLabels - need not be present
    if (!is.null(id)) {			# identify points in plot
      N <- object$dims$N
      id <-
        switch(mode(id),
               numeric = {
                 if ((id <= 0) || (id >= 1)) {
                   stop("'id' must be between 0 and 1")
                 }
                 aux <- ranef(object, standard = TRUE)
                 as.logical(abs(c(unlist(aux))) > -qnorm(id / 2))
               },
               call = eval(asOneSidedFormula(id)[[2]], data),
               stop("'id' can only be a formula or numeric")
               )
      if (length(id) == N) {
        ## id as a formula evaluated in data
        auxData[[".id"]] <- id
      }

      if (is.null(idLabels)) {
        idLabels <- rep(row.names(.x), nc)
      } else {
        if (mode(idLabels) == "call") {
          idLabels <-
            as.character(eval(asOneSidedFormula(idLabels)[[2]], data))
        } else if (is.vector(idLabels)) {
          if (length(idLabels <- unlist(idLabels)) != N) {
            stop("'idLabels' of incorrect length")
          }
          idLabels <- as.character(idLabels)
        } else {
          stop("'idLabels' can only be a formula or a vector")
        }
      }
      if (length(idLabels) == N) {
        ## idLabels as a formula evaluated in data
        auxData[[".Lid"]] <- idLabels
      }
    }
    data <-
        if (length(auxData)) {		# need collapsing
            auxData <- gsummary(data.frame(auxData),
                                groups = getGroups(object, level = level))
            auxData <- auxData[row.names(.x), , drop = FALSE]
            if (!is.null(auxData[[".id"]])) {
                id <- rep(auxData[[".id"]], nc)
            }
            if (!is.null(auxData[[".Lid"]])) {
                idLabels <- rep(auxData[[".Lid"]], nc)
            }
            cbind(fData, do.call("rbind", rep(list(auxData), nc)))
        } else {
            fData
        }
  }

  if(is.null(id)) id <- as.logical(as.character(id))
  idLabels <- as.character(idLabels)
  if (is.null(args$strip)) {
    args$strip <- function(...) strip.default(..., style = 1)
  }
  if (is.null(args$cex)) args$cex <- par("cex")
  if (is.null(args$adj)) args$adj <- par("adj")

  args <- c(list(eval(parse(text = dform)),
                 data = substitute(data),
                 panel = function(x, y, subscripts, ...){
                   x <- as.numeric(x)
                   y <- as.numeric(y)
                   dots <- list(...)
                   if (grid) panel.grid()
                   panel.xyplot(x, y, ...)
                   if (any(ids <- id[subscripts])) {
                       ltext(x[ids], y[ids], idLabels[subscripts][ids],
                             cex = dots$cex, adj = dots$adj)
                   }
                   if (!is.null(abline)) panel.abline(abline, ...)
                 }), args)
  if(type == "reff" && !std) {
    args[["scales"]] <- list(x = list(relation = "free"))
  }
  do.call(xyplot, args)
} ## {qqnorm.lmList}

ranef.lmList <-
  ##  Extracts the random effects from an lmList object.
  ##  If aug.frame is true, the returned data frame is augmented with a
  ##  values from the original data object, if available.  The variables
  ##  in the original data are collapsed over the cluster variable by the
  ##  function fun.
  function(object, augFrame = FALSE, data = NULL,
           which = NULL, FUN = mean, standard = FALSE,
           omitGroupingFactor = TRUE, ...)
{
  val <- coef(object, augFrame, data, which, FUN, omitGroupingFactor)
  effNames <- attr(val, "effectNames")
  effs <- val[, effNames, drop = FALSE]
  effs <-
    as.data.frame(lapply(effs, function(el) el - mean(el, na.rm = TRUE)))

  if(standard) {
    stdEff <- unlist(lapply(effs, function(el) sqrt(var(el[!is.na(el)]))))
    effs <- as.data.frame(as.matrix(effs) %*% diag(1/stdEff))
    attr(val, "label") <- "Standardized random effects"
  } else {
    attr(val, "label") <- "Random effects"
  }
  val[, effNames] <- effs
  attr(val, "standardized") <- standard
  class(val) <- unique(c("ranef.lmList", class(val)[-1]))
  val
}

residuals.lmList <-
  function(object, type = c("response", "pearson", "pooled.pearson"),
	   subset = NULL, asList = FALSE, ...)
{
  type <- match.arg(type)
  if(type == "pooled.pearson") {
    poolSD <- pooledSD(object)
  }
  if(!is.null(subset)) {
    if(is.character(subset)) {
      if (any(is.na(match(subset, names(object)))))
        stop("nonexistent groups requested in 'subset'")
    } else if (is.integer(subset)) {
      if (any(is.na(match(subset, seq_along(object)))))
        stop("nonexistent groups requested in 'subset'")
    } else
        stop("'subset' can only be character or integer")

    oclass <- class(object)
    oatt <- attr(object, "call")
    object <- object[subset]
    attr(object, "call") <- oatt
    class(object) <- oclass
  }
  val <-
    switch(type,
	   pooled.pearson = {
	     lapply(object, function(el, pSD)
	       if(!is.null(el)) resid(el)/pSD # else NULL
	     , pSD = poolSD)
	   },
	   pearson = lapply(object, function(el) {
	     if(!is.null(el)) {
	       aux <- resid(el)
	       aux/sqrt(sum(aux^2)/(length(aux) - length(coef(el))))
	     } # else NULL
	   }),
	   response = lapply(object, function(el)
               if(!is.null(el)) resid(el)) # else NULL
	   )
  if(!asList) {				#convert to array
    ngrps <- table(getGroups(object))[names(object)]
    if(any(aux <- vapply(object, is.null, NA))) {
      for(i in names(ngrps[aux])) {
	val[[i]] <- rep(NA, ngrps[i])
      }
    }
    val <- val[attr(object, "origOrder")] # putting in original order
    val <- setNames(unlist(val), rep(names(val), ngrps))
  }
  lab <-
    if (type == "response") {
      lab <- "Residuals"
      if(!is.null(aux <- attr(object, "units")$y)) paste(lab, aux) else lab
    }
    else "Standardized residuals"
  attr(val, "label") <- lab
  val
}

summary.lmList <-
  function(object, pool = attr(object, "pool"), ...)
{
    to.3d.array <-
        ## Convert the list to a 3d array watching for null elements
        function(lst, template) {
          if (!is.matrix(template))
            return(lst)

          ## Make empty array val[,,] and then fill it  -----
          dnames <- dimnames(template)
          use.i <- which(lengths(lst) > 0)
          ## TODO? just   identical(dnames[[1]], dnames[[2]]) :
          if (length(dnames[[1]]) == length(dnames[[2]]) &&
              all(dnames[[1]] == dnames[[2]])) { ## symmetric
            val <- array(NA_real_, dim=c(length(cfNms), length(cfNms), length(lst)),
                         dimnames=list(cfNms, cfNms, names(lst)))
            for (ii in use.i) {
              use <- dimnames(lst[[ii]])[[1]]
              val[use, use, ii] <- lst[[ii]]
              ##       ----
            }
          } else {
            val <- array(NA_real_, dim=c(length(cfNms), dim(template)[2], length(lst)),
                         dimnames=list(cfNms, dnames[[2]], names(lst)))
            for (ii in use.i) {
              use <- dimnames(lst[[ii]])[[1]]
              val[use, , ii] <- lst[[ii]]
              ##     ---
            }
          }
	  aperm(val, 3:1)
          ## val <- aperm(array(unlist(lapply(lst, function(el, template)
          ## 				 if(is.null(el)) { template }
          ## 				 else { el }, template = template)),
          ## 		   c(dim(template), length(lst)),
          ## 		   c(dnames, list(names(lst)))),
          ## 	     c(3, 2, 1))
          ## val[unlist(lapply(lst, is.null)), , ] <- NA
        }
    to.2d.array <-
        ## Convert the list to a 2d array watching for null elements
        function(lst, template)
        {
            if(is.null(template)) return(lst)
            template <- as.vector(template)
            val <- t(array(unlist(lapply(lst, function(el) if(is.null(el))
                                                             template else el)),
                           c(length(template), length(lst)),
                           list(names(template), names(lst))))
            val[vapply(lst, is.null, NA), ] <- NA
            val
        }
    ## Create a summary by applying summary to each component of the list
    sum.lst <- lapply(object, function(el) if(!is.null(el)) summary(el))
    nonNull <- !vapply(sum.lst, is.null, NA)
    if(!any(nonNull)) return(NULL)
    template <- sum.lst[[match(TRUE, nonNull)]] # the first one
    val <- as.list(setNames(nm = names(template)))
    for (i in names(template)) {
        val[[i]] <- lapply(sum.lst, `[[`, i)
        class(val[[i]]) <- "listof"
    }
    ## complete set of coefs [only used in to.3d.array()]
    cfNms <-
      unique(unlist(lapply(sum.lst[nonNull],
                              function(x) dimnames(x[['coefficients']])[[1]])))
    ## re-arrange the matrices into 3d arrays
    for(i in c("parameters", "cov.unscaled", "correlation", "coefficients"))
        if(length(val[[i]]))
            val[[i]] <- to.3d.array(val[[i]], template[[i]])
    ## re-arrange the vectors into 2d arrays
    for(i in c("df", "fstatistic"))
        val[[i]] <- to.2d.array(val[[i]], template[[i]])
    ## re-arrange the scalars into vectors
    for(i in c("sigma", "r.squared")) {
        ##    val[[i]] <- unlist(val[[i]]) - this deletes NULL components
        val[[i]] <- c(to.2d.array(val[[i]], template[[i]]))
    }
    ## select those attributes that do not vary with groups
    for(i in c("terms", "formula"))
        val[[i]] <- template[[i]]
    val[["call"]] <- attr(object, "call")
    if(inherits(object, "nlsList"))
        names(val[["call"]]["model"]) <- "object"
    val[["pool"]] <- pool
    if(pool) {
        poolSD <- pooledSD(object)
        dfRes <- attr(poolSD, "df")
        RSE <- c(poolSD)
        corRSE <- RSE/val$sigma
        pname <- if(inherits(object, "nlsList")) "parameters" else "coefficients"
        val[[pname]][,2,] <- val[[pname]][,2,] * corRSE
        val[[pname]][,3,] <- val[[pname]][,3,] / corRSE
        if(!inherits(object, "nlsList"))
            val[[pname]][,4,] <- 2*pt(abs(val[[pname]][,3,]), dfRes, lower.tail=FALSE)
        val[["df.residual"]] <- dfRes
        val[["RSE"]] <- RSE
    }
    attr(val, "groupsForm") <- attr(object, "groupsForm")
    class(val) <- "summary.lmList"
    val
}

# based on R's update.default
update.lmList <-
    function (object, formula., ..., evaluate = TRUE)
{
    call <- attr(object, "call")
    if (is.null(call))
	stop("need an object with call component")
    extras <- match.call(expand.dots = FALSE)$...
    if (!missing(formula.))
	call$object <- update.formula(formula(object), formula.)
    if(length(extras) > 0) {
	existing <- !is.na(match(names(extras), names(call)))
	## do these individually to allow NULL to remove entries.
	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
}

#update.lmList <-
#  function(object, formula, data, level, subset, na.action, pool, ...)
#{
#  thisCall <- as.list(match.call())[-(1:2)]
#  if (!missing(formula)) {
#    names(thisCall)[match(names(thisCall), "formula")] <- "object"
#  }
#  nextCall <- attr(object, "call")
#  nextCall[names(thisCall)] <- thisCall
#  if (!is.null(thisCall$object)) {
#    nextCall$object <- update(as.formula(nextCall$object), nextCall$object)
#  }
#  nextCall[[1]] <- quote(lmList)
#  eval(nextCall, envir = parent.frame(1))
#}