群NLS函数

Question

我有一个数据集，我想按组应用非线性最小二乘。这是我上一个问题的延续：NLS Function - Number of Iterations Exceeds max

数据集如下所示：

df
x        y    GRP
0        0      1
426   9.28      1
853   18.5      1
1279  27.8      1
1705  37.0      1
2131  46.2      1
0        0      2
450   7.28      2
800   16.5      2
1300  30.0      2
2000  40.0      2
2200  48.0      2  

如果我和一组人一起做这件事，就会是这样：

df1<-filter(df, GRP==1)

a.start <- max(df1$y)
b.start <- 1e-06
control1 <- nls.control(maxiter= 10000,tol=1e-02, warnOnly=TRUE)
nl.reg <- nls(y ~ a * (1-exp(-b * x)),data=df1,start= 
list(a=a.start,b=b.start),
           control= control1)
coef(nl.reg)[1]
coef(nl.reg)[2]

> coef(nl.reg)[1]
       a 
5599.075 
> coef(nl.reg)[2]
       b 
3.891744e-06 

然后，我也会为GRP2做同样的事情。我希望我的最后输出如下所示：

x        y    GRP                       a                       b
0        0      1                5599.075            3.891744e-06
426   9.28      1                5599.075            3.891744e-06
853   18.5      1                5599.075            3.891744e-06
1279  27.8      1                5599.075            3.891744e-06
1705  37.0      1                5599.075            3.891744e-06
2131  46.2      1                5599.075            3.891744e-06
0        0      2    New Value for a GRP2    New Value for b GRP2     
450   7.28      2    New Value for a GRP2    New Value for b GRP2
800   16.5      2    New Value for a GRP2    New Value for b GRP2
1300  30.0      2    New Value for a GRP2    New Value for b GRP2
2000  40.0      2    New Value for a GRP2    New Value for b GRP2
2200  48.0      2    New Value for a GRP2    New Value for b GRP2

理想情况下，我认为dplyr将是最好的方式，但我不知道如何去做。这就是我认为可能会是这样的：

control1 <- nls.control(maxiter= 10000,tol=1e-02, warnOnly=TRUE)
b.start <- 1e-06

df %>%
  group_by(GRP) %>%
  do(nlsfit = nls( form = y ~ a * (1-exp(-b * x)), data=., 
start= list( a=max(.$y), b=b.start),
      control= control1) ) %>%
  list(a = coef(nlsfit)[1], b = coef(nlsfit)[2])

错误：

 in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates

不过，不太清楚该如何做到这一点，任何帮助都是很好的。谢谢!

Answer 1

最初，我得到了相同的错误消息(re: NotFindObject‘y’in nls)，这与我最初尝试使用lapply-split-function范式时使用tidyverse表时所做的一样，并搜索："r使用nls内部函数“。我已经将attach的最初用法更改为list2env：

sapply(  split( df , df$GRP), function(d){ dat <- list2env(d)
    nlsfit <- nls( form = y ~ a * (1-exp(-b * x)), data=dat, start= list( a=max(y), b=b.start),
          control= control1) 

list(a = coef(nlsfit)[1], b = coef(nlsfit)[2])} )
#---

  1            2            
a 14.51827     441.5489     
b 2.139378e-06 -6.775562e-06

你也会收到你期待的警告。这些都可以被suppressWarnings( ... )抑制。

其中一个建议是使用attach。我当时非常不情愿地这样做，因为我经常警告新手不要使用attach。但在这里，它似乎迫使当地的环境建设。我更喜欢用list2env作为一种满足nls的机制。nls代码的顶部是我选择的原因：

if (!is.list(data) && !is.environment(data)) 
    stop("'data' must be a list or an environment")