R/varFunc.R

In nlme: Linear and Nonlinear Mixed Effects Models

###              Classes of variance functions
###
### Copyright 2007-2021  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/
#

##*## Generics that should be implemented for any varFunc class

varWeights <-
  ## Calculates the weights of the variance function
  function(object) UseMethod("varWeights")

##*## varFunc - a virtual class of variance functions

###*# Constructor

varFunc <-
  ## Can take as argument either a varFunc object, in which case it does
  ## nothing, a formula or a character string , in which case it
  ## calls varFixed
  function(object)
{
  if(is.null(object)) return(object)	# NULL object - no varFunc structure
  if (inherits(object, "varFunc")) {
    ## constructing from another varFunc object
    return(object)
  }
  if (inherits(object, "formula") || is.character(object)) {
    ## constructing from a formula of the form ~ x
    return(varFixed(asOneSidedFormula(object)))
  }

  stop("can only construct \"varFunc\" object from another \"varFunc\" object, a formula, or a character string")
}


###*# Methods for local generics

varWeights.varFunc <-
  function(object) attr(object, "weights")

###*# Methods for standard generics

coef.varFunc <-
  function(object, unconstrained = TRUE, allCoef = FALSE, ...)
{
  ### checking if initialized
  wPar <- attr(object, "whichFix")
  if (is.null(wPar) ||
      (length(object) != (length(wPar) - sum(wPar)))) {
    stop("cannot extract parameters of uninitialized object")
  }
  if (unconstrained) {
    if (allCoef) {
      val <- double(length(wPar))
      if (any(wPar)) {
        val[wPar] <- attr(object, "fixed")
      }
      if (any(!wPar)) {
        val[!wPar] <- as.vector(object)
      }
    } else {
      val <- as.vector(object)
    }
    val
  } else {
      stop(gettextf("do not know how to get coefficients for %s object",
                    dQuote(class(object)[1])), domain = NA)
  }
}

"covariate<-.varFunc" <-
  function(object, value)
{
  value <- as.numeric(value)
  if (!is.null(aux <- getCovariate(object))) {
    if (length(aux) != length(value)) {
      stop("cannot change the length of covariate in \"varFunc\" object")
    }
  }
  attr(object, "covariate") <- value
  object
}

formula.varFunc <-
  function(x, ...) eval(attr(x, "formula"))

getCovariate.varFunc <-
  function(object, form, data) attr(object, "covariate")

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

Initialize.varFunc <-
  function(object, data, ...)
{
  if (is.null(varWeights(object))) {
    attr(object, "weights") <- rep(1, dim(data)[1])
  }
  if (is.null(logLik(object))) {
    attr(object, "logLik") <- 0
  }
  object
}

## NB no nobs

logLik.varFunc <- function(object, data, ...)
{
    if (is.null(ll <- attr(object, "logLik")))
        ll
    else
        structure(ll, df = length(object), class = "logLik")
}

print.summary.varFunc <-
  function(x, header = TRUE, ...)
{
  ox <- x
  class(x) <- attr(x, "oClass")
  if (length(aux <- coef(x, unconstrained = FALSE, allCoef = TRUE)) > 0) {
    if (header) cat("Variance function:\n")
    cat(paste(" Structure: ", attr(x, "structName"), "\n", sep = ""))
    cat(paste(" Formula:", deparse(formula(x)),"\n"))
    cat(" Parameter estimates:\n")
    print(aux)
  } else {
    if (inherits(x, "varIdent")) {
      ## varIdent with no parameters - nothing to print
      return()
    }
    cat("Variance function structure of class", class(x)[1],
	"with no parameters, or uninitialized\n")
  }
  invisible(ox)
}

print.varFunc <-
  function(x, ...)
{
  if (length(aux <- coef(x, unconstrained = FALSE, allCoef = TRUE)) > 0) {
    cat("Variance function structure of class", class(x)[1],
	"representing\n")
    print(aux, ...)
  } else {
    cat("Variance function structure of class", class(x)[1],
	"with no parameters, or uninitialized\n")
  }
  invisible(x)
}

recalc.varFunc <-
  function(object, conLin, ...)
{
  conLin$Xy[] <- conLin$Xy * varWeights(object)
  conLin$logLik <- conLin$logLik + logLik(object)
  conLin
}

summary.varFunc <-
  function(object, structName = class(object)[1], ...)
{
  attr(object, "structName") <- structName
  attr(object, "oClass") <- class(object)
  class(object) <- "summary.varFunc"
  object
}

update.varFunc <-
  function(object, data, ...)
{
  if (needUpdate(object)) {
    covariate(object) <-
      eval(getCovariateFormula(object)[[2]], data)
  }
  object
}

##*## Classes that substitute for (i.e. inherit from) varFunc

###*# varFixed - fixed weights

####* Constructor

varFixed <-
  function(value = ~ 1)
{
  if (!inherits(value, "formula")) {
    stop("'value' must be a one sided formula")
  }
  form <- asOneSidedFormula(value)
  if (length(all.vars(getCovariateFormula(form))) == 0) {
    stop("'form' must have a covariate")
  }
  if (!is.null(getGroupsFormula(form))) {
    form <- getCovariateFormula(form)
    warning("ignoring 'group' in \"varFixed\" formula")
  }
  value <- numeric(0)
  attr(value, "formula") <- form
  class(value) <- c("varFixed", "varFunc")
  value
}

###*# Methods for standard generics

coef.varFixed <-
  function(object, ...) numeric(0)

"coef<-.varFixed" <-
  function(object, ..., value) object

Initialize.varFixed <-
  function(object, data, ...)
{
  form <- formula(object)
  if (any(is.na(match(all.vars(form), names(data))))) {
    ## cannot evaluate covariate on data
    stop("all variables used in 'formula' must be in 'data'")
  }
  attr(object, "needUpdate") <- FALSE
  attr(object, "covariate") <- getCovariate(data, form)
  attr(object, "logLik") <-
    sum(log(attr(object, "weights") <- 1/sqrt(abs(attr(object,"covariate")))))
  object
}

