R.version.string
package
sapply(package, packageVersion)[1] "R version 3.6.1 (2019-07-05)"
[1] "nlme"
$nlme
[1] 3 1 142list(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()
controlvals <- glsControl()
if (!missing(control))
controlvals[names(control)] <- control
if (!inherits(model, "formula") || length(model) != 3L) {
stop("\nmodel must be a formula of the form \"resp ~ pred\"")
}
method <- match.arg(method)
REML <- method == "REML"
groups <- if (!is.null(correlation))
getGroupsFormula(correlation)
glsSt <- glsStruct(corStruct = correlation, varStruct = varFunc(weights))
model <- terms(model, data = data)
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)
if (!is.null(groups)) {
groups <- eval(substitute(~1 | GR, list(GR = groups[[2L]])))
grps <- getGroups(dataMod, groups, level = length(getGroupsFormula(groups,
asList = TRUE)))
ord <- order(grps)
grps <- grps[ord]
dataMod <- dataMod[ord, , drop = FALSE]
revOrder <- match(origOrder, row.names(dataMod))
}
else grps <- NULL
X <- model.frame(model, dataMod)
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)
parAssign <- attr(X, "assign")
fTerms <- terms(as.formula(model), data = data)
namTerms <- attr(fTerms, "term.labels")
if (attr(fTerms, "intercept") > 0) {
namTerms <- c("(Intercept)", namTerms)
}
namTerms <- factor(parAssign, labels = namTerms)
parAssign <- split(order(parAssign), namTerms)
fixedSigma <- (controlvals$sigma > 0)
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,
sigma = controlvals$sigma, fixedSigma = fixedSigma)
glsEstControl <- controlvals["singular.ok"]
glsSt <- Initialize(glsSt, dataMod, glsEstControl)
parMap <- attr(glsSt, "pmap")
numIter <- numIter0 <- 0L
repeat {
oldPars <- c(attr(glsSt, "glsFit")[["beta"]], coef(glsSt))
if (length(coef(glsSt))) {
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)
if (!needUpdate(glsSt)) {
if (optRes$convergence)
stop(optRes$message)
break
}
numIter <- numIter + 1L
glsSt <- update(glsSt, dataMod)
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")
}
}
glsFit <- attr(glsSt, "glsFit")
namBeta <- names(glsFit$beta)
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 <- fitted(glsSt)
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
apVar <- if (controlvals$apVar)
glsApVar(glsSt, glsFit$sigma, .relStep = controlvals[[".relStep"]],
minAbsPar = controlvals[["minAbsParApVar"]], natural = controlvals[["natural"]])
else "Approximate variance-covariance matrix not available"
dims <- attr(glsSt, "conLin")[["dims"]]
dims[["p"]] <- p
attr(glsSt, "conLin") <- NULL
attr(glsSt, "glsFit") <- NULL
attr(glsSt, "fixedSigma") <- fixedSigma
grpDta <- inherits(data, "groupedData")
structure(class = "gls", list(modelStruct = glsSt, dims = dims,
contrasts = contr, coefficients = glsFit[["beta"]], varBeta = varBeta,
sigma = if (fixedSigma) controlvals$sigma else glsFit$sigma,
apVar = apVar, logLik = glsFit$logLik, numIter = if (needUpdate(glsSt)) numIter else numIter0,
groups = grps, call = Call, method = method, fitted = Fitted,
residuals = Resid, parAssign = parAssign, na.action = attr(dataMod,
"na.action")), namBetaFull = colnames(X), units = if (grpDta)
attr(data, "units"), labels = if (grpDta)
attr(data, "labels"))
})| gls | R Documentation |
This function fits a linear model using generalized least squares. The errors are allowed to be correlated and/or have unequal variances.
gls(model, data, correlation, weights, subset, method, na.action,
control, verbose)
## S3 method for class 'gls'
update(object, model., ..., evaluate = TRUE)
object
|
an object inheriting from class |
model
|
a two-sided linear formula object describing the model, with the response on the left of a |
model.
|
Changes to the model – see |
data
|
an optional data frame containing the variables named in |
correlation
|
an optional |
weights
|
an optional |
subset
|
an optional expression indicating which subset of the rows of |
method
|
a character string. If |
na.action
|
a function that indicates what should happen when the data contain |
control
|
a list of control values for the estimation algorithm to replace the default values returned by the function |
verbose
|
an optional logical value. If |
…
|
some methods for this generic require additional arguments. None are used in this method. |
evaluate
|
If |
an object of class “gls” representing the linear model fit. Generic functions such as print, plot, and summary have methods to show the results of the fit. See glsObject for the components of the fit. The functions resid, coef and fitted, can be used to extract some of its components.
Josﾃｩ Pinheiro and Douglas Bates bates@stat.wisc.edu
The different correlation structures available for the correlation argument are described in Box, G.E.P., Jenkins, G.M., and Reinsel G.C. (1994), Littel, R.C., Milliken, G.A., Stroup, W.W., and Wolfinger, R.D. (1996), and Venables, W.N. and Ripley, B.D. (2002). The use of variance functions for linear and nonlinear models is presented in detail in Carroll, R.J. and Ruppert, D. (1988) and Davidian, M. and Giltinan, D.M. (1995).
Box, G.E.P., Jenkins, G.M., and Reinsel G.C. (1994) “Time Series Analysis: Forecasting and Control”, 3rd Edition, Holden-Day.
Carroll, R.J. and Ruppert, D. (1988) “Transformation and Weighting in Regression”, Chapman and Hall.
Davidian, M. and Giltinan, D.M. (1995) “Nonlinear Mixed Effects Models for Repeated Measurement Data”, Chapman and Hall.
Littel, R.C., Milliken, G.A., Stroup, W.W., and Wolfinger, R.D. (1996) “SAS Systems for Mixed Models”, SAS Institute.
Pinheiro, J.C., and Bates, D.M. (2000) “Mixed-Effects Models in S and S-PLUS”, Springer, esp. pp. 100, 461.
Venables, W.N. and Ripley, B.D. (2002) “Modern Applied Statistics with S”, 4th Edition, Springer-Verlag.
corClasses, glsControl, glsObject, glsStruct, plot.gls, predict.gls, qqnorm.gls, residuals.gls, summary.gls, varClasses, varFunc
# AR(1) errors within each Mare
fm1 <- gls(follicles ~ sin(2*pi*Time) + cos(2*pi*Time), Ovary,
correlation = corAR1(form = ~ 1 | Mare))
# variance increases as a power of the absolute fitted values
fm2 <- update(fm1, weights = varPower())
list(function (corStruct = NULL, varStruct = NULL)
{
val <- list(corStruct = corStruct, varStruct = varStruct)
val <- val[!vapply(val, is.null, NA)]
class(val) <- c("glsStruct", "modelStruct")
val
})| glsStruct | R Documentation |
A generalized least squares structure is a list of model components representing different sets of parameters in the linear model. A glsStruct may contain corStruct and varFunc objects. NULL arguments are not included in the glsStruct list.
glsStruct(corStruct, varStruct)
corStruct
|
an optional |
varStruct
|
an optional |
a list of model variance-covariance components determining the parameters to be estimated for the associated linear model.
Josﾃｩ Pinheiro and Douglas Bates bates@stat.wisc.edu
corClasses, gls, residuals.glsStruct, varFunc
gls1 <- glsStruct(corAR1(), varPower())
list(function (maxIter = 50L, msMaxIter = 200L, tolerance = 0.000001,
msTol = 0.0000001, 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)
{
if (is.null(sigma))
sigma <- 0
else {
if (!is.finite(sigma) || length(sigma) != 1 || sigma <
0)
stop("Within-group std. dev. must be a positive numeric value")
}
list(maxIter = maxIter, 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)
})| glsControl | R Documentation |
The values supplied in the function call replace the defaults and a list with all possible arguments is returned. The returned list is used as the control argument to the gls function.
glsControl(maxIter, msMaxIter, tolerance, msTol, msVerbose,
singular.ok, returnObject = FALSE, apVar, .relStep,
opt = c("nlminb", "optim"), optimMethod,
minAbsParApVar, natural, sigma = NULL)
maxIter
|
maximum number of iterations for the |
msMaxIter
|
maximum number of iterations for the optimization step inside the |
tolerance
|
tolerance for the convergence criterion in the |
msTol
|
tolerance for the convergence criterion of the first outer iteration when |
msVerbose
|
a logical value passed as the |
singular.ok
|
a logical value indicating whether non-estimable coefficients (resulting from linear dependencies among the columns of the regression matrix) should be allowed. Default is |
returnObject
|
a logical value indicating whether the fitted object should be returned when the maximum number of iterations is reached without convergence of the algorithm. Default is |
apVar
|
a logical value indicating whether the approximate covariance matrix of the variance-covariance parameters should be calculated. Default is |
.relStep
|
relative step for numerical derivatives calculations. Default is |
opt
|
the optimizer to be used, either |
optimMethod
|
character - the optimization method to be used with the |
minAbsParApVar
|
numeric value - minimum absolute parameter value in the approximate variance calculation. The default is |
natural
|
logical. Should the natural parameterization be used for the approximate variance calculations? Default is |
sigma
|
optionally a positive number to fix the residual error at. If |
a list with components for each of the possible arguments.
Josﾃｩ Pinheiro and Douglas Bates bates@stat.wisc.edu; the sigma option: Siem Heisterkamp and Bert van Willigen.
gls
# decrease the maximum number iterations in the optimization call and # request that information on the evolution of the ms iterations be printed glsControl(msMaxIter = 20, msVerbose = TRUE)
list(function (glsSt, sigma, conLin = attr(glsSt, "conLin"),
.relStep = .Machine$double.eps^(1/3), minAbsPar = 0, natural = TRUE)
{
fixedSigma <- attr(glsSt, "fixedSigma")
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
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)) {
val <- solve(-val)
nP <- names(Pars)
dimnames(val) <- list(nP, nP)
attr(val, "Pars") <- Pars
attr(val, "natural") <- natural
val
}
else {
"Non-positive definite approximate variance-covariance"
}
}
})| gls-internal | R Documentation |
These are functions used by gls to call its compiled C code. They are exported to allow experimentation with modified versions.
glsApVar(glsSt, sigma, conLin = attr(glsSt, "conLin"),
.relStep = .Machine$double.eps^(1/3), minAbsPar = 0, natural = TRUE)
glsEstimate(object, conLin = attr(object, "conLin"),
control = list(singular.ok = FALSE))
glsSt, object
|
An object inheriting from class |
sigma
|
the estimated residual standard error: see |
conLin
|
A ‘condensed linear model’: see |
.relStep, minAbsPar, natural
|
Control values: see |
control
|
The relevant part of a |
list(function (pars, fun, ..., .relStep = .Machine$double.eps^(1/3),
minAbsPar = 0)
{
pars <- as.numeric(pars)
npar <- length(pars)
incr <- pmax(abs(pars), minAbsPar) * .relStep
baseInd <- diag(npar)
frac <- c(1, incr, incr^2)
cols <- list(0, baseInd, -baseInd)
for (i in seq_along(pars)[-npar]) {
cols <- c(cols, list(baseInd[, i] + baseInd[, -(1:i)]))
frac <- c(frac, incr[i] * incr[-(1:i)])
}
indMat <- do.call("cbind", cols)
shifted <- pars + incr * indMat
indMat <- t(indMat)
Xcols <- list(1, indMat, indMat^2)
for (i in seq_along(pars)[-npar]) {
Xcols <- c(Xcols, list(indMat[, i] * indMat[, -(1:i)]))
}
coefs <- solve(do.call("cbind", Xcols), apply(shifted, 2,
fun, ...))/frac
Hess <- diag(coefs[1 + npar + seq_along(pars)], ncol = npar)
Hess[row(Hess) > col(Hess)] <- coefs[-(1:(1 + 2 * npar))]
list(mean = coefs[1], gradient = coefs[1 + seq_along(pars)],
Hessian = (Hess + t(Hess)))
})| fdHess | R Documentation |
Evaluate an approximate Hessian and gradient of a scalar function using finite differences.
fdHess(pars, fun, ...,
.relStep = .Machine$double.eps^(1/3), minAbsPar = 0)
pars
|
the numeric values of the parameters at which to evaluate the function |
fun
|
a function depending on the parameters |
…
|
Optional additional arguments to |
.relStep
|
The relative step size to use in the finite differences. It defaults to the cube root of |
minAbsPar
|
The minimum magnitude of a parameter value that is considered non-zero. It defaults to zero meaning that any non-zero value will be considered different from zero. |
This function uses a second-order response surface design known as a “Koschal design” to determine the parameter values at which the function is evaluated.
A list with components
mean
|
the value of function |
gradient
|
an approximate gradient (of length |
Hessian
|
a matrix whose upper triangle contains an approximate Hessian. |
Josﾃｩ Pinheiro and Douglas Bates bates@stat.wisc.edu
(fdH <- fdHess(c(12.3, 2.34), function(x) x[1]*(1-exp(-0.4*x[2]))))
stopifnot(length(fdH$ mean) == 1,
length(fdH$ gradient) == 2,
identical(dim(fdH$ Hessian), c(2L, 2L)))