R/nlme.R

In nlme: Linear and Nonlinear Mixed Effects Models

###            Fit a general nonlinear mixed effects model
###
### Copyright 2006-2023  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/

nlme <-
  function(model,
	   data = sys.frame(sys.parent()),
	   fixed,
	   random = fixed,
           groups,
	   start,
           correlation = NULL,
           weights = NULL,
	   subset,
	   method = c("ML", "REML"),
	   na.action = na.fail,
	   naPattern,
	   control = list(),
	   verbose= FALSE)
{
  UseMethod("nlme")
}

nlme.nlsList <-
  function(model,
	   data = sys.frame(sys.parent()),
	   fixed,
	   random = fixed,
           groups,
	   start,
           correlation = NULL,
           weights = NULL,
	   subset,
	   method = c("ML", "REML"),
	   na.action = na.fail,
	   naPattern,
	   control = list(),
	   verbose= FALSE)
{
  ## control parameters
  controlvals <- nlmeControl()
  controlvals[names(control)] <- control

  thisCall <- as.list(match.call())[-1]
  ## checking the use of arguments defined within the function
  if (any(!is.na(match(names(thisCall),
		       c("fixed", "data", "start"))))) {
    warning("'nlme.nlsList' will redefine 'fixed', 'data', and 'start'")
  }
  method <- match.arg(method)
  REML <- method == "REML"
  ## add model, data, and optionally groups from the call that created model
  last.call <- as.list(attr(model, "call"))[-1]
  last.call$control <- NULL
  last.call$pool <- NULL
  thisCall[names(last.call)] <- last.call
  thisModel <- last.call[["model"]]
  thisCall[["model"]] <-
    eval(parse(text=paste( deparse (getResponseFormula (thisModel)[[2]]),
                          c_deparse(getCovariateFormula(thisModel)[[2]]),sep="~")))
  ## create "fixed" and "start"
  cf <- na.omit(coef(model))
  start <- list(fixed = unlist(lapply(cf, median, na.rm = TRUE)))
  pnames <- names(start$fixed) <- names(cf)
  thisCall[["fixed"]] <- lapply(as.list(pnames), function(el)
    eval(parse(text = paste(el, 1, sep = "~"))))
  if (missing(random)) {
    random <- thisCall[["fixed"]]
  }
  reSt <- reStruct(random, data = NULL)
  if (missing(groups)) {
    thisCall[["groups"]] <- groups <- getGroupsFormula(model)
  }
  if (length(reSt) > 1 || length(groups[[2]]) > 1) {
    stop("can only fit \"nlsList\" objects with single grouping variable")
  }
  ranForm <- formula(reSt)[[1]]
  if (!is.list(ranForm)) {
    ranForm <- list(ranForm)
  }
  mData <- thisCall[["data"]]
  if (is.null(mData)) {			# will try to construct
    allV <- unique(unlist(lapply(ranForm, function(el) all.vars(el[[3]]))))
    if (length(allV) > 0) {
      alist <- lapply(as.list(allV), as.name)
      names(alist) <- allV
      alist <- c(as.list(quote(data.frame)), alist)
      mode(alist) <- "call"
      mData <- eval(alist, sys.parent(1))
    }
  } else if (mode(mData) == "name" || mode(mData) == "call") {
    mData <- eval.parent(mData)
  }
  reSt <- reStruct(random, REML = REML, data = mData)
  names(reSt) <- deparse(groups[[2]])
  ## convert list of "name" objects to "character" vector
  rnames <- sapply(lapply(ranForm, `[[`, 2L), deparse)
  ## if the random effects are a subset of the coefficients,
  ## construct initial estimates for their var-cov matrix
  if (all(match(rnames, pnames, 0))) {
    madRes <- mad(resid(model), na.rm = TRUE)
    madRan <- unlist(lapply(cf, mad, na.rm = TRUE))
    madRan <- madRan[rnames]
    if (isInitialized(reSt)) {
      warning("initial value for 'reStruct' overwritten in 'nlme.nlsList'")
    }
    matrix(reSt) <- diag((madRan/madRes)^2, ncol = length(rnames))
  }
  thisCall[["start"]] <- start
  thisCall[["random"]] <- reSt
  val <- do.call(nlme.formula, thisCall, envir = parent.frame())
  val$origCall <- match.call()
  val
}


