robustbase (version 0.7-8)

nlrob: Robust Fitting of Nonlinear Regression Models

Description

nlrob fits a nonlinear regression model by robust M-estimators, using iterated reweighted least squares (IWLS).

Usage

nlrob(formula, data, start, weights = NULL, na.action = na.fail,
      psi = psi.huber, test.vec = c("resid", "coef", "w"), maxit = 20,
      acc = 1e-06, algorithm = "default", control = nls.control(),
      trace = FALSE, ...)

## S3 method for class 'nlrob':
fitted(object, ...)
## S3 method for class 'nlrob':
residuals(object, type = , ...)## S3 method for class 'nlrob':
predict(object, newdata, ...)

Arguments

formula

a nonlinear formula including variables and parameters of the model, such as y ~ f(x, theta) (cf. nls). (For some checks: if $f(.)$ is lin

data

an optional data frame containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which nlrob is called.

start

a named numeric vector of starting parameters estimates.

weights

an optional vector of weights to be used in the fitting process (for intrinsic weights, not the weights w used in the iterative (robust) fit). I.e., sum(w * e^2) is minimized with e = residuals, $

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) == le

psi

a function (possibly by name) of the form g(x, 'tuning constant(s)', deriv) that for deriv=0 returns psi(x)/x and for deriv=1 returns psi'(x). Tuning constants will be passed in via one or several arguments, depending on the ps

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 paramete

...

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

object

an Robject of class "nlrob", typically resulting from nlrob(..).

type

a string specifying the type of residuals desired. Currently, "response" and "working" are supported.

newdata

a data frame (or list) with the same names as the original data, see e.g., predict.nls.

Value

• nlrob() returns an object of S3 class "nlrob", also inheriting from class "nls", (see nls). It is a list with components
• FIXME
  ???
• There are methods (at least) for the generic accessor functions summary, coefficients, fitted.values and residuals. residuals(.), by default type = "response", returns the residuals $e_i$, defined above as $e_i = Y_i - f_(x_i, \theta)$. These differ from the standardized or weighted residuals which, e.g., are assumed to be normally distributed, and a version of which is returned in working.residuals component.

See Also

nls, rlm.

Examples

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

## note that selfstarting models don't work yet % <<< FIXME !!!

##--- without conditional linearity ---

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

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

##--- 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 <- nlrob(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 -- "moderate outlier":
DN2 <- DNase1
DN2[10,"density"] <- 2*DN2[10,"density"]

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

## robust
Rm3DN2 <- nlrob(density ~  Asym/(1 + exp(( xmid - log(conc) )/scal ) ),
                data = DN2, trace = TRUE,
                start = list( Asym = 3, xmid = 0, scal = 1 ))
Rm3DN2
summary(Rm3DN2)

## utility function sfsmisc::lseq() :
lseq <- function (from, to, length)
  2^seq(log2(from), log2(to), length.out = length)
## predict() {and plot}:
h.x <- lseq(min(DN2$conc), max(DN2$conc), length = 100)
nDat <- data.frame(conc = h.x)

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

plot(density ~ conc, data=DN2, log="x",
     main = format(formula(Rm3DN2)))
lines(h.x, h.p,  col="blue")
lines(h.x, h.rp, col="magenta")
legend("topleft", c("classical nls()", "robust nlrob()"),
       lwd = 1, col= c("blue", "magenta"), inset = 0.05)