使用 nls 或 nlsLM 来拟合全局和特定于组的参数

Question

我想使用 nls 来拟合全局参数和特定于组的参数。下面是我发现的最接近最小可重现示例的示例（在此处找到：https://stat.ethz.ch/pipermail/r-help/2015-September/432020.html）

#Generate some data
d <- transform(data.frame(x=seq(0,1,len=17),
     group=rep(c("A","B","B","C"),len=17)), y =
     round(1/(1.4+x^ifelse(group=="A", 2.3, ifelse(group=="B",3.1, 3.5))),2))

#Fit to model using nls
nls(y~1/(b+x^p[group]), data=d, start=list(b=1, p=rep(3,length(levels(d$group)))))

这给了我一个错误：

Error in numericDeriv(form[[3L]], names(ind), env, central = nDcentral) : Missing value or an infinity produced when evaluating the model

我无法弄清楚错误是否来自对起始值的错误猜测，或者此代码处理特定于组的参数的方式。似乎带有 p=rep(3,length(levels(d$group))) 的行是用于生成 c(3,3,3)，但是切换这部分代码并不能解决问题（与上述相同的错误）：

#Fit to model using nls
nls(y~1/(b+x^p[group]), data=d, start=list(b=1, p=c(3, 3, 3)))

切换到 nlsLM 会出现不同的错误，导致我认为我的组特定参数有问题：

#Generate some data
library(minpack.lm)
d <- transform(data.frame(x=seq(0,1,len=17),
                          group=rep(c("A","B","B","C"),len=17)), y =
                 round(1/(1.4+x^ifelse(group=="A", 2.3, ifelse(group=="B",3.1, 3.5))),2))

#Fit to model using nlsLM
nlsLM(y~1/(b+x^p[group]), data=d, start=list(b=1, p=c(3,3,3)))

错误：

Error in dimnames(x) <- dn : length of 'dimnames' [2] not equal to array extent

有什么想法吗？

Answer 1

我认为使用 nlme::gnls 可以更轻松地做到这一点：

fit2 <- nlme::gnls(y~1/(b+x^p),
           params = list(p~group-1, b~1), 
           data=d, 
           start = list(b=1, p = rep(3,3)))

结果：

Generalized nonlinear least squares fit
  Model: y ~ 1/(b + x^p) 
  Data: d 
  Log-likelihood: 62.05887

Coefficients:
p.groupA p.groupB p.groupC        b 
2.262383 2.895903 3.475324 1.407561 

Degrees of freedom: 17 total; 13 residual
Residual standard error: 0.007188101

params 参数允许您为每个非线性参数指定固定效应子模型。使用 p ~ b-1 通过对每个组的单独估计来参数化模型，而不是为第一组拟合基线（截距）值和连续组之间的差异。（在 R 的公式语言中，-1 或 +0 表示“将没有 intercept/set 截距的模型拟合为 0”，在这种情况下对应于分别拟合所有三个组。）

我很惊讶 gnls 和 nls 没有给出相同的结果（尽管两者都给出了合理的结果）；想进一步挖掘...

参数估计（代码如下）：

  term    nls  gnls
1 b      1.41  1.40
2 pA     2.28  2.28
3 pB     3.19  3.14
4 pC     3.60  3.51

par(las = 1, bty = "l")
plot(y~x, data = d, col = d$group, pch = 16)
xvec <- seq(0, 1, length = 21)
f <- function(x) factor(x, levels = c("A","B","C"))
## fit1 is nls() fit
ll <- function(g, c = 1) {
  lines(xvec, predict(fit1, newdata = data.frame(group=f(g), x = xvec)), col = c)
}
Map(ll, LETTERS[1:3], 1:3)
d2 <- expand.grid(x = xvec, group = f(c("A","B","C")))
pp <- predict(fit2, newdata = d2)
ll2 <- function(g, c = 1) {
  lines(xvec, pp[d2$group == g], lty = 2, col = c)
}
Map(ll2, LETTERS[1:3], 1:3)
legend("bottomleft", lty = 1:2, col = 1, legend = c("nls", "gnls"))

library(tidyverse)
library(broom)
library(broom.mixed)
(purrr::map_dfr(list(nls=fit1, gnls=fit2), tidy, .id = "pkg")
  %>% select(pkg, term, estimate)
  %>% group_by(pkg)
  ## force common parameter names
  %>% mutate(across(term, ~ c("b", paste0("p", LETTERS[1:3]))))
  %>% pivot_wider(names_from = pkg, values_from = estimate)
)

Answer 2

我可以通过将组的 class 从 chr 切换为 factor 来实现。注意在生成数据集时添加factor()。

> d <- transform(data.frame(
+       x=seq(0,1,len=17),
+       group=rep(factor(c("A","B","B","C")),len=17)),
+       y=round(1/(1.4+x^ifelse(group=="A", 2.3, ifelse(group=="B",3.1, 3.5))),2)
+     )
> str(d)
'data.frame':   17 obs. of  3 variables:
 $ x    : num  0 0.0625 0.125 0.1875 0.25 ...
 $ group: Factor w/ 3 levels "A","B","C": 1 2 2 3 1 2 2 3 1 2 ...
 $ y    : num  0.71 0.71 0.71 0.71 0.69 0.7 0.69 0.69 0.62 0.64 ...
> nls(y~1/(b+x^p[group]), data=d, start=list(b=1, p=c(3,3,3)))
Nonlinear regression model
  model: y ~ 1/(b + x^p[group])
   data: d
    b    p1    p2    p3 
1.406 2.276 3.186 3.601 
 residual sum-of-squares: 9.537e-05

Number of iterations to convergence: 5 
Achieved convergence tolerance: 4.536e-06

使用 nls 或 nlsLM 来拟合全局和特定于组的参数

Using nls or nlsLM to fit global and group-specific parameters

parameters

grouping

r

nls