nlme.formula <-
  function(model,
	   data = sys.frame(sys.parent()),
	   fixed,
	   random,
           groups,
	   start,
           correlation = NULL,
           weights = NULL,
	   subset,
	   method = c("ML", "REML"),
	   na.action = na.fail,
	   naPattern,
	   control = list(),
	   verbose= FALSE)
{
  ## This is the method that actually does the fitting
  finiteDiffGrad <- function(model, data, pars) {
    dframe <- data.frame(data, pars)
    base <- eval(model, dframe)
    nm <- colnames(pars)
    grad <- array(base, c(length(base), length(nm)), list(NULL, nm))
    ssize <- sqrt(.Machine$double.eps)
    for (i in nm) {
      diff <- pp <- pars[ , i]
      diff[pp == 0] <- ssize
      diff[pp != 0] <- pp[pp != 0] * ssize
      dframe[[i]] <- pp + diff
      grad[ , i] <- (base - eval(model, dframe))/diff
      dframe[[i]] <- pp
    }
    grad
  }

  ## keeping the call
  Call <- match.call()

  ## control parameters
  controlvals <- nlmeControl()
  if (!missing(control)) {
    controlvals[names(control)] <- control
  }
  ##
  ## checking arguments
  ##
  if (!inherits(model, "formula"))
    stop("'model' must be a formula")
  if (length(model)!=3)
    stop("model formula must be of the form \"resp ~ pred\"")
  ## if (length(attr(terms(model), "offset")))
  ##   stop("offset() terms are not supported")

  method <- match.arg(method)
  REML <- method == "REML"
  if (missing(random)) {
    random <- fixed
  }
  reSt <- reStruct(random, REML = REML, data = NULL)
  if (missing(groups)) {
    groups <- getGroupsFormula(reSt)
  }
  if (is.null(groups)) {
    if (inherits(data, "groupedData")) {
      groups <- getGroupsFormula(data)
      namGrp <- rev(names(getGroupsFormula(data, asList = TRUE)))
      Q <- length(namGrp)
      if (length(reSt) != Q) { # may need to repeat reSt
	if (length(reSt) != 1)
	  stop("incompatible lengths for 'random' and grouping factors")
        randL <- vector("list", Q)
        names(randL) <- rev(namGrp)
        for(i in 1:Q) randL[[i]] <- random
        reSt <- reStruct(randL, REML = REML, data = NULL)
      }
    } else {
      ## will assume single group
      groups <- ~ 1
      names(reSt) <- namGrp <- "1"
    }
  } else {
    g.exp <- eval(parse(text = paste0("~1 |", deparse(groups[[2]]))))
    namGrp <- rev(names(getGroupsFormula(g.exp, asList = TRUE)))
  }
  names(reSt) <- namGrp
  ##
  ## checking if self-starting formula is given
  ##
  if (missing(start) && !is.null(attr(eval(model[[3]][[1]]), "initial"))) {
    nlmeCall <- Call
    nlsLCall <- nlmeCall[c("","model","data")]
    nlsLCall[[1]] <- quote(nlme::nlsList)
    nm <- names(nlmeCall)
    for(i in c("fixed", "data", "groups", "start"))
      if(i %in% nm) nlmeCall[[i]] <- NULL
    nlmeCall[[1]] <- quote(nlme::nlme.nlsList)
    ## checking if "data" is not equal to sys.frame(sys.parent())
    if (is.null(dim(data))) {
      stop("'data' must be given explicitly to use 'nlsList'")
    }
    nlsLObj <- eval.parent(nlsLCall)
    nlmeCall[["model"]] <- nlsLObj
    val <- eval.parent(nlmeCall)
    val$origCall <- NULL
    return(val)
  }
  if (is.numeric(start)) {               # assume it is the fixed effects
    start <- list(fixed = start)
  }
  nlmeModel <- call("-", model[[2]], model[[3]])
  ##
  ## save writing list(...) when only one element
  ##

  if (!is.list(fixed))
    fixed <- list(fixed)

  fixed <- do.call(c, lapply(fixed, function(fix.i) {
      if (is.name(fix.i[[2]]))
        list(fix.i)
      else
        ## multiple parameters on left hand side
        eval(parse(text =
		     paste0("list(",
			    paste(paste(all.vars(fix.i[[2]]),
					deparse (fix.i[[3]]), sep = "~"),
				  collapse = ","), ")")))
  }))
  fnames <- lapply(fixed, function(fix.i) {
    this <- eval(fix.i)
    if (!inherits(this, "formula"))
      stop ("'fixed' must be a formula or list of formulae")
    if (length(this) != 3)
      stop ("formulae in 'fixed' must be of the form \"parameter ~ expr\"")
    if (!is.name(this[[2]]))
      stop ("formulae in 'fixed' must be of the form \"parameter ~ expr\"")
    as.character(this[[2]])
  })
  names(fixed) <- fnames

  ranForm <- formula(reSt)              # random effects formula(s)
  Q <- length(ranForm)                  # number of groups
  names(ranForm) <- namGrp
  rnames <- vector("list", Q)
  names(rnames) <- namGrp
  for(i in 1:Q) {
    rnames[[i]] <- character(length(ranForm[[i]]))
    for (j in seq_along(ranForm[[i]])) {
      this <- eval(ranForm[[i]][[j]])
      if (!inherits(this, "formula"))
        stop ("'random' must be a formula or list of formulae")
      if (length(this) != 3)
        stop ("formulae in 'random' must be of the form \"parameter ~ expr\"")
      if (!is.name(this[[2]]))
        stop ("formulae in 'random' must be of the form \"parameter ~ expr\"")
      rnames[[i]][j] <- deparse(this[[2]])
    }
    names(ranForm[[i]]) <- rnames[[i]]
  }
  ## all parameter names
  pnames <- unique(c(fnames, unlist(rnames)))
  ##
  ##  If data is a pframe, copy the parameters in the frame to frame 1
  ## Doesn't exist in R
  ##  if (inherits(data, "pframe")) {
  ##    pp <- parameters(data)
  ##    for (i in names(pp)) {
  ##      assign(i, pp[[i]])
  ##    }
  ##    attr(data,"parameters") <- NULL
  ##    class(data) <- "data.frame"
  ##  }

  ## check if corStruct is present and assign groups to its formula,
  ## if necessary
  if (!is.null(correlation)) {
    if(!is.null(corGrpsForm <- getGroupsFormula(correlation, asList = TRUE))) {
      corGrpsForm <- unlist(lapply(corGrpsForm,
                                   function(el) deparse(el[[2]])))
      corQ <- length(corGrpsForm)
      lmeGrpsForm <- unlist(lapply(splitFormula(groups),
                                   function(el) deparse(el[[2]])))
      lmeQ <- length(lmeGrpsForm)
      if (corQ <= lmeQ) {
        if (any(corGrpsForm != lmeGrpsForm[1:corQ])) {
          stop("incompatible formulas for groups in \"random\" and \"correlation\"")
        }
        if (corQ < lmeQ) {
          warning("cannot use smaller level of grouping for \"correlation\" than for \"random\". Replacing the former with the latter.")
          frm <- paste("~",
                       c_deparse(getCovariateFormula(formula(correlation))[[2]]),
                       "|", deparse(groups[[2]]))
          attr(correlation, "formula") <- eval(parse(text = frm))
        }
      } else {
        if (any(lmeGrpsForm != corGrpsForm[1:lmeQ])) {
          stop("incompatible formulas for groups in \"random\" and \"correlation\"")
        }
      }
    } else {
      ## using the same grouping as in random
      frm <- paste("~",
                   c_deparse(getCovariateFormula(formula(correlation))[[2]]),
                   "|", deparse(groups[[2]]))
      attr(correlation, "formula") <- eval(parse(text = frm))
      corQ <- lmeQ <- 1
    }
  } else {
    corQ <- lmeQ <- 1
  }
  ## create an nlme structure containing the random effects model and plug-ins
  nlmeSt <- nlmeStruct(reStruct = reSt, corStruct = correlation,
                       varStruct = varFunc(weights))

  ## extract a data frame with enough information to evaluate
  ## form, fixed, random, groups, correlation, and weights
  mfArgs <- list(formula = asOneFormula(formula(nlmeSt), model, fixed,
                                        groups, omit = c(pnames, "pi")),
		 data = data, na.action = na.action)
  if (!missing(subset)) {
    mfArgs[["subset"]] <- asOneSidedFormula(Call[["subset"]])[[2]]
  }
  mfArgs$drop.unused.levels <- TRUE
  dataMix <- do.call(model.frame, mfArgs)

  origOrder <- row.names(dataMix)	# preserve the original order
  ##
  ## Evaluating the groups expression
  ##
  grps <- getGroups(dataMix,
                    eval(parse(text = paste("~1", deparse(groups[[2]]), sep = "|"))))
  N <- dim(dataMix)[1]			# number of observations
  ##
  ## evaluating the naPattern expression, if any
  ##
  if (missing(naPattern)) naPat <- rep(TRUE, N)
  else naPat <- as.logical(eval(asOneSidedFormula(naPattern)[[2]], dataMix))
  origOrderShrunk <- origOrder[naPat]

  ## ordering data by groups
  if (inherits(grps, "factor")) {	# single level
    ord <- order(grps)	#"order" treats a single named argument peculiarly
    grps <- data.frame(grps)
    row.names(grps) <- origOrder
    names(grps) <- as.character(deparse((groups[[2]])))
  } else {
    ord <- do.call(order, grps)
    ## making group levels unique
    for(i in 2:ncol(grps)) {
      grps[, i] <-
        as.factor(paste(as.character(grps[, i-1]), as.character(grps[,i]),
                        sep = "/"))
    }
  }
  if (corQ > lmeQ) {
    ## may have to reorder by the correlation groups
    ord <- do.call(order, getGroups(dataMix,
                                    getGroupsFormula(correlation)))
  }
  grps <- grps[ord, , drop = FALSE]
  dataMix <- dataMix[ord, ,drop = FALSE]
  ##  revOrder <- match(origOrder, row.names(dataMix)) # putting in orig. order
  naPat <- naPat[ord]			# ordered naPat
  dataMixShrunk <- dataMix[naPat, , drop=FALSE]
  ##  ordShrunk <- ord[naPat]
  grpShrunk <- grps[naPat,, drop = FALSE]
  revOrderShrunk <- match(origOrderShrunk, row.names(dataMixShrunk))
  yShrunk <- eval(model[[2]], dataMixShrunk)

  ##
  ## defining list with parameter information
  ##
  contr <- list()
  plist <- vector("list", length(pnames))
  names(plist) <- pnames
  for (nm in pnames) {
    this <- list(fixed = !is.null(fixed[[nm]]),
                 random = as.list(lapply(ranForm, function(el, nm)
                   !is.null(el[[nm]]), nm = nm)))
    if (this[["fixed"]]) {
      ## Peter Dalgaard claims the next line should be this[["fixed"]][[3]][[1]] != "1"
      ## but the current version seems to work ok.
      if (fixed[[nm]][[3]] != "1") {
        as1F <- asOneSidedFormula(fixed[[nm]][[3]])
	this[["fixed"]] <-
          model.matrix(as1F,
                       model.frame(as1F, dataMix))
        auxContr <- attr(this[["fixed"]], "contrasts")
        contr <- c(contr, auxContr[is.na(match(names(auxContr), names(contr)))])
      }
    }
    if (any(unlist(this[["random"]]))) {
      for(i in 1:Q) {
        wch <- which(!is.na(match(rnames[[i]], nm)))
        if (length(wch) == 1) {           # only one formula for nm at level i
          if ((rF.i <- ranForm[[i]][[nm]][[3]]) != "1") {
            this[["random"]][[i]] <-
              model.matrix(asOneSidedFormula(rF.i),
                           model.frame(asOneSidedFormula(rF.i), dataMix))
            auxContr <- attr(this[["random"]][[i]], "contrasts")
            contr <-
              c(contr, auxContr[is.na(match(names(auxContr), names(contr)))])
          }
        } else if (length(wch) > 0) {    # multiple formulas
          this[["random"]][[i]] <- th.ran.i <-
            as.list(lapply(ranForm[[i]][wch], function(el, data) {
              if (el[[3]] == "1") TRUE
              else model.matrix(asOneSidedFormula(el[[3]]),
                                model.frame(asOneSidedFormula(el[[3]]), data))
            }, data = dataMix))
          for(j in seq_along(th.ran.i)) {
            if (is.matrix(th.ran.i[[j]])) {
              auxContr <- attr(th.ran.i[[j]], "contrasts")
              contr <-
                c(contr, auxContr[is.na(match(names(auxContr), names(contr)))])
            }
          }
        }
      }
    }
    plist[[nm]] <- this
  }
  ## Ensure that all elements of are matrices
  contrMat <- function(nm, contr, data)
  {
    ## nm is a term in a formula, and can be a call
    x <- eval(parse(text=nm), data)
    levs <- levels(x)
    val <- do.call(contr[[nm]], list(n = length(levs)))
    rownames(val) <- levs
    val
  }
  nms <- names(contr)[sapply(contr, is.character)]
  contr[nms] <- lapply(nms, contrMat, contr = contr, data = dataMix)

  if (is.null(sfix <- start$fixed))
    stop ("'start' must have a component called 'fixed'")
  ##
  ## Fixed effects names
  ##
  fn <- character(0)
  currPos <- 0
  fixAssign <- list()
  for(nm in fnames) {
    if (is.logical(f <- plist[[nm]]$fixed)) {
      currPos <- currPos + 1
      currVal <- list(currPos)
      if (all(unlist(lapply(plist[[nm]]$random, is.logical)))) {
        fn <- c(fn, nm)
        names(currVal) <- nm
      } else {
        aux <- paste(nm, "(Intercept)", sep=".")
        fn <- c(fn, aux)
        names(currVal) <- aux
      }
      fixAssign <- c(fixAssign, currVal)
    } else {
      currVal <- attr(f, "assign")
      fTerms <- terms(asOneSidedFormula(fixed[[nm]][[3]]), data=data)
      namTerms <- attr(fTerms, "term.labels")
      if (attr(fTerms, "intercept") > 0) {
        namTerms <- c("(Intercept)", namTerms)
      }
      namTerms <- factor(currVal, labels = namTerms)
      currVal <- split(order(currVal), namTerms)
      names(currVal) <- paste(nm, names(currVal), sep = ".")
      fixAssign <- c(fixAssign, lapply(currVal,
                                       function(el) el + currPos ))
      currPos <- currPos + length(unlist(currVal))
      fn <- c(fn, paste(nm, colnames(f), sep = "."))
    }
  }
  fLen <- length(fn)
  if (fLen == 0L || length(sfix) != fLen)
    stop ("starting values for the 'fixed' component are not the correct length")
  names(sfix) <- fn
  ##
  ## Random effects names
  ##
  rn <- wchRnames <- vector("list", Q)
  names(rn) <- names(wchRnames) <- namGrp
  for(i in 1:Q) {
    rn[[i]] <- character(0)
    uRnames <- unique(rnames[[i]])
    wchRnames[[i]] <- integer(length(uRnames))
    names(wchRnames[[i]]) <- uRnames
    for(j in seq_along(rnames[[i]])) {
      nm <- rnames[[i]][j]
      wchRnames[[i]][nm] <- wchRnames[[i]][nm] + 1
      r <- plist[[nm]]$random[[i]]
      if(is.list(r)) r <- r[[wchRnames[[i]][nm]]]
      if (is.logical(r)) {
        if (r) {
          rn[[i]] <- c(rn[[i]],
                       if (is.logical(plist[[nm]]$fixed)) nm
                       else paste(nm,"(Intercept)",sep="."))
        } # else: keep rn[[i]]
      } else {
        rn[[i]] <- c(rn[[i]], paste(nm, colnames(r), sep = "."))
      }
    }
  }
  Names(nlmeSt$reStruct) <- rn
  rNam <- unlist(rn)                    # unlisted names of random effects
  rlength <- lengths(rn)                # number of random effects per stratum
  rLen <- sum(rlength)                  # total number of random effects
  pLen <- rLen + fLen                   # total number of parameters
  ncols <- c(rlength, fLen, 1)
  Dims <- MEdims(grpShrunk, ncols)
  if (max(Dims$ZXlen[[1]]) < Dims$qvec[1]) {
    warning(gettextf("fewer observations than random effects in all level %s groups",
                     Q), domain = NA)
  }
  sran <- vector("list", Q)
  names(sran) <- namGrp
  if (!is.null(sran0 <- start$random)) {
    if (inherits(sran0, "data.frame")) {
      sran0 <- list(as.matrix(sran0))
    } else {
      if (!is.list(sran0)) {
        if (!is.matrix(sran0)) {
          stop("starting values for random effects should be a list, or a matrix")
        }
        sran0 <- list(as.matrix(sran0))
      }
    }
    if (is.null(namSran <- names(sran0))) {
      if (length(sran) != Q) {
        stop(gettextf("list with starting values for random effects must have names or be of length %d",
                      Q), domain = NA)
      }
      names(sran0) <- rev(namGrp)        # assume given in outer-inner order
    } else {
      if (any(noMatch <- is.na(match(namSran, namGrp)))) {
        stop(sprintf(ngettext(sum(noMatch),
                              "group name not matched in starting values for random effects: %s",
                              "group names not matched in starting values for random effects: %s"),
                     paste(namSran[noMatch], collapse=", ")), domain = NA)
      }
    }
  }
  for(i in 1:Q) {
    if (is.null(sran[[i]] <- sran0[[namGrp[i]]])) {
      sran[[i]] <- array(0, c(rlength[i], Dims$ngrps[i]),
                         list(rn[[i]], unique(as.character(grps[, Q-i+1]))))
    } else {
      if (!is.matrix(sran[[i]]))
        stop("starting values for the random components should be a list of matrices")
      dimsran <- dim(sran[[i]])
      if (dimsran[1] != Dims$ngrps[i]) {
        stop(gettextf("number of rows in starting values for random component at level %s should be %d",
                      namGrp[i], Dims$ngrps[i]), domain = NA)
      }
      if (dimsran[2] != rlength[i]) {
        stop(gettextf("number of columns in starting values for random component at level %s should be %d",
                      namGrp[i], rlength[i]), domain = NA)
      }
      dnamesran <- dimnames(sran[[i]])
      if (is.null(dnamesran[[1]])) {
        stop("starting values for random effects must include group levels")
      } else {
        levGrps <- unique(as.character(grps[, Q-i+1]))
        if(!all(sort(dnamesran[[1]]) == sort(levGrps))) {
          stop(gettextf("groups levels mismatch in 'random' and starting values for 'random' at level %s",
                        namGrp[i]), domain = NA)
        }
        sran[[i]] <- sran[[i]][levGrps, , drop = FALSE]
      }

      if (!is.null(dnamesran[[2]])) {
        if(!all(sort(dnamesran[[2]]) == sort(rn[[i]]))) {
          ## first try to resolve it
          for(j in seq_len(rlength[i])) {
            if (is.na(match(dnamesran[[2]][j], rn[[i]]))) {
              if (!is.na(mDn <- match(paste(dnamesran[[2]][j],
                                            "(Intercept)", sep="."), rn[[i]]))) {
                dnamesran[[2]][j] <- rn[[i]][mDn]
              } else {
                if (!is.na(mDn <- match(dnamesran[[2]][j],
                                        paste(rn[[i]], "(Intercept)", sep = ".")))) {
                  dnamesran[[2]][j] <- rn[[i]][mDn]
                } else {
                  stop (gettextf("names mismatch in 'random' and starting values  for 'random' at level %s",
                                 namGrp[i]), domain = NA)
                }
              }
            }
          }
          dimnames(sran[[i]]) <- dnamesran
        }
        sran[[i]] <- sran[[i]][, rn[[i]], drop = FALSE]
      } else {
        dimnames(sran[[i]])[[2]] <- rn[[i]]
      }
      sran[[i]] <- t(sran[[i]])
    }
  }
  names(sran) <- namGrp
  nPars <- length(unlist(sran)) + fLen  # total number of PNLS parameters
  ##
  ##   defining values of constants used in calculations
  ##
  NReal <- sum(naPat)
  ##
  ## Creating the fixed and random effects maps
  ##
  fmap <- list()
  n1 <- 1
  for(nm in fnames) {
    if (is.logical(f <- plist[[nm]]$fixed)) {
      fmap[[nm]] <- n1
      n1 <- n1 + 1
    } else {
      fmap[[nm]] <- n1:(n1+ncol(f) - 1)
      n1 <- n1 + ncol(f)
    }
  }
  rmap <- rmapRel <- vector("list", Q)
  names(rmap) <- names(rmapRel) <- namGrp
  n1 <- 1
  startRan <- 0
  for(i in 1:Q) {
    wchRnames[[i]][] <- 0
    rmap[[i]] <- rmapRel[[i]] <- list()
    for(nm in rnames[[i]]) {
      wchRnames[[i]][nm] <- wchRnames[[i]][nm] + 1
      r <- plist[[nm]]$random[[i]]
      if(is.list(r)) r <- r[[wchRnames[[i]][nm]]]
      if (is.logical(r)) {
        val <- n1
        n1 <- n1 + 1
      } else {
        val <- n1:(n1+ncol(r) - 1)
        n1 <- n1 + ncol(r)
      }
      if (is.null(rmap[[i]][[nm]])) {
        rmap[[i]][[nm]] <- val
        rmapRel[[i]][[nm]] <- val - startRan
      } else {
        rmap[[i]][[nm]] <- c(rmap[[i]][[nm]], list(val))
        rmapRel[[i]][[nm]] <- c(rmapRel[[i]][[nm]], list(val - startRan))
      }
    }
    startRan <- startRan + ncols[i]
  }

  ##
  ## defining the nlFrame, i.e., nlEnv, an environment in R :
  ##
  grpsRev <- rev(lapply(grps, as.character))
  bmap <- c(0, cumsum(unlist(lapply(sran, function(el) length(as.vector(el))))))
  nlEnv <- list2env(
    list(model = nlmeModel,
         data = dataMix,
         groups = grpsRev,
         plist = plist,
         beta = as.vector(sfix),
         bvec = unlist(sran),
         b = sran,
         X = array(0, c(N, fLen), list(NULL, fn)),
         Z = array(0, c(N, rLen), list(NULL, rNam)),
         fmap = fmap,
         rmap = rmap,
         rmapRel = rmapRel,
         bmap = bmap,
         level = Q,
         N = N,
         Q = Q,
         naPat = naPat,
         .parameters = c("bvec", "beta"),
         finiteDiffGrad = finiteDiffGrad))

  modelExpression <- ~ {
    pars <- getParsNlme(plist, fmap, rmapRel, bmap, groups, beta, bvec, b, level, N)
    res <- eval(model, data.frame(data, pars))
    if (!length(grad <- attr(res, "gradient"))) {
      grad <- finiteDiffGrad(model, data, pars)
    }
    for (nm in names(plist)) {
      gradnm <- grad[, nm]
      if (is.logical(f <- plist[[nm]]$fixed)) {
        if(f)
          X[, fmap[[nm]]] <- gradnm # else f == FALSE =^= 0
      } else
        X[, fmap[[nm]]] <- gradnm * f

      for(i in 1:Q) {
        if (is.logical(r <- plist[[nm]]$random[[i]])) {
          if (r)
            Z[, rmap[[i]][[nm]]] <- gradnm  # else r == FALSE =^= 0
        } else {
          rm.i <- rmap[[i]][[nm]]
          if (!is.list(rm.i)) {
            Z[, rm.i] <- gradnm * r
          } else {
            for(j in seq_along(rm.i)) {
              Z[, rm.i[[j]]] <-
                if (is.logical(rr <- r[[j]]))
                  gradnm
                else
                  gradnm * rr
            }
          }
        }
      }
    }
    result <- c(Z[naPat, ], X[naPat, ], res[naPat])
    result[is.na(result)] <- 0
    result
  }## {modelExpression}

  modelResid <-
    ~ eval(model,
           data.frame(data,
                      getParsNlme(plist, fmap, rmapRel, bmap, groups,
                                  beta, bvec, b, level, N)))[naPat]
  ww <- eval(modelExpression[[2]], envir = nlEnv)
  w <- ww[NReal * pLen + (1:NReal)]
  ZX <- array(ww[1:(NReal*pLen)], c(NReal, pLen),
              list(row.names(dataMixShrunk), c(rNam, fn)))
  w <- w + as.vector(ZX[, rLen + (1:fLen), drop = FALSE] %*% sfix)
  if (!is.null(start$random)) {
    startRan <- 0
    for(i in 1:Q) {
      w <- w + as.vector((ZX[, startRan + 1:ncols[i], drop = FALSE] *
                          t(sran[[i]])[as.character(grpShrunk[, Q-i+1]),,drop = FALSE]) %*%
                         rep(1, ncols[i]))
      startRan <- startRan + ncols[i]
    }
  }
  ## creating the condensed linear model
  attr(nlmeSt, "conLin") <-
    list(Xy = array(c(ZX, w), dim = c(NReal, sum(ncols)),
                    dimnames = list(row.names(dataMixShrunk),
                                    c(colnames(ZX), deparse(model[[2]])))),
	 dims = Dims, logLik = 0,
	 ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions
	 sigma = controlvals$sigma, auxSigma = 0)

  ## additional attributes of nlmeSt
  attr(nlmeSt, "resp") <- yShrunk
  attr(nlmeSt, "model") <- modelResid
  attr(nlmeSt, "local") <- nlEnv
  attr(nlmeSt, "NReal") <- NReal
  ## initialization
  nlmeSt <- Initialize(nlmeSt, dataMixShrunk, grpShrunk,
                       control = controlvals)
  parMap <- attr(nlmeSt, "pmap")

  decomp <- length(coef(nlmeSt)) == length(coef(nlmeSt$reStruct)) && !needUpdate(nlmeSt)
  if(decomp) { # can do one decomposition
    ## need to save conLin for calculating updating between steps
    oldConLin <- attr(nlmeSt, "conLin")
  }

  numIter <- 0				# number of iterations
  pnlsSettings <- c(controlvals$pnlsMaxIter, controlvals$minScale,
                    controlvals$pnlsTol, 0, 0, 0)
  nlModel <- nonlinModel(modelExpression, nlEnv)
  ##----------------------------------------------------------------------------
  repeat {                              ## alternating algorithm
    numIter <- numIter + 1
    ## LME step
    if (needUpdate(nlmeSt)) {             # updating varying weights
      nlmeSt <- update(nlmeSt, dataMixShrunk)
    }
    if (decomp) {
      attr(nlmeSt, "conLin") <- MEdecomp(oldConLin)
    }
    oldPars <- coef(nlmeSt)
    if (controlvals$opt == "nlminb") {
      control <-  list(trace = controlvals$msVerbose,
                       iter.max = controlvals$msMaxIter)
      keep <- c("eval.max", "abs.tol", "rel.tol", "x.tol", "xf.tol",
                "step.min", "step.max", "sing.tol", "scale.init", "diff.g")
      control <- c(control, controlvals[names(controlvals) %in% keep])
      optRes <- nlminb(c(coef(nlmeSt)),
                       function(nlmePars) -logLik(nlmeSt, nlmePars),
                       control = control)
      aConv <- coef(nlmeSt) <- optRes$par
      convIter <- optRes$iterations
      if(optRes$convergence && controlvals$msWarnNoConv)
        warning(paste(sprintf(
          "Iteration %d, LME step: nlminb() did not converge (code = %d).",
          numIter, optRes$convergence),
          if(convIter >= control$iter.max) "Do increase 'msMaxIter'!"
          else if(!is.null(msg <- optRes$message)) paste("PORT message:", msg)))
    } else { ## nlm(.)
      aNlm <- nlm(f = function(nlmePars) -logLik(nlmeSt, nlmePars),
                  p = c(coef(nlmeSt)), hessian = TRUE,
                  print.level = controlvals$msVerbose,
                  gradtol = if(numIter == 1) controlvals$msTol
                            else 100*.Machine$double.eps,
                  iterlim = if(numIter < 10) 10 else controlvals$msMaxIter,
                  check.analyticals = FALSE)
      aConv <- coef(nlmeSt) <- aNlm$estimate
      convIter <- aNlm$iterations
      if(aNlm$code > 2 && controlvals$msWarnNoConv)
        warning(paste(sprintf(
          "Iteration %d, LME step: nlm() did not converge (code = %d).",
          numIter, aNlm$code),
          if(aNlm$code == 4) "Do increase 'msMaxIter'!"))
    }
    nlmeFit <- attr(nlmeSt, "lmeFit") <- MEestimate(nlmeSt, grpShrunk)
    if (verbose) {
      cat("\n**Iteration", numIter)
      cat(sprintf("\nLME step: Loglik: %s, %s iterations: %d\n",
                  format(nlmeFit$logLik), controlvals$opt, convIter))
      print(nlmeSt)
    }

    ## PNLS step
    if (is.null(correlation)) {
      cF <- 1.0
      cD <- 1
    } else {
      cF <- corFactor(nlmeSt$corStruct)
      cD <- Dim(nlmeSt$corStruct)
    }
    vW <- if(is.null(weights)) 1.0 else varWeights(nlmeSt$varStruct)

    if (verbose) cat(" Beginning PNLS step: .. ")
    work <- .C(fit_nlme,
	       thetaPNLS = as.double(c(as.vector(unlist(sran)), sfix)),
               pdFactor = as.double(pdFactor(nlmeSt$reStruct)),
               as.integer(unlist(rev(grpShrunk))),
	       as.integer(unlist(Dims)),
               as.integer(attr(nlmeSt$reStruct, "settings"))[-(1:3)],
	       as.double(cF),
	       as.double(vW),
               as.integer(unlist(cD)),
	       settings = as.double(pnlsSettings),
	       additional = double(NReal * (pLen + 1)),
	       as.integer(!is.null(correlation)),
	       as.integer(!is.null(weights)),
               ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions
	       as.double(controlvals$sigma),
               nlModel,
	       NAOK = TRUE)
    if (verbose) cat(" completed fit_nlme() step.\n")
    if (work$settings[4] == 1) {
      ##      convResult <- 2
      msg <- "step halving factor reduced below minimum in PNLS step"
      if (controlvals$returnObject)
        warning(msg)
      else
        stop(msg)
    }
    ## dim(work$pdFactor) <- dim(pdMatrix(nlmeSt$reStruct[[1]]))
    ## matrix(nlmeSt$reStruct[[1]]) <- crossprod(work$pdFactor)
    ## fix from Setzer.Woodrow@epamail.epa.gov for nested grouping factors
    i.pdF <- 1
    for (i in seq_along(nlmeSt$reStruct)) {
      d.i <- dim(pdMatrix(nlmeSt$reStruct[[i]]))
      ni.pdF <- i.pdF + prod(d.i)
      pdF <- array(work$pdFactor[i.pdF:(ni.pdF -1)], dim = d.i)
      matrix(nlmeSt$reStruct[[i]]) <- crossprod(pdF)
      i.pdF <- ni.pdF
    }
    oldPars <- c(sfix, oldPars)
    for(i in 1:Q) sran[[i]][] <- work$thetaPNLS[(bmap[i]+1):bmap[i+1]]
    sfix[] <- work$thetaPNLS[nPars + 1 - (fLen:1)]
    if (verbose) {
      cat("PNLS step: RSS = ", format(work$settings[6]), "\n fixed effects: ")
      for (i in 1:fLen) cat(format(sfix[i])," ")
      cat("\n iterations:", work$settings[5],"\n")
    }
    aConv <- c(sfix, coef(nlmeSt)) # 2nd part added by SDR 04/19/2002
    w[] <- work$additional[(NReal * pLen) + 1:NReal]
    ZX[] <- work$additional[1:(NReal * pLen)]
    w <- w + as.vector(ZX[, rLen + (1:fLen), drop = FALSE] %*% sfix)
    startRan <- 0
    for(i in 1:Q) {
      gr.i <- as.character(grpShrunk[, Q-i+1])
      w <- w + as.vector((ZX[, startRan + 1:ncols[i], drop = FALSE] *
                          t(sran[[i]])[gr.i,, drop = FALSE]) %*%
                         rep(1, ncols[i]))
      startRan <- startRan + ncols[i]
    }
    if (decomp) {
      oldConLin$Xy[] <- c(ZX, w)
      oldConLin$logLik <- 0
    } else {
      attr(nlmeSt, "conLin")$Xy[] <- c(ZX, w)
      attr(nlmeSt, "conLin")$logLik <- 0
    }

    conv <- abs((oldPars - aConv)/
                ifelse(abs(aConv) < controlvals$tolerance, 1, aConv))
    aConv <- c(fixed = max(conv[1:fLen]))
    conv <- conv[-(1:fLen)]
    for(i in names(nlmeSt)) {
      if (any(parMap[,i])) {
	aConv <- c(aConv, max(conv[parMap[,i]]))
	names(aConv)[length(aConv)] <- i
      }
    }

    if (verbose) {
      cat(sprintf("Convergence crit. (must all become <= tolerance = %g):\n",
                  controlvals$tolerance))
      print(aConv)
    }

    if ((max(aConv) <= controlvals$tolerance) ||
        (aConv["fixed"] <= controlvals$tolerance && convIter == 1)) {
      ##      convResult <- 0
      break
    }
    if (numIter >= controlvals$maxIter) {
      ##      convResult <- 1
      msg <- gettextf(
	"maximum number of iterations (maxIter = %d) reached without convergence",
	controlvals$maxIter)
      if (controlvals$returnObject) {
	warning(msg, domain=NA) ; break
      } else
	stop(msg, domain=NA)
    }
  } ## end{ repeat } (nlme steps) ----------------------------------------------

  ## wrapping up
  nlmeFit <-
    if (decomp)
      MEestimate(nlmeSt, grpShrunk, oldConLin)
    else
      MEestimate(nlmeSt, grpShrunk)
  ## degrees of freedom for fixed effects tests
  fixDF <- getFixDF(ZX[, rLen + (1:fLen), drop = FALSE],
                    grpShrunk, attr(nlmeSt, "conLin")$dims$ngrps, fixAssign)

  ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions
  attr(fixDF, "varFixFact") <- varFix <- nlmeFit$sigma * nlmeFit$varFix
  varFix <- crossprod(varFix)

  dimnames(varFix) <- list(fn, fn)
  ##
  ## fitted.values and residuals (in original order)
  ##
  Resid <-
    if (decomp)
      resid(nlmeSt, level = 0:Q, oldConLin)[revOrderShrunk, ]
    else
      resid(nlmeSt, level = 0:Q)[revOrderShrunk, ]
  Fitted <- yShrunk[revOrderShrunk] - Resid
  rownames(Resid) <- rownames(Fitted) <- origOrderShrunk
  grpShrunk <- grpShrunk[revOrderShrunk, , drop = FALSE]
  attr(Resid, "std") <- nlmeFit$sigma/(varWeights(nlmeSt)[revOrderShrunk])
  ## inverting back reStruct
  nlmeSt$reStruct <- solve(nlmeSt$reStruct)
  attr(nlmeSt, "fixedSigma") <- (controlvals$sigma > 0)
  ## saving part of dims
  dims <- attr(nlmeSt, "conLin")$dims[c("N", "Q", "qvec", "ngrps", "ncol")]
  ## getting the approximate var-cov of the parameters
  apVar <-
    if (controlvals$apVar)
      lmeApVar(nlmeSt, nlmeFit$sigma,
               .relStep = controlvals[[".relStep"]],
               minAbsPar = controlvals[["minAbsParApVar"]],
               natural = controlvals[["natural"]])
    else
      "Approximate variance-covariance matrix not available"
  ## putting sran in the right format
  sran <- lapply(sran, t)
  ## getting rid of condensed linear model, fit, and other attributes
###- oClass <- class(nlmeSt)
  attributes(nlmeSt) <- attributes(nlmeSt)[
    c("names", "class", "pmap", "fixedSigma")]
##- class(nlmeSt) <- oClass
  ##
  ## creating the  nlme object
  ##
  isGrpd <- inherits(data, "groupedData")
  structure(class = c("nlme", "lme"),
            list(modelStruct = nlmeSt,
		 dims = dims,
                 contrasts = contr,
		 coefficients = list(fixed = sfix, random = rev(sran)),
		 varFix = varFix,
		 sigma = nlmeFit$sigma,
		 apVar = apVar,
		 logLik = nlmeFit$logLik,
		 numIter = numIter,
		 groups = grpShrunk,
		 call = Call,
		 method = method,
		 fitted = Fitted,
		 residuals = Resid,
		 plist = plist,
                 map = list(fmap=fmap,rmap=rmap,rmapRel=rmapRel,bmap=bmap),
                 fixDF = fixDF
               , formula = model # => reliable formula.nlme() [PR#18599]
                 ),
            ## saving labels and units for plots
            units  = if(isGrpd) attr(data, "units"),
            labels = if(isGrpd) attr(data, "labels"))
} ## {nlme.formula}

