2

I have been trying to fit non linear quantile regression model with nlrq() (from quantreg package). But what I found that the function is very sensitive to starting values. So I tried to use self starter SSexp() with it. What I read that self starter functions made for nls() also works with it. But all the time I am getting the following error message :

"Error in getInitial(formula, mf) : unused argument (mf)"

my input code was:

fit12 <- nlrq(visit.rate ~ SSexp(het.total, y0, b), 
                               data = flower_mat,tau = 0.75)

I have kind of vague idea that getInitial is somehow connected to create self starter function. But frankly speaking I have no idea where exactly the problem is. Can anybody help me regarding this?

Additional Information

# R version 3.1.2 -- "Pumpkin Helmet"
# quantreg version: 5.05
# nlrwr version : 1.1-0

# toy dataset and code
x <- c(0.3,0.6,0.9,1,1.5,2, 2.1, 2.5,3, 3,5,10,11,12,14,13,17,21,23,27,30, 50)
y <- c(0,0.1,0.3,0.4,0.6,0.2,0.27,0.2,0.25,0.4,0.15,0.05,0.25,0.2,0.3,0.35,0.1,0.15,0.1,0.14,0,0)
dat <- data.frame(x,y,stringsAsFactors=FALSE)
plot(y~x)
library(quantreg)
library(nlrwr)
fit.1 <- nlrq(y ~ SSexp(x, y0, b), 
              data = dat,tau = 0.75)
# above mentioned error message
with(dat, SSexp(x, 0.2, 4))
# this is working
getInitial((y ~ SSexp(x, y0, b)), dat)
# showing following error:
# Error in getInitial((y ~ SSexp(x, y0, b)), data = dat) : unused argument (dat)
lmo
  • 37,904
  • 9
  • 56
  • 69

2 Answers2

4

Both getInitial() and SSexp() from nlrwr are causing problems here. So, I would not use nlrwr in your case.

Plus, it attaches a lot of packages to the search path, which is generally considered as a bad practice.

  • nlrwr::getInitial() accepts only one argument, while stats::getInitial() accepts a varying number of arguments with two at least. So, the nlrwr version is very prone to break several things in R!

  • nlrwr::SSexp() is a poor self-start function. It does not calculate gradient numerically, it does not provide a pnames attribute, and moreover, it fails to calculate initial estimates for parameters in case the dependent variable contains values <= 0, as in your toy example.

I would rather use this form of the SSexp() function (take care: arguments are different and more in the spirit of the other self-start functions):

SSexp <- structure(function (input, A, rc) {
    .expr2 <- exp(rc * input)
    .value <- A * .expr2
    .actualArgs <- as.list(match.call()[c("A", "rc")])
    if (all(unlist(lapply(.actualArgs, is.name)))) {
        .grad <- array(0, c(length(.value), 2), list(NULL, c("A", "rc")))
        .grad[, "A"] <- .expr2
        .grad[, "rc"] <- A * (.expr2 * input)
        dimnames(.grad) <- list(NULL, .actualArgs)
        attr(.value, "gradient") <- .grad
    }
    .value
}
, initial = function (mCall, data, LHS) {
    xy <- data.frame(sortedXyData(mCall[["input"]], LHS, data))
    if (nrow(xy) < 3)
        stop("Too few distinct input values to fit an exponential")
    xy$logy <- log(xy$y)
    ## Keep only finite cases (if there are y <= 0)
    xyfinite <- xy[is.finite(xy$logy), ]
    if (nrow(xyfinite) < 2)
        stop("Too few distinct LHS values > 0 to fit an exponential")
    res <- lsfit(xyfinite$x, xyfinite$logy)$coef
    value <- c(exp(res[1]), res[2])
    setNames(value, mCall[c("A", "rc")])
}
, pnames = c("A", "rc"), class = "selfStart")

dat <- data.frame(
    x = c(0.3, 0.6, 0.9, 1, 1.5, 2, 2.1, 2.5, 3, 3, 5, 10, 11, 12, 14, 13, 17,
          21, 23, 27, 30, 50),
    y = c(0, 0.1, 0.3, 0.4, 0.6, 0.2, 0.27, 0.2, 0.25, 0.4, 0.15, 0.05, 0.25,
          0.2, 0.3, 0.35, 0.1, 0.15, 0.1, 0.14, 0,0))
plot(y ~ x, data = dat)
library(quantreg)
dat.nlrq <- nlrq(y ~ SSexp(x, A, rc), data = dat, tau = 0.75)
summary(dat.nlrq)
lines(0:50, predict(dat.nlrq, newdata = list(x = 0:50)))
Carl Witthoft
  • 20,573
  • 9
  • 43
  • 73
0

I suspect the problem is that one of the dependency packages loaded when you load nlrwr is dlc, which has its own version of getInitial.

Try stats::getInitial(etc) -- tho' I will warn you that when I tried it with your data and code, I get "Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : NA/NaN/Inf in 'y'"

Carl Witthoft
  • 20,573
  • 9
  • 43
  • 73