R/gls.R

In nlme: Linear and Nonlinear Mixed Effects Models

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

gls <-
    ## fits linear model with serial correlation and variance functions,
    ## by maximum likelihood using a Newton-Raphson algorithm.
    function(model,
             data = sys.frame(sys.parent()),
             correlation = NULL,
             weights = NULL,
             subset,
             method = c("REML", "ML"),
             na.action = na.fail,
             control = list(),
             verbose = FALSE)
{
    Call <- match.call()
    ## control parameters
    controlvals <- glsControl()
    if (!missing(control)) controlvals[names(control)] <- control

    ##
    ## checking arguments
    ##
    if (!inherits(model, "formula") || length(model) != 3L) {
        stop("model must be a formula of the form \"resp ~ pred\"")
    }
    method <- match.arg(method)
    REML <- method == "REML"
    ## check if correlation is present and has groups
    groups <- if (!is.null(correlation)) getGroupsFormula(correlation) ## else NULL
    ## create a gls structure containing the plug-ins
    glsSt <-
        glsStruct(corStruct = correlation, varStruct = varFunc(weights))

    ## we need to resolve '.' in the formula here
    model <- terms(model, data=data)
    ## extract a data frame with enough information to evaluate
    ## formula, groups, corStruct, and varStruct
    mfArgs <- list(formula = asOneFormula(formula(glsSt), model, groups),
                   data = data, na.action = na.action)
    if (!missing(subset)) {
        mfArgs[["subset"]] <- asOneSidedFormula(Call[["subset"]])[[2L]]
    }
    mfArgs$drop.unused.levels <- TRUE
    dataMod <- do.call(model.frame, mfArgs)
    origOrder <- row.names(dataMod)	# preserve the original order
    if (!is.null(groups)) {
        ## sort the model.frame by groups and get the matrices and parameters
        ## used in the estimation procedures
        ## always use innermost level of grouping
        groups <- eval(substitute(~ 1 | GR, list(GR = groups[[2L]])))
        grps <- getGroups(dataMod, groups,
                          level = length(getGroupsFormula(groups, asList = TRUE)))
        ## ordering data by groups
        ord <- order(grps)
        grps <- grps[ord]
        dataMod <- dataMod[ord, ,drop = FALSE]
        revOrder <- match(origOrder, row.names(dataMod)) # putting in orig. order
    } else grps <- NULL

    ## obtaining basic model matrices
    X <- model.frame(model, dataMod)
    Terms <- attr(X, "terms")
    if (length(attr(Terms, "offset")))
        stop("offset() terms are not supported")
    ## keeping the contrasts for later use in predict
    contr <- lapply(X, function(el)
                    if (inherits(el, "factor")) contrasts(el))
    contr <- contr[!unlist(lapply(contr, is.null))]
    X <- model.matrix(model, X)
    if(ncol(X) == 0L) stop("no coefficients to fit")
    y <- eval(model[[2L]], dataMod)
    N <- nrow(X)
    p <- ncol(X)				# number of coefficients
    parAssign <- attr(X, "assign")
    namTerms <- attr(Terms, "term.labels")
    if (attr(Terms, "intercept") > 0) {
        namTerms <- c("(Intercept)", namTerms)
    }
    namTerms <- factor(parAssign, labels = namTerms)
    parAssign <- split(order(parAssign), namTerms)
    fixedSigma <- (controlvals$sigma > 0) ## 17-11-2015; Fixed sigma patch
    ## creating the condensed linear model
    attr(glsSt, "conLin") <-
        list(Xy = array(c(X, y), c(N, ncol(X) + 1L), list(row.names(dataMod),
	     c(colnames(X), deparse(model[[2]])))),
             dims = list(N = N, p = p, REML = as.integer(REML)), logLik = 0,
             ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions
             sigma = controlvals$sigma, fixedSigma = fixedSigma)

    ## initialization
    glsEstControl <- controlvals["singular.ok"]
    glsSt <- Initialize(glsSt, dataMod, glsEstControl)
    parMap <- attr(glsSt, "pmap")

    ##
    ## getting the fitted object, possibly iterating for variance functions
    ##
    numIter <- numIter0 <- 0L
    repeat {
        oldPars <- c(attr(glsSt, "glsFit")[["beta"]], coef(glsSt))
        if (length(coef(glsSt))) {		# needs ms()
            optRes <- if (controlvals$opt == "nlminb") {
                nlminb(c(coef(glsSt)),
                       function(glsPars) -logLik(glsSt, glsPars),
                       control = list(trace = controlvals$msVerbose,
                       iter.max = controlvals$msMaxIter))
            } else {
                optim(c(coef(glsSt)),
                      function(glsPars) -logLik(glsSt, glsPars),
                      method = controlvals$optimMethod,
                      control = list(trace = controlvals$msVerbose,
                      maxit = controlvals$msMaxIter,
                      reltol = if(numIter == 0L) controlvals$msTol
                      else 100*.Machine$double.eps))
            }
            coef(glsSt) <- optRes$par
        } else {
            optRes <- list(convergence = 0)
        }
        attr(glsSt, "glsFit") <- glsEstimate(glsSt, control = glsEstControl)
        ## checking if any updating is needed
        if (!needUpdate(glsSt)) {
            if (optRes$convergence)
                stop(optRes$message)
            break
        }
        ## updating the fit information
        numIter <- numIter + 1L
        glsSt <- update(glsSt, dataMod)
        ## calculating the convergence criterion
        aConv <- c(attr(glsSt, "glsFit")[["beta"]], coef(glsSt))
        conv <- abs((oldPars - aConv)/ifelse(aConv == 0, 1, aConv))
        aConv <- c("beta" = max(conv[1:p]))
        conv <- conv[-(1:p)]
        for(i in names(glsSt)) {
            if (any(parMap[,i])) {
                aConv <- c(aConv, max(conv[parMap[,i]]))
                names(aConv)[length(aConv)] <- i
            }
        }
        if (verbose) {
            cat("\nIteration:",numIter)
            cat("\nObjective:", format(optRes$value), "\n")
            print(glsSt)
            cat("\nConvergence:\n")
            print(aConv)
        }
        if (max(aConv) <= controlvals$tolerance) {
            break
        }
        if (numIter > controlvals$maxIter) {
            stop("maximum number of iterations reached without convergence")
        }
    }
    ## wrapping up
    glsFit <- attr(glsSt, "glsFit")
    namBeta <- names(glsFit$beta)
    ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions
    attr(glsSt, "fixedSigma") <- fixedSigma
    attr(parAssign, "varBetaFact") <- varBeta <-
        glsFit$sigma * glsFit$varBeta * sqrt((N - REML * p)/(N - p))
    varBeta <- crossprod(varBeta)
    dimnames(varBeta) <- list(namBeta, namBeta)
    ##
    ## fitted.values and residuals (in original order)
    ##
    Fitted <- fitted(glsSt)
    ## putting groups back in original order, if present
    if (!is.null(grps)) {
        grps <- grps[revOrder]
        Fitted <- Fitted[revOrder]
        Resid <- y[revOrder] - Fitted
        attr(Resid, "std") <- glsFit$sigma/varWeights(glsSt)[revOrder]
    } else {
        Resid <- y - Fitted
        attr(Resid, "std") <- glsFit$sigma/varWeights(glsSt)
    }
    names(Resid) <- names(Fitted) <- origOrder

    ## getting the approximate var-cov of the parameters
    apVar <-
        if (controlvals$apVar)
            glsApVar(glsSt, glsFit$sigma,
                     .relStep = controlvals[[".relStep"]],
                     minAbsPar = controlvals[["minAbsParApVar"]],
                     natural = controlvals[["natural"]])
        else
            "Approximate variance-covariance matrix not available"
    ## getting rid of condensed linear model and fit
    dims <- attr(glsSt, "conLin")[["dims"]]
    dims[["p"]] <- p
    attr(glsSt, "conLin") <- NULL
    attr(glsSt, "glsFit") <- NULL
    attr(glsSt, "fixedSigma") <- fixedSigma ## 17-11-2015; Fixed sigma patch; ..
    grpDta <- inherits(data, "groupedData")
    ##
    ## creating the  gls object
    ##
    structure(class = "gls",
              list(modelStruct = glsSt,
                   dims = dims,
                   contrasts = contr,
                   coefficients = glsFit[["beta"]],
                   varBeta = varBeta,
	## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions
                   sigma = if(fixedSigma) controlvals$sigma else glsFit$sigma,
                   apVar = apVar,
                   logLik = glsFit$logLik,
                   numIter = if (needUpdate(glsSt)) numIter else numIter0,
                   groups = grps,
                   call = Call,
                   terms = Terms,
                   method = method,
                   fitted = Fitted,
                   residuals = Resid,
                   parAssign = parAssign,
		   na.action = attr(dataMod, "na.action")),
	      namBetaFull = colnames(X),
	      ## saving labels and units for plots
	      units = if(grpDta) attr(data, "units"),
	      labels= if(grpDta) attr(data, "labels"))
}