###
##' Calculate the parameters from the fixed and random effects
getParsNlme <-
  function(plist, fmap, rmapRel, bmap, groups, beta, bvec, b, level, N)
{
  pars <- array(0, c(N, length(plist)), list(NULL, names(plist)))
  ## for random effects below
  iQ <- if (level > 0) {
    Q <- length(groups)
    (Q - level + 1L):Q
  } else integer() # empty

  for (nm in names(plist)) {
    ## 1) Fixed effects
    if (is.logical(f <- plist[[nm]]$fixed)) {
      if (f)
        pars[, nm] <- beta[fmap[[nm]]]
      ## else pars[, nm] <- 0  (as f == FALSE)
    } else
      pars[, nm] <- f %*% beta[fmap[[nm]]]

    ## 2) Random effects
    for(i in iQ)
      if(!is.null(rm.i. <- rmapRel[[i]][[nm]])) {
        b.i <- b[[i]]
        b.i[] <- bvec[(bmap[i] + 1):bmap[i+1]]
        ## NB: some groups[[i]] may be *new* levels, i.e. non-matching:
        gr.i <- match(groups[[i]], colnames(b.i)) # column numbers + NA
        if (is.logical(r <- plist[[nm]]$random[[i]])) {
          if (r)
            pars[, nm] <- pars[, nm] + b.i[rm.i., gr.i]
          ## else r == FALSE =^= 0
        } else if (!is.list(r)) {
          pars[, nm] <- pars[, nm] +
            (r * t(b.i)[gr.i, rm.i., drop=FALSE]) %*% rep(1, ncol(r))
        } else {
          b.i.gi <- b.i[, gr.i, drop=FALSE]
          for(j in seq_along(rm.i.)) {
            if (is.logical(rr <- r[[j]])) {
              if(rr)
                pars[, nm] <- pars[, nm] + b.i.gi[rm.i.[[j]], ]
              ## else rr == FALSE =^= 0
            }
            else
              pars[, nm] <- pars[, nm] +
                (rr * t(b.i.gi[rm.i.[[j]], , drop=FALSE])) %*% rep(1, ncol(rr))
          }
        }
      } # for( i ) if(!is.null(rm.i. ..))
  }
  pars
}

