.packageName <- "lme4"
# Ensure that the stats4 and lattice packages are available
.onLoad <- function(lib, pkg)
{
    require("methods")
#    require("stats4")
    require("lattice")
}

setClass("pdMat",      # parameterized positive-definite matrices
         representation(form="formula",    # a model-matrix formula
                        Names="character", # column (and row) names
                        param="numeric",   # parameter vector
                        Ncol="integer",    # number of columns
                        factor="matrix",   # factor of the pos-def matrix
                        logDet="numeric"   # logarithm of the absolute value
                        ## of the determinant of the factor (i.e. half
                        ## the logarithm of the determinant of the matrix)
                        ),
         prototype(form=formula(NULL),
                   Names=character(0),
                   param=numeric(0),
                   Ncol=as.integer(0),
                   factor=matrix(numeric(0),0,0),
                   logDet=numeric(0))
         )

#setClass("pdSymm", contains="pdMat")    # general symmetric pd matrices

#setClass("pdScalar", contains="pdMat") # special case of positive scalars
setClass("pdLogChol", contains="pdMat") # default parameterization
setClass("pdNatural", contains="pdMat") # log sd and logistic of correlation
#setClass("pdMatrixLog", contains="pdSymm") # matrix logarithm parameterization

setClass("pdDiag", contains="pdMat")    # diagonal pd matrices

setClass("pdIdent", contains="pdMat")   # positive multiple of the identity

setClass("pdCompSymm", contains="pdMat") # compound symmetric pd matrices

setClass("pdBlocked",                   # block-diagonal pd matrices
         representation("pdMat", components = "list"))

                       # positive-definite symmetric matrices as matrices
setClass("pdmatrix", contains="matrix")

                       # factors of positive-definite symmetric matrices
setClass("pdfactor", representation("matrix", logDet = "numeric"))

                       # correlation matrices and standard deviations
setClass("corrmatrix", representation("matrix", stdDev = "numeric"))

                       # a single level in the random effects structure
setClass("lmeLevel",
         representation(precision="pdMat", # the relative precision matrix
                        groups="factor",   # grouping factor for the level
                        columns="integer", # columns in model-matrix for level
                        parsInd="integer", # indices into the parameter vector
                        originalRows="list",
                        decomposedRows="list",
                        storedRows="list",
                        nrow="integer",
                        updateFactor="matrix", # used for EM update
                                               # and gradient calculation
                        hessianArray="array",
                        nlev="integer"))

                       # basic LME object representation
setClass("reStruct",
         representation(random="list", fixed="formula",
                        offset="numeric",
                        dirtyDecomposed="logical", useWeighted="logical",
                        dirtyStored="logical", dirtyBbetas="logical",
                        logLik="numeric",
                        analyticHessian="logical",
                        REML="logical",
                        reverseOrder="integer",
                        origOrder="integer",
                        original="matrix",
                        weighted="matrix", stored="matrix",
                        decomposed="matrix", bbetas="numeric",
                        dontCopy = "logical", assign.X = "ANY"),
         prototype=list(fixed = formula(NULL), dirtyBbetas = TRUE,
                        dirtyDecomposed=TRUE, REML=FALSE, dirtyStored=TRUE,
                        useWeighted=FALSE, logLik=as.numeric(NA),
                        dontCopy = FALSE, analyticHessian=FALSE))

setClass("lmeLevelList", contains="list")

setClass("summary.pdMat", representation(cor = "corrmatrix",
                                         structName = "character",
                                         noCorrelation = "logical",
                                         formula = "formula"),
         prototype=list(structName="", formula=formula(NULL)))

setClass("summary.reStruct",
         representation(fixed="formula",
                        coefficients="matrix",
                        scale="numeric",
                        denomDF="integer",
                        REML="logical",
                        ngrps="integer",
                        nobs="integer",
                        corFixed="corrmatrix",
                        reSumry="list",
                        useScale="logical",
                        showCorrelation="logical"
                        ))

setClass("summary.lme",
         representation(call = "call",
                        logLik = "logLik",
                        AIC = "numeric",
                        BIC = "numeric",
                        re ="summary.reStruct",
                        residuals = "numeric")
         )

setClass("summary.glmm", representation(method="character",
                                        family="character",
                                        link="character"),
         contains="summary.lme")

## Current versions of the methods package define

setClass("lme",
         representation(reStruct = "reStruct",
                        frame = "data.frame",
                        na.action = "ANY",
                        fitted = "numeric",
                        call = "call"),
         prototype(frame=data.frame(), fitted = numeric(0)))

## This is needed for the family slot of glmmStruct
setOldClass("family")

## Structure for fitting glmm classes
setClass("glmm",
         representation(family="family", # The glm family
                        origy="numeric",
                        n="numeric",
                        prior.weights="numeric",
                        init.weights="numeric",
                        init.y="numeric",
                        method="character"),
         contains = "lme")

setClass("lmList",
         representation(call = "call",
                        pool = "logical"),
         contains = "list")

setClass("VarCorr",
         representation(scale="numeric",
                        reSumry="list",
                        useScale="logical"),
         prototype = list(scale = 1.0, useScale = TRUE))










## deepayan experimenting with possible groupedData structures:



setClass("groupedData",
         representation(data = "data.frame",
                        formula = "formula",
                        outer = "formula",
                        inner = "formula",
                        labels = "list",
                        units = "list"))


## these are now obsolete: (nlme data are now in lme4 as S4 groupedData objects)

# ## Temporarily added so that groupedData objects are also data.frame objects
# setOldClass(c("nfnGroupedData", "nfGroupedData", "groupedData",
#               "data.frame"))

# setOldClass(c("nffGroupedData", "nfGroupedData", "groupedData",
#               "data.frame"))

# setOldClass(c("nmGroupedData", "groupedData", "data.frame"))

## The generics Names, Names<-, pdFactor, pdMatrix, and corMatrix
## will be deprecated in nlme_4.0
if (!isGeneric("Names")) {
    setGeneric("Names",
               function(object, ...)
           {
               .Deprecated("names")
               standardGeneric("Names")
           })
}

if (!isGeneric("Names<-")) {
    setGeneric("Names<-",
               function(object, value)
           {
               .Deprecated("names<-")
               standardGeneric("Names<-")
           })
}

if (!isGeneric("pdFactor")) {
    setGeneric("pdFactor",
               function(object) standardGeneric("pdFactor"))
}

if (!isGeneric("pdMatrix")) {
    setGeneric("pdMatrix",
               function(object) standardGeneric("pdMatrix"))
}

if (!isGeneric("corFactor")) {
    ## extractor for transpose inverse square root factor of corr matrix
    setGeneric("corFactor",
               function(object, ...) standardGeneric("corFactor"))
}

if (!isGeneric("corMatrix")) {
    ## extractor for correlation matrix or the transpose inverse
    ## square root matrix
    setGeneric("corMatrix",
               function(object, ...) standardGeneric("corMatrix"))
}

if (!isGeneric("isInitialized")) {
    setGeneric("isInitialized",
               function(object)
               #standardGeneric("isInitialized"), valueClass = "logical")
               standardGeneric("isInitialized"))
}

if (!isGeneric("logDet")) {
    setGeneric("logDet",
               function(object, covariate = getCovariate(object), ...)
               standardGeneric("logDet"))#,
               #valueClass = "numeric")
}

if (!isGeneric("matrix<-")) {
    setGeneric("matrix<-",
               function(object, value)
               standardGeneric("matrix<-"))
}

if (!isGeneric("coef<-")) {
    setGeneric("coef<-",
               function(object, value)
               standardGeneric("coef<-"))
}

if (!isGeneric("weighted<-")) {
    setGeneric("weighted<-", function(x, value) standardGeneric("weighted<-"))
}

if (!isGeneric("model.matrix<-")) {
    setGeneric("model.matrix<-", function(x, value)
               standardGeneric("model.matrix<-"))
}

if (!isGeneric("getGroups")) {
    ## Return the groups associated with object according to form.
    setGeneric("getGroups",
               function(object, form, level, data, sep)
               standardGeneric("getGroups"))
}

if (!isGeneric("getGroupsFormula")) {
    ## Return the formula(s) for the groups associated with object.
    ## The result is a one-sided formula unless asList is TRUE in which case
    ## it is a list of formulas, one for each level.
    setGeneric("getGroupsFormula",
               function(object, asList = FALSE, sep = "/")
               standardGeneric("getGroupsFormula"))
}

if (!isGeneric("getCovariate")) {
    ## Return the primary covariate associated with object
    setGeneric("getCovariate",
               function(object, form = formula(object), data = list())
               standardGeneric("getCovariate"))
}

if (!isGeneric("getResponse")) {
    ## Return the primary covariate associated with object
    setGeneric("getResponse",
               function(object, form = formula(object))
               standardGeneric("getResponse"))
}

if (!isGeneric("LMEgradient")) {
    setGeneric("LMEgradient",
               function(x, A, nlev) standardGeneric("LMEgradient"))
}

if (!isGeneric("LMEhessian")) {
    setGeneric("LMEhessian",
               function(x, A, H, nlev)
               standardGeneric("LMEhessian"))
}

setGeneric("lme",
           function(formula, data, random, correlation, weights, subset,
                    method, na.action, control, model, x)
           standardGeneric("lme"))

if (!isGeneric("EMupdate<-")) {
    setGeneric("EMupdate<-",
               function(x, nlev, value) standardGeneric("EMupdate<-"))
}

if (!isGeneric("reStruct")) {
    setGeneric("reStruct",
               function(fixed, random, data, weights, REML,
                        nextraCols=0, analyticHessian=FALSE)
               standardGeneric("reStruct"))
}

if (!isGeneric("glmmStruct")) {
    setGeneric("glmmStruct",
               function(formula, random, family, data, nextraCols=0,
                        method="PQL", ...)
               standardGeneric("glmmStruct"))
}

if (!isGeneric("pdgradient")) {
    setGeneric("pdgradient", function(x) standardGeneric("pdgradient"))
}

if (!isGeneric("EMsteps<-")) {
    setGeneric("EMsteps<-", function(x, value) standardGeneric("EMsteps<-"))
}

if (!isGeneric("LMEoptimize<-")) {
    setGeneric("LMEoptimize<-", function(x, value)
               standardGeneric("LMEoptimize<-"))
}

if (!isGeneric("response<-")) {
    setGeneric("response<-", function(x, value)
               standardGeneric("response<-"))
}

if (!isGeneric("fixef")) {
    setGeneric("fixef", function(object, ...)
               standardGeneric("fixef"))
}

if (!isGeneric("fixef<-")) {
    setGeneric("fixef<-",
               function(object, value) standardGeneric("fixef<-"))
}

## fixed.effects was an alternative name
fixed.effects = function(object, ...) {
    .Deprecated("fixef")
    mCall = match.call()
    mCall[[1]] = as.name("fixef")
    eval(mCall, parent.frame())
}

if (!isGeneric("ranef")) {
    setGeneric("ranef", function(object, ...)
               standardGeneric("ranef"))
}

## random.effects was an alternative name for ranef
random.effects = function(object, ...) {
    .Deprecated("ranef")
    mCall = match.call()
    mCall[[1]] = as.name("ranef")
    eval(mCall, parent.frame())
}

if (!isGeneric("BIC")) {
    setGeneric("BIC", function(object, ...)
               standardGeneric("BIC"))
}

setMethod("BIC", "logLik",
          function(object, ...)
          -2 * (c(object) - attr(object, "df") * log(attr(object, "nobs"))/2)
          )

if (!isGeneric("getFixDF")) {
    ## Return the formula(s) for the groups associated with object.
    ## The result is a one-sided formula unless asList is TRUE in which case
    ## it is a list of formulas, one for each level.
    setGeneric("getFixDF", function(object) standardGeneric("getFixDF"))
}

## FIXME: Can this be replaced by confint?
if (!isGeneric("intervals")) {
    setGeneric("intervals",
               function(object, level = 0.95, ...)
               standardGeneric("intervals"))
}

if (!isGeneric("lmList")) {
    setGeneric("lmList",
               function(formula, data, level, subset, na.action, pool)
               standardGeneric("lmList"))
}

if (!isGeneric("GLMM")) {
    setGeneric("GLMM",
               function(formula, family, data, random, control, niter,
                        method, verbose, ...)
               standardGeneric("GLMM"))
}

if (!isGeneric("pooledSD")) {
    setGeneric("pooledSD", function(object) standardGeneric("pooledSD"))
}

if (!isGeneric("VarCorr")) {
    setGeneric("VarCorr", function(x)
               standardGeneric("VarCorr"))
}
setMethod("VarCorr", signature(x="reStruct"),
          function(x)
      {
          nobs = dim(x@original)[1]
          nfixed = length(x@random[['*fixed*']]@columns)
          denomDF <- nobs - ifelse(x@REML, nfixed, 0)
          sigma <-
              abs(x@bbetas[x@random[['*response*']]@storedRows[[1]]])/denomDF
          new("VarCorr",
              scale=sigma,
              reSumry=lapply(rev(x@random)[-c(1, 2)],
                            function(x) summary(solve(x@precision))))
      })

setMethod("VarCorr", signature(x="lme"),
          function(x)
      {
          x <- x@reStruct
          callGeneric()
      })

setMethod("VarCorr", signature(x="glmm"),
          function(x)
      {
          useScale <- !(x@family$family %in% c("binomial", "poisson"))
          x <- x@reStruct
          x <- callGeneric()
          x@useScale <- useScale
          x
      })

setMethod("show", signature(object="VarCorr"),
          function(object)
      {
          digits = max(3, getOption("digits") - 2)
          useScale = length(object@useScale) > 0 && object@useScale[1]
          sc = ifelse(useScale, object@scale,  1.)
          reStdDev <- lapply(object@reSumry, function(x, sc) sc*x@cor@stdDev,
                             sc = sc)
          reLens = unlist(lapply(reStdDev, length))
          reMat <- array('', c(sum(reLens), 4),
                         list(rep('', sum(reLens)),
                              c("Groups", "Name", "Variance", "Std.Dev.")))
          reMat[1+cumsum(reLens)-reLens, 1] = names(reLens)
          reMat[,2] = unlist(lapply(reStdDev, names))
          reMat[,3] = format(unlist(reStdDev)^2, digits = digits)
          reMat[,4] = format(unlist(reStdDev), digits = digits)
          if (any(reLens > 1) &&
              !all(sapply(object@reSumry,
                          function(x) x@noCorrelation))) {
              maxlen = max(reLens)
              corr =
                  do.call("rbind",
                          lapply(object@reSumry,
                                 function(x, maxlen) {
                                     if (x@noCorrelation) {
                                         matrix("", dim(x@cor)[1], maxlen)
                                     } else {
                                         cc = format(round(x@cor, 3),
                                                     nsmall = 3)
                                         cc[!lower.tri(cc)] = ""
                                         nr = dim(cc)[1]
                                         cbind(cc, matrix("",
                                                          nr, maxlen-nr))
                                     }
                                 }, maxlen))
              colnames(corr) = c("Corr", rep("", maxlen - 1))
              reMat = cbind(reMat, corr)
          }
          if (useScale) {
              reMat <- rbind(reMat, c("Residual", "",
                                      format(sc^2, digits = digits),
                                      format(sc, digits = digits),
                                      rep('', ncol(reMat) - 4)))
          }
          print(reMat, quote = FALSE)
      })
