sfsmisc (version 0.95-2)

rnls: Robust Nonlinear Regression

Description

Compute robustified nonlinear regressions, i.e., robustly estimate parameters in the nonlinear model. The fitting is done by iterated reweighted least squares (IWLS) as in rlm() of package MASS. Additionally, see nls.

Usage

rnls(formula, data, start, weights = NULL,
     na.action = na.fail, psi = MASS::psi.huber,
     test.vec = c("resid", "coef", "w"),
     maxit = 20, acc = 1e-06, algorithm = "default",
     control, trace = FALSE, ...)
## S3 method for class 'rnls':
fitted(object, \dots)
## S3 method for class 'rnls':
predict(object, newdata, \dots)
## S3 method for class 'rnls':
residuals(object, \dots)
## S3 method for class 'rnls':
summary(object, \dots)
## S3 method for class 'rnls':
print(x, \dots)

Arguments

formula

a nonlinear formula including variables and parameters and which results in the residuals of the model, such as y ~ f(x, alpha). (For some checks: if $f(.)$ is linear, then we need parentheses, e.g., y ~ (a + b * x)

data

a data frame in which to do the computations.

start

a named numeric vector of starting parameters estimates.

weights

an optional weighting vector (for intrinsic weights, not the weights w used in the iterative (robust) fit). I.e., sum(w * e^2) is minimized with e = residuals, $e_i = y_{i} - f(xreg_{i}, alpha's)$, an

na.action

a function which indicates what should happen when the data contain NAs. The default action is for the procedure to fail. If NAs are present, use na.exclude to have residuals with length == nrow(data) == lengt

psi

a function (possibly by name) of the form g(x, ..., deriv) that for deriv=0 returns $\psi(x)/x$ and for deriv=1 returns $\psi'(x)$. Tuning constants will be passed in via .... Psi func

test.vec

character string specifying the convergence criterion. The relative change is tested for residuals with a value of "resid" (the default), for coefficients with "coef", and for weights with "w".

maxit

maximum number of iterations in the robust loop.

acc

convergence tolerance for the robust fit.

algorithm

character string specifying the algorithm to use for nls. The default algorithm is a Gauss-Newton algorithm. The other alternative is "plinear", the Golub-Pereyra algorithm for partially linear least-squares models.

control

an optional list of control settings for nls. See nls.control for the names of the settable control values and their effect.

trace

logical value indicating if a trace of the nls iteration progress should be printed. Default is FALSE. If TRUE, in each robust iteration, the residual sum-of-squares and the parameter v

...

potentially arguments to be passed to the psi function (see above).

object, x

an Robject of class "rnls", as returned from rnls.

newdata

a data frame or list containing variables identically named as in data, with (new) values at which to evaluate (predict) the fitted model.

Value

• An object of S3 class "rnls", inheriting from class "nls", (see nls). It is a list with components
• comp1
  Description of 'comp1'
• comp2
  Description of 'comp2'
• ...

See Also

nls, rlm.

Examples

DNase1 <- DNase[ DNase$Run == 1, ]

## note that selfstarting models doesn't work yes % <<< FIXME !!!

##--- without conditional linearity ---

## classical
fm3DNase1 <- nls( density ~ Asym/(1 + exp(( xmid - log(conc) )/scal ) ),
                  data = DNase1,
                  start = list( Asym = 3, xmid = 0, scal = 1 ),
                  trace = TRUE )
summary( fm3DNase1 )

## robust
frm3DNase1 <- rnls(density ~ Asym/(1 + exp(( xmid - log(conc) )/scal ) ),
                  data = DNase1, trace = TRUE,
                  start = list( Asym = 3, xmid = 0, scal = 1 ))
## FAILS (summary.nls): summary( frm3DNase1 )%% <<< FIXME

##--- using conditional linearity ---

## classical
fm2DNase1 <- nls( density ~ 1/(1 + exp(( xmid - log(conc) )/scal ) ),
                  data = DNase1,
                  start = c( xmid = 0, scal = 1 ),
                  alg = "plinear", trace = TRUE )
summary( fm2DNase1 )

## robust
if(FALSE) { # currently fails %% FIXME error in nls's nlsModel.plinear()
frm2DNase1 <- rnls(density ~ 1/(1 + exp(( xmid - log(conc) )/scal ) ),
                  data = DNase1, start = c( xmid = 0, scal = 1 ),
                  alg = "plinear", trace = TRUE )
summary( frm2DNase1 )
} # not yet

### -- new examples
DNase1[10,"density"] <- 2*DNase1[10,"density"]

fm3DNase1 <- nls(density ~  Asym/(1 + exp(( xmid - log(conc) )/scal ) ),
                       data = DNase1, trace = TRUE,
                       start = list( Asym = 3, xmid = 0, scal = 1 ))

## robust
frm3DNase1 <- rnls(density ~  Asym/(1 + exp(( xmid - log(conc) )/scal ) ),
                   data = DNase1, trace = TRUE,
                   start = list( Asym = 3, xmid = 0, scal = 1 ))
frm3DNase1coef(frm3DNase1)

## predict() {and plot}:
h.x <- lseq(min(DNase1$conc), max(DNase1$conc), length = 100)
nDat <- data.frame(conc = h.x)

h.p  <- predict(fm3DNase1, newdata = nDat)# classical
h.rp <- predict(frm3DNase1, newdata= nDat)# robust

plot(density ~ conc, data=DNase1, log="x",
     main = deparse(frm3DNase1$call$formula))
lines(h.x, h.p,  col="blue")
lines(h.x, h.rp, col="magenta")
legend("topleft", c("classical nls()", "robust    rnls()"),
       lwd = 1, col= c("blue", "magenta"))