###
###  Methods for standard generics
###

formula.nlme <- function(x, ...) x$formula %||% eval(x$call[["model"]])

predict.nlme <-
  function(object, newdata, level = Q, asList = FALSE, na.action = na.fail,
	   naPattern = NULL, ...)
{
  ##
  ## method for predict() designed for objects inheriting from class nlme
  ##
  Q <- object$dims$Q
  if (missing(newdata)) {		# will return fitted values
    val <- fitted(object, level, asList)
    if (length(level) == 1) return(val)
    return(data.frame(object[["groups"]][,level[level != 0], drop = FALSE],
		      predict = val))
  }
  maxQ <- max(level)			# maximum level for predictions
  nlev <- length(level)
  newdata <- as.data.frame(newdata)
  if (maxQ > 0) {			# predictions with random effects
    whichQ <- Q - (maxQ-1):0
    reSt <- object$modelStruct$reStruct[whichQ]
    ##    nlmeSt <- nlmeStruct(reStruct = reSt)
    groups <- getGroupsFormula(reSt)
    if (any(is.na(match(all.vars(groups), names(newdata))))) {
      ## groups cannot be evaluated in newdata
      stop("cannot evaluate groups for desired levels on 'newdata'")
    }
  } else {
    reSt <- NULL
  }

  ## use xlev to make sure factor levels are the same as in contrasts
  ## and to support character-type 'newdata' for factors
  contr <- object$contrasts
  dataMix <- model.frame(formula = asOneFormula(
                   formula(object),
                   object$call$fixed, formula(reSt), naPattern,
                   omit = c(names(object$plist), "pi",
                            deparse(getResponseFormula(object)[[2]]))),
                 data = newdata, na.action = na.action,
                 drop.unused.levels = TRUE,
                 xlev = lapply(contr, rownames))
  origOrder <- row.names(dataMix)	# preserve the original order
  whichRows <- match(origOrder, row.names(newdata))

  if (maxQ > 0) {
    ## sort the model.frame by groups and get the matrices and parameters
    ## used in the estimation procedures
    grps <- getGroups(newdata,          # (unused levels are dropped here)
                      eval(substitute(~ 1 | GRPS,
                                      list(GRPS = groups[[2]]))))
    ## ordering data by groups
    if (inherits(grps, "factor")) {	# single level
      grps <- grps[whichRows]
      oGrps <- data.frame(grps)
      ## checking if there are missing groups
      if (any(naGrps <- is.na(grps))) {
	grps[naGrps] <- levels(grps)[1L]	# input with existing level
      }
      ord <- order(grps)     #"order" treats a single named argument peculiarly
      grps <- data.frame(grps)
      row.names(grps) <- origOrder
      names(grps) <- names(oGrps) <- as.character(deparse((groups[[2L]])))
    } else {
      grps <- oGrps <- grps[whichRows, ]
      ## checking for missing groups
      if (any(naGrps <- is.na(grps))) {
	## need to input missing groups
	for(i in names(grps)) {
	  grps[naGrps[, i], i] <- levels(grps[,i])[1L]
	}
	naGrps <- t(apply(naGrps, 1, cumsum)) # propagating NAs
      }
      ord <- do.call(order, grps)
      ## making group levels unique
      for(i in 2:ncol(grps)) {
	grps[, i] <-
          as.factor(paste(as.character(grps[, i-1]),
                          as.character(grps[, i  ]), sep = "/"))
      }
    }
    ## if (match(0, level, nomatch = 0)) {
    ##   naGrps <- cbind(FALSE, naGrps)
    ## }
    ## naGrps <- as.matrix(naGrps)[ord, , drop = FALSE]
    naGrps <- cbind(FALSE, naGrps)[ord, , drop = FALSE]
    grps <- grps[ord, , drop = FALSE]
    dataMix <- dataMix[ord, ,drop = FALSE]
    revOrder <- match(origOrder, row.names(dataMix)) # putting in orig. order
  }
  ## restore contrasts for the below model.matrix() calls (which may use
  ## different factor variables, so we avoid passing the whole contr list)
  for(i in intersect(names(dataMix), names(contr))) {
    attr(dataMix[[i]], "contrasts") <- contr[[i]]
  }

  N <- nrow(dataMix)
  ##
  ## evaluating the naPattern expression, if any
  ##
  naPat <- if(is.null(naPattern)) rep(TRUE, N)
           else
             as.logical(eval(asOneSidedFormula(naPattern)[[2]], dataMix))
  ##
  ## Getting  the plist for the new data frame
  ##
  plist <- object$plist
  fixed <- eval(object$call$fixed)
  if (!is.list(fixed))
    fixed <- list(fixed)

  fixed <- do.call(c, lapply(fixed, function(fix.i) {
      if (is.name(fix.i[[2]]))
        list(fix.i)
      else
        ## multiple parameters on left hand side
        eval(parse(text = paste0("list(", paste(paste(all.vars(fix.i[[2]]),
                                                      deparse (fix.i[[3]]),
                                                      sep = "~"),
                                                collapse = ","), ")")))
  }))
  fnames <- unlist(lapply(fixed, function(el) deparse(el[[2]])))
  names(fixed) <- fnames
  fix <- fixef(object)
  ##  fn <- names(fix)
  for(nm in fnames) {
    if (!is.logical(plist[[nm]]$fixed)) {
      oSform <- asOneSidedFormula(fixed[[nm]][[3]])
      plist[[nm]]$fixed <- model.matrix(oSform, model.frame(oSform, dataMix))
    }
  }

  if (maxQ > 0) {
    grpsRev <- lapply(rev(grps), as.character)
    ranForm <- formula(reSt)[whichQ]
    namGrp <- names(ranForm)
    rnames <- lapply(ranForm, function(el)
      unlist(lapply(el, function(el1) deparse(el1[[2]]))))
    for(i in seq_along(ranForm)) {
      names(ranForm[[i]]) <- rnames[[i]]
    }
    ran <- ranef(object)
    ran <- if(is.data.frame(ran)) list(ran) else rev(ran)
    ##    rn <- lapply(ran[whichQ], names)
    ran <- lapply(ran, t)
    ranVec <- unlist(ran)
    for(nm in names(plist)) {
      for(i in namGrp) {
        if (!is.logical(plist[[nm]]$random[[i]])) {
          wch <- which(!is.na(match(rnames[[i]], nm)))
          plist[[nm]]$random[[i]] <-
            if (length(wch) == 1) {         # only one formula for nm
              oSform <- asOneSidedFormula(ranForm[[i]][[nm]][[3]])
              model.matrix(oSform, model.frame(oSform, dataMix))
            } else {                        # multiple formulae
              lapply(ranForm[[i]][wch], function(el) {
                if (el[[3]] == "1") {
                  TRUE
                } else {
                  oSform <- asOneSidedFormula(el[[3]])
                  model.matrix(oSform, model.frame(oSform, dataMix))
                } })
            }
        }
      }
    }
  } else {
    namGrp <- ""
    grpsRev <- ranVec <- ran <- NULL
  }
  val <- vector("list", nlev)
  modForm <- getCovariateFormula(object)[[2]]
  omap <- object$map
  for(i in 1:nlev) {
    val[[i]] <- eval(modForm,
                     data.frame(dataMix,
                                getParsNlme(plist, omap$fmap, omap$rmapRel,
                                            omap$bmap, grpsRev, fix, ranVec, ran,
                                            level[i], N)))[naPat]
  }
  names(val) <- c("fixed", rev(namGrp))[level + 1]
  val <- as.data.frame(val)

  if (maxQ > 0) {
    val <- val[revOrder, , drop = FALSE]
    if (any(naGrps)) {
      val[naGrps] <- NA
    }
  }
  ## putting back in original order

  if (maxQ > 1) {                      # making groups unique
    for(i in 2:maxQ)
      oGrps[, i] <-
        as.factor(paste(as.character(oGrps[,i-1]),
                        as.character(oGrps[,i  ]), sep = "/"))
  }
  if (length(level) == 1) {
    val <- val[,1] ## ?? not in predict.lme()
    if (level > 0) { # otherwise 'oGrps' are typically undefined
      grps <- as.character(oGrps[, level])
      if (asList) {
        val <- split(val, ordered(grps, levels = unique(grps)))
      } else {
        names(val) <- grps
      }
    }
    lab <- "Predicted values"
    if (!is.null(aux <- attr(object, "units")$y)) {
      lab <- paste(lab, aux)
    }
    attr(val, "label") <- lab
    val
  } else {
    data.frame(oGrps, predict = val)
  }
}