print.summary.varFixed <-
  function(x, header = TRUE, ...)
{
  cat("Variance function:\n")
  cat(" Structure: fixed weights\n")
  cat(paste(" Formula:", deparse(formula(x)),"\n"))
  invisible(x)
}

summary.varFixed <-
  function(object, structName, ...)
{
  class(object) <- "summary.varFixed"
  object
}

###*# varIdent - equal variances per stratum structure

####* Constructor

varIdent <-
  function(value = numeric(0), form = ~ 1, fixed = NULL)
{
  if (is.null(getGroupsFormula(form))) { # constant value
    value <- numeric(0)
    attr(value, "fixed") <- NULL	# nothing to estimate
  } else {
    if ((lv <- length(value)) > 0) {		# initialized
      if (is.null(grpNames <- names(value)) && (lv > 1)) {
	stop("initial values must have group names in 'varIdent'")
      }
      value <- unlist(value)		# may be a list with names
      if (any(value <= 0)) {
	stop("initial values for 'varIdent' must be > 0")
      }
      value <- log(value)               # unconstrained form
    } else grpNames <- NULL
    attr(value, "groupNames") <- grpNames
    if (!is.null(fixed)) {
       fix <- attr(value, "fixed") <- log(unlist(fixed))
      if (is.null(fixNames <- names(fix))) {
	stop("fixed parameters must have names in 'varIdent'")
      }
      if (!is.null(attr(value, "groupNames"))) {
	attr(value, "groupNames") <- c(attr(value, "groupNames"), fixNames)
      }
    }
  }
  attr(value, "formula") <- asOneSidedFormula(form)
  class(value) <- c("varIdent", "varFunc")
  value
}

###*# Methods for standard generics

coef.varIdent <-
  function(object, unconstrained = TRUE, allCoef = FALSE, ...)
{
  if (!is.null(getGroupsFormula(object)) &&
      !is.null( wPar <- attr(object, "whichFix"))) {
    ## different groups variances
    if (unconstrained && !allCoef) {
      return(as.vector(object))
    }
    val <- double(length(wPar))
    if (any(wPar)) {
      val[wPar] <- attr(object, "fixed")
    }
    if (any(!wPar)) {
      val[!wPar] <- as.vector(object)
    }
    if (!unconstrained) {
      val <- c(1, exp(val))
      names(val) <- attr(object, "groupNames")
      if (!allCoef) {
	val <- val[c(FALSE, !attr(object, "whichFix"))]
      }
    }
    val
  } else {
    numeric(0)
  }
}

"coef<-.varIdent" <-
  function(object, ..., value)
{
  if (!(is.null(grps <- getGroups(object)) ||
       all(attr(object, "whichFix")))) {
    ## different group variances & varying parameters
    value <- as.numeric(value)
    nGroups <- length(attr(object, "groupNames"))
#    if (nGroups == 0) {
#      stop("Cannot assign parameters of uninitialized varIdent object")
#    }
    ## fix from PR#9831
    nFixed  <- sum(as.numeric(attr(object,"whichFix")))
    if (length(value) != nGroups - nFixed - 1) {
      stop("cannot change the length of the \"varIdent\" parameter after initialization")
    }
    object[] <- value
    natPar <- coef(object, FALSE, allCoef = TRUE)
    attr(object, "logLik") <-
      sum(log(attr(object, "weights") <- 1/natPar[grps]))
  }
  object
}

Initialize.varIdent <-
  function(object, data, ...)
{
  if (!is.null(form <- formula(object)) &&
      !is.null(grpForm <- getGroupsFormula(form))) {
    if (length(coef(object)) > 0) { # initialized - nothing to do
      return(object)
    }
    strat <- attr(object, "groups") <-
      as.character(getGroups(data, form,
                             level = length(splitFormula(grpForm, sep = "*")),
                             sep = "*"))
    if (length((uStrat <- unique(strat))) == 1) {
      ## equal variances structure
      return(Initialize(varIdent(), data))
    }
    if (!is.null(fix <- attr(object, "fixed"))) {
      fixNames <- names(fix)
      if (any(is.na(match(fixNames, uStrat)))) {
	stop("fixed parameter names in 'varIdent' must be a subset of group names")
      }
      uStratVar <- uStrat[is.na(match(uStrat, fixNames))] # varying strata
      uStrat <- c(uStratVar, fixNames)
    } else {				# nothing fixed
      uStratVar <- uStrat
    }
    if ((nStratVar <- length(uStratVar)) == 0) {
      stop("cannot fix variances in all groups")
    }
    if (nStratVar > 1) {
      if (length(object) <= 1) {
	## repeat for all groups
	oldAttr <- attributes(object)
	if (length(object) > 0) {		# initialized
	  object <- rep(as.vector(object), nStratVar - 1)
	} else {			# uninitialized
	  object <- rep(0, nStratVar - 1)
	}
	attributes(object) <- oldAttr
	attr(object, "groupNames") <- uStrat
      } else {
	if (length(as.vector(object)) != (len <- (nStratVar - 1))) {
            stop(gettextf("initial value for \"varIdent\" should be of length %d",
                          len), domain = NA)
	}
	if (!is.null(stN <- attr(object, "groupNames"))) {
	  missStrat <- uStrat[is.na(match(uStrat, stN))]
	  if (length(missStrat) != 1) {
	    stop("names of starting value for \"varIdent\" object must contain all but one of the stratum levels")
	  }
	  stN <-  c(missStrat, stN)
	  if ((length(stN) != length(uStrat)) ||
	      any(sort(stN) != sort(uStrat))) {
	    stop("nonexistent group names for initial values in 'varIdent'")
	  }
	  attr(object, "groupNames") <- stN
	} else {
	  attr(object, "groupNames") <- uStrat
	}
      }
    } else {				# fixed for all but one strata
      oldAttr <- attributes(object)
      object <- numeric(0)
      attributes(object) <- oldAttr
      attr(object, "groupNames") <- uStrat
    }
    attr(object, "whichFix") <-
      !is.na(match(attr(object, "groupNames")[-1], names(fix)))
    if (all(attr(object, "whichFix"))) {
      if (all(attr(object, "fixed") == 0)) {
	## equal variances structure
	return(Initialize(varIdent(), data))
      } else {
	oldAttr <- attributes(object)
	object <- numeric(0)
	attributes(object) <- oldAttr
      }
    }
    ## initializing weights and logDet
    attr(object, "logLik") <-
      sum(log(attr(object, "weights") <-
              1/coef(object, FALSE, allCoef = TRUE)[strat]))
    object
  } else {				# no strata
    attr(object, "whichFix") <- TRUE
    NextMethod()
  }
}