## Methods for dealing with formulas

splitFormula <-
    function(form, sep = "/")
{
    ## split, on the sep call, the rhs of a formula into a list of subformulas
    if (inherits(form, "formula") ||
        mode(form) == "call" && form[[1]] == as.name("~"))
        return(splitFormula(form[[length(form)]], sep = sep))
    if (mode(form) == "call" && form[[1]] == as.name(sep))
        return(do.call("c", lapply(as.list(form[-1]), splitFormula, sep = sep)))
    if (mode(form) == "(") return(splitFormula(form[[2]], sep = sep))
    if (length(form) < 1) return(NULL)
    list(asOneSidedFormula(form))
}

subFormula <- function(form, pos = 2)
{
    ## extract component pos of form as a formula preserving the environment
    comp = form[[pos]]
    val = eval(substitute(~comp))
    environment(val) = environment(form)
    val
}

getCovariateFormula <- function(object)
{
    ## Return the primary covariate formula as a one sided formula
    form <- formula(object)
    form <- form[[length(form)]]
    if (length(form) == 3 && form[[1]] == as.name("|")){ # conditional expression
        form <- form[[2]]
    }
    eval(substitute(~form))
}

getResponseFormula <- function(object)
{
    ## Return the response formula as a one sided formula
    form <- formula(object)
    if (!(inherits(form, "formula") && (length(form) == 3)))
        stop("object must yield a two-sided formula")
    subFormula(form, 2)
}

setMethod("getGroupsFormula", signature(object = "ANY"),
          function(object, asList = FALSE, sep = "/")
      {
          form = formula(object)
          if (!inherits(form, "formula")) stop("object must yield a formula")
          rhs = form[[length(form)]]
          if (length(rhs) < 2 || rhs[[1]] != as.name("|")) return(NULL)
          if (asList) {
              val = splitFormula(asOneSidedFormula(rhs[[3]]), sep = sep)
              names(val) = unlist(lapply(val, function(el) deparse(el[[2]])))
              return(val)
          }
          asOneSidedFormula(rhs[[3]])
      })


setMethod("getGroups", signature(object="data.frame", form="formula"),
          function(object, form, level, data, sep)
              eval(getGroupsFormula(form)[[2]], object))
setMethod("GLMM", signature(data = "missing"),
          function(formula, family, data, random, control, niter,
                   method, verbose, ...)
      {
          nCall = mCall = match.call()
          nCall$data = list()
          .Call("nlme_replaceSlot", eval(nCall, parent.frame()), "call",
                mCall, PACKAGE = "lme4")
      })

setMethod("GLMM", signature(formula = "missing", data = "groupedData"),
          function(formula, family, data, random, control, niter,
                   method, verbose, ...)
      {
          nCall = mCall = match.call()
          resp = getResponseFormula(data)[[2]]
          cov = getCovariateFormula(data)[[2]]
          nCall$formula = eval(substitute(resp ~ cov))
          .Call("nlme_replaceSlot", eval(nCall, parent.frame()), "call",
                mCall, PACKAGE = "lme4")
      })

setMethod("GLMM", signature(formula = "formula", data = "groupedData",
                           random = "missing"),
          function(formula, family, data, random, control, niter,
                   method, verbose, ...)
      {
          nCall = mCall = match.call()
          cov = formula[[3]]
          grps = getGroupsFormula(data)[[2]]
          nCall$random = eval(substitute(~ cov | grps))
          .Call("nlme_replaceSlot", eval(nCall, parent.frame()), "call",
                mCall, PACKAGE = "lme4")
      })


setMethod("GLMM", signature(formula = "formula", random = "formula"),
          function(formula, family, data, random, control, niter,
                   method, verbose, ...)
      {
          nCall = mCall = match.call()
          nCall$random = lapply(getGroupsFormula(random, asList = TRUE),
                                function(x, form) form,
                                form = pdLogChol(getCovariateFormula(random)))
          .Call("nlme_replaceSlot", eval(nCall, parent.frame()), "call",
                mCall, PACKAGE = "lme4")
      })

setMethod("GLMM", signature(formula = "formula", random = "list"),
          function(formula, family, data, random, control, niter,
                   method, verbose, nEM.IRLS, model, x, ...)
      {
          if (missing(nEM.IRLS))
              nEM.IRLS <- 1
          if (missing(verbose)) verbose = FALSE
          m <- Call <- match.call()
          method <- if(missing(method)) "PQL" else
                    match.arg(method, c("PQL", "Laplace"))
          nm <- names(m)[-1]
          dontkeep <-
              is.element(nm, c("correlation", "control", "niter",
                               "verbose", "nEM.IRLS", "method"))
          for(i in nm[dontkeep]) m[[i]] <- NULL
          m[[1]] <- as.name("glmmStruct")
          m$nextraCols <- 1
          m$method <- method
          fit <- eval(m, parent.frame())

          off <- fit@reStruct@offset
          w <-  fit@prior.weights
          origy <- fit@origy
          fam <- fit@family

          ## We always do the PQL steps. For that, we extract the reStruct
          ## slot from fit
          fit <- as(fit, "reStruct")
          control = if (missing(control)) lmeControl() else
                    do.call("lmeControl", control)
          EMsteps(fit) <- control

          control$niterEM <- nEM.IRLS
          converged <- FALSE
          eta <- .Call("nlme_reStruct_fitted", fit, NULL, PACKAGE="lme4")

          for(i in seq(length=if(missing(niter)) 20 else niter)) {
              ##update zz and wz
              mu <- fam$linkinv(eta)
              mu.eta.val <- fam$mu.eta(eta)
              zz <- eta + (origy - mu)/mu.eta.val  - off
              wz <- w * mu.eta.val^2 / fam$variance(mu)

              response(fit) <- zz
              weighted(fit) <- sqrt(abs(wz))
              EMsteps(fit) <- control
              LMEoptimize(fit) <- control
              if(verbose) {
                  cat("iteration", i, "\n")
###             class(fit) <- "glmmStruct"
###             fit@logLik <- as.numeric(NA)
###             cat("Approximate logLik:",
###                 .Call("nlme_glmmLa2_logLikelihood",
###                     fit, NULL, PACKAGE="lme4"), "\n")
###            class(fit) <- "reStruct"
###            fit@logLik <- as.numeric(NA)
                  cat("Parameters:", coef(fit), "\n")
                  cat("Fixed Effects:", fixef(fit), "\n")
              }
              etaold <- eta
              eta <- .Call("nlme_reStruct_fitted", fit, NULL, PACKAGE="lme4")
              if(sum((eta-etaold)^2) < 1e-6*sum(eta^2)) {
                  converged <- TRUE
                  break
              }
          }
          if (control$msMaxIter > 0 && !converged)
              stop("IRLS iterations in GLMM failed to converge")
          ## We recreate the glmm object and set the reStruct slot
          .Call("nlme_replaceSlot", fit, "logLik",
                as.numeric(NA), PACKAGE = "lme4")
          .Call("nlme_replaceSlot", fit, "dontCopy", TRUE, PACKAGE = "lme4")
          fit <- .Call("nlme_replaceSlot", eval(m, parent.frame()),
                       "reStruct", fit, PACKAGE = "lme4")
          .Call("nlme_replaceSlot", fit, c("reStruct", "dontCopy"),
                TRUE, PACKAGE = "lme4")
          fit <- .Call("nlme_glmmLaplace_solveOnly", fit,
                       500, 1, PACKAGE="lme4")
          ## dirtyStored is FALSE but should be TRUE
#          .Call("nlme_replaceSlot", fit, c("reStruct", "dirtyStored"),
#                TRUE, PACKAGE = "lme4")
#          .Call("nlme_replaceSlot", fit, "reStruct",
#                .Call("nlme_commonDecompose", fit@reStruct,
#                      NULL, PACKAGE="lme4"), PACKAGE = "lme4")
#          .Call("nlme_replaceSlot", fit, c("reStruct", "dontCopy"),
#                FALSE, PACKAGE = "lme4")

          if (method != "PQL") {
              ## Do the 2nd order Laplace fit here
              LMEoptimize(fit) <- control
          }
          ## dirtyStored is FALSE but should be TRUE
          .Call("nlme_replaceSlot", fit, c("reStruct", "dirtyStored"),
                TRUE, PACKAGE = "lme4")
          .Call("nlme_replaceSlot", fit, "reStruct",
                .Call("nlme_commonDecompose", fit@reStruct,
                      NULL, PACKAGE="lme4"), PACKAGE = "lme4")
          .Call("nlme_replaceSlot", fit, c("reStruct", "dontCopy"),
                FALSE, PACKAGE = "lme4")
          ## zero some of the matrix slots
          if (!missing(x) && x == FALSE)
              .Call("nlme_replaceSlot", fit, c("reStruct", "original"),
                    matrix(0.0, nrow = 0, ncol = 0), PACKAGE = "lme4")
          .Call("nlme_replaceSlot", fit, c("reStruct", "decomposed"),
                matrix(0.0, nrow = 0, ncol = 0), PACKAGE = "lme4")
          .Call("nlme_replaceSlot", fit, c("reStruct", "weighted"),
                matrix(0.0, nrow = 0, ncol = 0), PACKAGE = "lme4")

          if (!missing(model) && model == FALSE)
              .Call("nlme_replaceSlot", fit, "frame",
                    data.frame(), PACKAGE = "lme4")
          .Call("nlme_replaceSlot", fit, "fitted",
                fam$linkinv(if (is.null(fit@na.action)) {
                    fitted(fit@reStruct)[fit@reStruct@reverseOrder]
                } else {
                    napredict(attr(data, "na.action"),
                              fitted(fit@reStruct)[fit@reStruct@reverseOrder])
                }), PACKAGE = "lme4")
          fit
      })

setMethod("summary", signature(object="glmm"),
          function(object, ...) {
              llik <- logLik(object)    # has an oldClass
              resd <- residuals(object, type="pearson")
              if (length(resd) > 5) {
                  resd <- quantile(resd)
                  names(resd) <- c("Min","Q1","Med","Q3","Max")
              }
              ans <- new("summary.glmm",
                         call = object@call,
                         logLik = llik,
                         AIC = AIC(llik),
                         BIC = BIC(llik),
                         re = summary(as(object, "reStruct")),
                         residuals = resd,
                         method = object@method,
                         family = object@family$family,
                         link = object@family$link)
              ans@re@useScale = !(ans@family %in% c("binomial", "poisson"))
              ans
          })

setMethod("show", "summary.glmm",
          function(object) {
              rdig <- 5
              cat("Generalized linear mixed-effects model fit by ")
              cat(switch(object@method, PQL="PQL\n",
                         Laplace="2nd order Laplace\n"))
              cat(" Family:", object@family, "with",
                  object@link, "link\n")
              cat(" Data:", deparse( object@call$data ), "\n")
              if (!is.null(object@call$subset)) {
                  cat("  Subset:",
                      deparse(asOneSidedFormula(object@call$subset)[[2]]),"\n")
              }
              print(data.frame(AIC = object@AIC, BIC = object@BIC,
                               logLik = c(object@logLik), row.names = ""))
              cat("\n")
              object@re@showCorrelation = TRUE
              show(object@re)
              ## Should this be part of the show method for summary.reStruct?
              cat("\nNumber of Observations:", object@re@nobs)
              cat("\nNumber of Groups: ")
              ngrps <- object@re@ngrps
              if ((length(ngrps)) == 1) {
                  cat(ngrps,"\n")
              } else {				# multiple nesting
                  cat("\n")
                  print(ngrps)
              }
              invisible(object)
          })

setMethod("show", "glmm",
          function(object)
      {
          sumry = summary(object)
          rdig <- 5
          cat("Generalized linear mixed-effects model\n")
          cat(" Family:", sumry@family, "with",
              sumry@link, "link\n")
          cat(" Data:", deparse( sumry@call$data ), "\n")
          if (!is.null(sumry@call$subset)) {
              cat("  Subset:",
                  deparse(asOneSidedFormula(sumry@call$subset)[[2]]),"\n")
          }
          cat(paste(" log-", ifelse(sumry@re@REML, "restricted-", ""),
                    "likelihood: ", sep = ''), sumry@logLik, "\n")
          sumry@re@showCorrelation = FALSE
          saveopt = options(show.signif.stars=FALSE)
          on.exit(saveopt)
          show(sumry@re)
          options(saveopt)
          on.exit()
          cat("\nNumber of Observations:", sumry@re@nobs, "\n")
          invisible(object)
      })

setMethod("getResponse", signature(object="glmm"),
          function(object, form)
      {
          object@origy
      })

setMethod("logLik", signature(object="glmm"),
          function(object)
      {
          value = .Call("nlme_glmmLaplace_logLikelihood", object,
                         NULL,          # do not pass new parameter value
                         50, 1, PACKAGE="lme4")
          p = length(object@reStruct@random[["*fixed*"]]@columns)
          # df calculated from sigma + fixed effects + random effects pars
          attr(value, "df") = 1 + p +
              sum(unlist(lapply(object@reStruct@random,
                                function(x)length(coef(x)))))
          attr(value, "nall") = length(object@fitted)
          attr(value, "nobs") = length(object@fitted) - p
          class(value) = "logLik"
          value
      })

