selfStart
Construct Self-starting Nonlinear Models

Description

Constructs an objective function for nls that can compute its own starting parameters.

Usage

# Generic function
selfStart(model, initial, parameters, template)
## Default S3 method:
selfStart(model, initial, parameters, template)
## S3 method for class 'formula':
selfStart(model, initial, parameters, template)

Arguments

model either a function of data and parameters defining a nonlinear model, or a formula whose right hand side defines a nonlinear model.
initial a function taking three arguments:
mCall a call to the function model with all argument names spelled out completely, converted to a list with as.list.
data a named list of data vectors in which to interpret the expressions in mCall and LHS.
LHS the expression from the left hand side of the model formula in the call to nls.
This function should return values for the parameters in model that nls can use as starting values when it tries to minimize sum((LHS - model(data,parameters)^2)).
parameters a character vector specifying the arguments to the function model or the variables on the right hand side of the formula model for which initial estimates should be calculated.
template an optional prototype for the calling sequence of the returned object, passed as the function.arg argument to the deriv function. By default, a template is generated, with the covariates in model coming first and the parameters in model coming last in the calling sequence. This argument is not used in the default method.
Value
returns a function of class selfStart. The returned model function has the attributes initial and pnames, copied from the arguments initial and parameters, respectively.
Note
nls can fail to find an answer if it is given poor starting values for the parameters. You can use selfStart to attach a method for computing starting parameters to a model function so nls has a greater chance of converging. The initial estimates often come from linearizing the model or fitting a simplified model to different parts of the data.
See Also
nls, deriv
Examples
# broken-stick regression: two lines, one in the left part of the
# range of x with slope1 and the other in the right part with slope2.
# The lines meet at the point (breakPointX,breakPointY).
# Compute initial estimates by fitting lines to the right- and
# left-most 4 points.
SSbsr <- selfStart(
    model = function(x, breakPointX, breakPointY, slope1, slope2)
            {
                ifelse(x > breakPointX,
                    (x - breakPointX) * slope2 + breakPointY,
                    (x - breakPointX) * slope1 + breakPointY)
            },
    initial = function(mCall, data, LHS)
            {
                xy <- sortedXyData(mCall[["x"]], LHS, data)
                n <- nrow(xy)
                if (n < 7) {
                    stop("need at least 7 points to compute initial estimates")
                }
                first <- seq(from=1, length.out=4)
                last <- seq(length.out=4, to=n)
                fit1 <- coefficients(lm.fit(cbind(1,xy[first,"x"]), xy[first,"y"]))
                fit2 <- coefficients(lm.fit(cbind(1,xy[last,"x"]), xy[last,"y"]))
                slope1 <- fit1[[2]]
                slope2 <- fit2[[2]]
                bpXY <- solve(cbind(-c(slope1,slope2), 1),
                            as.matrix(c(fit1[[1]], fit2[[1]])))
                breakPointX <- bpXY[1]
                breakPointY <- bpXY[2]
                # the names of the output list must match the
                # names the user gave in the call
                structure(list(breakPointX, breakPointY, slope1, slope2),
                          names=as.character(mCall[c("breakPointX",
                                "breakPointY", "slope1", "slope2")]))
            },
    parameters = c("breakPointX", "breakPointY", "slope1", "slope2"))

bsrData <- data.frame( x = c(8.4, 2.8, 6.9, 0, 4.1, 0.3, 7.2, 0.6, 8.1, 8.8, 5.9, 5.3), y = c(0.62, 0.39, 0.52, 0.33, 0.45, 0.36, 0.56, 0.33, 0.63, 0.65, 0.49, 0.48))

nls(y ~ SSbsr(x, bpX, bpY, s1, s2), data=bsrData)

## self-starting logistic model SSlogis_1 <- selfStart(~ Asym/(1 + exp((xmid - x)/scal)), initial = function(mCall, data, LHS) { xy <- sortedXyData(mCall[["x"]], LHS, data) if(nrow(xy) < 4) { stop("Too few distinct x values to fit a logistic") } z <- xy[["y"]] if (min(z) <= 0) { z <- z + 0.05 * max(z) } # avoid zeroes z <- z/(1.05 * max(z)) # scale to within unit height xy[["z"]] <- log(z/(1 - z)) # logit transformation aux <- coef(lm(x ~ z, xy)) pars <- as.vector(coef(nls(y ~ 1/(1 + exp((xmid - x)/scal)), start = list(xmid = aux[1], scal = aux[2]), data = xy, algorithm = "plinear"))) value <- list(pars[3], pars[1], pars[2]) names(value) <- mCall[c("Asym", "xmid", "scal")] value }, parameters = c("Asym", "xmid", "scal")) logisData <- data.frame( x = 1:10, y = c(0.48, 0.95, 0.98, 1.07, 1.15, 1.85, 3.32, 3.76, 4.08, 4.22)) getInitial(y ~ SSlogis_1(x, A, m, s), data=logisData)

Package stats version 6.0.0-69
Package Index