如何将 nlfb 对象转换为 nls 对象
How to convert an nlfb object in a nls object
我使用 nlmrt 包中的 nlfb 执行了 NLS 回归。和
theta_scaled 是我的初始参数向量
> theta_scaled
[1] 0.60000 0.73624 -0.77962
和residuals_function一个函数，它为每个参数向量计算我的数据的残差向量，我使用了调用
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.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
我想绘制一些诊断图、计算置信区间、预测等。但是，"nlmrt" 对象没有这些奇特的方法。我想将该对象转换为 nls 对象，但我不能使用 wrapnls，因为我使用 nlfb 而不是 nlxb 进行回归。有什么办法吗？
PS 如果你想知道为什么我不能使用 nlxb，原因是我使用 NLS 回归来校准复杂的流体动力学代码。因此，我的模型没有简单的分析公式。但是，由于我可以 运行 （或多或少）任意输入和参数的代码，我可以编写一个残差函数并使用 nlfb.
EDIT 正如 G. Grothendieck 在评论中正确指出的那样，我应该提供一个工作示例。就我而言，他的例子还可以。但是，他的回答不起作用：
#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"
在运行宁nlfb后尝试运行宁nls。
例如，使用从 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
已添加
请注意，这实际上并没有将对象从 nlfb 转换为 nls，而是从 nlfb 找到的值开始执行第二次优化，产生不同的值，但也许这已经足够好了。
如果不是，这会将 ans1n 转换为 "nls" class 对象。我们先用nls2计算出一个"nls"对象。如此生成的 "nls" 对象将适用于大多数用途，但不适用于 confint。为了让它工作，我们需要插入一个 call 组件，如下所示（直到这个功能被添加到 nls2）。现在 confint 应该 运行.
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
我使用 nlmrt 包中的 nlfb 执行了 NLS 回归。和
theta_scaled 是我的初始参数向量
> theta_scaled
[1] 0.60000 0.73624 -0.77962
和residuals_function一个函数，它为每个参数向量计算我的数据的残差向量，我使用了调用
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.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
我想绘制一些诊断图、计算置信区间、预测等。但是，"nlmrt" 对象没有这些奇特的方法。我想将该对象转换为 nls 对象，但我不能使用 wrapnls，因为我使用 nlfb 而不是 nlxb 进行回归。有什么办法吗？
PS 如果你想知道为什么我不能使用 nlxb，原因是我使用 NLS 回归来校准复杂的流体动力学代码。因此，我的模型没有简单的分析公式。但是，由于我可以 运行 （或多或少）任意输入和参数的代码，我可以编写一个残差函数并使用 nlfb.
EDIT 正如 G. Grothendieck 在评论中正确指出的那样，我应该提供一个工作示例。就我而言，他的例子还可以。但是，他的回答不起作用：
#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"
在运行宁nlfb后尝试运行宁nls。
例如，使用从 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
已添加
请注意，这实际上并没有将对象从 nlfb 转换为 nls，而是从 nlfb 找到的值开始执行第二次优化，产生不同的值，但也许这已经足够好了。
如果不是，这会将 ans1n 转换为 "nls" class 对象。我们先用nls2计算出一个"nls"对象。如此生成的 "nls" 对象将适用于大多数用途，但不适用于 confint。为了让它工作，我们需要插入一个 call 组件，如下所示（直到这个功能被添加到 nls2）。现在 confint 应该 运行.
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