如何在 R 中使用共享参数执行非线性最小二乘？答案

【问题标题】：How to perform nonlinear least squares with shared parameters in R?如何在 R 中使用共享参数执行非线性最小二乘？
【发布时间】：2020-08-21 22:37:16
【问题描述】：

我想在 R 中执行非线性最小二乘回归，同时最小化三个模型的平方残差（见下文）。现在，三个模型共享一些参数，在我的示例中，参数b 和d。

有没有办法使用nls() 或包minpack.lm 或nlsr 来做到这一点？

所以，理想情况下，我想生成目标函数（所有模型的最小二乘之和）并一次回归所有参数：a1、a2、a3、b、@ 987654331@、c2、c3 和 d。

（我试图避免运行三个独立的回归，然后对 b 和 d 进行一些平均。）

my_model <- function(x, a, b, c, d) {
  a * b ^ (x - c) + d
}

# x values
x <- seq(0, 10, 0.2)

# Shared parameters
b <- 2
d <- 10

a1 <- 1
c1 <- 1
y1 <- my_model(x,
               a = a1,
               b = b,
               c = c1,
               d = d) + rnorm(length(x))

a2 <- 2
c2 <- 5
y2 <- my_model(x,
               a = a2,
               b = b,
               c = c2,
               d = d) + rnorm(length(x))

a3 <- -2
c3 <- 3
y3 <- my_model(x,
               a = a3,
               b = b,
               c = c3,
               d = d) + rnorm(length(x))

plot(
  y1 ~ x,
  xlim = range(x),
  ylim = d + c(-50, 50),
  type = 'b',
  col = 'red',
  ylab = 'y'
)
lines(y2 ~ x, type = 'b', col = 'green')
lines(y3 ~ x, type = 'b', col = 'blue')

【问题讨论】：

我对统计方法很有信心。我的问题实际上是关于 R 中的实现。

标签： r non-linear-regression nls

【解决方案1】：

下面我们运行nls（使用稍微修改的模型）和nlxb（来自nlsr）但nlxb在收敛之前停止。尽管存在这些问题，但这两个问题确实给出了在视觉上很好地拟合数据的结果。这些问题表明模型本身存在问题，因此在 Other 部分中，在nlxb 输出的指导下，我们展示了如何修复模型，给出适合数据的原始模型的子模型轻松使用nls 和nlxb，并且非常适合。在备注部分的最后，我们以可重复的形式提供数据。

nls

假设最后的注释中显示的设置可重现，通过定义右侧矩阵，其列分别乘以每个线性参数 a1、a2、a3 和 d，重新表述 nls plinear 算法的问题。 plinear 不需要那些简化设置的起始值。它将分别报告为 .lin1、.lin2、.lin3 和 .lin4。

为了获得起始值，我们使用了一个没有分组的更简单的模型，并使用同名包中的nls2 对 b 从 1 到 10 和 c 从 1 到 10 进行网格搜索。我们还发现nls仍然会产生错误，但是通过在公式中使用abs，如图所示，它运行完成。

模型的问题表明它存在根本问题，在“其他”部分中，我们将讨论如何解决它。

xx <- c(x, x, x)
yy <- c(y1, y2, y3)

# startingi values using nls2
library(nls2)
fo0 <- yy ~ cbind(b ^ abs(xx - c), 1)
st0 <- data.frame(b = c(1, 10), c = c(1, 10))
fm0 <- nls2(fo0, start = st0, alg = "plinear-brute")

# run nls using starting values from above
g <- rep(1:3, each = length(x))   
fo <- yy ~ cbind((g==1) * b ^ abs(xx - c[g]), 
                 (g==2) * b ^ abs(xx - c[g]),  
                 (g==3) * b ^ abs(xx - c[g]), 
                 1) 
st <- with(as.list(coef(fm0)), list(b = b, c = c(c, c, c)))
fm <- nls(fo, start = st, alg = "plinear")

plot(yy ~ xx, col = g)
for(i in unique(g)) lines(predict(fm) ~ xx, col = i, subset = g == i)

fm

给予：

Nonlinear regression model
  model: yy ~ cbind((g == 1) * b^abs(xx - c[g]), (g == 2) * b^abs(xx -     c[g]), (g == 3) * b^abs(xx - c[g]), 1)
   data: parent.frame()
     b     c1     c2     c3  .lin1  .lin2  .lin3  .lin4 
 1.997  0.424  1.622  1.074  0.680  0.196 -0.532  9.922 
 residual sum-of-squares: 133

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

（剧情后续）

nlsr