setReplaceMethod("LMEoptimize", signature(x="glmm",
                                          value="list"),
                 function(x, value)
             {
                 value$analyticGradient <- FALSE
                 if (value$msMaxIter < 1)
                     return(x)
                 xval <- -.Call("nlme_glmmLaplace_logLikelihood", x,
                                NULL,
                                50, 1, PACKAGE="lme4")
                 xval =
                     if (xval > 0)
                         xval+1
                     else abs(min(xval/2, xval+1))
                 fixlen <- length(fixef(x))
                 if (value$optimizer == "optim") {
                     optimRes =
                         if (value$analyticGradient) {
                             optim(fn = function(params)
                               {
                                   fixef(x) <- params[seq(length=fixlen)]
                                   .Call("nlme_glmmLaplace_logLikelihood",
                                         x,
                                         params[-seq(length=fixlen)],
                                         50, 1, PACKAGE="lme4")
                               },
                                   gr = function(params)
                                   LMEgradient(.Call("nlme_commonDecompose",
                                                     x, params,
                                                     PACKAGE="lme4")),
                                   par = c(fixef(x), coef(x)),
                                   ##hessian = TRUE,
                                   method = "BFGS",
                                   control = list(trace = value$msVerbose,
                                   reltol = value$msTol,
                                   fnscale = -1,
#                                   fnscale = -xval,
#                                   parscale = 1/value$msScale(coef(x)),
                                   maxit = value$msMaxIter))
                         } else {
                             optim(fn = function(params)
                               {
                                   fixef(x) <- params[seq(length=fixlen)]
                                   .Call("nlme_glmmLaplace_logLikelihood",
                                         x,
                                         params[-seq(length=fixlen)],
                                         50, 1, PACKAGE="lme4")
                               },
                                   par = c(fixef(x), coef(x)),
                                   #hessian = TRUE,
                                   method = "BFGS",
                                   control = list(trace = value$msVerbose,
                                   reltol = value$msTol,
                                   fnscale = -1,
#                                   fnscale = -xval,
#                                   parscale = 1/value$msScale(coef(x)),
                                   maxit = value$msMaxIter))
                         }
                     if (optimRes$convergence != 0) {
                         warning("optim failed to converge")
                     }
                     fixef(x) <- optimRes$par[seq(length=fixlen)]
                     coef(x@reStruct) <- optimRes$par[-seq(length=fixlen)]
                     x@reStruct@logLik <- as.numeric(NA)
                     .Call("nlme_glmmLaplace_solveOnly", x,
                           500, 1, PACKAGE="lme4")
#                  } else if (value$optimizer == "ms") {
#                      pars <- coef(x)
#                      .Call("nlme_msOptimize", value$msMaxIter,
#                            value$msTol, rep(1.0, length(pars)),
#                            value$msVerbose, x, pars,
#                            value$analyticGradient,
#                            PACKAGE = "lme4")
                 } else {
#                     typsize <- 1/value$msScale(coef(x))
                     typsize <- rep(1.0, length(coef(x))+fixlen)
                     if (is.null(value$nlmStepMax))
                         value$nlmStepMax <-
                             max(100 * sqrt(sum((c(fixef(x),
                                                   coef(x))/typsize)^2)), 100)
                     nlmRes =
                         nlm(f = if (value$analyticGradient) {
                             function(params)
                             {
                                 fixef(x) <- params[seq(length=fixlen)]
                                 -.Call("nlme_glmmLaplace_logLikelihood",
                                        x,
                                        params[-seq(length=fixlen)],
                                        50, 1, PACKAGE="lme4")
#                                  x = .Call("nlme_commonDecompose",
#                                             x, params,
#                                             PACKAGE="lme4")
#                                  grad = -LMEgradient(x)
#                                  ans = -x@logLik
#                                  attr(ans, "gradient") = grad
#                                  ans
                             }
                         } else {
                             function(params)
                             {
                                 fixef(x) <- params[seq(length=fixlen)]
                                 -.Call("nlme_glmmLaplace_logLikelihood",
                                        x,
                                        params[-seq(length=fixlen)],
                                        50, 1, PACKAGE="lme4")
                             }
                         },
                             p = c(fixef(x), coef(x)),
                             #hessian = TRUE,
                             print.level = if (value$msVerbose) 2 else 0,
                             steptol = value$msTol,
                             gradtol = value$msTol,
                             stepmax = value$nlmStepMax,
                             typsize=typsize,
#                             fscale=xval,
                             iterlim = value$msMaxIter)
                     fixef(x) <- nlmRes$estimate[seq(length=fixlen)]
                     coef(x@reStruct) <-
                         nlmRes$estimate[-seq(length=fixlen)]
                     x@reStruct@logLik <- as.numeric(NA)
                     .Call("nlme_glmmLaplace_solveOnly", x,
                           500, 1, PACKAGE="lme4")
                 }
             })

setReplaceMethod("fixef", signature(object="glmm", value="numeric"),
          function(object, value) {
              fixef(object@reStruct) <- value
              object
          })

### Local variables:
### mode: R
### End:
setMethod("glmmStruct", signature(formula = "formula",
                                  random = "list"),
          function(formula, random, family, data, nextraCols, method, ...)
      {

          addToFamily <- function(fam)
          {
              if (method == "Laplace")
                  fam$mu2.eta2 <-
                      if (fam$family %in% c("binomial",
                                            "quasibinomial") &&
                          fam$link == "logit") {
                          function(eta, mu, mu.eta) mu.eta*(1-2*mu)
                      } else if (fam$family %in% c("poisson",
                                                   "quasipoisson") &&
                                 fam$link == "log") {
                          function(eta, mu, mu.eta) mu.eta
                      } else {
                          stop("can not handle link ",
                               fam$link , " for family ",
                               fam$family)
                      }
              fam
          }


          mcall1 <- mcall2 <- match.call()

          mcall1$random <- mcall1$nextraCols <- mcall1$method <- NULL
          mcall1[[1]] <- as.name("glm")
          glmFit <- eval(mcall1, parent.frame())
          rm(mcall1)

          nm <- names(mcall2)[-1]
          keep <-
              is.element(nm, c("data", "subset", "na.action",
                               "xlev", "offset", "weights"))
          for(i in nm[!keep]) mcall2[[i]] <- NULL
          allvars <- c(unlist(lapply(random, all.vars)), names(random))
          mcall2$formula <-
              as.formula(paste(paste(deparse(formula), collapse = ''),
                               paste(allvars, collapse="+"), sep = "+"))
          environment(mcall2$formula) <- environment(formula)
          mcall2$drop.unused.levels <- TRUE
          mcall2[[1]] <- as.name("model.frame")
          data <-  eval(mcall2, parent.frame())
          if (is.null(glmFit$offset))
              glmFit$offset <- 0.0
          data[, attr(attr(data, 'terms'), 'response')] <-
              glmFit$linear.predictor + glmFit$residuals - glmFit$offset

          ans <- new("glmm",
                     reStruct = reStruct(fixed=formula,
                                         random=random,
                                         data=data,
                                         weights=sqrt(abs(glmFit$weights)),
                                         REML=FALSE,
                                         nextraCols=nextraCols),
                     frame = data, call = match.call(),
                     na.action = attr(data, "na.action"),
                     family = addToFamily(glmFit$family),
                     method = method)
          rm(data)
          origOrder <- ans@reStruct@origOrder

          ans@origy <- as.numeric(glmFit$y[origOrder])
          ans@n <- local({
              y <- model.response(glmFit$model, "numeric")
              mt <- attr(glmFit$model, "terms")
              x <- if (!is.empty.model(mt))
                  model.matrix(mt, glmFit$model, glmFit$contrasts)
              x <- as.matrix(x)
              nobs <- NROW(y)
              weights <- glmFit$prior.weights
              ## calculates mustart and may change y and weights and set n (!)
              eval(ans@family$initialize)
              n
          })[origOrder]
          if (length(glmFit$offset) == 1)
              ans@reStruct@offset <- glmFit$offset
          else ans@reStruct@offset <- glmFit$offset[origOrder]
          if (length(glmFit$prior.weights) == 1)
              ans@prior.weights <- glmFit$prior.weights
          else ans@prior.weights <- glmFit$prior.weights[origOrder]
          if (length(glmFit$weights) == 1)
              ans@init.weights <- glmFit$weights
          else ans@init.weights <- glmFit$weights[origOrder]
          ans@init.y <- getResponse(ans@reStruct)
          ans@reStruct@logLik <- as.numeric(NA)
          ans
      })

glmmLa2RespWt <-
    function(fam, eta, origy, w, off)
{
    mu <- fam$linkinv(eta)
    mu.eta.val <- fam$mu.eta(eta)
    list(eta + (origy - mu)/mu.eta.val  - off,
         sqrt(abs(w * mu.eta.val)))
}

glmmLa2Wt2 <-
    function(fam, eta, w)
{
    abs(w)*fam$mu2.eta2(eta, fam$linkinv(eta), fam$mu.eta(eta))
}

glmmLa2LogLikComp <-
    function(x, eta)
{
    fam <- x@family
    origy <- x@origy
    w <- x@prior.weights

    mu <- fam$linkinv(eta)
    nobs <- length(origy)
    aic <- fam$aic(origy, x@n, mu, w, sum(fam$dev.resids(origy, mu, w)))
    ## allow for estimated dispersion
    if (fam$family %in% c("gaussian", "Gamma", "inverse.gaussian",
                          "quasibinomial", "quasipoisson"))
        1 - aic/2
    else -aic/2
###    m <- if (any(x@n > 1)) x@n else w
###    sum(dbinom(round(m*origy), round(m), mu, log=TRUE))
}
# Provisional method.
# The definition of the method should change when groupedData is defined
#  as an S4 class. (update: which has now happened)

# setMethod("formula",
#           "groupedData",
#           function(x, ...) attr(x, "formula"),
#           valueClass = "formula")



## experimental groupedData class -- deepayan


groupedData <-
    function(formula, data = sys.parent(1),  ## parent.frame() ? is this a data.frame ?
             order.groups = TRUE,
             FUN = function(x) max(x, na.rm = TRUE),
             outer = ~0, inner = ~0,
             labels = list(), units = list())
{
    new("groupedData",
        data = data,
        formula = formula,
        outer = outer,
        inner = inner,
        labels = labels,
        units = units)
}



setMethod("show", "groupedData",
          function(object) {
              cat("groupedData with Formula: ")
              show(formula(object))
              cat("\n")
              show(object@data)
          })
          


convertFromS3groupedData <- function(from)
{
    data <- from
    class(data) <- "data.frame"
    attr(data, "labels") <- NULL
    attr(data, "units") <- NULL
    attr(data, "inner") <- NULL
    attr(data, "outer") <- NULL
    attr(data, "formula") <- NULL
    attr(data, "order.groups") <- NULL

    ans <- groupedData(data = data,
                       formula = attr(from, "formula"),
                       outer = if (is.null(attr(from, "outer")))
                       formula(~0) else attr(from, "outer"),
                       inner = if (is.null(attr(from, "inner")))
                       formula(~0) else attr(from, "inner"),
                       labels = if (is.null(attr(from, "labels")))
                       list() else attr(from, "labels"),
                       units = if (is.null(attr(from, "units")))
                       list() else attr(from, "units"))
    ans
}

## Method for x[i, j]

setMethod("[",
          signature(x = "groupedData", i = "missing", j = "ANY",
                    drop = "logical"),

          function(x, i, j, drop) {
              if (drop) warning("drop=TRUE will be ignored")
              x@data <- x@data[, j, drop = FALSE]
              x
          })


setMethod("[",
          signature(x = "groupedData", i = "ANY", j = "missing",
                    drop = "logical"),

          function(x, i, j, drop) {
              if (drop) warning("drop=TRUE will be ignored")
              x@data <-
                  if (nargs() == 4) x@data[i, , drop = FALSE]
                  else x@data[i, drop = FALSE]
              x
          })


setMethod("[",
          signature(x = "groupedData", i = "missing", j = "ANY",
                    drop = "missing"),

          function(x, i, j, drop) {
              x@data <- x@data[, j, drop = FALSE]
              x
          })


setMethod("[",
          signature(x = "groupedData", i = "ANY", j = "missing",
                    drop = "missing"),

          function(x, i, j, drop) {
              x@data <-
                  if (nargs() == 3) x@data[i, ]
                  else x@data[i]
              x
          })


setMethod("[",
          signature(x = "groupedData", i = "ANY", j = "ANY",
                    drop = "missing"),

          function(x, i, j, drop) {
              x@data <- x@data[i, j, drop = FALSE]
              x
          })


setMethod("[",
          signature(x = "groupedData", i = "ANY", j = "ANY",
                    drop = "logical"),
          
          function(x, i, j, drop) {
              if (drop) warning("drop=TRUE will be ignored")
              x@data <- x@data[i, j, drop = FALSE]
              x
          })


## Method for x[[i, j]]


setMethod("[[",
          signature(x = "groupedData", i = "ANY", j = "missing"),
          function(x, i, j) {
              x@data[[i]]
          })

setMethod("[[",
          signature(x = "groupedData", i = "ANY", j = "ANY"),
          function(x, i, j) {
              x@data[[i, j]]
          })


## Method for x$name

setMethod("$",
          signature(x = "groupedData", name = "ANY"),
          function(x, name) {
              x@data[[name]]
          })



## Method for x$name <-


setMethod("$<-",
          signature(x = "groupedData", name = "ANY", value = "ANY"),
          function(x, name, value) {
              x@data[[name]] <- value
              x
          })





## Method for x[...] <-


setMethod("[<-",
          signature(x = "groupedData", i = "ANY", j = "ANY", value = "ANY"),
          function(x, i, j, value) {
              x@data[i, j] <- value
              x
          })

setMethod("[<-",
          signature(x = "groupedData", i = "ANY", j = "missing", value = "ANY"),
          function(x, i, j, value) {
              if (nargs() == 3) x@data[i] <- value
              else x@data[i,] <- value
              x
          })

setMethod("[<-",
          signature(x = "groupedData", i = "missing", j = "ANY", value = "ANY"),
          function(x, i, j, value) {
              x@data[, j] <- value
              x
          })
















setMethod("formula",
          "groupedData",
          function(x, ...) x@formula)


#setMethod("getData",
#          "groupedData",
#          function(x, ...) x@data)


setAs("groupedData", "data.frame",
      function(from) from@data)




collectorPlot <-
    function(formula, data, displayfunction,
             labels, xargs,
             ...)
{
    ## It's probably sufficient to have displayfunction a quoted fn name (like "xyplot")
    ## Think about it

    dotargs <- list(...)


    if (!is.function(displayfunction)) {
        ## needs more sophistication/checking ?
        displayfunction <- get(displayfunction)
    }

    ## adjust arguments - dotargs and xargs

    # add labels
    xargs$xlab <- labels$x
    xargs$ylab <- labels$y

    ## Combine args:
    dotargs <- c(dotargs, xargs[!(names(xargs) %in% names(dotargs))])

    do.call("displayfunction",
            c(list(formula = formula, data = data), dotargs))
}

setMethod("summary", signature(object = "groupedData"),
          function(object, ...) summary(as(object, "data.frame"), ...))

