R 中的非线性回归,nls:奇异梯度

Dav*_*.F. 0 regression r nls ggplot2 ggplotly

我想将我的数据拟合到已使用 Matlab 优化的特定函数中。

我收到以下错误:“警告消息:计算失败stat_smooth():奇异梯度”

请帮忙!这是我的 R 代码:

tibble
       x     y     SEM
 1     1 0.0342 0.00532
 2     3 0.0502 0.00639
 3     5 0.0700 0.0118 
 4    10 0.123  0.0269 
 5    20 0.154  0.0125 
 6    30 0.203  0.0190 
 7    40 0.257  0.0255 
 8    50 0.287  0.0266 
 9    60 0.345  0.0347 
10    90 0.442  0.0398 
11   120 0.569  0.0570 
12   180 0.726  0.0406 
13   240 0.824  0.0150 
14   360 0.868  0.00821
15  1440 0.890  0.0246 

tibble %>% 
  ggplot(aes(x, y)) +
  geom_point()+
  geom_errorbar(aes(ymin=y-SEM, ymax=y+SEM), width=25)+
  geom_ribbon(aes(ymin = y-2.575*SEM, ymax = y+2.575*SEM), alpha = 0.1)+
  geom_smooth(method="nls", 
              formula= y ~ (1-((k2/(k2-k1))*exp(-k1*x))+((k1/(k2-k1))*exp(-k2*x))),
              se=F,
              method.args = list(start=list(k1=0.006999, k2=849.6)))
Run Code Online (Sandbox Code Playgroud)

All*_*ron 6

就目前情况而言,nls无法确定 k2 的最佳值,因为随着 k2 趋于无穷大,平方和逐渐减小到大约 0.0225 的值。因此,不存在使平方和最小化的 k2 的有限值。由于 k2 趋于无穷大,因此它有效地抵消了,使公式等于:

y ~ 1 - exp(-k1*x)
Run Code Online (Sandbox Code Playgroud)

这意味着对于任何有限的 k2 值,该公式始终比原始公式更好地拟合数据。

简而言之,k2 是一个冗余参数,只会使您的拟合效果恶化。

因此你的情节可以是这样的:

ggplot(df, aes(x, y)) +
  geom_point() +
  geom_errorbar(aes(ymin = y - SEM, ymax = y + SEM), width = 25) +
  geom_ribbon(aes(ymin = y - 2.575 * SEM, ymax = y + 2.575 * SEM), alpha = 0.1) +
  geom_smooth(method = "nls", 
              formula = y ~ 1 - exp(-k1 * x),
              se = FALSE,
              method.args = list(start = list(k1 = 0.006999)))
Run Code Online (Sandbox Code Playgroud)

在此输入图像描述

或者,正如 G. Grothendieck 建议的那样,使用额外的参数来优化像这样的拟合(假设它在您的用例中具有物理意义)

ggplot(df, aes(x, y)) +
  geom_point() +
  geom_errorbar(aes(ymin = y - SEM, ymax = y + SEM), width = 25) +
  geom_ribbon(aes(ymin = y - 2.575 * SEM, ymax = y + 2.575 * SEM), alpha = 0.1) +
  geom_smooth(method = "nls", 
              formula = y ~ k2 * (1 - exp(-k1 * x)),
              se = FALSE,
              method.args = list(start = list(k1 = 0.006999, k2 = 1)))
Run Code Online (Sandbox Code Playgroud)

在此输入图像描述