### Auxiliary functions used internally in gls and its methods

glsApVar.fullGlsLogLik <- function(Pars, object, conLin, dims, N)
{
    fixedSigma <- attr(object, "fixedSigma")
    ## logLik as a function of sigma and coef(glsSt)
    npar <- length(Pars)
    if (!fixedSigma) {
        lsigma <- Pars[npar]              # within-group std. dev.
        Pars <- Pars[-npar]
        sigma <- 0
    } else {
        sigma <- conLin$sigma
    }
    coef(object) <- Pars
    conLin <- recalc(object, conLin)
    val <- .C(gls_loglik,
              as.double(conLin$Xy),
              as.integer(unlist(dims)),
              logLik = double(1L),
              lRSS = double(1L), sigma = as.double(sigma), NAOK = TRUE)[c("logLik", "lRSS")]
    if (!fixedSigma) {
        aux <- 2 * (val[["lRSS"]] - lsigma)
        conLin[["logLik"]] + val[["logLik"]] + (N * aux - exp(aux))/2
    } else {
        val[["logLik"]]
    }
}

glsApVar <-
    function(glsSt, sigma, conLin = attr(glsSt, "conLin"),
             .relStep = .Machine$double.eps^(1/3), minAbsPar = 0,
             natural = TRUE)
{
    fixedSigma <- attr(glsSt, "fixedSigma")
    ## calculate approximate variance-covariance matrix of all parameters
    ## except the coefficients
    if (length(glsCoef <- coef(glsSt)) > 0L) {
        cSt <- glsSt[["corStruct"]]
        if (natural && !is.null(cSt) && inherits(cSt, "corSymm")) {
            cStNatPar <- coef(cSt, unconstrained = FALSE)
            class(cSt) <- c("corNatural", "corStruct")
            coef(cSt) <- log((1 + cStNatPar)/(1 - cStNatPar))
            glsSt[["corStruct"]] <- cSt
            glsCoef <- coef(glsSt)
        }
        dims <- conLin$dims
        N <- dims$N - dims$REML * dims$p
        conLin[["logLik"]] <- 0               # making sure
        Pars <- if(fixedSigma) glsCoef else c(glsCoef, lSigma = log(sigma))
        val <- fdHess(Pars, glsApVar.fullGlsLogLik, glsSt, conLin, dims, N,
                      .relStep = .relStep, minAbsPar = minAbsPar)[["Hessian"]]
        if (all(eigen(val, only.values=TRUE)$values < 0)) {
            ## negative definite - OK
            val <- solve(-val)
            nP <- names(Pars)
            dimnames(val) <- list(nP, nP)
            attr(val, "Pars") <- Pars
            attr(val, "natural") <- natural
            val
        } else {
            ## problem - solution is not a maximum
            "Non-positive definite approximate variance-covariance"
        }
    } # else NULL
}