needUpdate.varIdent <-
  function(object) FALSE

recalc.varIdent <-
  function(object, conLin, ...)
{
  if (is.null(formula(object))) conLin else NextMethod()
}

summary.varIdent <-
  function(object,
	   structName = if (is.null(formula(object))) "Constant variance"
	                else "Different standard deviations per stratum",
           ...)
{
    summary.varFunc(object, structName)
}


###*# varPower - power of variance covariate variance structure

####* Constructor

varPower <-
  function(value = numeric(0), form = ~ fitted(.), fixed = NULL)
{
  value <- unlist(value)		# may be given as a list
  fixed <- attr(value, "fixed") <- unlist(fixed)
  attr(value, "formula") <- form <- asOneSidedFormula(form)
  if (length(all.vars(getCovariateFormula(form))) == 0) {
    stop("'form' must have a covariate")
  }
  if (!is.null(getGroupsFormula(form))) {
    if (is.null(grpNames <- names(value)) && (length(value) > 1)) {
      stop("initial values must have group names in 'varPower'")
    }
    if (!is.null(fixed)) {
      if (is.null(names(fixed))) {
	stop("fixed parameters must have group names in 'varPower'")
      }
    }
    attr(value, "groupNames") <- c(grpNames, names(fixed))
  } else {                              # single parameter
    attr(value, "whichFix") <- !is.null(fixed)
  }
  class(value) <- c("varPower", "varFunc")
  value
}

###*# Methods for standard generics

coef.varPower <-
  function(object, unconstrained = TRUE, allCoef = FALSE, ...)
{
  if (((length(object) == 0) &&
       (!allCoef || is.null(attr(object, "fixed")))) ||
      is.null(wPar <- attr(object, "whichFix"))) {
    ## uninitialized
    return(numeric(0))
  }
  val <- double(length(wPar))
  if (any(wPar)) {
    val[wPar] <- attr(object, "fixed")
  }
  if (any(!wPar)) {
    val[!wPar] <- as.vector(object)
  }
  if (!is.null(getGroupsFormula(object))) {
    ##different values per group
    names(val) <- attr(object, "groupNames")
  } else {
    names(val) <- "power"
  }
  if (!allCoef) {
    val <- val[!wPar]
  }
  val
}

"coef<-.varPower" <-
  function(object, ..., value)
{
  if ((len <- length(object)) > 0) {		# varying parameters
    value <- as.numeric(value)
    if (length(value) != len) {
      stop("cannot change the length of the \"varStruct\" parameter after initialization")
    }
    object[] <- value
    aux <- coef(object, FALSE, allCoef = TRUE)
    if (!is.null(grps <- getGroups(object))) {
      aux <- aux[grps]
    }
    covariateObj <- getCovariate(object)
    if(is.null(covariateObj)) covariateObj <- NA
    attr(object, "logLik") <-
      sum(log(attr(object, "weights") <- abs(covariateObj)^(-aux)))
  } else {
    stop("cannot change coefficients before initialization or when all parameters are fixed")
  }
  object
}

Initialize.varPower <-
  function(object, data, ...)
{
  form <- formula(object)
  if (all(!is.na(match(all.vars(getCovariateFormula(form)), names(data))))) {
    ## can evaluate covariate on data
    attr(object, "needUpdate") <- FALSE
    attr(object, "covariate") <- getCovariate(data, form)
  } else {
    attr(object, "needUpdate") <- TRUE
  }
  if (!is.null(grpForm <- getGroupsFormula(form))) {
    strat <- as.character(getGroups(data, form,
                            level = length(splitFormula(grpForm, sep = "*")),
                            sep = "*"))
    uStrat <- unique(strat)
    if (length(uStrat) > 1) {		# multi-groups
      attr(object, "groups") <- strat
      if (!is.null(attr(object, "fixed"))) {
	fixNames <- names(attr(object, "fixed"))
	if (is.null(fixNames)) {
	  stop("fixed parameters must have group names")
	}
	if (any(is.na(match(fixNames, uStrat)))) {
	  stop("mismatch between group names and fixed values names")
	}
      } else {
	fixNames <- NULL
      }
      uStratVar <- uStrat[is.na(match(uStrat, fixNames))]
      nStratVar <- length(uStratVar)
      attr(object, "whichFix") <- !is.na(match(uStrat, fixNames))
      if (nStratVar > 0) {
	if (length(object) <= 1) {
	  ## repeat for all groups
	  names(object) <- NULL
	  oldAttr <- attributes(object)
	  if (length(object) > 0) {
	    object <- rep(as.vector(object), nStratVar)
	  } else {
	    object <- rep(0, nStratVar)
	  }
	  attributes(object) <- oldAttr
	  attr(object, "groupNames") <- uStrat
	  names(object) <- uStratVar
	} else {
	  if (length(as.vector(object)) != nStratVar) {
              stop(gettextf("initial value for \"varPower\" should be of length %d",
                            nStratVar), domain = NA)
	  }
	  stN <- attr(object, "groupNames") # must have names
	  if (length(stN) != length(uStrat) ||
	      any(sort(stN) != sort(uStrat))) {
	    stop("nonexistent group names for initial values in \"varPower\"")
	  }
	}
      } else {				# all parameters are fixed
	if (all(attr(object, "fixed") == 0)) {
	  ## equal variances structure
	  return(Initialize(varIdent(), data))
	} else {
	  oldAttr <- attributes(object)
	  object <- numeric(0)
	  attributes(object) <- oldAttr
	  attr(object, "groupNames") <- uStrat
	}
      }
    } else {                            # single stratum
      attr(object, "formula") <- getCovariateFormula(formula(object))
      attr(object, "whichFix") <- !is.null(attr(object, "fixed"))
    }
  }
  if (is.null(getGroupsFormula(object))) {
    ## single stratum
    if (attr(object, "whichFix")) {
      if (attr(object, "fixed") == 0) {
        ## equal variances structure
        return(Initialize(varIdent(), data))
      } else {				# fixed power
        oldAttr <- attributes(object)
        object <- numeric(0)
        attributes(object) <- oldAttr
      }
    } else {
      len <- length(as.vector(object))
      if (len == 0) {			# uninitialized
        oldAttr <- attributes(object)
        object <- 0
        attributes(object) <- oldAttr
      } else if (len > 1) {
        stop("initial value for \"varPower\" should be of length 1")
      }
    }
  }
  if (!is.null(covar <- getCovariate(object))) {
    natPar <- coef(object, allCoef = TRUE)
    if (!is.null(grps <- getGroups(object))) {
      natPar <- natPar[grps]
    }
    attr(object, "logLik") <-
      sum(log(attr(object, "weights") <- abs(covar^(-natPar))))
    object
  } else {
    NextMethod()
  }
}

