R/varFunc.R

Defines functions varPower summary.varIdent recalc.varIdent needUpdate.varIdent Initialize.varIdent coef.varIdent varIdent summary.varFixed print.summary.varFixed Initialize.varFixed coef.varFixed varFixed update.varFunc summary.varFunc recalc.varFunc print.varFunc print.summary.varFunc logLik.varFunc Initialize.varFunc getGroups.varFunc getCovariate.varFunc formula.varFunc coef.varFunc varWeights.varFunc

Documented in coef.varFixed coef.varFunc coef.varIdent getCovariate.varFunc getGroups.varFunc Initialize.varFixed Initialize.varFunc Initialize.varIdent logLik.varFunc needUpdate.varIdent print.varFunc recalc.varFunc recalc.varIdent summary.varFixed summary.varFunc summary.varIdent update.varFunc varFixed varIdent varPower varWeights.varFunc

###              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
}

Try the nlme package in your browser

Any scripts or data that you put into this service are public.

nlme documentation built on Aug. 9, 2022, 5:10 p.m.