glsEstimate <-
    function(object, conLin = attr(object, "conLin"),
             control = list(singular.ok = FALSE))
{
    dd <- conLin$dims
    p <- dd$p
    oXy <- conLin$Xy
    fixSig <- conLin$fixedSigma ## 17-11-2015; Fixed sigma patch; ..
    sigma <- conLin$sigma
    conLin <- recalc(object, conLin)	# updating for corStruct and varFunc
    val <- .C(gls_estimate,
              as.double(conLin$Xy),
              as.integer(unlist(dd)),
              beta = double(p),
              sigma = as.double(sigma), ## 17-11-2015; Fixed sigma patch; ..
              logLik = double(1L),
              varBeta = double(p * p),
              rank = integer(1),
              pivot = as.integer(1:(p + 1L)),
              NAOK = TRUE)[c("beta","sigma","logLik","varBeta", "rank", "pivot")]
    rnk <- val[["rank"]]
    rnkm1 <- rnk - 1
    if (!control$singular.ok && rnkm1 < p) {
        stop(gettextf("computed \"gls\" fit is singular, rank %s", rnk),
             domain = NA)
    }
    N <- dd$N - dd$REML * p
    namCoef <- colnames(oXy)[val[["pivot"]][1:rnkm1] + 1L]	# coef names
    varBeta <- t(array(val[["varBeta"]], c(rnkm1, rnkm1),
                       list(namCoef, namCoef)))
    beta <- val[["beta"]][1:rnkm1]
    names(beta) <- namCoef
    fitted <- c(oXy[, namCoef, drop = FALSE] %*% beta)
    resid <- oXy[, p + 1] - fitted
    ll <- conLin$logLik + val[["logLik"]]
    logLik <-
        if (!fixSig) {
            (N * (logb(N) - (1 + logb(2 * pi))))/2 + ll
            ## formula 2.21 on page 70  if sigma is estimated ML formula or 2.23 page 76 with REML
        } else {
            (-N/2) * logb(2*pi)  + ll
        }
    list(logLik = logLik, beta = beta,
         sigma = val[["sigma"]], varBeta = varBeta, fitted = fitted,
         resid = resid, auxSigma = sqrt(sum((resid)^2))/sqrt(N))
}

### Methods for standard generics

ACF.gls <-
    function(object, maxLag, resType = c("pearson", "response", "normalized"),
             form = ~1, na.action = na.fail, ...)
{
    resType <- match.arg(resType)
    res <- resid(object, type = resType)
    wchRows <- NULL
    if (is.null(grps <- getGroups(object))) {
        ## check if formula defines groups
        if (!is.null(grpsF <- getGroupsFormula(form))) {
            if (is.null(data <- getData(object))) {
                ## will try to construct
                allV <- all.vars(grpsF)
                if (length(allV) > 0L) {
                    alist <- lapply(as.list(allV), as.name)
                    names(alist) <- allV
                    alist <- c(as.list(quote(data.frame)), alist)
                    mode(alist) <- "call"
                    data <- eval(alist, sys.parent(1))
                }
            }
            grps <- model.frame(grpsF, data, na.action = na.action)
            wchRows <- !is.na(match(row.names(data), row.names(grps)))
            grps <- getGroups(grps, grpsF)
        }
    }
    if (!is.null(grps)) {
        if (!is.null(wchRows)) {
            res <- res[wchRows]
        }
        res <- split(res, grps)
    } else {
        res <- list(res)
    }
    if(missing(maxLag)) {
        maxLag <- min(c(maxL <- max(lengths(res)) - 1L,
                        as.integer(10 * log10(maxL + 1))))
    }
    val <- lapply(res,
                  function(el, maxLag) {
                      N <- maxLag + 1L
                      tt <- double(N)
                      nn <- integer(N)
                      N <- min(c(N, n <- length(el)))
                      nn[1:N] <- n + 1L - 1:N
                      ## el <- el - mean(el)
                      for(i in 1:N) {
                          el1 <- el[1:(n-i+1L)]
                          el2 <- el[i:n]
                          tt[i] <- sum(el1 * el2)
                      }
                      array(c(tt,nn), c(length(tt), 2L))
                  }, maxLag = maxLag)
    val0 <- apply(sapply(val, function(x) x[,2L]), 1, sum)
    val1 <- apply(sapply(val, function(x) x[,1L]), 1, sum)/val0
    val2 <- val1/val1[1L]
    structure(data.frame(lag = 0:maxLag, ACF = val2),
	      n.used = val0,
	      class = c("ACF", "data.frame"))
}