summary.varPower <-
  function(object, structName = "Power of variance covariate", ...)
{
  if (!is.null(getGroupsFormula(object))) {
    structName <- paste(structName, " different strata", sep = ",")
  }
  summary.varFunc(object, structName)
}

update.varPower <-
  function(object, data, ...)
{
  val <- NextMethod()
  if (length(val) == 0) {		# chance to update weights
    aux <- coef(val, allCoef = TRUE)
    if (!is.null(grps <- getGroups(val))) {
      aux <- aux[grps]
    }
    attr(val, "logLik") <-
      sum(log(attr(val, "weights") <- abs(getCovariate(val))^(-aux)))
  }
  val
}

###*# varExp - exponential of variance covariate variance structure

####* Constructor

varExp <-
  function(value = numeric(0), form = ~ fitted(.), fixed = NULL)
{
  value <- unlist(value)		# may be given as a list
  fixed <- attr(value, "fixed") <- unlist(fixed)
  attr(value, "formula") <- form <- asOneSidedFormula(form)
  if (length(all.vars(getCovariateFormula(form))) == 0) {
    stop("'form' must have a covariate")
  }
  if (!is.null(getGroupsFormula(form))) {
    if (is.null(grpNames <- names(value)) && (length(value) > 1)) {
      stop("initial values must have group names in 'varExp'")
    }
    if (!is.null(fixed)) {
      if (is.null(names(fixed))) {
	stop("fixed parameters must have group names in 'varExp'")
      }
    }
    attr(value, "groupNames") <- c(grpNames, names(fixed))
  } else {
    attr(value, "whichFix") <- !is.null(fixed)
  }
  class(value) <- c("varExp", "varFunc")
  value
}

###*# Methods for standard generics

coef.varExp <-
  function(object, unconstrained = TRUE, allCoef = FALSE, ...)
{
  if (((length(object) == 0) &&
       (!allCoef || is.null(attr(object, "fixed")))) ||
      is.null( wPar <- attr(object, "whichFix"))) {
    return(numeric(0))
  }
  val <- double(length(wPar))
  if (any(wPar)) {
    val[wPar] <- attr(object, "fixed")
  }
  if (any(!wPar)) {
    val[!wPar] <- as.vector(object)
  }
  if (!is.null(getGroupsFormula(object))) {
    ##different values per group
    names(val) <- attr(object, "groupNames")
  } else {
    names(val) <- "expon"
  }
  if (!allCoef) {
    val <- val[!wPar]
  }
  val
}

"coef<-.varExp" <-
  function(object, ..., value)
{
  if (length(object) > 0) {		# varying parameters
    value <- as.numeric(value)
    if (length(value) != length(object)) {
      stop("cannot change the length of the \"varExp\" parameter after initialization")
    }
    object[] <- value
    aux <- coef(object, FALSE, allCoef = TRUE)
    if (!is.null(grps <- getGroups(object))) {
      aux <- aux[grps]
    }
    attr(object, "logLik") <-
      sum(log(attr(object, "weights") <- exp(-aux * getCovariate(object))))
  } else {
    stop("cannot change coefficients before initialization or when all parameters are fixed")
  }
  object
}

Initialize.varExp <-
  function(object, data, ...)
{
  form <- formula(object)
  if (all(!is.na(match(all.vars(getCovariateFormula(form)), names(data))))) {
    ## can evaluate covariate on data
    attr(object, "needUpdate") <- FALSE
    attr(object, "covariate") <- getCovariate(data, form)
  } else {
    attr(object, "needUpdate") <- TRUE
  }
  if (!is.null(grpForm <- getGroupsFormula(form))) {
    strat <- as.character(getGroups(data, form,
                            level = length(splitFormula(grpForm, sep = "*")),
                            sep = "*"))
    uStrat <- unique(strat)
    if (length(uStrat) > 1) {		# multi-groups
      attr(object, "groups") <- strat
      if (!is.null(attr(object, "fixed"))) {
	fixNames <- names(attr(object, "fixed"))
	if (is.null(fixNames)) {
	  stop("fixed parameters must have group names")
	}
	if (any(is.na(match(fixNames, uStrat)))) {
	  stop("mismatch between group names and fixed values names")
	}
      } else {
	fixNames <- NULL
      }
      uStratVar <- uStrat[is.na(match(uStrat, fixNames))]
      nStratVar <- length(uStratVar)
      attr(object, "whichFix") <- !is.na(match(uStrat, fixNames))
      if (nStratVar > 0) {
	if (length(object) <= 1) {
	  ## repeat for all groups
	  names(object) <- NULL
	  oldAttr <- attributes(object)
	  if (length(object) > 0) {
	    object <- rep(as.vector(object), nStratVar)
	  } else {
	    object <- rep(0, nStratVar)
	  }
	  attributes(object) <- oldAttr
	  attr(object, "groupNames") <- uStrat
	  names(object) <- uStratVar
	} else {
	  if (length(as.vector(object)) != nStratVar) {
              stop(gettextf("initial value for \"varExp\" should be of length %d",
                            nStratVar), domain = NA)
	  }
	  stN <- attr(object, "groupNames") #must have names
	  if ((length(stN) != length(uStrat)) ||
	      any(sort(stN) != sort(uStrat))) {
	    stop("nonexistent group names for initial values in \"varExp\"")
	  }
	}
      } else {
	if (all(attr(object, "fixed") == 0)) {
	  ## equal variances structure
	  return(Initialize(varIdent(), data))
	} else {
	  oldAttr <- attributes(object)
	  object <- numeric(0)
	  attributes(object) <- oldAttr
	  attr(object, "groupNames") <- uStrat
	}
      }
    } else {                            # single stratum
      attr(object, "formula") <- getCovariateFormula(formula(object))
      attr(object, "whichFix") <- !is.null(attr(object, "fixed"))
    }
  }
  if (is.null(getGroupsFormula(object))) {
    ## single stratum
    if (attr(object, "whichFix")) {
      if (!attr(object, "fixed")) {
        ## equal variances structure
        return(Initialize(varIdent(), data))
      } else {
        oldAttr <- attributes(object)
        object <- numeric(0)
        attributes(object) <- oldAttr
      }
    } else {
      len <- length(as.vector(object))
      if (len == 0) {			# uninitialized
        oldAttr <- attributes(object)
        object <- 0
        attributes(object) <- oldAttr
      } else if (len > 1) {
        stop("initial value for \"varExp\" should be of length 1")
      }
    }
  }
  if (!is.null(covar <- getCovariate(object))) {
    natPar <- coef(object, allCoef = TRUE)
    if (!is.null(grps <- getGroups(object))) {
      natPar <- natPar[grps]
    }
    attr(object, "logLik") <-
      sum(log(attr(object, "weights") <- exp(-natPar * covar)))
    object
  } else {
    NextMethod()
  }
}