## based on R's update.default
update.nlme <-
  function (object, model., ..., evaluate = TRUE)
{
  if (is.null(call <- object$call))
    stop("need an object with call component")
  extras <- match.call(expand.dots = FALSE)$...
  if (!missing(model.))
    call$model <- update.formula(formula(object), model.)
  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.nlme <-
##  function(object, model, data, fixed, random, groups, start, correlation,
##           weights, subset, method, na.action, naPattern, control,
##           verbose, ...)
## {
##  thisCall <- as.list(match.call())[-(1:2)]
##  if (!is.null(thisCall$start) && is.numeric(start)) {
##    thisCall$start <- list(fixed = start)
##  }
##  if (!is.null(nextCall <- object$origCall) &&
##      (is.null(thisCall$fixed) && !is.null(thisCall$random))) {
##    nextCall <- as.list(nextCall)[-1]
##  } else {
##    nextCall <- as.list(object$call)[-1]
##    if (is.null(thisCall$fixed)) {        # no changes in fixef model
##      if (is.null(thisCall$start)) {
##        thisCall$start <- list(fixed = fixef(object))
##      } else {
##        if (is.null(thisCall$start$fixed)) {
##          thisCall$start$fixed <- fixef(object)
##        }
##      }
##    }
##    if (!is.null(thisCall$start$random)) {  # making start random NULL
##      thisCall$start$random <- NULL
##    }
##    if (is.null(thisCall$random) && is.null(thisCall$subset)) {
##      ## no changes in ranef model and no subsetting
##      thisCall$random <- object$modelStruct$reStruct
##    }
##  }
##  if (!is.null(thisCall$model)) {
##    thisCall$model <- update(formula(object), model)
##  }
##  if (is.na(match("correlation", names(thisCall))) &&
##      !is.null(thCor <- object$modelStruct$corStruct)) {
##    thisCall$correlation <- thCor
##  }
##  if (is.na(match("weights", names(thisCall))) &&
##      !is.null(thWgt <- object$modelStruct$varStruct)) {
##    thisCall$weights <- thWgt
##  }
##  nextCall[names(thisCall)] <- thisCall
##  do.call(nlme, nextCall)
## }

###*### nlmeStruct - a model structure for nlme fits

nlmeStruct <-
  ## constructor for nlmeStruct objects
  function(reStruct, corStruct = NULL, varStruct = NULL)#, resp = NULL,
    ## model = NULL, local = NULL, N = NULL, naPat = NULL)
{

  val <- list(reStruct = reStruct, corStruct = corStruct,
              varStruct = varStruct)
  structure(val[!vapply(val, is.null, NA)], # removing NULL components
            settings = attr(val$reStruct, "settings"),
            ## attr(val, "resp") <- resp
            ## attr(val, "model") <- model
            ## attr(val, "local") <- local
            ## attr(val, "N") <- N
            ## attr(val, "naPat") <- naPat
            class = c("nlmeStruct", "lmeStruct", "modelStruct"))
}

##*## nlmeStruct methods for standard generics

fitted.nlmeStruct <-
  function(object, level = Q,  conLin = attr(object, "conLin"), ...)
{
  Q <- attr(object, "conLin")$dims[["Q"]]
  attr(object, "resp") - resid(object, level, conLin)
}


residuals.nlmeStruct <-
  function(object, level = Q, conLin = attr(object, "conLin"), ...)
{
  Q <- conLin$dims[["Q"]]
  stopifnot(length(level) >= 1)
  loc <- attr(object, "local")
  oLev <- get("level", envir = loc)
  on.exit(assign("level", oLev, envir = loc))
  dn <- c("fixed", rev(names(object$reStruct)))[level + 1]
  val <- array(0, c(attr(object, "NReal"), length(level)),
               list(dimnames(conLin$Xy)[[1]], dn))
  for(i in seq_along(level)) {
    assign("level", level[i], envir = loc, immediate = TRUE)
    val[, i] <- c(eval(attr(object, "model")[[2]], envir = loc))
  }
  val
}

nlmeControl <-
  ## Set control values for iterations within nlme
  function(maxIter = 50, pnlsMaxIter = 7, msMaxIter = 50,
	   minScale = 0.001, tolerance = 1e-5, niterEM = 25,
           pnlsTol = 0.001, msTol = 0.000001,
           returnObject = FALSE, msVerbose = FALSE, msWarnNoConv = TRUE,
           gradHess = TRUE, apVar = TRUE,
           .relStep = .Machine$double.eps^(1/3), minAbsParApVar = 0.05,
	   opt = c("nlminb", "nlm"), natural = TRUE, sigma = NULL, ...)
{
  if(is.null(sigma))
    sigma <- 0
  else if(!is.finite(sigma) || length(sigma) != 1 || sigma < 0)
    stop("Within-group std. dev. must be a positive numeric value")
  list(maxIter = maxIter, pnlsMaxIter = pnlsMaxIter, msMaxIter = msMaxIter,
       minScale = minScale, tolerance = tolerance, niterEM = niterEM,
       pnlsTol = pnlsTol, msTol = msTol,
       returnObject = returnObject, msVerbose = msVerbose,
       msWarnNoConv = msWarnNoConv,
       gradHess = gradHess, apVar = apVar, .relStep = .relStep,
       minAbsParApVar = minAbsParApVar,
       opt = match.arg(opt), natural = natural, sigma=sigma, ...)
}

nonlinModel <- function(modelExpression, env,
                        paramNames = get(".parameters", envir = env)) {
  modelExpression <- modelExpression[[2]]
  thisEnv <- environment()
  offset <- 0
  ind <- vector("list", length(paramNames))
  names(ind) <- paramNames
  for( i in paramNames ) {
    v.i <- get(i, envir = env)
    ind[[ i ]] <- offset + seq_along(v.i)
    offset <- offset + length(v.i)
  }
  modelValue <- eval(modelExpression, env)
  rm(i, offset, paramNames)
  function (newPars) {
    if(!missing(newPars)) {
      for(i in names(ind))
	assign( i, unname(newPars[ ind[[i]] ]), envir = env)
      assign("modelValue", eval(modelExpression, env), envir = thisEnv)
    }
    modelValue
  }
}


## Local Variables:
## ess-indent-offset: 2
## End: