nlrob               package:robustbase               R Documentation

_R_o_b_u_s_t _F_i_t_t_i_n_g _o_f _N_o_n_l_i_n_e_a_r _R_e_g_r_e_s_s_i_o_n _M_o_d_e_l_s

_D_e_s_c_r_i_p_t_i_o_n:

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

_U_s_a_g_e:

     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, ...)
     ## S3 method for class 'nlrob':
     predict(object, newdata, ...)

_A_r_g_u_m_e_n_t_s:

 formula: a nonlinear 'formula' including variables and parameters of
          the model, such as 'y ~ f(x, alpha)' (cf. 'nls'). (For some
          checks: if f(.) is linear, then we need parentheses, e.g., 'y
          ~ (a + b * x)'.) Do not use as 'w' as variable or parameter
          name! 

    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, e[i] = y[i] - f(xreg[i],
          alpha's), and 'w' are the robust weights from 'resid *
          weights'.

na.action: a function which indicates what should happen when the data
          contain 'NA's.  The default action is for the procedure to
          fail.  If NAs are present, use 'na.exclude' to have residuals
          with 'length == nrow(data) == length(w)', where 'w' are the
          weights used in the iterative robust loop.  This is better if
          the explanatory variables in 'formula' are time series (and
          so the NA location is important). For this reason, 'na.omit',
          which leads to omission of cases with missing values on any
          required variable, is not suitable here since the residuals
          length is different from 'nrow(data) == length(w)'. 

     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 psi
          function. (see also 'rlm').

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 values are printed at the
          conclusion of each 'nls' iteration. When the '"plinear"'
          algorithm is used, the conditional estimates of the linear
          parameters are printed after the nonlinear parameters.

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

  object: an R object of class '"nlrob"', typically resulting from
          'nlrob(..)'.

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

_V_a_l_u_e:

     An object of S3 class '"nlrob"', also inheriting from class "nls",
     (see 'nls').

     There methods (at least) for the generic accessor functions
     'summary', 'coefficients', 'fitted.values' and 'residuals'.

     It is a list with components 

  FIXME : ???

_N_o_t_e:

     This function used to be named 'rnls' and has been in package
     'sfsmisc', but will be deprecated and dropped there, eventually.

_A_u_t_h_o_r(_s):

     Andreas Ruckstuhl (inspired by 'rlm'() and 'nls'()), in July 1994
     for S-plus.
      Christian Sangiorgio did the update to R and corrected some
     errors, from June 2002 to January 2005, and Andreas contributed
     slight changes and the first methods in August 2005.
      Help page, testing, more cleanup, methods: Martin Maechler.

_S_e_e _A_l_s_o:

     'nls', 'rlm'.

_E_x_a_m_p_l_e_s:

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

     ## note that selfstarting models don't work yet 

     ##--- 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
     Rm3DNase1 <- nlrob(density ~ Asym/(1 + exp(( xmid - log(conc) )/scal ) ),
                        data = DNase1, trace = TRUE,
                        start = list( Asym = 3, xmid = 0, scal = 1 ))
     summary( Rm3DNase1 )

     ##--- 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 
     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
     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
     Rm3DNase1 <- nlrob(density ~  Asym/(1 + exp(( xmid - log(conc) )/scal ) ),
                        data = DNase1, trace = TRUE,
                        start = list( Asym = 3, xmid = 0, scal = 1 ))
     Rm3DNase1 ## summary() is not yet there {FIXME}

     ## 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(DNase1$conc), max(DNase1$conc), length = 100)
     nDat <- data.frame(conc = h.x)

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

     plot(density ~ conc, data=DNase1, log="x",
          main = deparse(Rm3DNase1$call$formula))
     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)