summary.varExp <-
  function(object, structName = "Exponential of variance covariate", ...)
{
  if (!is.null(getGroupsFormula(object))) {
    structName <- paste(structName, " different strata", sep = ",")
  }
  summary.varFunc(object, structName)
}

update.varExp <-
  function(object, data, ...)
{
  val <- NextMethod()
  if (length(val) == 0) {		# chance to update weights
    aux <- coef(val, allCoef = TRUE)
    if (!is.null(grps <- getGroups(val))) {
      aux <- aux[grps]
    }
    attr(val, "logLik") <-
      sum(log(attr(val, "weights") <- exp(-aux * getCovariate(val))))
  }
  val
}

###*# varConstPower - Constant plus power of covariance function
###*#               variance structure

####* Constructor for the varConstPower class
varConstPower <-
  function(const = numeric(0), power = numeric(0),
	   form = ~ fitted(.), fixed = NULL)
{
  value <- list(const = const, power = power)
  form <- asOneSidedFormula(form)
  if (length(all.vars(getCovariateFormula(form))) == 0)
    stop("'form' must have a covariate")
  CPconstr <- function(val, form, nam) {
    if ((lv <- length(val)) == 0) return(val)
    if (lv > 2) {
      stop(gettextf("%s can have at most two components", nam), domain = NA)
    }
    if (is.null(nv <- names(val))) {
      names(val) <- c("const", "power")[1:lv]
    } else {
      if (any(is.na(match(nv, c("const", "power"))))) {
        stop(gettextf("%s can only have names \"const\" and \"power\"", nam),
             domain = NA)
      }
    }
    nv <- names(val)
    if (data.class(val) == "list") {
      val <- lapply(val, unlist)
      grpNames <- unique(unlist(lapply(val, names)))
    } else {				# must be a vector or a scalar
      if (!is.numeric(val)) {
        stop(gettextf("%s can only be a list or numeric", nam), domain = NA)
     }
      val <- as.list(val)
      names(val) <- nv
      grpNames <- NULL
    }
    if (!is.null(getGroupsFormula(form))) {
      if (any(vapply(val, function(el) length(el) > 1 && is.null(names(el)), NA)))
        stop(gettextf("%s must have group names in 'varConstPower'", nam),
             domain = NA)
      attr(val, "groupNames") <- grpNames
    }
    if (length(const <- val$const) > 0) {
      if (any(const <= 0))
	stop("constant in \"varConstProp\" structure must be > 0")
      const <- log(const)
    }
    list(const = const, power = val$power)
  }
  ## initial value may be given as a vector or list. If groups are
  ## present and different initial values are given for each group, then
  ## it must be a list with components "const" and/or "power"
  value <- CPconstr(value, form, "Value")
  fixed <- CPconstr(fixed, form, "Fixed")
  attr(value, "formula") <- form
  attr(value, "groupNames") <-
    unique(c(attr(value, "groupNames"),
	   attr(attr(value[["const"]], "fixed"), "groupNames"),
	   attr(attr(value[["power"]], "fixed"), "groupNames")))
  for (i in names(fixed)) {
    attr(value[[i]], "fixed") <- c(fixed[[i]])
  }
  if (is.null(getGroupsFormula(form))) {   # no groups
    whichFix <- array(FALSE, c(2,1), list(c("const", "power"), NULL))
    whichFix[,1] <- vapply(value, function(el) !is.null(attr(el, "fixed")), NA)
    attr(value, "whichFix") <- whichFix
  }
  class(value) <- c("varConstPower", "varFunc")
  value
}

###*# Methods for standard generics

coef.varConstPower <-
    function (object, unconstrained = TRUE, allCoef = FALSE, ...)
{
    wPar <- attr(object, "whichFix")
    nonInit <- !lengths(object)
    nonInit <- is.null(wPar) || (any(nonInit) && !all(c(wPar[nonInit, ])))
    if (nonInit || (!allCoef && (length(unlist(object)) == 0))) {
        return(numeric(0))
    }
    val <- array(0, dim(wPar), dimnames(wPar))
    for (i in names(object)) {
        if (any(wPar[i, ])) {
            val[i, wPar[i, ]] <- attr(object[[i]], "fixed")
        }
        if (any(!wPar[i, ])) {
            val[i, !wPar[i, ]] <- c(object[[i]])
        }
    }
    if (!unconstrained) {
        val[1, ] <- exp(val[1, ])
    }
    if (!allCoef) {
        val <- list(const = if (!all(wPar[1, ])) val[1, !wPar[1, ]], # else NULL
                    power = if (!all(wPar[2, ])) val[2, !wPar[2, ]])
        val <- lapply(val, function(el) if(length(el) == 1) as.vector(el) else el)
        unlist(val[!vapply(val, is.null, NA)])
    }
    else {
        val[, 1:ncol(val)]
    }
}