anova.gls <-
    function(object, ..., test = TRUE, type = c("sequential", "marginal"),
             adjustSigma = NA, Terms, L, verbose = FALSE)
{
    fixSig <- attr(object$modelStruct, "fixedSigma")
    fixSig <- !is.null(fixSig) && fixSig
    Lmiss <- missing(L)
    ## returns the likelihood ratio statistics, the AIC, and the BIC
    dots <- list(...)
    if ((rt <- length(dots) + 1L) == 1L) {    ## just one object
        if (!inherits(object,"gls"))
            stop("object must inherit from class \"gls\"")
	if(is.na(adjustSigma))
	    ## REML correction already applied to gnls objects
	    adjustSigma <- inherits(object, "gnls")
        dims <- object$dims
        N <- dims$N
        p <- dims$p
        REML <- dims$REML
        assign <- object$parAssign
        vBeta <- attr(assign, "varBetaFact")
	if (!REML && adjustSigma)
	    ## using REML-like estimate of sigma under ML
	    vBeta <- sqrt((N - p)/N) * vBeta
        c0 <- solve(t(vBeta), coef(object))
        nTerms <- length(assign)
        dDF <- N - p
        lab <- paste("Denom. DF:", dDF,"\n")
        if (missing(Terms) && Lmiss) {
            ## returns the F.table (Wald) for the fixed effects
            type <- match.arg(type)
            Fval <- Pval <- double(nTerms)
            nDF <- integer(nTerms)
            for(i in 1:nTerms) {
                nDF[i] <- length(assign[[i]])
		c0i <-
		    if (type == "sequential") # type I SS
			c0[assign[[i]]]
		    else ## "marginal"
			qr.qty(qr(vBeta[, assign[[i]], drop = FALSE]), c0)[1:nDF[i]]
                Fval[i] <- sum(c0i^2)/nDF[i]
                Pval[i] <- pf(Fval[i], nDF[i], dDF, lower.tail=FALSE)
            }
            ##
            ## fixed effects F-values, df, and p-values
            ##
	    aod <- data.frame(numDF = nDF, "F-value" = Fval, "p-value" = Pval,
			check.names=FALSE)
            rownames(aod) <- names(assign)
        } else {
            if (Lmiss) {                 # terms is given
                if (is.numeric(Terms) && all(Terms == as.integer(Terms))) {
                    if (min(Terms) < 1 || max(Terms) > nTerms) {
                        stop(gettextf("'Terms' must be between 1 and %d",
                                      nTerms), domain = NA)
                    }
                } else {
                    if (is.character(Terms)) {
                        if (any(noMatch <- is.na(match(Terms, names(assign))))) {
                            stop(sprintf(ngettext(sum(noMatch),
                                                  "term %s not matched",
                                                  "terms %s not matched"),
                                         paste(Terms[noMatch], collapse = ", ")),
                                 domain = NA)
                        }
                    } else {
                        stop("terms can only be integers or characters")
                    }
                }
                lab <-
                    paste(lab, "F-test for:",
                          paste(names(assign[Terms]),collapse=", "),"\n")
                L <- diag(p)[unlist(assign[Terms]),,drop=FALSE]
            } else {
                L <- as.matrix(L)
                if (ncol(L) == 1L) L <- t(L)     # single linear combination
                nrowL <- nrow(L)
                ncolL <- ncol(L)
                if (ncol(L) > p) {
                    stop(sprintf(ngettext(ncol(L),
                                          "'L' must have at most %d column",
                                          "'L' must have at most %d columns"),
                                 ncol(L)), domain = NA)
                }
                dmsL1 <- rownames(L)
                L0 <- array(0, c(nrowL, p), list(NULL, names(coef(object))))
                if (is.null(dmsL2 <- colnames(L))) {
                    ## assume same order as effects
                    L0[, 1:ncolL] <- L
                } else {
                    if (any(noMatch <- is.na(match(dmsL2, colnames(L0))))) {
                        stop(sprintf(ngettext(sum(noMatch),
                                              "effect %s not matched",
                                              "effects %s not matched"),
                                     paste(dmsL2[noMatch],collapse=", ")),
                             domain = NA)
                    }
                    L0[, dmsL2] <- L
                }
                L <- L0[noZeroRowL <- as.logical((L0 != 0) %*% rep(1, p)), , drop = FALSE]
                nrowL <- nrow(L)
                noZeroColL <- as.logical(c(rep(1,nrowL) %*% (L != 0)))
                rownames(L) <- if(is.null(dmsL1)) 1:nrowL else dmsL1[noZeroRowL]
                lab <- paste(lab, "F-test for linear combination(s)\n")
            }
            nDF <- sum(svd.d(L) > 0)
            c0 <- c(qr.qty(qr(vBeta %*% t(L)), c0))[1:nDF]
            Fval <- sum(c0^2)/nDF
            Pval <- pf(Fval, nDF, dDF, lower.tail=FALSE)
            aod <- data.frame(numDF = nDF, "F-value" = Fval, "p-value" = Pval,
                              check.names=FALSE)
            if (!Lmiss) {
                attr(aod, "L") <-
                    if(nrow(L) > 1) L[, noZeroColL, drop = FALSE] else L[, noZeroColL]
            }
        }
        attr(aod, "label") <- lab
        attr(aod,"rt") <- rt
        class(aod) <- c("anova.lme", "data.frame")
        aod
    }
    ##
    ## Otherwise construct the likelihood ratio and information table
    ## objects in ... may inherit from gls, lm, lmList, and lme (for now)
    ##
    else {
        Call <- match.call()
        Call[[1]] <- quote(nlme::anova.lme)
        eval.parent(Call)
    }
}

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

coef.gls <-
    function(object, allCoef = FALSE, ...)
{
    val <- object$coefficients
    if (allCoef) {
        namFull <- attr(object, "namBetaFull")
        if (length(val) != (lF <- length(namFull))) {
            aux <- rep(NA, lF)
            names(aux) <- namFull
            aux[names(val)] <- val
            val <- aux
        }
    }
    val
}

comparePred.gls <-
    function(object1, object2, primary = NULL,
             minimum = min(primary), maximum = max(primary),
             length.out = 51L, level = NULL, ...)
{
    if (length(level) > 1L) {
        stop("only one level allowed for predictions")
    }
    args <- list(object = object1,
                 primary = primary,
                 level = level,
                 length.out = length.out)
    if (!is.null(primary)) {
	  ## need to do this before forcing the evaluations
	  primary <- eval(asOneSidedFormula(primary)[[2]],
                         eval(object1$call$data))
        args[["minimum"]] <- minimum
        args[["maximum"]] <- maximum
    }
    val1 <- do.call(augPred, args)
    dm1 <- dim(val1)
    c1 <- deparse(substitute(object1))
    levels(val1[,4])[1] <- c1
    args[["object"]] <- object2
    val2 <- do.call(augPred, args)
    dm2 <- dim(val2)
    c2 <- deparse(substitute(object2))
    levels(val2[, 4L])[1] <- c2
    val2 <- val2[val2[, 4L] != "original", , drop = FALSE]
    names(val2) <- names(val1)

    if (dm1[1L] == dm2[1L]) {
        lv1 <- sort(levels(val1[, 2L]))
        lv2 <- sort(levels(val2[, 2L]))
        if ((length(lv1) != length(lv2)) || any(lv1 != lv2)) {
            stop(gettextf("%s and %s must have the same group levels", c1, c2),
                 domain = NA)
        }
        val <- rbind(val1[, -4L], val2[, -4L])
        val[, ".type"] <-
            ordered(c(as.character(val1[,4L]), as.character(val2[,4L])),
                    levels = c(c1, c2, "original"))
        attr(val, "formula") <- attr(val1, "formula")
    } else {				# one may have just "fixed"
        if (dm1[1L] > dm2[1L]) {
            mult <- dm1[1L] %/% dm2[1L]
            if ((length(levels(val2[, 2])) != 1L) ||
                (length(levels(val1[, 2])) != mult))
            {
                stop("wrong group levels")
            }
            val <-
                data.frame(c(val1[,1L], rep(val2[,1L], mult)), rep(val1[,1L], 2L),
                           c(val1[,3L], rep(val2[,3L], mult)),
                           ordered(c(as.character(val1[,4L]),
                                     rep(as.character(val2[,4L]), mult)),
                                   levels = c(c1, c2, "original")))
            attr(val, "formula") <- attr(val1, "formula")
        } else {
            mult <- dm2[1L] %/% dm1[1L]
            if ((length(levels(val1[, 2])) != 1L) ||
                (length(levels(val2[, 2])) != mult))
            {
                stop("wrong group levels")
            }
            val <-
                data.frame(c(rep(val1[,1L], mult), val2[,1L]), rep(val2[,1L], 2L),
                           c(rep(val1[,3L], mult), val2[,3L]),
                           ordered(c(rep(as.character(val1[,4L]), mult),
                                     as.character(val1[,4L])),
                                   levels = c(c1, c2, "original")))
            attr(val, "formula") <- attr(val2, "formula")
        }
    }
    class(val) <- c("comparePred", "augPred", class(val))
    attr(val, "labels") <- attr(val1, "labels")
    attr(val, "units") <- attr(val1, "units")
    val
}

