Data engineering and BI courses are free!
nlrob fits a nonlinear regression model by robust methods.
Per default, by an M-estimator, using iterated reweighted least
squares (called nlrob(formula, data, start, lower, upper,
weights = NULL, na.action = na.fail,
method = c("M", "MM", "tau", "CM", "mtl"),
psi = .Mwgt.psi1("huber", cc=1.345), scale = NULL,
test.vec = c("resid", "coef", "w"), maxit = 20,
tol = 1e-06, acc, algorithm = "default", doCov = FALSE,
control = if(method == "M") nls.control() else
nlrob.control(method, optArgs = list(trace=trace), ...),
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, ...)
data, the variables are taken
from environment(formula), typically the environment from
which nlrob is called.method = "M".start,
lower or upper. For (the default) method = "M",
if the bounds are unspecified w used in the
iterative (robust) fit). I.e.,
sum(w * e^2) is minimized with e = residuals,
$e_i 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"M", for historical and back-compatibility
reasons. For the other methods, primarily see
nlrob.algorithmsg(x, 'tuning
constant(s)', deriv) that for deriv=0 returns
$\psi(x)/x$ and for deriv=1 returns
$\psi'(x)$. Note that tuning constants can not
be pNULL (default), a positive number
specifying a scale kept fixed during the iterations (and
returned as Scale component)."resid" (the default), for coefficients with
"coef", and for weights with "w".tol, now deprecated.nls, see there, only when method = "M". The
default algorithm is a Gauss-Newton algorithm.nlrob() should compute the
asymptotic variance-covariance matrix (see vcov)
already. This used to be hard-wired to TRUE; however, the
default has nls iteration progress should be printed. Default is
FALSE.
If TRUE, in each robust iteration, the residual
sum-of-squares and the parameter v"nlrob", typically
resulting from nlrob(..).nlrob: only when method is not "M",
optional arguments for nlrob.control;
for other functions: potentially optional arguments passed to the"response" and "working" are supported.data, see e.g., predict.nls.nlrob() returns an object of S3 class "nlrob",
for method = "M" also inheriting from class "nls", (see
nls). It is a list with several components; they are not documented yet,
as some of them will probably change.
Instead, rather use summary(), coefficients() (aka coef())
fitted.values(), residuals(), sigma() and
vcov(), the latter for the variance-covariance matrix of
the estimated parameters, as returned by coef(), i.e., not
including the variance of the errors.
For nlrob() results, estimethod() returns the
method
argument used.
residuals(.), by default type = "response", returns
the residuals $e_i$, defined above as
$e_i = Y_i - f_(x_i, \hat\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.
method = "M", iterated reweighted least squares
(nls(*,
weights= .) where weights $w_i$ are proportional to
$\psi(r_i/ \hat{\sigma})$. All other methods minimize differently, and work without
nls. See nlrob.algorithms for details.
nls, rlm.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
frm2DNase1 <- nlrob(density ~ 1/(1 + exp(( xmid - log(conc) )/scal ) ),
data = DNase1, start = c( xmid = 0, scal = 1 ),
alg = "plinear", trace = TRUE )
summary( frm2DNase1 )
## Confidence for linear parameter is quite smaller than "Asym" above
c1 <- coef(summary(RmN1))
c2 <- coef(summary(frm2DNase1))
rownames(c2)[rownames(c2) == ".lin"] <- "Asym"
stopifnot(all.equal(c1[,1:2], c2[rownames(c1), 1:2], tol = 0.09)) # 0.07315
### -- 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) # -> robustness weight of obs. 10 ~= 0.037
confint(Rm3DN2, method = "Wald")
## 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)
## See ?nlrob.algorithms for examples
DNase1 <- DNase[DNase$Run == 1,]
form <- density ~ Asym/(1 + exp(( xmid -log(conc) )/scal ))
gMM <- nlrob(form, data = DNase1, method = "MM",
lower = c(Asym = 0, xmid = 0, scal = 0),
upper = 3, trace = TRUE)
## "CM" (and "mtl") additionally need bounds for "sigma" :
gCM <- nlrob(form, data = DNase1, method = "CM",
lower = c(Asym = 0, xmid = 0, scal = 0, sigma = 0),
upper = c(3,3,3, sigma = 0.8))Run the code above in your browser using DataLab