"coef<-.varConstPower" <-
  function(object, ..., value)
{
  if (length(unlist(object)) > 0) {	# varying parameters
    value <- as.numeric(value)
    if (length(value) != length(unlist(object))) {
      stop("cannot change the length of the parameter after initialization")
    }
    start <- 0
    for(i in names(object)) {
      if (aux <- length(object[[i]])) {
	object[[i]][] <- value[start + (1:aux)]
	start <- start + aux
      }
    }
    natPar <- as.matrix(coef(object, FALSE, allCoef = TRUE))
    if (!is.null(grps <- getGroups(object))) {
      natPar <- natPar[, grps]
    }
    attr(object, "logLik") <-
      sum(log(attr(object, "weights") <-
	      1/(natPar[1,] + abs(getCovariate(object))^natPar[2,])))
  } else {
    stop("cannot change coefficients before initialization or when all parameters are fixed")
  }
  object
}

Initialize.varConstPower <-
  function(object, data, ...)
{
  form <- formula(object)
  if (all(!is.na(match(all.vars(getCovariateFormula(form)), names(data))))) {
    ## can evaluate covariate on data
    attr(object, "needUpdate") <- FALSE
    attr(object, "covariate") <- getCovariate(data, form)
  } else {
    attr(object, "needUpdate") <- TRUE
  }
  dfltCoef <- c(const = log(0.1), power = 0)
  if (!is.null(grpForm <- getGroupsFormula(form))) {
    strat <- as.character(getGroups(data, form,
                            level = length(splitFormula(grpForm, sep = "*")),
                            sep = "*"))
    uStrat <- unique(strat)
    whichFix <- array(FALSE, c(2, length(uStrat)),
		      list(c("const", "power"), uStrat))
    if (length(uStrat) > 1) {		# multi-groups
      attr(object, "groups") <- strat
      for(i in names(object)) {
	if (!is.null(attr(object[[i]], "fixed"))) {
	  fixNames <- names(attr(object[[i]], "fixed"))
	  if (is.null(fixNames)) {
	    stop("fixed parameters must have group names")
	  }
	  if (any(is.na(match(fixNames, uStrat)))) {
	    stop("mismatch between group names and fixed values names")
	  }
	} else {
	  fixNames <- NULL
	}
	uStratVar <- uStrat[is.na(match(uStrat, fixNames))]
	nStratVar <- length(uStratVar)
	whichFix[i,] <- !is.na(match(uStrat, fixNames))
	if (nStratVar > 0) {
	  if (length(object[[i]]) <= 1) {
	    ## repeat for all groups
	    names(object[[i]]) <- NULL
	    oldAttr <- attributes(object[[i]])
	    if (length(object[[i]]) > 0) {
	      object[[i]] <- rep(as.vector(object[[i]]), nStratVar)
	    } else {
	      object[[i]] <- rep(dfltCoef[i], nStratVar)
	    }
	    attributes(object[[i]]) <- oldAttr
	    names(object[[i]]) <- uStratVar
	  } else {
	    if (length(as.vector(object[[i]])) != nStratVar) {
                stop(gettext("initial value should be of length %d",
                             nStratVar), domain = NA)
	    }
	    stN <- names(object[[i]]) # must have names
	    if ((length(stN) != length(uStratVar)) ||
		any(sort(stN) != sort(uStratVar))) {
	      stop("nonexistent group names for initial values")
	    }
	  }
	}
      }
      if (all(whichFix) &&
	  all(attr(object[["const"]], "fixed") == 0) &&
	  all(attr(object[["power"]], "fixed") == 0)) {
	## equal variances structure
	return(Initialize(varIdent(), data))
      }
      for(i in names(object)) {
	if (all(whichFix[i,])) {
	  oldAttr <- attributes(object[[i]])
	  object[[i]] <- numeric(0)
	  attributes(object[[i]]) <- oldAttr
	}
      }
      attr(object, "whichFix") <- whichFix
      attr(object, "groupNames") <- uStrat
      return(NextMethod())
    }
  }
  ## single stratum
  whichFix <- attr(object, "whichFix")
  if (all(whichFix) &&
      !any(unlist(lapply(object, function(el) attr(el, "fixed"))))) {
    ## equal variances structure
    return(Initialize(varIdent(), data))
  }
  for(i in names(object)) {
    if (all(whichFix[i,])) {
      oldAttr <- attributes(object[[i]])
      object[[i]] <- numeric(0)
      attributes(object[[i]]) <- oldAttr
    } else {
      if (length(object[[i]]) == 0) {
	object[[i]] <- dfltCoef[i]
      }
    }
  }
  aux <- 2 - sum(whichFix[,1])
  if (length(as.vector(unlist(object))) != aux) {
      stop(gettext("initial value should be of length %d", aux), domain = NA)
  }
  NextMethod()
}

summary.varConstPower <-
  function(object, structName = "Constant plus power of variance covariate",
           ...)
{
  if (!is.null(getGroupsFormula(object))) {
    structName <- paste(structName, " different strata", sep = ",")
  }
  summary.varFunc(object, structName)
}

update.varConstPower <-
  function(object, data, ...)
{
  val <- NextMethod()
  if (length(unlist(val)) == 0) {	# chance to update weights
    aux <- as.matrix(coef(val, FALSE, allCoef = TRUE))
    if (!is.null(grps <- getGroups(val))) {
      aux <- aux[, grps]
    }
    attr(val, "logLik") <-
      sum(log(attr(val, "weights") <-
	      1/(aux[1,] + abs(getCovariate(val))^aux[2,])))
  }
  val
}


## NB "varConstProp" was derived from "varConstPower" with minimal code changes.

###*# varConstProp - constant plus proportional variance function structure