使用 nlsr 会像这样完成。不需要网格搜索起始值，也不需要添加abs。 b 和 d 值似乎与 nls 解决方案相似，但其他系数不同。从视觉上看，这两种解决方案似乎都适合数据。

另一方面，从 JSingval 列中，我们看到 jacobian 秩不足，导致它停止并且不产生 SE 值，并且收敛性存在疑问（尽管考虑到在视觉上没有显示的情节，这可能就足够了，看起来很合适）。我们将在“其他”部分讨论如何解决此问题。

g1 <- g == 1; g2 <- g == 2; g3 <- g == 3
fo2 <- yy ~ g1 * (a1 * b ^ (xx - c1) + d) + 
            g2 * (a2 * b ^ (xx - c2) + d) + 
            g3 * (a3 * b ^ (xx - c3) + d)
st2 <- list(a1 = 1, a2 = 1, a3 = 1, b = 1, c1 = 1, c2 = 1, c3 = 1, d = 1)
fm2 <- nlxb(fo2, start = st2)
fm2

给予：

vn: [1] "yy" "g1" "a1" "b"  "xx" "c1" "d"  "g2" "a2" "c2" "g3" "a3" "c3"
no weights
nlsr object: x 
residual sumsquares =  133.45  on  153 observations
    after  16    Jacobian and  22 function evaluations
  name            coeff          SE       tstat      pval      gradient    JSingval   
a1               3.19575            NA         NA         NA    9.68e-10        4097  
a2               0.64157            NA         NA         NA   8.914e-11       662.5  
a3              -1.03096            NA         NA         NA  -1.002e-09       234.9  
b                1.99713            NA         NA         NA   -2.28e-08       72.57  
c1               2.66146            NA         NA         NA   -2.14e-09       10.25  
c2               3.33564            NA         NA         NA  -3.955e-11   1.585e-13  
c3                2.0297            NA         NA         NA  -7.144e-10   1.292e-13  
d                9.92363            NA         NA         NA  -2.603e-12   3.271e-14

我们可以使用 nls2 作为第二阶段来计算 SE，但这仍然不能解决奇异值所暗示的整个问题。

summary(nls2(fo2, start = coef(fm2), algorithm = "brute-force"))

给予：

Formula: yy ~ g1 * (a1 * b^(xx - c1) + d) + g2 * (a2 * b^(xx - c2) + d) + 
    g3 * (a3 * b^(xx - c3) + d)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
a1  3.20e+00   5.38e+05     0.0        1    
a2  6.42e-01   3.55e+05     0.0        1    
a3 -1.03e+00   3.16e+05     0.0        1    
b   2.00e+00   2.49e-03   803.4   <2e-16 ***
c1  2.66e+00   9.42e-02    28.2   <2e-16 ***
c2  3.34e+00   2.43e+05     0.0        1    
c3  2.03e+00   8.00e+05     0.0        1    
d   9.92e+00   4.42e+05     0.0        1    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.959 on 145 degrees of freedom

Number of iterations to convergence: 8 
Achieved convergence tolerance: NA

其他

当nls 无法拟合模型时，通常表明模型本身存在问题。在上面 nlsr 输出中的 JSingval 列的指导下尝试一下，这表明 c 参数或 d 可能是问题所在，我们发现如果我们将所有 c 参数值修复为 0，那么模型是在给定足够好的起始值的情况下很容易拟合，并且仍然给出较低的残差平方和。

library(nls2)

fo3 <- yy ~ cbind((g==1) * b ^ xx, (g==2) * b ^ xx, (g==3) * b ^ xx, 1) 
st3 <-  coef(fm0)["b"]
fm3 <- nls(fo3, start = st3, alg = "plinear")

给予：

Nonlinear regression model
  model: yy ~ cbind((g == 1) * b^xx, (g == 2) * b^xx, (g == 3) * b^xx,     1)
   data: parent.frame()
      b   .lin1   .lin2   .lin3   .lin4 
 1.9971  0.5071  0.0639 -0.2532  9.9236 
 residual sum-of-squares: 133

Number of iterations to convergence: 4 
Achieved convergence tolerance: 1.67e-09

尽管参数少了 3 个，但以下 anova 表示的参数与上面的 fm 相当：

anova(fm3, fm)

给予：

Analysis of Variance Table

Model 1: yy ~ cbind((g == 1) * b^xx, (g == 2) * b^xx, (g == 3) * b^xx, 1)
Model 2: yy ~ cbind((g == 1) * b^abs(xx - c[g]), (g == 2) * b^abs(xx - c[g]), (g == 3) * b^abs(xx - c[g]), 1)
  Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