fitted.gls <-
    function(object, ...)
{
    val <- napredict(object$na.action, object$fitted)
    lab <- "Fitted values"
    if (!is.null(aux <- attr(object, "units")$y)) {
        lab <- paste(lab, aux)
    }
    attr(val, "label") <- lab
    val
}


formula.gls <- function(x, ...)
{
    if (is.null(x$terms)) { # gls objects from nlme <= 3.1-155
        eval(x$call$model)
    } else formula(x$terms)
}

getGroups.gls <- function(object, form, level, data, sep) object$groups

getGroupsFormula.gls <- function(object, asList = FALSE, sep)
{
    if (!is.null(cSt <- object$modelStruct$corStruct))
        getGroupsFormula(cSt, asList)
    ## else NULL
}

getResponse.gls <-
    function(object, form)
{
    val <- resid(object) + fitted(object)
    if (is.null(lab <- attr(object, "labels")$y)) {
        lab <- deparse(getResponseFormula(object)[[2]])
    }
    if (!is.null(aux <- attr(object, "units")$y)) {
        lab <- paste(lab, aux)
    }
    attr(val, "label") <- lab
    val
}

intervals.gls <-
    function(object, level = 0.95, which = c("all", "var-cov", "coef"), ...)
{
    which <- match.arg(which)
    val <- list()
    dims <- object$dims
    ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions
    fixSig <- attr(object$modelStruct, "fixedSigma")
    fixSig <- !is.null(fixSig) && fixSig

    if (which != "var-cov") {		# coefficients included
        len <- -qt((1-level)/2, dims$N - dims$p) * sqrt(diag(object$varBeta))
        est <- coef(object)
        val[["coef"]] <-
            array(c(est - len, est, est + len), c(length(est), 3L),
                  list(names(est), c("lower", "est.", "upper")))
        attr(val[["coef"]], "label") <- "Coefficients:"
    }

    if (which != "coef") {		# variance-covariance included
        if (is.null(aV <- object$apVar)) {	# only sigma
            Nr <- if (inherits(object, "gnls")) {   #always REML-like sigma
                      dims$N - dims$p
                  } else {
                      dims$N - dims$REML * dims$p
                  }
            ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions
            if(!fixSig) {
                est <- object$sigma * sqrt(Nr)
                val[["sigma"]] <-
                    structure(c(lower = est/sqrt(qchisq((1+level)/2, Nr)),
                                "est."= object$sigma,
                                upper = est/sqrt(qchisq((1-level)/2, Nr))),
                              label = "Residual standard error:")
            } else {
                est <- 1
		val[["sigma"]] <- structure(setNames(rep(object$sigma, 3),
                                                     c("lower", "est.", "upper")),
                                            label = "Fixed Residual standard error:")
            }
        } else {
            if (is.character(aV)) {
                stop(gettextf("cannot get confidence intervals on var-cov components: %s",
                              aV), domain = NA)
            }
            len <- -qnorm((1-level)/2) * sqrt(diag(aV))
            est <- attr(aV, "Pars")
            nP <- length(est)
            glsSt <- object[["modelStruct"]]
            pmap <- attr(glsSt, "pmap")
            if (!all(whichKeep <- apply(pmap, 2, any))) {
                ## need to deleted components with fixed coefficients
                aux <- glsSt[whichKeep]
                class(aux) <- class(glsSt)
                attr(aux, "settings") <- attr(glsSt, "settings")
                attr(aux, "pmap") <- pmap[, whichKeep, drop = FALSE]
                glsSt <- aux
            }
            cSt <- glsSt[["corStruct"]]
            if (!is.null(cSt) && inherits(cSt, "corSymm") && attr(aV, "natural")) {
                ## converting to corNatural
                class(cSt) <- c("corNatural", "corStruct")
                glsSt[["corStruct"]] <- cSt
            }
            namG <- names(glsSt)
            auxVal <- vector("list", length(namG) + 1L)
            names(auxVal) <- c(namG, "sigma")
            aux <-
                array(c(est - len, est, est + len),
                      c(nP, 3), list(NULL, c("lower", "est.", "upper")))
	    auxVal[["sigma"]] <-
		structure(if(!fixSig) { # last param. = log(sigma)
			      s <- exp(aux[nP, ])
			      aux <- aux[-nP,, drop = FALSE]
			      s
			  } else
			      c(object$sigma, object$sigma, object$sigma),
			  label = "Residual standard error:")
            rownames(aux) <- names(coef(glsSt, FALSE))
            for(i in 1:3) {
                coef(glsSt) <- aux[,i]
                aux[,i] <- coef(glsSt, unconstrained = FALSE)
            }
            for(i in namG) {
                au.i <- aux[pmap[,i], , drop = FALSE]
                dimnames(au.i)[[1]] <- substring(dimnames(au.i)[[1]], nchar(i, "c") + 2L)
                attr(au.i, "label") <-
                    switch(i,
                           corStruct = "Correlation structure:",
                           varStruct = "Variance function:",
                           paste0(i, ":"))
                auxVal[[i]] <- au.i
            }
            val <- c(val, auxVal)
        }
    }
    structure(val, level = level, class = "intervals.gls")
}