####* Constructor  for the "varConstProp" class
varConstProp <-
  function(const = numeric(0), prop = numeric(0),
	   form = ~ fitted(.), fixed = NULL)
{
  value <- list(const = const, prop = prop)
  form <- asOneSidedFormula(form)
  if (length(all.vars(getCovariateFormula(form))) == 0)
    stop("'form' must have a covariate")
  CPconstr <- function(val, form, nam) {
    if ((lv <- length(val)) == 0) return(val)
    if (lv > 2) {
      stop(gettextf("%s can have at most two components", nam), domain = NA)
    }
    if (is.null(nv <- names(val))) {
      names(val) <- c("const", "prop")[1:lv]
    } else {
      if (any(is.na(match(nv, c("const", "prop"))))) {
        stop(gettextf("%s can only have names \"const\" and \"prop\"", nam),
             domain = NA)
      }
    }
    nv <- names(val)
    if (data.class(val) == "list") {
      val <- lapply(val, unlist)
      grpNames <- unique(unlist(lapply(val, names)))
    } else {				# must be a vector or a scalar
      if (!is.numeric(val)) {
        stop(gettextf("%s can only be a list or numeric", nam), domain = NA)
     }
      val <- as.list(val)
      names(val) <- nv
      grpNames <- NULL
    }
    if (!is.null(getGroupsFormula(form))) {
      if (any(vapply(val, function(el) length(el) > 1 && is.null(names(el)), NA)))
        stop(gettextf("%s must have group names in 'varConstProp'", nam),
             domain = NA)
      attr(val, "groupNames") <- grpNames
    }
    if (length(const <- val$const) > 0) {
      if (any(const <= 0))
	stop("constant in \"varConstProp\" structure must be > 0")
      const <- log(const)
    }
    list(const = const, prop = val$prop)
  }
  ## initial value may be given as a vector or list. If groups are
  ## present and different initial values are given for each group, then
  ## it must be a list with components "const" and/or "prop"
  value <- CPconstr(value, form, "Value")
  fixed <- CPconstr(fixed, form, "Fixed")
  attr(value, "formula") <- form
  attr(value, "groupNames") <-
    unique(c(attr(value, "groupNames"),
	   attr(attr(value[["const"]], "fixed"), "groupNames"),
	   attr(attr(value[["prop"]],  "fixed"), "groupNames")))
  for (i in names(fixed)) {
    attr(value[[i]], "fixed") <- c(fixed[[i]])
  }
  if (is.null(getGroupsFormula(form))) {   # no groups
    whichFix <- array(FALSE, c(2,1), list(c("const", "prop"), NULL))
    whichFix[,1] <- vapply(value, function(el) !is.null(attr(el, "fixed")), NA)
    attr(value, "whichFix") <- whichFix
  }
  class(value) <- c("varConstProp", "varFunc")
  value
}

###*# Methods for standard generics

coef.varConstProp <-
    function (object, unconstrained = TRUE, allCoef = FALSE, ...)
{
    wPar <- attr(object, "whichFix")
    nonInit <- !lengths(object)
    nonInit <- is.null(wPar) || (any(nonInit) && !all(c(wPar[nonInit, ])))
    if (nonInit || (!allCoef && (length(unlist(object)) == 0))) {
        return(numeric(0))
    }
    val <- array(0, dim(wPar), dimnames(wPar))
    for (i in names(object)) {
        if (any(wPar[i, ])) {
            val[i, wPar[i, ]] <- attr(object[[i]], "fixed")
        }
        if (any(!wPar[i, ])) {
            val[i, !wPar[i, ]] <- c(object[[i]])
        }
    }
    if (!unconstrained) {
        val[1, ] <- exp(val[1, ])
    }
    if (!allCoef) {
        val <- list(const = if (!all(wPar[1, ])) val[1, !wPar[1, ]], # else NULL
                    prop  = if (!all(wPar[2, ])) val[2, !wPar[2, ]])
        val <- lapply(val, function(el) if(length(el) == 1L) as.vector(el) else el)
        unlist(val[!vapply(val, is.null, NA)])
    }
    else {
        val[, 1:ncol(val)]
    }
}

"coef<-.varConstProp" <-
  function(object, ..., value)
{
  if (length(unlist(object)) > 0) {	# varying parameters
    value <- as.numeric(value)
    if (length(value) != length(unlist(object))) {
      stop("cannot change the length of the parameter after initialization")
    }
    start <- 0
    for(i in names(object)) {
      if (aux <- length(object[[i]])) {
	object[[i]][] <- value[start + (1:aux)]
	start <- start + aux
      }
    }
    natPar <- as.matrix(coef(object, FALSE, allCoef = TRUE))
    if (!is.null(grps <- getGroups(object))) {
      natPar <- natPar[, grps]
    }
    attr(object, "logLik") <-
      sum(log(attr(object, "weights") <-
	      1/sqrt(natPar[1,]^2 + natPar[2, ]^2 * getCovariate(object)^2)))
  } else {
    stop("cannot change coefficients before initialization or when all parameters are fixed")
  }
  object
}