1    148        134                         
2    145        133  3  0.385    0.14   0.94

我们可以像这样使用nlxb 重做fm3：

fo4 <- yy ~ g1 * (a1 * b ^ xx + d) + 
            g2 * (a2 * b ^ xx + d) + 
            g3 * (a3 * b ^ xx + d)
st4 <- list(a1 = 1, a2 = 1, a3 = 1, b = 1, d = 1)
fm4 <- nlxb(fo4, start = st4)
fm4

给予：

nlsr object: x 
residual sumsquares =  133.45  on  153 observations
    after  24    Jacobian and  33 function evaluations
  name            coeff          SE       tstat      pval      gradient    JSingval   
a1              0.507053      0.005515      91.94  1.83e-132   8.274e-08        5880  
a2             0.0638554     0.0008735      73.11  4.774e-118    1.26e-08        2053  
a3             -0.253225      0.002737     -92.54  7.154e-133  -4.181e-08        2053  
b                1.99713      0.002294      870.6  2.073e-276   -2.55e-07       147.5  
d                9.92363       0.09256      107.2  3.367e-142  -1.219e-11       10.26

注意

下面假设的输入与问题中的相同，除了我们另外设置种子以使其可重复。

set.seed(123)

my_model <- function(x, a, b, c, d) a * b ^ (x - c) + d

x <- seq(0, 10, 0.2)

b <- 2; d <- 10 # shared

a1 <- 1; c1 <- 1
y1 <- my_model(x, a = a1, b = b, c = c1, d = d) + rnorm(length(x))

a2 <- 2; c2 <- 5
y2 <- my_model(x, a = a2, b = b, c = c2, d = d) + rnorm(length(x))

a3 <- -2; c3 <- 3
y3 <- my_model(x, a = a3, b = b, c = c3, d = d) + rnorm(length(x))

【讨论】：

为什么SE、tstat和pval都等于NA？
在 nlsr 部分添加了讨论，解释了问题及其解决方案，并多次修改了重要的 Other 部分。

【解决方案2】：

我不确定这是否真的是最好的方法，但您可以使用 optim() 最小化残差平方和。

#start values
params <- c(a1=1, a2=1, a3=1, b=1, c1=1, c2=1, c3=1,d=1)
# minimize total sum of squares of residuals
fun <- function(p) {
  sum(
    (y1-my_model(x, p["a1"], p["b"], p["c1"], p["d"]))^2 + 
    (y2-my_model(x, p["a2"], p["b"], p["c2"], p["d"]))^2 +
    (y3-my_model(x, p["a3"], p["b"], p["c3"], p["d"]))^2
  )
}
out <- optim(params, fun, method="BFGS")
out$par
#         a1         a2         a3          b         c1         c2         c3 
#  0.8807542  1.0241804 -2.8805848  1.9974615  0.7998103  4.0030597  3.5184600 
#          d 
#  9.8764917

我们可以在图像顶部添加绘图

curve(my_model(x, out$par["a1"], out$par["b"], out$par["c1"], out$par["d"]), col="red", add=T)
curve(my_model(x, out$par["a2"], out$par["b"], out$par["c2"], out$par["d"]), col="green", add=T)
curve(my_model(x, out$par["a3"], out$par["b"], out$par["c3"], out$par["d"]), col="blue", add=T)

【讨论】：

赞赏。但是通过使用optim，它并没有真正进行回归......所以我无法获得结果统计信息，例如估计的置信区间等。
@rmango 我不确定你为什么不认为这是回归。线性模型也只是最小二乘拟合。您只需要做出额外的建模假设来估计错误和间隔。目前尚不清楚您需要什么类型的模型。如果您对自己的统计方法有信心，那么如果您能更清楚地描述它会很有帮助。
如果您指的是误差项：我假设一个均值为 0 的正态分布误差和一些固定（未知）方差，正如 rnorm(length(x)) 所暗示的那样。
nls2 可用作获取 SE 的第二阶段。使用fo2 和我的回答中定义的变量：library(nls2); fit <- nls2(fo2, start = out$par, algorithm = "brute-force") 但是，问题是雅可比在我的回答中讨论的等级不足。注意svd(fit$m$Rmat())$d 中的零奇异值对收敛性产生怀疑。
@G.Grothendieck 你不能使用引导程序来获取标准错误吗？但即便如此，解释 b 之类的 p 值似乎也很棘手。对 0 的测试真的意味着它很重要吗？如果它是 1，那意味着它并没有真正做任何事情。