logLik.gls <-
    function(object, REML, ...)
{
    ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions
    fixSig <- attr(object[["modelStruct"]], "fixedSigma")
    fixSig <- !is.null(fixSig) && fixSig
    p <- object$dims$p
    N <- object$dims$N
    Np <- N - p
    estM <- object$method
    if (missing(REML)) REML <- estM == "REML"
    val <- object[["logLik"]]
    if (REML && estM == "ML") { # have to correct logLik
        val <- val + (p * (log(2 * pi) + 1) + Np * log(1 - p/N) +
                      sum(log(abs(svd.d(object$varBeta))))) / 2
    }
    else if (!REML && (estM == "REML")) { # have to correct logLik
        val <- val - (p * (log(2*pi) + 1) + N * log(1 - p/N) +
                      sum(log(abs(svd.d(object$varBeta))))) / 2
    }
    structure(val,
              nall = N,
              nobs = N - REML * p,
              ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions
              df = p + length(coef(object[["modelStruct"]])) + as.integer(!fixSig),
              class = "logLik")
}

nobs.gls <- function(object, ...) object$dims$N


plot.gls <-
    function(x, form = resid(., type = "pearson") ~ fitted(.), abline,
             id = NULL, idLabels = NULL, idResType = c("pearson", "normalized"),
             grid, ...)
    ## Diagnostic plots based on residuals and/or fitted values
{
    plot.lme(x, form=form, abline=abline, id=id, idLabels=idLabels, idResType=idResType, grid=grid, ...)
}

predict.gls <-
    function(object, newdata, na.action = na.fail, ...)
{
    ##
    ## method for predict() designed for objects inheriting from class gls
    ##
    if (missing(newdata)) {		# will return fitted values
        return(fitted(object))
    }
    form <- delete.response(object[["terms"]])
    ## making sure factor levels are the same as in contrasts
    ## and supporting character-type 'newdata' for factors (all via xlev)
    contr <- object$contrasts           # these are in matrix form
    dataMod <- model.frame(formula = form, data = newdata,
                           na.action = na.action, drop.unused.levels = TRUE,
                           xlev = lapply(contr, rownames))
    N <- nrow(dataMod)
    if (length(all.vars(form)) > 0) {
        X <- model.matrix(form, dataMod, contr)
    } else {
        X <- array(1, c(N, 1), list(row.names(dataMod), "(Intercept)"))
    }
    cf <- coef(object)
    val <- c(X[, names(cf), drop = FALSE] %*% cf)
    lab <- "Predicted values"
    if (!is.null(aux <- attr(object, "units")$y)) {
        lab <- paste(lab, aux)
    }
    structure(val, label = lab)
}

print.intervals.gls <-
    function(x, ...)
{
    cat(paste0("Approximate ", attr(x, "level") * 100,
               "% confidence intervals\n"))
    for(i in names(x)) {
        aux <- x[[i]]
        cat("\n ",attr(aux, "label"), "\n", sep = "")
        attr(aux, "label") <- NULL
        print(if(i == "sigma") c(aux) else as.matrix(aux), ...)
    }
    invisible(x)
}

print.gls <-
    ## method for print() used for gls objects
    function(x, ...)
{
    ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions
    fixSig <- attr(x[["modelStruct"]], "fixedSigma")
    fixSig <- !is.null(fixSig) && fixSig
    dd <- x$dims
    mCall <- x$call
    if (inherits(x, "gnls")) {
        cat("Generalized nonlinear least squares fit\n")
    } else {
        cat("Generalized least squares fit by ")
        cat(if(x$method == "REML") "REML\n" else "maximum likelihood\n")
    }
    cat("  Model:", deparse(mCall$model), "\n")
    cat("  Data:", deparse( mCall$data ), "\n")
    if (!is.null(mCall$subset)) {
        cat("  Subset:", deparse(asOneSidedFormula(mCall$subset)[[2]]),"\n")
    }
    if (inherits(x, "gnls")) {
        cat("  Log-likelihood: ", format(x$logLik), "\n", sep = "")
    } else {
        cat("  Log-", if(x$method == "REML") "restricted-" else "",
            "likelihood: ", format(x$logLik), "\n", sep = "")
    }
    cat("\nCoefficients:\n")
    print(coef(x), ...)
    cat("\n")
    if (length(x$modelStruct) > 0L) {
        print(summary(x$modelStruct), ...)
    }
    cat("Degrees of freedom:", dd[["N"]],"total;",dd[["N"]] - dd[["p"]],
        "residual\n")
    cat("Residual standard error:", format(x$sigma),"\n")
    invisible(x)
}

print.summary.gls <-
    function(x, verbose = FALSE, digits = .Options$digits, ...)
{
    dd <- x$dims
    fixSig <- attr(x[["modelStruct"]], "fixedSigma")
    fixSig <- !is.null(fixSig) && fixSig
    verbose <- verbose || attr(x, "verbose")
    mCall <- x$call
    if (inherits(x, "gnls")) {
        cat("Generalized nonlinear least squares fit\n")
    } else {
        cat("Generalized least squares fit by ")
        cat(if(x$method == "REML") "REML\n" else "maximum likelihood\n")
    }
    cat("  Model:", deparse(mCall$model), "\n")
    cat("  Data:", deparse( mCall$data ), "\n")
    if (!is.null(mCall$subset)) {
        cat("  Subset:", deparse(asOneSidedFormula(mCall$subset)[[2]]),"\n")
    }
    print(data.frame(AIC=x$AIC, BIC=x$BIC, logLik=as.vector(x$logLik),
		     row.names = " "),
	  ...)
    if (verbose) { cat("Convergence at iteration:",x$numIter,"\n") }
    if (length(x$modelStruct)) {
        cat("\n")
        print(summary(x$modelStruct), ...)
    }
    cat("\nCoefficients:\n")
    xtTab <- as.data.frame(x$tTable)
    wchPval <- match("p-value", names(xtTab))
    for(i in names(xtTab)[-wchPval]) {
        xtTab[, i] <- format(zapsmall(xtTab[, i]))
    }
    xtTab[,wchPval] <- format(round(xtTab[,wchPval], 4L))
    if (any(wchLv <- (as.double(levels(xtTab[, wchPval])) == 0))) {
        levels(xtTab[, wchPval])[wchLv] <- "<.0001"
    }
    row.names(xtTab) <- dimnames(x$tTable)[[1]]
    print(xtTab, ...)
    if (nrow(x$tTable) > 1L) {
        corr <- x$corBeta
        class(corr) <- "correlation"
        print(corr, title = "\n Correlation:", ...)
    }
    cat("\nStandardized residuals:\n")
    print(x$residuals, ...)
    cat("\n")
    cat("Residual standard error:", format(x$sigma),"\n")
    cat("Degrees of freedom:", dd[["N"]],"total;", dd[["N"]] - dd[["p"]],
        "residual\n")
    invisible(x)
}