setMethod("plot",
          signature(x = "groupedData", y="missing"),
          function(x,
                   formula,
                   inner = NULL, outer = NULL,
                   ...) {


              respVar <- getResponseFormula(x)[[2]]
              isNumericResponse <- is.numeric(eval(respVar, x@data))

              covVar <- getCovariateFormula(x)[[2]]
              if (!(isTrivialCovariate <- covVar == 1))
                  isNumericCovariate <- is.numeric(eval(covVar, x@data))

              groupVar <- getGroupsFormula(x)[[2]]
              isMultipleGroups <- length(groupVar) > 2




              if (isMultipleGroups) return("multiple groups not implemented yet")


              if (!isNumericResponse) return("non-numeric response plots not implemented yet")


              if (isTrivialCovariate || !isNumericCovariate)
              {
                  xargs <- list()
                  displayfunction = "dotplot"
                  form = as.formula(paste(groupVar, respVar, sep = "~"))
                  groups = if (isTrivialCovariate) NULL else covVar

                 xlab <-
                      if ("y" %in% names(x@labels)) x@labels$y
                      else as.character(respVar)
                  if ("y" %in% names(x@units)) xlab <- paste(xlab, x@units$y)
                  ylab <- as.character(groupVar)
                  labels <- list(x = xlab, y = ylab)



                  if (!is.null(groups))
                      xargs$groups <- groups



                  ## groups needs to be combined with inner, perhaps
                  ##if (!isTrivialCovariate && is.null(inner)) inner = covVar
              }

              else if (isNumericCovariate)
              {
                  ## Custom panel function
                  panel.nfn = function(x, y, ...) {
                      panel.xyplot(x, y, ...)
                      y.avg <- tapply(y, x, mean) # lines through average y
                      y.avg <- y.avg[!is.na(y.avg)]
                      if (length(y.avg) > 0) {
                          xvals <- as.numeric(names(y.avg))
                          ord <- order(xvals)
                          panel.xyplot(xvals[ord], y.avg[ord], type = "l", ...)
                      }
                  }

                  xlab <-
                      if ("x" %in% names(x@labels)) x@labels$x
                      else as.character(covVar)
                  if ("x" %in% names(x@units)) xlab <- paste(xlab, x@units$x)
                  ylab <-
                      if ("y" %in% names(x@labels)) x@labels$y
                      else as.character(respVar)
                  if ("y" %in% names(x@units)) ylab <- paste(ylab, x@units$y)
                  labels <- list(x = xlab, y = ylab)


                  xargs <- list(aspect = "xy", grid = TRUE)

                  displayfunction = "xyplot"
                  form = as.formula( paste(   paste(respVar, covVar, sep = "~"), groupVar, sep = "|") )

                  groups = NULL
                  if (is.null(inner)) inner <- x@inner
                  if (inherits(inner, "formula") && !inner[[2]]==0) groups <- inner[[2]]
                  ## what's innerGroups ?

                  if (!is.null(groups))
                  {
                      xargs$groups <- groups
                      xargs$panel <- function(x, y, grid = TRUE, ...)
                      {
                          if (grid) panel.grid()
                          panel.superpose(x, y, ...)
                      }
                      xargs$panel.groups <- panel.nfn
                  }
                  else
                  {
                      xargs$panel <- function(x, y, grid = TRUE, ...)
                      {
                          if (grid) panel.grid()
                          panel.nfn(x, y, ...)
                      }
                  }


              }

              else return("didn't expect this to happen!")


              collectorPlot(formula = if (missing(formula)) form else formula,
                            data = x@data,
                            displayfunction = displayfunction,
                            labels = labels,
                            xargs = xargs,
                            ...)

          })





## Tools for plot methods



## was old nlme plot method (???)
plot.nmGroupedData <-
      function(x, collapseLevel = Q, displayLevel = collapseLevel,
                  outer = NULL, inner = NULL, preserve = NULL, FUN = mean,
                          subset = NULL, key = TRUE, grid = TRUE, ...)
{
      args <- list(outer = outer, inner = inner, key = key, grid = grid, ...)
        Q <- length(getGroupsFormula(x, asList = TRUE))
        if (is.null(preserve) && (collapseLevel < Q) && (!is.null(inner))) {
                if (is.logical(inner)) {
                          preserve <- attr(x, "inner")[[displayLevel]]
                      } else {
                                preserve <- inner
                            }
            }
        x <- collapse(x, collapseLevel, displayLevel, outer, inner,
                      preserve, FUN, subset)
        args[["innerGroups"]] <- attr(x, "innerGroups")
        args[["x"]] <- x
        do.call("plot", args)
  }








## The next 2 methods added for groupedData objects, but
## this may not be the best way to do this

setMethod("lmList", signature(formula = "formula", data = "groupedData"),
          function(formula, data, level, subset, na.action, pool)
      {
          data <- data@data
          standardGeneric("lmList")
      })

setMethod("lmList", signature(formula = "missing", data = "groupedData"),
          function(data, level, subset, na.action, pool)
      {
          formula <- data@formula ## formula(data) doesn't work !
          data <- data@data
          standardGeneric("lmList")
      })

          

setMethod("lmList", signature(formula = "formula", data = "data.frame"),
          function(formula, data, level, subset, na.action, pool)
      {
          mCall = frmCall = match.call()
          resp = getResponseFormula(formula)[[2]]
          cov = getCovariateFormula(formula)[[2]]
          lmForm = eval(substitute(resp ~ cov))
          gfm = getGroupsFormula(formula)
          if (is.null(gfm)) gfm = getGroupsFormula(data)
          if (is.null(gfm))
              stop("Unable to determine a grouping formula from either the formula or the data")
          val <- lapply(split(data, eval(gfm[[2]], data)),
                        function(dat, formula)
                    {
                        ans <- try({
                            data <- as.data.frame(dat)
                            lm(formula = formula, data = data)
                        })
                        if (inherits(ans, "try-error"))
                            NULL
                        else ans
                    }, formula = lmForm)
          if (missing(pool)) pool = TRUE
          new("lmList", val, call = mCall, pool = pool)
      })


