# 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)
| 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:
| ||||||
| 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. |
# 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)