residuals.gls <-
    function(object, type = c("response", "pearson", "normalized"), ...)
{
    type <- match.arg(type)
    val <- object$residuals
    if (type != "response") {
        val <- val/attr(val, "std")
        lab <- "Standardized residuals"
        if (type == "normalized") {
            if (!is.null(cSt <- object$modelStruct$corStruct)) {
                ## normalize according to inv-trans factor
                val <- recalc(cSt, list(Xy = as.matrix(val)))$Xy[, 1]
                lab <- "Normalized residuals"
            }
        }
    } else {
        lab <- "Residuals"
        if (!is.null(aux <- attr(object, "units")$y))
            lab <- paste(lab, aux)
    }
    attr(val, "label") <- lab
    if (!is.null(object$na.action)) {
        res <- naresid(object$na.action, val)
        attr(res, "std") <- naresid(object$na.action, attr(val, "std"))
        attr(res, "label") <- attr(val, "label")
        res
    } else val
}

summary.gls <- function(object, verbose = FALSE, ...) {
    ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions
    fixSig <- attr(object[["modelStruct"]], "fixedSigma")
    fixSig <- !is.null(fixSig) && fixSig
    ##
    ## generates an object used in the print.summary method for lme
    ##
    ##  variance-covariance estimates for coefficients
    ##
    stdBeta <- sqrt(diag(as.matrix(object$varBeta)))
    corBeta <- t(object$varBeta/stdBeta)/stdBeta
    ##
    ## coefficients, std. deviations and z-ratios
    ##
    beta <- coef(object)
    dims <- object$dims
    dimnames(corBeta) <- list(names(beta),names(beta))
    object$corBeta <- corBeta
    tTable <- data.frame(beta, stdBeta, beta/stdBeta, beta)
    dimnames(tTable)<-
        list(names(beta),c("Value","Std.Error","t-value","p-value"))
    tTable[, "p-value"] <- 2 * pt(-abs(tTable[,"t-value"]), dims$N - dims$p)
    object$tTable <- as.matrix(tTable)
    ##
    ## residuals
    ##
    resd <- resid(object, type = "pearson")
    if (length(resd) > 5) {
        resd <- quantile(resd, na.rm = TRUE)
        names(resd) <- c("Min","Q1","Med","Q3","Max")
    }
    object$residuals <- resd
    ##
    ## generating the final object
    ##
    aux <- logLik(object)
    structure(c(object, list(BIC = BIC(aux), AIC = AIC(aux))),
	      verbose = verbose,
	      class = c("summary.gls", class(object)))
}

update.gls <-
    function (object, model., ..., evaluate = TRUE)
{
    call <- object$call
    if (is.null(call))
	stop("need an object with call component")
    extras <- match.call(expand.dots = FALSE)$...
    if (!missing(model.))
	call$model <- update.formula(formula(object), model.)
    if(length(extras) > 0L) {
        ## the next two lines allow partial matching of argument names
        ## in the update.  This is nonstandard but required for examples
        ## in chapter 5 of Pinheiro and Bates (2000).
        glsa <- names(as.list(args(gls)))
        names(extras) <- glsa[pmatch(names(extras), glsa[-length(glsa)])]
	existing <- !is.na(match(names(extras), names(call)))
	## do these individually to allow NULL to remove entries.
	for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
	if(any(!existing)) {
	    call <- c(as.list(call), extras[!existing])
	    call <- as.call(call)
	}
    }
    if(evaluate) eval(call, parent.frame())
    else call
}