Initialize.varConstProp <-
  function(object, data, ...)
{
  form <- formula(object)
  if (all(!is.na(match(all.vars(getCovariateFormula(form)), names(data))))) {
    ## can evaluate covariate on data
    attr(object, "needUpdate") <- FALSE
    attr(object, "covariate") <- getCovariate(data, form)
  } else {
    attr(object, "needUpdate") <- TRUE
  }
  dfltCoef <- c(const = 0.1, prop = 0.1)
  if (!is.null(grpForm <- getGroupsFormula(form))) {
    strat <- as.character(getGroups(data, form,
                            level = length(splitFormula(grpForm, sep = "*")),
                            sep = "*"))
    uStrat <- unique(strat)
    whichFix <- array(FALSE, c(2, length(uStrat)),
		      list(c("const", "prop"), uStrat))
    if (length(uStrat) > 1) {		# multi-groups
      attr(object, "groups") <- strat
      for(i in names(object)) {
	if (!is.null(attr(object[[i]], "fixed"))) {
	  fixNames <- names(attr(object[[i]], "fixed"))
	  if (is.null(fixNames)) {
	    stop("fixed parameters must have group names")
	  }
	  if (any(is.na(match(fixNames, uStrat)))) {
	    stop("mismatch between group names and fixed values names")
	  }
	} else {
	  fixNames <- NULL
	}
	uStratVar <- uStrat[is.na(match(uStrat, fixNames))]
	nStratVar <- length(uStratVar)
	whichFix[i,] <- !is.na(match(uStrat, fixNames))
	if (nStratVar > 0) {
	  if (length(object[[i]]) <= 1) {
	    ## repeat for all groups
	    names(object[[i]]) <- NULL
	    oldAttr <- attributes(object[[i]])
	    if (length(object[[i]]) > 0) {
	      object[[i]] <- rep(as.vector(object[[i]]), nStratVar)
	    } else {
	      object[[i]] <- rep(dfltCoef[i], nStratVar)
	    }
	    attributes(object[[i]]) <- oldAttr
	    names(object[[i]]) <- uStratVar
	  } else {
	    if (length(as.vector(object[[i]])) != nStratVar) {
                stop(gettext("initial value should be of length %d",
                             nStratVar), domain = NA)
	    }
	    stN <- names(object[[i]]) # must have names
	    if ((length(stN) != length(uStratVar)) ||
		any(sort(stN) != sort(uStratVar))) {
	      stop("nonexistent group names for initial values")
	    }
	  }
	}
      }
      if (all(whichFix) &&
	  all(attr(object[["const"]], "fixed") == 0) &&
	  all(attr(object[["prop"]], "fixed") == 0)) {
	## equal variances structure
	return(Initialize(varIdent(), data))
      }
      for(i in names(object)) {
	if (all(whichFix[i,])) {
	  oldAttr <- attributes(object[[i]])
	  object[[i]] <- numeric(0)
	  attributes(object[[i]]) <- oldAttr
	}
      }
      attr(object, "whichFix") <- whichFix
      attr(object, "groupNames") <- uStrat
      return(NextMethod())
    }
  }
  ## single stratum
  whichFix <- attr(object, "whichFix")
  if (all(whichFix) &&
      !any(unlist(lapply(object, function(el) attr(el, "fixed"))))) {
    ## equal variances structure
    return(Initialize(varIdent(), data))
  }
  for(i in names(object)) {
    if (all(whichFix[i,])) {
      oldAttr <- attributes(object[[i]])
      object[[i]] <- numeric(0)
      attributes(object[[i]]) <- oldAttr
    } else {
      if (length(object[[i]]) == 0) {
	object[[i]] <- dfltCoef[i]
      }
    }
  }
  aux <- 2 - sum(whichFix[,1])
  if (length(as.vector(unlist(object))) != aux) {
      stop(gettext("initial value should be of length %d", aux), domain = NA)
  }
  NextMethod()
}

summary.varConstProp <-
  function(object, structName = "Constant plus proportion of variance covariate",
           ...)
{
  if (!is.null(getGroupsFormula(object))) {
    structName <- paste(structName, " different strata", sep = ",")
  }
  summary.varFunc(object, structName)
}

update.varConstProp <-
  function(object, data, ...)
{
  val <- NextMethod()
  if (length(unlist(val)) == 0) {	# chance to update weights
    aux <- as.matrix(coef(val, FALSE, allCoef = TRUE))
    if (!is.null(grps <- getGroups(val))) {
      aux <- aux[, grps]
    }
    attr(val, "logLik") <-
      sum(log(attr(val, "weights") <-
	      1/sqrt((aux[1, ]^2 + aux[2, ]^2 * getCovariate(val)^2))))
  }
  val
}

###*# varComb - combination of variance function structures

####* Constructor for the  "varComb" class
varComb <-
  function(...)
{
  val <- list(...)
  if (!all(vapply(val, inherits, NA, "varFunc"))) {
    stop("all arguments to 'varComb' must be of class \"varFunc\".")
  }
  if (is.null(names(val))) {
    names(val) <- LETTERS[seq_along(val)]
  }
  class(val) <- c("varComb", "varFunc")
  val
}

####* Methods for local generics


varWeights.varComb <-
  function(object)
{
  apply(as.data.frame(lapply(object, varWeights)), 1, prod)
}

###*# Methods for standard generics

coef.varComb <-
  function(object, unconstrained = TRUE, allCoef = FALSE, ...)
{
  unlist(lapply(object, coef, unconstrained, allCoef))
}

"coef<-.varComb" <-
  function(object, ..., value)
{
  plen <- attr(object, "plen")
  if ((len <- sum(plen)) > 0) {		# varying parameters
    if (length(value) != len) {
      stop("cannot change parameter length of initialized \"varComb\" object")
    }
    start <- 0
    for (i in seq_along(object)) {
      if (plen[i] > 0) {
	coef(object[[i]]) <- value[start + (1:plen[i])]
	start <- start + plen[i]
      }
    }
  }
  object
}

formula.varComb <-
  function(x, ...) lapply(x, formula)

Initialize.varComb <- function(object, data, ...)
{
    val <- lapply(object, Initialize, data)
    structure(val,
              plen = vapply(val, function(el) length(coef(el)), 1L),
              class = c("varComb", "varFunc"))
}

logLik.varComb <- function(object, ...)
{
    lls <- lapply(object, logLik)
    structure(sum(unlist(lls)),
              df = sum(vapply(lls, attr, 1, "df")),
              class = "logLik")
}

needUpdate.varComb <- function(object) any(vapply(object, needUpdate, NA))

print.varComb <-
  function(x, ...)
{
  cat("Combination of:\n")
  lapply(x, print)
  invisible(x)
}

print.summary.varComb <-
  function(x, ...)
{
  cat(attr(x, "structName"),"\n")
  lapply(x, print, FALSE)
  invisible(x)
}

summary.varComb <-
  function(object, structName = "Combination of variance functions:", ...)
{
  object[] <- lapply(object, summary)
  attr(object, "structName") <- structName
  class(object) <- c("summary.varComb", class(object))
  object
}

update.varComb <-
  function(object, data, ...)
{
  object[] <- lapply(object, update, data)
  object
}