setMethod("coef", signature(object = "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 = !unlist(lapply(coefs, is.null))
          if (sum(non.null) > 0) {
              template = coefs[non.null][[1]]
              if (is.numeric(template)) {
                  co <- matrix(template,
                               ncol = length(template),
                               nrow = length(coefs),
                               byrow = TRUE,
                               dimnames = list(names(object), names(template)))
                  for (i in names(object)) {
                      co[i,] = if (is.null(coefs[[i]])) { NA } else 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
      })

setMethod("show", signature(object = "lmList"),
          function(object)
      {
          mCall = object@call
          cat("Call:", deparse(mCall), "\n")
          cat("Coefficients:\n")
          invisible(print(coef(object)))
          if (object@pool) {
              cat("\n")
              poolSD = pooledSD(object)
              dfRes = attr(poolSD, "df")
              RSE = c(poolSD)
              cat("Degrees of freedom: ", length(unlist(lapply(object, fitted))),
                  " total; ", dfRes, " residual\n", sep = "")
              cat("Residual standard error:", format(RSE))
              cat("\n")
          }
      })

setMethod("pooledSD", signature(object = "lmList"),
          function(object)
      {
          sumsqr <- apply(sapply(object,
                                 function(el) {
                                     if (is.null(el)) {
                                         c(0,0)
                                     } else {
                                         res = resid(el)
                                         c(sum(res^2), length(res) - length(coef(el)))
                                  }
                              }), 1, sum)
          if (sumsqr[2] == 0) {
              stop("No degrees of freedom for estimating std. dev.")
          }
          val <- sqrt(sumsqr[1]/sumsqr[2])
          attr(val, "df") <- sumsqr[2]
          val
      })

setMethod("intervals", signature(object = "lmList", level = "ANY"),
          function(object, level = 0.95, ...)
          cat("intervals method for lmList not yet implemented\n"))

setMethod("plot", signature(x = "lmList"),
          function(x, y, ...)
          cat("plot method for lmList not yet implemented\n"))

setMethod("update", signature(object = "lmList"),
          function(object, formula., ..., evaluate = TRUE)
      {
          call <- object@call
          if (is.null(call))
              stop("need an object with call component")
          extras <- match.call(expand.dots = FALSE)$...
          if (!missing(formula.))
              call$formula <- update.formula(formula(object), formula.)
          if (length(extras) > 0) {
              existing <- !is.na(match(names(extras), names(call)))
              for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
              if (any(!existing)) {
                  call <- c(as.list(call), extras[!existing])
                  call <- as.call(call)
              }
          }
          if (evaluate)
              eval(call, parent.frame())
          else call
      })

setMethod("formula", signature(x = "lmList"),
          function(x, ...) x@call[["formula"]])
lmeControl <-
  ## Control parameters for lme
  function(maxIter = 50, msMaxIter = 50, tolerance =
           sqrt((.Machine$double.eps)), niterEM = 35,
           msTol = sqrt(.Machine$double.eps), msScale, msVerbose = FALSE,
           returnObject = FALSE, gradHess = TRUE, apVar = TRUE,
           .relStep = (.Machine$double.eps)^(1/3), minAbsParApVar = 0.05,
           nlmStepMax = NULL,
           natural = TRUE, optimizer="nlm", EMverbose=FALSE,
           analyticGradient = TRUE,
           analyticHessian=FALSE)
{
    if (missing(msScale)) msScale = function(start) {
        scale <- abs(start)
        nonzero <- scale > 0
        if (any(nonzero)) {
            scale[nonzero] <- 1/scale[nonzero]
            scale[!nonzero] <- median(scale[nonzero])
        }
        else {
            scale <- rep(1, length(scale))
        }
        scale
    }
    list(maxIter = maxIter, msMaxIter = msMaxIter, tolerance = tolerance,
         niterEM = niterEM, msTol = msTol, msScale = msScale,
         msVerbose = msVerbose, returnObject = returnObject,
         gradHess = gradHess , apVar = apVar, .relStep = .relStep,
         nlmStepMax = nlmStepMax,
         minAbsParApVar = minAbsParApVar, natural = natural,
         optimizer=optimizer, EMverbose=EMverbose,
         analyticHessian=analyticHessian,
         analyticGradient=analyticGradient)
}

setMethod("lme", signature(data = "missing"),
          function(formula, data, random, correlation, weights, subset,
                   method, na.action, control, model, x)
      {
          nCall = mCall = match.call()
          nCall$data = list()
          .Call("nlme_replaceSlot", eval(nCall, parent.frame()), "call",
                mCall, PACKAGE = "lme4")
      })

setMethod("lme", signature(formula = "missing", data = "groupedData"),
          function(formula, data, random, correlation, weights, subset,
                   method, na.action, control, model, x)
      {
          nCall = mCall = match.call()
          resp = getResponseFormula(data)[[2]]
          cov = getCovariateFormula(data)[[2]]
          nCall$formula = eval(substitute(resp ~ cov))
          .Call("nlme_replaceSlot", eval(nCall, parent.frame()), "call",
                mCall, PACKAGE = "lme4")
      })

setMethod("lme", signature(formula = "formula", data = "groupedData",
                           random = "missing"),
          function(formula, data, random, correlation, weights, subset,
                   method, na.action, control, model, x)
      {
          nCall = mCall = match.call()
          cov = formula[[3]]
          grps = getGroupsFormula(data)[[2]]
          nCall$random = eval(substitute(~ cov | grps))
          .Call("nlme_replaceSlot", eval(nCall, parent.frame()), "call",
                mCall, PACKAGE = "lme4")
      })


setMethod("lme", signature(formula = "formula", random = "formula"),
          function(formula, data, random, correlation, weights, subset,
                   method, na.action, control, model, x)
      {
          nCall = mCall = match.call()
          nCall$random = lapply(getGroupsFormula(random, asList = TRUE),
                                function(x, form) form,
                                form = pdLogChol(getCovariateFormula(random)))

          nCall$data <- as(data, "data.frame")

          .Call("nlme_replaceSlot", eval(nCall, parent.frame()), "call",
                mCall, PACKAGE = "lme4")
      })

setMethod("lme", signature(formula = "formula", random = "list"),
          function(formula, data, random, correlation, weights, subset,
                   method, na.action, control, model, x)
      {
          if (missing(model))
              model = TRUE
          if (missing(x))
              x = TRUE
          random = lapply(random, function(x)
                          if(inherits(x, "formula")) pdLogChol(x) else x)
          method = if (missing(method)) "REML" else
                   match.arg(method, c("REML", "ML"))
          controlvals <- if (missing(control)) lmeControl() else
                            do.call("lmeControl", control)
          mCall <- match.call(expand.dots = FALSE)
          mCall[[1]] <- as.name("model.frame")
          names(mCall)[2] <- "formula"
          mCall$random <- mCall$correlation <- mCall$method <-
              mCall$control <- NULL
          form <- formula
          form[[3]] <- (~a+b)[[2]]
          form[[3]][[2]] <- formula[[3]]
          form[[3]][[3]] <-
              as.formula((parse(text=paste("~",
                                paste(names(random),
                                      collapse = "+")))[[1]]))[[2]]
          for (pdm in random) {
              tmp <- form
              tmp[[3]] <- (~a+b)[[2]]
              tmp[[3]][[2]] <- form[[3]]
              tmp[[3]][[3]] <- formula(pdm)[[2]]
              form <- tmp
          }
          environment(form) <- environment(formula)
          mCall$formula <- form
          mCall$drop.unused.levels <- TRUE
          data <- eval(mCall, parent.frame())
          re <- reStruct(fixed = formula, random = random,
                         data = data,
                         REML = method != "ML",
                         analyticHessian=controlvals$analyticHessian)
          .Call("nlme_replaceSlot", re, "dontCopy", TRUE, PACKAGE = "lme4")
          .Call("nlme_replaceSlot", re, "analyticHessian",
                FALSE, PACKAGE = "lme4")
          EMsteps(re) <- controlvals
          .Call("nlme_replaceSlot", re, "analyticHessian",
                controlvals$analyticHessian, PACKAGE = "lme4")
          LMEoptimize(re) <- controlvals
          .Call("nlme_replaceSlot", re, "dontCopy", FALSE, PACKAGE = "lme4")
          ## zero some of the matrix slots
          if (x == FALSE)
              .Call("nlme_replaceSlot", re, "original",
                    matrix(0.0, nrow = 0, ncol = 0), PACKAGE = "lme4")
          .Call("nlme_replaceSlot", re, "decomposed",
                matrix(0.0, nrow = 0, ncol = 0), PACKAGE = "lme4")
          .Call("nlme_replaceSlot", re, "weighted",
                matrix(0.0, nrow = 0, ncol = 0), PACKAGE = "lme4")
          if (model == FALSE)
              data = data.frame()
          new("lme", reStruct = re, call = match.call(),
              fitted = if (is.null(attr(data, "na.action"))) {
                  fitted(re)[re@reverseOrder]
              } else {
                  napredict(attr(data, "na.action"), fitted(re)[re@reverseOrder])
              },
              frame = data, na.action = attr(data, "na.action"))
      })

setAs("lme", "reStruct",
      function(from) from@reStruct,
      function(from, value) {
          if (nrow(from@frame) != nrow(value@original))
              stop("Dimension mismatch between model.frame and original matrix")
          from@reStruct <- value
          from
      })

setMethod("fitted", signature=c(object="lme"),
          function(object, ...)
      {
          object@fitted
      })


setMethod("residuals", signature=c(object="lme"),
          function(object, ...)
      {
          re <- as(object, "reStruct")
          if (is.null(object@na.action)) {
              (getResponse(object@reStruct) - fitted(object@reStruct))[re@reverseOrder]
          } else {
              napredict(object@na.action,
                        (getResponse(object@reStruct) -
                         fitted(object@reStruct))[object@reStruct@reverseOrder])
          }
      })


setMethod("logLik", signature(object="lme"),
          function(object) logLik(object@reStruct))

setMethod("summary", signature(object="lme"),
          function(object, ...) {
              llik <- logLik(object)    # has an oldClass
              resd <- residuals(object, type="pearson")
              if (length(resd) > 5) {
                  resd <- quantile(resd)
                  names(resd) <- c("Min","Q1","Med","Q3","Max")
              }
              new("summary.lme",
                  call = object@call,
                  logLik = llik,
                  AIC = AIC(llik),
                  BIC = BIC(llik),
                  re = summary(as(object, "reStruct")),
                  residuals = resd)
          })

setMethod("show", "summary.lme",
          function(object) {
              rdig <- 5
              cat("Linear mixed-effects model fit by ")
              cat(ifelse(object@re@REML, "REML\n", "maximum likelihood\n") )
              cat(" Data:", deparse( object@call$data ), "\n")
              if (!is.null(object@call$subset)) {
                  cat("  Subset:",
                      deparse(asOneSidedFormula(object@call$subset)[[2]]),"\n")
              }
              print(data.frame(AIC = object@AIC, BIC = object@BIC,
                               logLik = c(object@logLik), row.names = ""))
              cat("\n")
              object@re@useScale = TRUE
              object@re@showCorrelation = TRUE
              show(object@re)
              ## Should this be part of the show method for summary.reStruct?
              cat("\nNumber of Observations:", object@re@nobs)
              cat("\nNumber of Groups: ")
              ngrps <- object@re@ngrps
              if ((length(ngrps)) == 1) {
                  cat(ngrps,"\n")
              } else {				# multiple nesting
                  cat("\n")
                  print(ngrps)
              }
              invisible(object)
          })

setMethod("show", "lme",
          function(object)
      {
          sumry = summary(object)
          rdig <- 5
          cat("Linear mixed-effects model\n")
          cat(" Data:", deparse( sumry@call$data ), "\n")
          if (!is.null(sumry@call$subset)) {
              cat("  Subset:",
                  deparse(asOneSidedFormula(sumry@call$subset)[[2]]),"\n")
          }
          cat(paste(" log-", ifelse(sumry@re@REML, "restricted-", ""),
                    "likelihood: ", sep = ''), sumry@logLik, "\n")
          sumry@re@useScale = TRUE
          sumry@re@showCorrelation = FALSE
          saveopt = options(show.signif.stars=FALSE)
          on.exit(saveopt)
          show(sumry@re)
          options(saveopt)
          on.exit()
          cat("\nNumber of Observations:", sumry@re@nobs, "\n")
          invisible(object)
      })


setMethod("isInitialized", "lmeLevelList",
          function(object) all(sapply(object[seq(length=length(object)-2)],
                                      function(x) isInitialized(x@precision))),
          valueClass = "logical")

setMethod("anova", signature(object = "lme"),
          function(object, ...)
          cat("anova method for lme not yet implemented\n"))

setMethod("fixef", signature(object = "lme"),
          function(object, ...)
      {
          object = object@reStruct
          callGeneric()
      })

setMethod("formula", "lme", function(x, ...) x@call$formula)

setMethod("intervals", signature(object = "lme", level = "ANY"),
          function(object, level = 0.95, ...)
          cat("intervals method for lme not yet implemented\n"))

setMethod("plot", signature(x = "lme"),
          function(x, y, ...)
          cat("plot method for lme not yet implemented\n"))

setMethod("ranef", signature(object = "lme"),
          function(object, ...)
      {
          object = object@reStruct
          callGeneric()
      })

setMethod("coef", signature(object = "lme"),
          function(object, ...)
      {
          object = object@reStruct
          callGeneric()
      })

setMethod("update", signature(object = "lme"),
          function(object, formula., ..., evaluate = TRUE)
      {
          call <- object@call
          if (is.null(call))
              stop("need an object with call component")
          extras <- match.call(expand.dots = FALSE)$...
          if (!missing(formula.))
              call$formula <- update.formula(formula(object), formula.)
          if (length(extras) > 0) {
              existing <- !is.na(match(names(extras), names(call)))
              for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
              if (any(!existing)) {
                  call <- c(as.list(call), extras[!existing])
                  call <- as.call(call)
              }
          }
          if (evaluate)
              eval(call, parent.frame())
          else call
      })

setMethod("getGroups", signature(object="reStruct",
                                 form="missing",
                                 data="missing",
                                 sep="missing"),
          function(object, form, level, data, sep)
      {
          object <- object@reStruct
          callGeneric()
      })

setMethod("getResponse", signature(object="lme"),
          function(object, form)
      {
          object <- object@reStruct
          callGeneric()
      })

### Local variables:
### mode: R
### End:
### pdCompSymm: Compound symmetry structure

setGeneric('pdCompSymm',
           function(value, form, nam, data, ...)
           standardGeneric('pdCompSymm'))

setMethod("pdCompSymm",
          signature(value = 'formula', form = 'missing',
                    nam = 'missing', data = 'missing'),
          function(value, form, nam, data, ...) {
              new('pdCompSymm', form = value)
          })

setMethod("isInitialized", "pdCompSymm",
          function(object) (length(object@ncol) != 0))

setReplaceMethod("coef",
                 signature(object = "pdCompSymm", value = "numeric"),
                 function(object, value) {
                     .Call("pdCompSymm_coefGets", object, value,
                           PACKAGE = "lme4")
                 })

setAs('pdCompSymm', 'pdmatrix',
      function(from) {
          if (length(from@Ncol) == 1 && from@Ncol >= 1 &&
              length(from@param) == 2) {
              value <- exp(2 * from@param) * diag(from@Ncol)
              nam <- from@Names
              if (length(nam) == ncol(value)) {
                  dimnames(value) <- list(nam, nam)
              }
              return(value)
          } else {
              stop("Uninitialized pdCompSymm object")
          }
      })
### pdDiag - diagonal structure parameterized by the logarithm of
###   the square root of the diagonal terms.

setGeneric('pdDiag',
           function(value, form, nam, data, ...)
           standardGeneric('pdDiag'))

setMethod("pdDiag",
          signature(value = 'formula', form = 'missing',
                    nam = 'missing', data = 'missing'),
          function(value, form, nam, data, ...) {
              new('pdDiag', form = value)
          })

## Methods for the pdDiag class

setReplaceMethod("coef",
                 signature(object = "pdDiag", value = "numeric"),
                 function(object, value) {
                     .Call("pdDiag_coefGets", object, value, PACKAGE = "lme4")
#                     lenVal <- length(value)
#                     if (lenVal <= 0)
#                         stop('coef for a pdDiag object must have length > 0')
#                     lenPar <- length(object@param)
#                     if (lenPar == 0)
#                         lenPar <- length(object@Names)
#                     if (lenPar != 0 && lenPar != lenVal)
#                         stop("coef for a pdDiag object must be same as its number of rows")
#                     object@param <- value
#                     object@Ncol <- lenVal
#                     object@factor <- diag(exp(value), ncol = lenVal)
#                     object@logDet <- sum(value)
#                     object
                 })

setAs("pdDiag", "pdmatrix",
      function(from) {
          if (!isInitialized(from))
              stop(paste("Uninitialized", class(from), "object"))
          value <- diag(exp(2 * from@param), ncol = from@Ncol)
          nam <- from@Names
          if (length(nam) == length(from@param)) {
              dimnames(value) <- list(nam, nam)
          }
          new("pdmatrix", value)
      },
      function(from, value) {
          nc <- ncol(value)
          if (!identical(nc, dim(value)[2]))
              stop("value must be a square matrix")
          if (length(from@param) < 1) {
              from@Ncol <- nc
          }
          if (from@Ncol != nc)
              stop("can not change length of an initialized pdMat object")
          Names <- dimnames(value)[[2]]
          if (!is.null(Names))
              from@Names <- Names
          coef(from) <- log(diag(value))/2
          from
      })

setMethod("solve", signature(a="pdDiag", b="missing"),
          function(a, b) {
              if (!isInitialized(a))
                  stop(paste("Uninitialized", class(a), "object"))
              coef(a) <- -a@param
              a
          })

setMethod("LMEgradient",
          signature(x="pdDiag", A="matrix", nlev="numeric"),
          function(x, A, nlev) {
              .Call("pdDiag_LMEgradient", x, A, nlev, PACKAGE="lme4")
          })

setMethod("LMEhessian",
          signature(x="pdDiag", A="matrix", H="array",
                    nlev="numeric"),
          function(x, A, H, nlev)
      {
          theta <- x@param
          q <- length(theta)

          ## the part involving D_iD_j
          ans <-
              if (q > 1)
                  diag(-exp(2*theta)*colSums(A*A))
              else as.matrix(-exp(2*theta)*A*A)

          ## add the part not involving D_iD_j
          tmp <- matrix(0, nrow=q, ncol=q)
          for (j in 1:q)
              tmp[, j] <- H[seq(from = 1+(j-1)*q*q*(q+1), length=q, by = q+1)]
          ans <- ans + tmp*outer(theta, theta,
                                 function(v, w) exp(2*(v+w)))
          nm <- names(x)
          ans <- ans+t(ans)
          if (!is.null(nm))
              dimnames(ans) <- list(nm, nm)
          ans
      })

setReplaceMethod("EMupdate",
                 signature(x="pdDiag", nlev="numeric", value="matrix"),
                 function(x, nlev, value) {
                     .Call("pdDiag_EMupdate", x, nlev, value, PACKAGE="lme4")
                 })
### pdIdent: multiple of the identity matrix - the parameter is
###  half the logarithm of the multiple.

## constructor for pdIdent

setGeneric('pdIdent',
           function(value, form, nam, data, ...)
           standardGeneric('pdIdent'))

setMethod("pdIdent",
          signature(value = 'formula', form = 'missing',
                    nam = 'missing', data = 'missing'),
          function(value, form, nam, data, ...) {
              new('pdIdent', form = value)
          })

## methods for the pdIdent class

setReplaceMethod("coef",
                 signature(object = "pdIdent", value = "numeric"),
                 function(object, value) {
                     if (length(value) != 1)
                         stop("coef for pdIdent class must be of length 1")
                     if (length(object@Ncol) < 1 || object@Ncol < 1)
                         stop(paste("Changing parameter of uninitialized",
                                    class(object), "object"))
                     object@param <- value
                     diag(object@factor) <- exp(value)
                     object@logDet <- object@Ncol * object@param
                     object
                 })

setAs('pdIdent', 'pdmatrix',
      function(from) {
          if (!isInitialized(from))
              stop("Uninitialized pdIdent object")
          Ncol <- from@Ncol
          value <- diag(exp(2 * from@param), Ncol)
          nam <- from@Names
          if (length(nam) == Ncol) {
              dimnames(value) <- list(nam, nam)
          }
          new("pdmatrix", value)
      },
      function(from, value) {
          nc <- ncol(value)
          if (!identical(nc, dim(value)[1]))
              stop("value must be a square matrix")
          if (length(from@param) < 1) {
              from@Ncol <- nc
              from@factor <- diag(nrow = nc, ncol = nc)
          }
          if (from@Ncol != nc)
              stop("can not change length of an initialized pdMat object")
          if (length(colnames(value)) != nc) {
              if (length(from@Names) == 0)
                  from@Names <- paste("V", 1:nc, sep = "")
          } else {
              from@Names <- colnames(value)
          }
          from@param <- mean(log(diag(value)))/2
          diag(from@factor) <-  exp(from@param)
          from@logDet <- nc * from@param
          from
      })

setAs("pdIdent", "corrmatrix",
      function(from) {
          if (!isInitialized(from))
              stop(paste("Uninitialized", class(from), "object"))
          Ncol <- from@Ncol
          val <- diag(Ncol)
          stdDev <- rep(exp(from@param), Ncol)
          if (length(nm <- from@Names) == 0) {
              nm <- paste("V", 1:Ncol, sep = "")
          }
          dimnames(val) <- list(nm, nm)
          names(stdDev) <- nm
          new("corrmatrix", val, stdDev = stdDev)
      })

setMethod("solve", signature(a="pdIdent", b="missing"),
          function(a, b)
      {
          if (!isInitialized(a))
              stop(paste("Uninitialized", class(a), "object"))
          a@param <- -a@param
          diag(a@factor) <- exp(a@param)
          a@logDet <- a@Ncol * a@param
          a
      })

setMethod("LMEgradient",
          signature(x="pdIdent", A="matrix", nlev="numeric"),
          function(x, A, nlev) {
              .Call("pdIdent_gradient", x, A, nlev, PACKAGE="lme4")
          })

setReplaceMethod("EMupdate",
                 signature(x="pdIdent", nlev="numeric", value="matrix"),
                 function(x, nlev, value) {
                     .Call("pdIdent_EMupdate", x, nlev, value, PACKAGE="lme4")
                 })

setMethod("pdgradient", "pdIdent",
          function(x) {
              mat <- as(x, "pdmatrix")
              dn <- dimnames(mat)
              if (!is.null(dn)) dn <- c(list(dimnames(mat)), NULL)
              array(2. * c(mat), dim = c(dim(mat), 1), dimnames = dn)
          })
### pdLogChol - a general positive definite structure parameterized
###   by the non-zero elements of the Cholesky factor.  The logarithms of
###   the diagonal elements are the first Ncol elements of the parameter
###   vector

setGeneric('pdLogChol',
           function(value, form=formula(NULL), nam = character(), data=list(),
                    ...)
           standardGeneric('pdLogChol'))

setMethod('pdLogChol',
          signature(value = 'formula', form = 'missing',
                    nam = 'missing', data = 'missing'),
          function(value, form, nam, data, ...) {
              new('pdLogChol', form = value)
          })

## Methods for the pdLogChol class

setReplaceMethod("coef",
                 signature(object = "pdLogChol", value = "numeric"),
                 function(object, value) {
                     .Call("pdLogChol_coefGets", object, value, PACKAGE = "lme4")
#                     npar <- length(value)
#                     if (npar < 1)
#                         stop('coef for a pdLogChol object must have length > 0')
#                     if (npar != length(object@param)) {
#                         Ncol <- round((sqrt(8*length(value) + 1) - 1)/2)
#                         np <- (Ncol * (Ncol + 1))/2
#                         if (np != npar)
#                             stop(paste("coef for a pdLogChol object cannot have",
#                                        "length", npar))
#                         lenPar <- length(object@param)
#                         if (lenPar <= 0 && length(object@Names) > 0) {
#                             lenPar <- length(object@Names)
#                             lenPar <- (lenPar * (lenPar+1))/2
#                         }
#                         if (lenPar && lenPar != npar)
#                             stop("coef for a pdLogChol object has inconsistent length")
#                         object@Ncol <- as.integer(Ncol)
#                         object@factor <- matrix(0., Ncol, Ncol)
#                     }
#                     Ncol <- object@Ncol
#                     fact <- object@factor
#                     diag(fact) <- exp(value[1:Ncol])
#                     fact[row(fact) < col(fact)] <- value[-(1:Ncol)]
#                     object@param <- value
#                     object@factor <- fact
#                     object@logDet <- sum(value[1:Ncol])
#                     object
                 })

setAs("pdLogChol", "pdmatrix",
      function(from) new("pdmatrix", crossprod(from@factor)),
      function(from, value) {
          nc <- ncol(value)
          if (!identical(nc, dim(value)[2]))
              stop("value must be a square matrix")
          if (length(from@param) < 1) {
              from@Ncol <- nc
          }
          if (from@Ncol != nc)
              stop("can not change length of an initialized pdMat object")
          Names <- dimnames(value)[[2]]
          if (!is.null(Names))
              from@Names <- Names
          fact <- .Call("nlme_Chol", as(value, "matrix"), PACKAGE="lme4")
          from@factor <- fact
          from@logDet = sum(log(diag(fact)))
          from@param <- c(log(diag(fact)), fact[col(fact) > row(fact)])
          from
      })

setMethod("solve", signature(a="pdLogChol", b="missing"),
          function(a, b) {
              if (!isInitialized(a))
                  stop(paste("Uninitialized", class(a), "object"))
              as(a, "pdmatrix") <- crossprod(t(solve(a@factor)))
              a
          })

#setMethod("summary", signature(object="pdLogChol"),
#          function(object, structName, noCorrelation, ...) {
#              if (missing(structName)) structName =
#                   "General positive-definite, Log-Cholesky parametrization"
#              if (missing(noCorrelation)) noCorrelation = FALSE
#              callNextMethod()
#          })

setMethod("LMEgradient",
          signature(x="pdLogChol", A="matrix", nlev="numeric"),
          function(x, A, nlev)
          .Call("pdLogChol_LMEgradient", x, A, nlev, PACKAGE="lme4")
          )

setMethod("LMEhessian",
          signature(x="pdLogChol", A="matrix", H="array",
                    nlev="numeric"),
          function(x, A, H, nlev)
      {
          .Call("pdLogChol_LMEhessian", x, A, H, nlev, PACKAGE="lme4")
      })

setReplaceMethod("EMupdate",
                 signature(x="pdLogChol", nlev="numeric", value="matrix"),
                 function(x, nlev, value) {
                     .Call("pdLogChol_EMupdate", x, nlev, value, PACKAGE="lme4")
#                     if (!isInitialized(x))
#                         stop(paste("Uninitialized", class(x), "object"))
#                     if (any(dim(value) != dim(x)))
#                         stop(paste("value must be a matrix of dimension",
#                                    dim(x)))
#                     if (length(nlev) != 1 || nlev <= 0)
#                         stop("nlev must be > 0")
#                     as(x, "pdmatrix") <- nlev*crossprod(t(solve(value)))
#                     x
                 })

setMethod("pdgradient", "pdLogChol",
          function(x) {
              .Call("pdLogChol_pdgradient", x, PACKAGE="lme4")
#              fact <- as(x, "pdfactor")
#              pars <- x@param
#              Ncol <- ncol(fact)
#              dn <- dimnames(fact)
#              if (!is.null(dn)) dn <- c(list(dimnames(fact)), NULL)
#              val <- array(0., dim = c(dim(fact), length(pars)),
#                           dimnames = dn)
#              nc <- 1:Ncol
#              val[cbind(nc, nc, nc)] <- exp(pars[nc])
#              offdiag <- row(fact) < col(fact)
#              val[cbind(row(fact)[offdiag], col(fact)[offdiag],
#                  (Ncol+1):length(pars))] <- 1
#              for (i in seq(along = pars)) {
#                  pr <- crossprod(fact, val[,,i])
#                  val[,,i] <- pr + t(pr)
#              }
#              val
          })

## General methods for the virtual class pdMat

setMethod("formula", "pdMat", function(x, ...) x@form,
          valueClass = "formula")

setMethod("dim", "pdMat", function(x) c(x@Ncol, x@Ncol),
          valueClass = "integer")

setMethod("names", "pdMat", function(x) x@Names,
          valueClass = "character")

setMethod("coef", signature(object="pdMat"),
          function(object, ...) object@param,
          valueClass = "numeric")

setMethod("isInitialized", "pdMat",
          function(object) length(object@param) != 0,
          valueClass = "logical")

setReplaceMethod("names", "pdMat",
                 function(x, value) {
                     x@Names <- value
                     x
                 })

## The generics Names and Names<- will be deprecated in nlme_4.0
setMethod("Names", "ANY", function(object, ...) names(object))
setReplaceMethod("Names", "ANY",
                 function(object, value) "names<-"(object, value))

## In general we use coersion to the pdmatrix class as the pdMatrix method

setMethod("pdMatrix", signature(object="pdMat"),
          function(object) as(object, "pdmatrix"))

setAs("pdMat", "matrix",
      function(from) as(as(from, "pdmatrix"), "matrix"),
      function(from, value) {
          as(from, "pdmatrix") <- new("pdmatrix", value)
          from
      })

## In general we use coersion to the pdfactor class as the pdFactor method

setMethod("pdFactor", signature(object="pdMat"),
          function(object) as(object, "pdfactor"))

## This coersion should be overridden in classes with a simpler way of
## obtaining the factor

setAs("pdMat", "pdfactor",
      function(from) {
          if (!isInitialized(from))
              stop(paste("Uninitialized", class(from), "object"))
          new("pdfactor", from@factor, logDet = from@logDet)
      },
      function(from, value) {
          as(from, "pdmatrix") <- crossprod(value)
          from
      })

setAs("matrix", "pdmatrix",
      function(from) new("pdmatrix", from))

setAs("pdfactor", "pdmatrix",
      function(from) new("pdmatrix", crossprod(from)))

setAs("pdmatrix", "pdfactor",
      function(from) {
          val <- new("pdfactor",
                     .Call("nlme_Chol", as(from, "pdmatrix"), PACKAGE="lme4"))
          val@logDet <- sum(log(diag(val)))
          val
      },
      function(from, value) {
          as(from, "pdmatrix") <- crossprod(value)
          from
      })

## We prefer
##  as(pdm, "pdmatrix") <- val
## to
##  matrix(pdm) <- val
## The latter is included for backward compatibility and will be deprecated.

setReplaceMethod("matrix", signature(object="pdMat", value="matrix"),
                 function(object, value) {
                     as(object, "pdmatrix") <- value
                 })

## In general we use coersion to the pdfactor class as the pdFactor method

setMethod("pdFactor", signature(object="pdMat"),
          function(object) as(object, "pdfactor"),
          valueClass = "pdfactor")

setMethod("logDet", signature(object="pdDiag", covariate="missing"),
          function(object, covariate, ...) {
              if (!isInitialized(object))
                  stop(paste("Uninitialized", class(object), "object"))
              object@logDet
          },
          valueClass = "numeric")

## In general we use coersion to the corrmatrix class as the corMatrix method

setMethod("corMatrix", signature(object="pdMat"),
          function(object) as(object, "corrmatrix"),
          valueClass = "corrmatrix")

## This should be overridden in classes with a simpler way of obtaining
## the correlation matrix

setAs("pdMat", "corrmatrix",
      function(from) {
          Var  <- as(from, "pdmatrix")
          stdDev <- sqrt(diag(Var))
          colNames <- colnames(Var)
          if (is.null(colNames))
              colNames <- rownames(Var, do.NULL = FALSE, prefix="V")
          names(stdDev) <- colNames
          value <- array(t(Var/stdDev)/stdDev, dim(Var),
                         list(colNames, colNames))
          diag(value) <- rep(1.0, from@Ncol)
          new("corrmatrix", value, stdDev = stdDev)
      })

## gradient of the positive-definite matrix with respect to the parameters
## This method uses finite differences.  It should be overridden in
## explicit pdMat classes for which the gradient can be calculated explicitly

setMethod("pdgradient", "pdMat",
          function(x) {
              m0 <- as(x, "pdmatrix")
              pars <- coef(x)
              val <- array(0., c(dim(m0), length(pars)))
              eps <- sqrt(.Machine$double.eps)
              ind <- seq(along = pars)
              for (i in ind) {
                  coef(x) <- pars + eps * pars * (ind == i)
                  val[,,i] <- (as(x, "pdmatrix") - m0)/(pars[i] * eps)
              }
              val
          })

setMethod("show", "pdMat",
          function(object) {
              if (isInitialized(object)) {
                  cat("positive definite matrix of class",
                      class(object), "\n")
                  print(as(object, "pdmatrix"))
              } else {
                  cat(paste("Uninitialized", class(object), "object\n"))
              }
          })

setMethod("solve", signature(a="pdMat", b="missing"),
          function(a, b, ...) {
              if (!isInitialized(a))
                  stop(paste("Uninitialized", class(a),
                             "object", deparse(substitute(a))))
              as(a, "matrix") <- solve(as(a, "matrix"))
              a
          })

setMethod("summary", signature(object="pdMat"),
          function(object, structName, noCorrelation , ...) {
              cl = as.character(class(object))
              if (missing(structName))
                  structName = cl
              if (missing(noCorrelation))
                  noCorrelation = switch(cl, pdDiag = , pdIdent = TRUE, FALSE)
              if (isInitialized(object)) {
                  new("summary.pdMat", cor = as(object, "corrmatrix"),
                      structName = structName,
                      noCorrelation = noCorrelation,
                      formula = formula(object))
              } else {
                  object
              }
          })

setAs("pdmatrix", "corrmatrix",
      function(from) {
          ss = sqrt(diag(from))
          new("corrmatrix", t(from/ss)/ss, stdDev = ss)
      })

## This function is for testing only.  It will eventually be removed.
##
pdgradNumeric <- function(x)
{
    m0 <- as(x, "pdmatrix")
    pars <- coef(x)
    val <- array(0., c(dim(m0), length(pars)))
    eps <- sqrt(.Machine$double.eps)
    ind <- seq(along = pars)
    for (i in ind) {
        coef(x) <- pars + eps * pars * (ind == i)
        val[,,i] <- (as(x, "pdmatrix") - m0)/(pars[i] * eps)
    }
    val
}

### pdNatural - a general positive definite structure parameterized
###   by the log of the square root of the diagonal elements and the
###   generalized logit of the correlations. This is NOT an unrestricted
###   parametrization

setGeneric('pdNatural',
           function(value, form=formula(NULL), nam = character(),
                    data=list(), ...)
           standardGeneric('pdNatural'))

setMethod('pdNatural',
          signature(value = 'formula', form = 'missing',
                    nam = 'missing', data = 'missing'),
          function(value, form, nam, data, ...) {
              new('pdNatural', form = value)
          })

setMethod('pdNatural',
          signature(value = 'pdMat', form = 'missing',
                    nam = 'missing', data = 'missing'),
          function(value, form, nam, data) {
              val <- new('pdNatural', form = value@form, Names = value@Names)
              as(val, "corrmatrix") <- as(value, "corrmatrix")
              val
          })


setReplaceMethod("coef",
                 signature(object = "pdNatural", value = "numeric"),
                 function(object, value) {
                     npar <- length(value)
                     lenPar <- length(object@param)
                     if (npar != lenPar) {
                         Ncol <- round((sqrt(8*length(value) + 1) - 1)/2)
                         np <- (Ncol * (Ncol + 1))/2
                         if (np != npar)
                             stop(paste("coef for a pdNatural object cannot have",
                                        "length", npar))
                         if (lenPar <= 0 && length(object@Names) > 0) {
                             lenPar <- length(object@Names)
                             lenPar <- (lenPar * (lenPar+1))/2
                         }
                         if (lenPar && lenPar != npar)
                             stop("coef for a pdNatural object has inconsistent length")
                         object@Ncol <- as.integer(Ncol)
                         object@factor <- matrix(0., Ncol, Ncol)
                     }
                     Ncol <- object@Ncol
                     corr <- diag(nrow = Ncol, ncol = Ncol)
                     expz <- exp(value[-(1:Ncol)])
                     corr[upper.tri(corr)] <- (expz - 1)/(1 + expz)
                     corr[lower.tri(corr)] <- t(corr)[lower.tri(corr)]
                     corrFact <- La.chol(corr)
                     object@factor <- t(exp(value[1:Ncol]) * t(corrFact))
                     object@param <- value
                     object@logDet <-
                         sum(log(diag(corrFact)))+sum(value[1:Ncol])
                     object
                 })

setAs("pdNatural", "corrmatrix",
      function(from) {
          if (!isInitialized(from))
              stop(paste("Uninitialized", class(from), "object"))
          value <- .Call("pdNatural_corrmatrix", from, PACKAGE = "lme4")
          if (length(from@Names) == from@Ncol)
              dimnames(value) <- list(from@Names, from@Names)
          new("corrmatrix", value, stdDev = exp(from@param[1:from@Ncol]))
      },
      function(from, value) {
          nc <- ncol(value)
          if (!identical(nc, dim(value)[1]))
              stop("value must be a square matrix")
          if (length(from@param) < 1) {
              from@Ncol <- nc
          }
          if (from@Ncol != nc)
              stop("can not change length of an initialized pdMat object")
          Names <- dimnames(value)[[2]]
          if (!is.null(Names))
              from@Names <- Names
          rho <- value[upper.tri(value)]
          coef(from) <- c(log(value@stdDev), log((1.+rho)/(1.-rho)))
          from
      })

setAs("pdNatural", "pdmatrix",
      function(from) {
          if (!isInitialized(from))
              stop(paste("Uninitialized", class(from), "object"))
          new("pdmatrix", .Call("pdNatural_pdmatrix", from, PACKAGE="lme4"))
      },
      function(from, value) {
          stdDev <- sqrt(diag(value))
          value <- array(t(value/stdDev)/stdDev, dim(value))
          colNames <- colnames(value)
          diag(value) <- rep(1.0, from@Ncol)
          if (!is.null(colNames)) {
              dimnames(value) <- list(colNames, colNames)
              names(stdDev) <- colNames
          }
          as(from, "corrmatrix") <-
              new("corrmatrix", value, stdDev = stdDev)
          from
      })

setMethod("solve", signature(a = "pdNatural", b = "missing"),
          function(a, b) {
              if (!isInitialized(a))
                  stop(paste("Unitialized", class(a), "object"))
              Ncol <- a@Ncol
              if (Ncol > 1) {
                  as(a, "pdmatrix") <- solve(as(a, "pdmatrix"))
              } else {
                  coef(a) <- -coef(a)
              }
              a
          })

setMethod("LMEgradient",
          signature(x="pdNatural", A="matrix", nlev="numeric"),
          function(x, A, nlev)
          .Call("pdNatural_LMEgradient", x, A, nlev, PACKAGE="lme4")
          )

#setMethod("summary", signature(object="pdNatural"),
#          function(object,
#                   structName = "General positive-definite, Natural parametrization",
#                   ...) {
#              callNextMethod()
#          })
if (!isGeneric("lmeLevel")) {
    setGeneric("lmeLevel",
               function(precision, groups, columns, modelMatrix)
               standardGeneric("lmeLevel"))
}

setMethod("lmeLevel", signature(precision="pdMat",
                               groups="factor",
                               columns="integer",
                               modelMatrix="matrix"),
          function(precision, groups, columns, modelMatrix)
      {
          groups <- factor(groups, levels=unique(groups))
          nlev <- length(levels(groups))
          n <- length(columns)
          mat <- matrix(0.0, nrow=n, ncol=n)
          diag(mat) <- (9/64) * diag(crossprod(modelMatrix[, columns]))/nlev
          col.names <- dimnames(modelMatrix)[[2]][columns]
          dimnames(mat) <- list(col.names, col.names)
          as(precision, "matrix") <- mat
          ans <- new("lmeLevel",
                     precision=precision,
                     groups=groups,
                     columns=columns,
                     originalRows=split(as.integer(seq(length=length(groups))),
                     groups),
                     decomposedRows=vector("list", nlev),
                     storedRows=vector("list", nlev),
                     nrow=as.integer(n),
                     updateFactor=matrix(0.0, nrow=n, ncol=n),
                     nlev=as.integer(nlev))
          ans
      })

setMethod("lmeLevel", signature(precision="missing",
                               groups="missing",
                               columns="integer",
                               modelMatrix="matrix"),
          function(precision, groups, columns, modelMatrix)
      {
          ans <- new("lmeLevel",
                     precision=new("pdLogChol"),
                     groups=factor("a")[-1,drop=TRUE],
                     columns=columns,
                     originalRows=list(as.integer(seq(1, nrow(modelMatrix)))),
                     decomposedRows=vector("list", 1),
                     storedRows=vector("list", 1),
                     nrow=as.integer(0),
                     nlev=as.integer(1))
          ans
      })

setMethod("coef", signature(object="lmeLevel"),
          function(object, ...) coef(object@precision))

setReplaceMethod("coef", signature(object="lmeLevel", value="numeric"),
                 function(object, value)
             {
                 coef(object@precision) <- value
                 object
             })

setMethod("LMEhessian", signature(x="lmeLevel", A="missing",
                                  H="missing",
                                  nlev="missing"),
          function(x, A, H, nlev)
      {
          LMEhessian(x@precision, x@updateFactor, x@hessianArray,
                     nlev=x@nlev)
      })

setMethod("LMEgradient", signature(x="lmeLevel", A="missing", nlev="missing"),
          function(x, A, nlev)
      {
          LMEgradient(x@precision, x@updateFactor, x@nlev)
      })

setMethod("LMEgradient", signature(x="lmeLevel", A="matrix", nlev="missing"),
          function(x, A, nlev)
      {
          LMEgradient(x@precision, A, x@nlev)
      })

setReplaceMethod("EMupdate", signature(x="lmeLevel", nlev="missing", value="matrix"),
                 function(x, nlev, value)
             {
                 EMupdate(x@precision, x@nlev) <- value
                 x
             })
###      Methods for the class of random-effects structures.

setMethod("reStruct", signature(fixed = "formula",
                                random = "list",
                                data = "data.frame",
                                weights = "missing",
                                REML = "logical"),
          function(fixed, random, data, weights, REML, nextraCols,
                   analyticHessian)
      {
          reStruct(fixed, random, data, weights=numeric(0), REML,
                                nextraCols, analyticHessian)
      })

setMethod("reStruct", signature(fixed = "formula",
                                random = "list",
                                data = "data.frame",
                                weights = "numeric",
                                REML = "logical"),
          function(fixed, random, data, weights, REML, nextraCols,
                   analyticHessian)
      {
          ## given a matrix of no. of rows needed in the last row of a
          ## level for each column, create the storedRows or
          ## decomposedRows list
          createRows <-
              function(rowMatrixEnd)
              {
                  rowMatrixBeg <- as.vector(rowMatrixEnd)
                  nonzero <- rowMatrixBeg[rowMatrixBeg != 0]
                  rowMatrixBeg[rowMatrixBeg != 0] <-
                      c(1, nonzero[-length(nonzero)]+1)
                  dim(rowMatrixBeg) <- dim(rowMatrixEnd)
                  value <- lapply(seq(length=nrow(rowMatrixBeg)),
                                  function(i)
                              {
                                  ind <- rowMatrixBeg[i, ] != 0
                                  beg <- rowMatrixBeg[i, ind]
                                  end <- rowMatrixEnd[i, ind]
                                  lapply(seq(along=beg),
                                         function(i) seq(from=beg[i],
                                                         to=end[i]))
                              })
                  names(value) <- dimnames(rowMatrixEnd)[[1]]
                  value
              }
          ## order from inner to outer groups
          random <- rev(random)
          analyticHessian <- as.logical(analyticHessian)
          self <- new("reStruct",
                      fixed=fixed,
                      REML=REML,
                      analyticHessian=analyticHessian)

          ## Get the grouping factors
          groups <- as.data.frame(lapply(names(random),
                                         function(expr)
                                         eval(parse(text=expr),
                                              data,
                                              environment(fixed))))
          for (i in seq(to=1, length=length(groups)-1, by = -1)) {
              groups[, i] <-
                  as.factor(paste(as.character(groups[[i+1]]),
                                  as.character(groups[[i]]),
                                  sep = "/"))
          }
          names(groups) <- names(random)
          ## save the old ordering
          row.names(groups) <- row.names(data)
          ## generate the new ordering for the model.frame
          self@origOrder <- do.call("order", groups)
          groups[, 1] <- as.factor(groups[,1])
          self@reverseOrder <- order(self@origOrder)
          ## Reorder the model.frame
          data[,] <- data[self@origOrder, , drop = FALSE]
          row.names(data) <- row.names(groups)[self@origOrder]
          groups[,] <- groups[self@origOrder, , drop = FALSE]
          row.names(groups) <- row.names(data)
          ## Create the model.matrix for fixed and response
          self@original <- model.matrix(fixed, data)
          self@assign.X <- attr(self@original, "assign")
          self@original <- cbind(self@original, model.response(data))
          nCol <- ncol(self@original)
          p <- nCol-1
          ## Set the column name for response
          colnames(self@original)[nCol] <-
              colnames(data)[attr(attr(data, "terms"),
                                        "response")]

          qVector <- integer(length(random))
          for (i in seq(length=length(random),
                        by = -1, to = 1)) {
              self@original <- cbind(model.matrix(formula(random[[i]]),
                                                  data),
                                     self@original)
              qVector[i] <- ncol(self@original)-nCol
              nCol <- ncol(self@original)
          }

          nextraCols <- as.integer(nextraCols)
          if (nextraCols > 0) {
              self@original <- cbind(self@original,
                                     matrix(as.numeric(NA),
                                            nrow=nrow(self@original),
                                            ncol=nextraCols))
              nCol <- ncol(self@original)
          }

          ncols <- qVector
          qVector <- c(qVector, p, 1)
          indx <- seq(along=random)
          names(indx) <- names(random)
          reStructColumns <- lapply(seq(along=qVector),
                                    function(i, end=cumsum(qVector))
                                    as.integer(seq(length=qVector[i],
                                                   to=end[i])))
          if (length(weights) > 0) {
              weighted(self) <- weights[self@origOrder]
              rm(weights)
          }
          Q <- ncol(groups)
          random <- lapply(indx, function(i)
                           lmeLevel(random[[i]],
                                    groups[[i]],
                                    reStructColumns[[i]],
                                    if (self@useWeighted)
                                    self@weighted
                                    else self@original))
          if (analyticHessian) {
              for (i in seq(along=random)) {
                  q <- length(random[[i]]@columns)
                  random[[i]]@hessianArray <-
                      array(0.0, rep(q, 4))
              }
          }
          class(random) <- "lmeLevelList"
          self@random <- random
          rm(random)
          ## FIXME - there should be a less grubby way of doing this
          pars <- lapply(self@random, coef)
          maxInd <- 0
          for (i in seq(along = pars)) {
              self@random[[i]]@parsInd <- as.integer(maxInd + seq(along = pars[[i]]))
              maxInd <- maxInd + length(pars[[i]])
          }

          self@random[["*fixed*"]] <- lmeLevel(columns=reStructColumns[[Q+1]],
                                               modelMatrix=self@original)
          self@random[["*response*"]] <- lmeLevel(columns=reStructColumns[[Q+2]],
                                                  modelMatrix=self@original)
          N <- nrow(groups)

          self <- .Call("nlme_reStructDims", self, PACKAGE="lme4")

          self@decomposed <- self@stored <-
              matrix(0, nrow = self@random[[Q+2]]@storedRows[[1]], ncol=nCol)
          self@bbetas <- numeric(nrow(self@stored))

          self
      })

setMethod("coef", "reStruct",
          function(object, ...)
          unlist(lapply(object@random, coef)))

setReplaceMethod("coef", signature(object="reStruct",
                                   value="numeric"),
                 function(object, value)
             {
                 names(value) = NULL
                 for (i in seq(length=length(object@random)-2)) {
                     coef(object@random[[i]]) =
                         value[object@random[[i]]@parsInd]
                 }
                 object@logLik = as.numeric(NA)
                 object
             })

setMethod("getResponse", signature(object="reStruct"),
          function(object, form)
      {
          object@original[, object@random[["*response*"]]@columns]
      })

setReplaceMethod("response", signature(x="reStruct", value="numeric"),
                 function(x, value)
             {
                 if (length(value) != nrow(x@original))
                     stop("Dimension mismatch in model.matrix")
                 x@original[, x@random[["*response*"]]@columns] = value
                 x@dirtyDecomposed = TRUE
                 x@useWeighted = FALSE
                 x
             })

setReplaceMethod("weighted", signature(x="reStruct",
                                       value="numeric"),
                 function(x, value)
             {
                 if (all(value == 1.0)) {
                     if (x@useWeighted) {
                         x@useWeighted = FALSE
                         x@dirtyDecomposed = TRUE
                     }
                 } else {
                     ## make sure the row number of original
                     ## model.matrix match length of weights
                     if (nrow(x@original) != length(value) &&
                         length(value) != 1)
                         stop("Dimension mismatch in setting weighted")
                     if (identical(dim(x@weighted), x@original))
                         x@weighted[] = x@original * value
                     else x@weighted = x@original * value
                     x@useWeighted = TRUE
                     x@dirtyDecomposed = TRUE
                 }
                 x
             })
setMethod("model.matrix", "reStruct",
          function(object, ...) object@original)

setReplaceMethod("model.matrix", signature(x="reStruct",
                                           value="matrix"),
                 function(x, value)
             {
                 ## FIXME: check that value is a model.matrix

                 ## make sure the dimensions of old and new
                 ## model.matrix match
                 if (!identical(dim(value), dim(x@original)))
                     stop("Dimension mismatch in model.matrix")
                 x@original = value
                 x@dirtyDecomposed = TRUE
                 x@useWeighted = FALSE
                 x
             })

setMethod("getGroups", signature(object="reStruct",
                                 form="missing",
                                 data="missing",
                                 sep="missing"),
          function(object, form, level, data, sep)
      {
          Q <- length(object@random)-2
          if (missing(level))
              level <- Q
          else if (length(level) > Q || max(level) > Q)
              stop ("Invalid value for level")
          val <- lapply(object@random[seq(length=Q)[level]],
                        function(ranlev) ranlev@groups)
          if (length(val) == 1) { # single group
              valname <- names(val)
              val <- val[[1]]
              attr(val, "label") <- valname
              val
          } else {
              as.data.frame(val)
          }
      })

setMethod("logLik", signature(object="reStruct"),
          function(object)
      {
          value = .Call("nlme_logLikelihood", object,
                         NULL,          # do not pass new parameter value
                         PACKAGE="lme4")
          p = length(object@random[["*fixed*"]]@columns)
          # df calculated from sigma + fixed effects + random effects pars
          attr(value, "df") = 1 + p +
              sum(unlist(lapply(object@random, function(x)length(coef(x)))))
          attr(value, "nall") = nrow(object@original)
          attr(value, "nobs") = nrow(object@original) - p
          class(value) = "logLik"
          value
      })

setReplaceMethod("EMsteps", signature(x="reStruct", value="list"),
                 function(x, value)
             {
                 x <- .Call("nlme_reStructEMsteps", x, value$niterEM,
                            value$EMverbose,
                            PACKAGE="lme4")
#                  verbose = value$EMverbose
#                  randIndx = seq(length=length(x@random)-2)
#                  for (i in seq(length=value$niterEM)) {
#                      x = .Call("nlme_commonDecompose", x, NULL,
#                                 PACKAGE="lme4")
#                      for (j in randIndx)
#                          EMupdate(x@random[[j]]) =
#                              x@random[[j]]@updateFactor
#                      if (verbose) {
#                          cat("\n**EM Iteration", i,
#                              .Call("nlme_logLikelihood", x,
#                                    NULL, # do not pass new parameter value
#                                    PACKAGE="lme4"), coef(x), "\n")
#                      }
#                      x@logLik = as.numeric(NA)
#                  }
                 .Call("nlme_commonDecompose", x, NULL,
                       PACKAGE="lme4")
             })

setReplaceMethod("LMEoptimize", signature(x="reStruct",
                                          value="list"),
                 function(x, value)
             {
                 if (value$msMaxIter < 1)
                     return(x)
                 xval = -.Call("nlme_logLikelihood",
                                x,
                                NULL,  # use already set parameters
                                PACKAGE="lme4")
                 xval =
                     if (xval > 0)
                         xval+1
                     else abs(min(xval/2, xval+1))
                 if (value$optimizer == "optim") {
                     optimRes =
                         if (value$analyticGradient) {
                             optim(fn = function(params)
                                   .Call("nlme_logLikelihood",
                                         x,
                                         params,  # new parameter value
                                         PACKAGE="lme4"),
                                   gr = function(params)
                                   LMEgradient(.Call("nlme_commonDecompose",
                                                     x, params,
                                                     PACKAGE="lme4")),
                                   par = c(coef(x)), #hessian = TRUE,
                                   method = "BFGS",
                                   control = list(trace = value$msVerbose,
                                   reltol = value$msTol,
                                   fnscale = -1,
#                                   fnscale = -xval,
#                                   parscale = 1/value$msScale(coef(x)),
                                   maxit = value$msMaxIter))
                         } else {
                             optim(fn = function(params)
                                   .Call("nlme_logLikelihood",
                                         x,
                                         params,  # new parameter value
                                         PACKAGE="lme4"),
                                   par = c(coef(x)), #hessian = TRUE,
                                   method = "BFGS",
                                   control = list(trace = value$msVerbose,
                                   reltol = value$msTol,
                                   fnscale = -1,
#                                   fnscale = -xval,
#                                   parscale = 1/value$msScale(coef(x)),
                                   maxit = value$msMaxIter))
                         }
                     if (optimRes$convergence != 0) {
                         warning("optim failed to converge")
                     }
                     .Call("nlme_commonDecompose", x, optimRes$par,
                           PACKAGE="lme4")
#                  } else if (value$optimizer == "ms") {
#                      pars <- coef(x)
#                      .Call("nlme_msOptimize", value$msMaxIter,
#                            value$msTol, rep(1.0, length(pars)),
#                            value$msVerbose, x, pars,
#                            value$analyticGradient,
#                            PACKAGE = "lme4")
                 } else {
#                     typsize <- 1/value$msScale(coef(x))
                     typsize <- rep(1.0, length(coef(x)))
                     if (is.null(value$nlmStepMax))
                         value$nlmStepMax <-
                             max(100 * sqrt(sum((coef(x)/typsize)^2)), 100)
                     nlmRes =
                         nlm(f = if (value$analyticGradient) {
                             function(params)
                             {
                                 x = .Call("nlme_commonDecompose",
                                            x, params,
                                            PACKAGE="lme4")
                                 grad = -LMEgradient(x)
                                 ans = -x@logLik
                                 attr(ans, "gradient") = grad
                                 ans
                             }
                         } else {
                             function(params)
                             {
                                 -.Call("nlme_logLikelihood",
                                        x,
                                        params,
                                        PACKAGE="lme4")
                             }
                         },
                             p = c(coef(x)), #hessian = TRUE,
                             print.level = if (value$msVerbose) 2 else 0,
                             steptol = value$msTol,
                             gradtol = value$msTol,
                             stepmax = value$nlmStepMax,
                             typsize=typsize,
#                             fscale=xval,
                             iterlim = value$msMaxIter)
                     .Call("nlme_commonDecompose", x, nlmRes$estimate,
                           PACKAGE="lme4")
                 }
             })


setMethod("LMEgradient", signature(x="reStruct", A="missing", nlev="missing"),
          function(x, A, nlev)
      {
          unlist(lapply(x@random[seq(length=length(x@random)-2)],
                        LMEgradient))
      })

setMethod("LMEhessian", signature(x="reStruct", A="missing",
                                  H="missing",
                                  nlev="missing"),
          function(x, A, H, nlev)
      {
          if (!x@analyticHessian)
              stop("Can not calculate analytic hessian")
          x <- .Call("nlme_commonDecompose", x, NULL, PACKAGE="lme4")
          lapply(x@random[seq(length=length(x@random)-2)],
                 LMEhessian)
      })

setMethod("fitted", signature(object="reStruct"),
          function(object, ...)
      {
          .Call("nlme_reStruct_fitted", object, NULL, PACKAGE="lme4")
      })

setMethod("fixef", signature(object="reStruct"),
          function(object, ...) {
              fixd = object@random[['*fixed*']]
              val = object@bbetas[fixd@storedRows[[1]]]
              nn = dimnames(object@original)[[2]][fixd@columns]
              if (length(nn) == length(val)) names(val) = nn
              val
          })

setReplaceMethod("fixef", signature(object="reStruct", value="numeric"),
          function(object, value) {
              fixdRows = object@random[['*fixed*']]@storedRows[[1]]
              object@bbetas[fixdRows] <- value
              object
          })

setMethod("ranef", signature(object="reStruct"),
          function(object, ...) {
              lapply(object@random[-(length(object@random) - c(1,0))],
                     function(x) matrix(object@bbetas[unlist(x@storedRows)],
                                        ncol = length(x@columns),
                                        byrow = TRUE))
          })

setMethod("summary", "reStruct",
          function(object, ...) {
              fixd = object@random[['*fixed*']]
              fstrRows = fixd@storedRows[[1]]
              fcols = fixd@columns
              rsp = object@random[['*response*']]
              sigma = abs(object@bbetas[rsp@storedRows[[1]]])
              nobs = dim(object@original)[1]
              nfixed = length(fixd@columns)
              denomDF = nobs - ifelse(object@REML, nfixed, 0)
              sigma = sigma/sqrt(denomDF)
              fcoef = object@bbetas[fstrRows]
              DF = getFixDF(object)$X
              rinv =
                  .Call("nlme_commonDecompose", object, NULL,
                        PACKAGE="lme4")@stored[fstrRows, fcols, drop=FALSE]
              se = sqrt(rowSums(rinv * rinv))
              corF = new("corrmatrix", crossprod(t(rinv/se)), stdDev = se)
              coefs = cbind(fcoef, se, DF)
              nn = dimnames(object@original)[[2]][fcols]
              dimnames(coefs) =
                  list(nn, c("Estimate", "Std. Error", "DF"))
              rnd = rev(object@random)[-(1:2)]
              ngrps = sapply(rnd, function(x) x@nlev)
              new("summary.reStruct",
                  fixed = object@fixed,
                  coefficients = as.matrix(coefs),
                  scale = sigma,
                  denomDF = as.integer(denomDF),
                  REML = object@REML,
                  ngrps = ngrps,
                  nobs = nobs,
                  corFixed = corF,
                  reSumry = lapply(rnd,
                      function(x) summary(solve(x@precision))))
          })

setMethod("getFixDF", signature(object="reStruct"),
          function(object)
      {
          ## calculates degrees of freedom for fixed effects Wald tests
#          Q <- length(object@random)-2
#          columns = object@random[["*fixed*"]]@columns
#          X = object@original[, columns, drop = FALSE]
#          ngrps = unlist(lapply(object@random, function(lmeLevel)
#                                 lmeLevel@nlev))
#          names(ngrps) = names(object@random)
          val = .Call("nlme_getFixDF", object, PACKAGE = "lme4")
          names(val$X) =
              colnames(object@original)[object@random[["*fixed*"]]@columns]
                                        # Convert R's assign to S-PLUS style
          assign = object@assign.X
          terms = terms(object@fixed)
          namTerms = attr(terms, "term.labels")
          if (attr(terms, "intercept") > 0) {
              namTerms = c("(Intercept)", namTerms)
          }
          names(val$terms) = namTerms
          namTerms = factor(assign, labels = namTerms)
          attr(val, "assign") = split(order(assign), namTerms)
          val
      })
#          N <- nrow(X)
#          p <- ncol(X)
#          Qp1 <- Q + 1
#          namX <- colnames(X)
#          ngrps <- rev(ngrps)[-(1:2)]
#          stratNam <- c(names(ngrps), "Residual")
#          dfX <- dfTerms <- c(ngrps, N) - c(0, ngrps)
#          names(dfX) <- names(dfTerms) <- stratNam
#          valX <- double(p)
#          names(valX) <- namX
#          namTerms <- names(assign)
#          valTerms <- double(length(assign))
#          names(valTerms) <- namTerms
#          if (any(notIntX <- apply(X, 2, function(el) any(el != el[1])))) {
#              ## percentage of groups for which columns of X are inner
#              innP <- .Call("nlme_inner_perc_table",
#                            object,
#                            PACKAGE = "lme4")
#              dimnames(innP) <- list(namX, stratNam)
#              ## strata in which columns of X are estimated
#              ## ignoring fractional inner percentages for now
#              stratX <- stratNam[apply(innP, 1, function(el, index) max(index[el > 0]),
#                                       index = 1:Qp1)]
#              ## strata in which terms are estimated
#              notIntTerms <- unlist(lapply(assign,
#                                           function(el, notIntX) {
#                                               any(notIntX[el])
#                                           }, notIntX = notIntX))
#              stratTerms <- stratNam[unlist(lapply(assign,
#                                                   function(el, stratX, stratNam) {
#                                                       max(match(stratX[el], stratNam))
#                                                   },
#                                                   stratX = stratX, stratNam = stratNam))][notIntTerms]
#              stratX <- stratX[notIntX]
#              xDF <- table(stratX)
#              dfX[names(xDF)] <- dfX[names(xDF)] - xDF
#              if (!all(notIntX)) {                # correcting df for intercept
#                  dfX[1] <- dfX[1] - 1
#              } else {
#                  dfX[-1] <- dfX[-1] + 1
#              }
#              valX[notIntX] <- dfX[stratX]
#              ## number of parameters in each term
#              pTerms <- unlist(lapply(assign, length))[notIntTerms]
#              tDF <- tapply(pTerms, stratTerms, sum)
#              dfTerms[names(tDF)] <- dfTerms[names(tDF)] - tDF
#              if (!all(notIntTerms)) {
#                  dfTerms[1] <- dfTerms[1] - 1
#              } else {
#                  dfTerms[-1] <- dfTerms[-1] + 1
#              }
#              valTerms[notIntTerms] <- dfTerms[stratTerms]
#          } else {
#              notIntTerms <- unlist(lapply(assign,
#                                           function(el, notIntX) {
#                                               any(notIntX[el])
#                                           }, notIntX = notIntX))
#          }
#          if (!all(notIntX)) {  #intercept included
#              valX[!notIntX] <- max(dfX)
#              if (!all(notIntTerms)) {
#                  valTerms[!notIntTerms] <- max(dfTerms)
#              }
#          }
#          val <- list(X = valX, terms = valTerms)
#          attr(val, "assign") <- assign
#          val
#      })

setMethod("show", signature(object="summary.reStruct"),
          function(object) {
              digits = max(3, getOption("digits") - 2)
              useScale = length(object@useScale) > 0 && object@useScale[1]
              sc = ifelse(useScale, object@scale,  1.)
              reStdDev = lapply(object@reSumry, function(x, sc) sc*x@cor@stdDev,
                                sc = sc)
              reLens = unlist(lapply(reStdDev, length))
              reMat = array('', c(sum(reLens), 4),
                             list(rep('', sum(reLens)),
                                  c("Groups", "Name", "Variance", "Std.Dev.")))
              reMat[1+cumsum(reLens)-reLens, 1] = names(reLens)
              reMat[,2] = unlist(lapply(reStdDev, names))
              reMat[,3] = format(unlist(reStdDev)^2, digits = digits)
              reMat[,4] = format(unlist(reStdDev), digits = digits)
              if (any(reLens > 1) &&
                  !all(sapply(object@reSumry,
                              function(x) x@noCorrelation))) {
                  maxlen = max(reLens)
                  corr =
                      do.call("rbind",
                              lapply(object@reSumry,
                                     function(x, maxlen) {
                                         if (x@noCorrelation) {
                                             matrix("", dim(x@cor)[1], maxlen)
                                         } else {
                                             cc = format(round(x@cor, 3),
                                                         nsmall = 3)
                                             cc[!lower.tri(cc)] = ""
                                             nr = dim(cc)[1]
                                             cbind(cc, matrix("",
                                                              nr, maxlen-nr))
                                         }
                                     }, maxlen))
                  colnames(corr) = c("Corr", rep("", maxlen - 1))
                  reMat = cbind(reMat, corr)
              }
              if (useScale) {
                  reMat = rbind(reMat, c("Residual", "",
                                          format(sc^2, digits = digits),
                                          format(sc, digits = digits),
                                          rep('', ncol(reMat) - 4)))
              }
              cat("Random effects:\n")
              print(reMat, quote = FALSE)
              cm = object@coefficients
              if (useScale) {
                  cm[,"Std. Error"] = sc * cm[, "Std. Error"]
                  stat = cm[,1]/cm[,2]
                  pval = 2*pt(abs(stat), cm[,3], lower = FALSE)
                  nms = colnames(cm)
                  cm = cbind(cm, stat, pval)
                  colnames(cm) = c(nms, "t value", "Pr(>|t|)")
              } else {
                  cat("\nEstimated scale (compare to 1) ", object@scale, "\n")
                  stat = cm[,1]/cm[,2]
                  pval = 2*pnorm(abs(stat), lower = FALSE)
                  nms = colnames(cm)
                  cm = cbind(cm, stat, pval)
                  colnames(cm) = c(nms, "z value", "Pr(>|z|)")
              }
              cat("\nFixed effects:",
                  paste(deparse(object@fixed),
                        sep = '\n', collapse = '\n'), "\n")
              if (R.version$major > 1 || R.version$major == "1" &&
                  R.version$minor >= 8) {
                  printCoefmat(cm, tst.ind = 4, zap.ind = 3)
              } else {
                  print.coefmat(cm, tst.ind = 4, zap.ind = 3)
              }
              if (length(object@showCorrelation) > 0 && object@showCorrelation[1]) {
                  correl = object@corFixed
                  rn = rownames(cm)
                  dimnames(correl) = list(
                          abbreviate(rn, minlen=11), abbreviate(rn, minlen=6))
                  if (!is.null(correl)) {
                      p = NCOL(correl)
                      if (p > 1) {
                          cat("\nCorrelation of Fixed Effects:\n")
                          correl = format(round(correl, 3), nsmall = 3)
                          correl[!lower.tri(correl)] = ""
                          print(correl[-1, -p, drop=FALSE], quote = FALSE)
                      }
                  }
              }
          })


### Local variables:
### mode: R
### End:
.First.lib <- function(lib, pkg) {
    if ("package:nlme" %in% search()) {
        stop(paste("Package lme4 conflicts with package nlme.\n",
                   "To attach lme4 you must restart R without package nlme."))
    }
    library.dynam(pkg, pkg, lib)
    require(lattice)
    cat(paste(" This package is in development.  For production work use\n",
              "lme from package nlme or glmmPQL from package MASS.\n"))
}