Variogram.gls <-
    function(object, distance, form = ~1,
             resType = c("pearson", "response", "normalized"),
             data, na.action = na.fail, maxDist, length.out = 50,
             collapse = c("quantiles", "fixed", "none"), nint = 20, breaks,
             robust = FALSE, metric = c("euclidean", "maximum", "manhattan"),
             ...)
{
    resType <- match.arg(resType)
    ## checking if object has a corSpatial element
    csT <- object$modelStruct$corStruct
    wchRows <- NULL
    if (missing(distance)) {
        if (missing(form) && inherits(csT, "corSpatial")) {
            distance <- getCovariate(csT)
            grps <- getGroups(object)
        } else {
            metric <- match.arg(metric)
            if (missing(data)) {
                data <- getData(object)
            }
            if (is.null(data)) {			# will try to construct
                allV <- all.vars(form)
                if (length(allV) > 0L) {
                    alist <- lapply(as.list(allV), as.name)
                    names(alist) <- allV
                    alist <- c(as.list(quote(data.frame)), alist)
                    mode(alist) <- "call"
                    data <- eval(alist, sys.parent(1))
                }
            }
            grpsF <- getGroupsFormula(form)
            grps <- NULL
            if (is.null(grpsF) || is.null(grps <- getGroups(data, grpsF))) {
                ## try to get from object
                grps <- getGroups(object)
            }
            covForm <- getCovariateFormula(form)
            covar <-
                if (length(all.vars(covForm)) > 0L) {
                    if (attr(terms(covForm), "intercept") == 1L) {
                        covForm <- eval(substitute(~ CV - 1, list(CV = covForm[[2]])))
                    }
                    covar <- model.frame(covForm, data, na.action = na.action)
                    ## making sure grps is consistent
                    wchRows <- !is.na(match(row.names(data), row.names(covar)))
                    if (!is.null(grps)) {
                        grps <- grps[wchRows, drop = TRUE]
                    }
                    as.data.frame(unclass(model.matrix(covForm, covar)))
                } else if (is.null(grps))
                    1:nrow(data)
                else
                    data.frame(dist = unlist(tapply(rep(1, nrow(data)), grps, cumsum)))

            distance <-
                if (is.null(grps))
                    dist(as.matrix(covar), method = metric)
                else {
                    covar <- split(covar, grps)
                    ## getting rid of 1-observation groups
                    covar <- covar[sapply(covar, function(el) nrow(as.matrix(el))) > 1]
                    lapply(covar,
                           function(el, metric) dist(as.matrix(el), method=metric),
                           metric = metric)
                }
        }
    }
    res <- resid(object, type = resType)
    if (!is.null(wchRows)) {
        res <- res[wchRows]
    }
    if (is.null(grps)) {
        val <- Variogram(res, distance)
    } else {
        res <- split(res, grps)
        res <- res[lengths(res) > 1L] # no 1-observation groups
        levGrps <- levels(grps)
        val <- structure(vector("list", length(levGrps)), names = levGrps)
        for(i in levGrps) {
            val[[i]] <- Variogram(res[[i]], distance[[i]])
        }
        val <- do.call(rbind, val)
    }
    if (!missing(maxDist)) {
        val <- val[val$dist <= maxDist, ]
    }
    collapse <- match.arg(collapse)
    if (collapse != "none") {             # will collapse values
        dst <- val$dist
        udist <- sort(unique(dst))
        ludist <- length(udist)
        if (!missing(breaks)) {
            if (min(breaks) > udist[1L]) {
                breaks <- c(udist[1L], breaks)
            }
            if (max(breaks) < udist[2L]) {
                breaks <- c(breaks, udist[2L])
            }
            if (!missing(nint) && nint != (length(breaks) - 1L)) {
                stop("'nint' is not consistent with 'breaks'")
            }
            nint <- length(breaks) - 1L
        }
        if (nint < ludist) {
            if (missing(breaks)) {
                if (collapse == "quantiles") {    # break into equal groups
                    breaks <- unique(quantile(dst, seq(0, 1, 1/nint)))
                } else {                          # fixed length intervals
                    breaks <- seq(udist[1L], udist[length(udist)], length = nint + 1)
                }
            }
            cutDist <- cut(dst, breaks)
        } else {
            cutDist <- dst
        }
        val <- lapply(split(val, cutDist),
                      function(el, robust) {
                          nh <- nrow(el)
                          vrg <- el$variog
                          if (robust) {
                              vrg <- ((mean(vrg^0.25))^4)/(0.457+0.494/nh)
                          } else {
                              vrg <- mean(vrg)
                          }
                          dst <- median(el$dist)
                          data.frame(variog = vrg, dist = dst)
                      }, robust = robust)
        val <- do.call(rbind, as.list(val))
        val$n.pairs <- unclass(table(na.omit(cutDist)))
    }
    row.names(val) <- 1:nrow(val)
    if (inherits(csT, "corSpatial") && resType != "normalized") {
        ## will keep model variogram
        sig2 <- if (resType == "pearson") 1 else object$sigma^2
        attr(val, "modelVariog") <-
            Variogram(csT, sig2 = sig2, length.out = length.out)
    }
    structure(val, collapse = collapse != "none",
	      class = c("Variogram", "data.frame"))
}

###*### glsStruct - a model structure for gls fits

glsStruct <-
    ## constructor for glsStruct objects
    function(corStruct = NULL, varStruct = NULL)
{
    val <- list(corStruct = corStruct, varStruct = varStruct)
    val <- val[!vapply(val, is.null, NA)] # removing NULL components
    class(val) <- c("glsStruct", "modelStruct")
    val
}

##*## glsStruct methods for standard generics

fitted.glsStruct <-
    function(object, glsFit = attr(object, "glsFit"), ...)
{
    glsFit[["fitted"]]
}

Initialize.glsStruct <-
    function(object, data, control = list(singular.ok = FALSE), ...)
{
    if (length(object)) {
        object[] <- lapply(object, Initialize, data)
        theta <- lapply(object, coef)
        len <- lengths(theta)
        num <- seq_along(len)
        pmap <-
            if (sum(len) > 0)
                outer(rep(num, len), num, "==")
            else
                array(FALSE, c(1, length(len)))
        dimnames(pmap) <- list(NULL, names(object))
        attr(object, "pmap") <- pmap
        attr(object, "glsFit") <-
            glsEstimate(object, control = control)
        if (needUpdate(object)) {
            object <- update(object, data)
        }
    }
    object
}

logLik.glsStruct <-
    function(object, Pars, conLin = attr(object, "conLin"), ...)
{
    coef(object) <- Pars			# updating parameter values
    conLin <- recalc(object, conLin)	# updating conLin
    val <- .C(gls_loglik,
              as.double(conLin[["Xy"]]),
              as.integer(unlist(conLin[["dims"]])),
              logLik = as.double(conLin[["logLik"]]),
              ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions
              double(1L), as.double(conLin$sigma), NAOK = TRUE)
    val[["logLik"]]
}

residuals.glsStruct <-
    function(object, glsFit = attr(object, "glsFit"), ...)
{
    glsFit[["resid"]]
}

varWeights.glsStruct <-
    function(object)
{
    if (is.null(object$varStruct)) rep(1, attr(object, "conLin")$dims$N)
    else varWeights(object$varStruct)
}

## Auxiliary control functions

glsControl <-
    ## Control parameters for gls
    function(maxIter = 50L, msMaxIter = 200L, tolerance = 1e-6, msTol = 1e-7,
             msVerbose = FALSE, singular.ok = FALSE,
             returnObject = FALSE,
             apVar = TRUE, .relStep = .Machine$double.eps^(1/3),
	     opt = c("nlminb", "optim"),  optimMethod = "BFGS",
             minAbsParApVar = 0.05, natural = TRUE, sigma = NULL)
{
    ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions
    if(is.null(sigma))
	sigma <- 0
    else {
	if(!is.finite(sigma) || length(sigma) != 1 || sigma < 0)
	    stop("Within-group std. dev. must be a positive numeric value")
	## if(missing(apVar)) apVar <- FALSE # not yet implemented
    }
    list(maxIter = maxIter, msMaxIter = msMaxIter, tolerance = tolerance,
         msTol = msTol, msVerbose = msVerbose,
         singular.ok = singular.ok,
         returnObject = returnObject, apVar = apVar,
         minAbsParApVar = minAbsParApVar, .relStep = .relStep,
         opt = match.arg(opt),
	 optimMethod = optimMethod, natural = natural, sigma = sigma)
}