I performed an NLS regression using nlfb from the nlmrt package. With
theta_scaled being my initial parameter vector
> theta_scaled
[1] 0.60000 0.73624 -0.77962
and residuals_function a function which for each parameter vector would compute the residuals vector over my data, I used the call
results.nlmrt <-nlfb(start = theta_scaled, resfn = residuals_function, jacfn = NULL, trace = TRUE, lower = rep(-1, 3), upper = rep(1, 3), control = list(ndstep = 1e-5))
Results:
> results.nlmrt
nlmrt class object: x
residual sumsquares = 7929.7 on 292 observations
after 5 Jacobian and 19 function evaluations
name coeff SE tstat pval gradient JSingval
a 0.999997 0.5262 1.9 0.05838 -6.343 403.5
b 0.707415 0.07837 9.027 0 323.8 67.68
c -0.631532 0.01454 -43.44 0 894.8 9.951
I would like to plot some diagnostic plots, compute confidence intervals, predictions, etc. However, the "nlmrt" object doesn't have these fancy methods. I would like to convert the object to an nls object, but I can't use wrapnls because I used nlfb instead of nlxb for my regression. Is there any way out?
PS if you're wondering why I cannot use nlxb, the reason is that I used NLS regression to calibrateg a complex fluid dynamics code. Thus I don't have a simple analytical formula for my model. However, since I can run the code for (more or less) arbitrary inputs and parameters, I could write a residuals function and use nlfb.
EDIT as G. Grothendieck rightly notes in the comments, I should have provided a working example. As far as I'm concerned, his example is ok. However, his answer doesn't work:
#G. Grothendieck's answer
library(nlmrt)
library(nls2)
# setup and run nlfb example
shobbs.res <- function(x) {
if(length(x) != 3) stop("hobbs.res -- parameter vector n!=3")
tt <- 1:12
res <- 100.0*x[1]/(1+x[2]*10.*exp(-0.1*x[3]*tt)) - y
}
y <- c(5.308, 7.24, 9.638, 12.866, 17.069, 23.192, 31.443,
38.558, 50.156, 62.948, 75.995, 91.972)
st <- c(b1=1, b2=1, b3=1)
low <- -Inf
up <- Inf
ans1n <- nlfb(st, shobbs.res)
# get nls object
ans <- nls2(y ~ shobbs.res(c(b1, b2, b3)) + y, start = coef(ans1n), alg = "brute")
class(ans)
# my additions
# confidence intervals don't seem to work
> confint(ans)
Waiting for profiling to be done...
Error in prof$getProfile() : 'control$maxiter' absent
# apparently, there were convergence issues:
> ans
Nonlinear regression model
model: y ~ shobbs.res(c(b1, b2, b3)) + y
data: NULL
b1 b2 b3
1.962 4.909 3.136
residual sum-of-squares: 2.587
Number of iterations to convergence: 3
Achieved convergence tolerance: NA
# even if I try to access the object's attributes, I can't find what I'm looking for
> str(ans)
List of 3
$ m :List of 16
..$ resid :function ()
..$ fitted :function ()
..$ formula :function ()
..$ deviance :function ()
..$ lhs :function ()
..$ gradient :function ()
..$ conv :function ()
..$ incr :function ()
..$ setVarying:function (vary = rep_len(TRUE, length(useParams)))
..$ setPars :function (newPars)
..$ getPars :function ()
..$ getAllPars:function ()
..$ getEnv :function ()
..$ trace :function ()
..$ Rmat :function ()
..$ predict :function (newdata = list(), qr = FALSE)
..- attr(*, "class")= chr "nlsModel"
$ call : language nls2(formula = y ~ shobbs.res(c(b1, b2, b3)) + y, start = coef(ans1n), algorithm = "brute")
$ convInfo:List of 3
..$ isConv : logi TRUE
..$ finIter: int 3
..$ finTol : logi NA
- attr(*, "class")= chr "nls"
Try running nls after running nlfb.
As an example, using the following modified from the examples section of help("nlmrt-package") :
library(nlmrt)
# setup and run nlfb example
shobbs.res <- function(x) 100.0*x[1]/(1+x[2]*10.*exp(-0.1*x[3]*seq(12))) - y
y <- c(5.308, 7.24, 9.638, 12.866, 17.069, 23.192, 31.443,
38.558, 50.156, 62.948, 75.995, 91.972)
st <- c(b1=1, b2=1, b3=1)
ans1n <- nlfb(st, shobbs.res)
print(coef(ans1n)) ##
## b1 b2 b3
## 1.9619 4.9092 3.1357
## attr(,"pkgname")
## [1] "nlmrt"
# get nls object
ans <- nls(y ~ shobbs.res(c(b1, b2, b3)) + y, start = coef(ans1n))
confint(ans)
## Waiting for profiling to be done...
## 2.5% 97.5%
## b1 1.742980 2.272051
## b2 4.563383 5.357094
## b3 2.981855 3.292911
Added
Note that this does not actually convert an object from nlfb to nls but rather performs a second optimization starting from the value found by nlfb producing a different value but perhaps that is good enough.
If not, this will convert ans1n to an "nls" class object. We first use nls2 to calculate an "nls" object. The "nls" object so produced will work for most purposes but not with confint. To get that to work we need to insert a call component into it as shown below (until this functionality is added to nls2). Now confint should run.
library(nls2)
ans_nls2 <- nls2(y ~ shobbs.res(c(b1, b2, b3)) + y, start = coef(ans1n), alg = "brute")
coef(ans_nls2) # same as ans1n above but as nls object
## b1 b2 b3
## 1.9619 4.9092 3.1357
ans_nls2$call <- ans$call # insert call component into nls object
confint(ans_nls2)
## Waiting for profiling to be done...
## 2.5% 97.5%
## b1 1.7430 2.2721
## b2 4.5634 5.3571
## b3 2.9819 3.2929