我想在 R 中执行非线性最小二乘回归，同时最小化三个模型的平方残差(见下文)。现在，这三个模型共享一些参数，在我的示例中，参数 b 和 d。
有没有一种方法可以使用 nls()，或者包 minpack.lm 或 nlsr 来做到这一点？
因此，理想情况下，我想生成目标函数(所有模型的最小二乘总和)并一次回归所有参数:a1、a2 , a3, b, c1, 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')
最佳答案
下面我们运行 nls(使用稍微修改的模型)和 nlxb(来自 nlsr)，但 nlxb 在收敛之前停止。尽管存在这些问题，但这两种方法确实给出了在视觉上很好地符合数据的结果。这些问题表明模型本身存在问题，因此在 Other 部分，在 nlxb 输出的指导下，我们展示了如何修复模型，给出原始模型的子模型使用 nls 和 nlxb 轻松拟合数据的模型，也给出了很好的拟合。在注释部分的末尾，我们以可重现的形式提供了数据。
nls
假设最后的注释中显示的设置可重复，通过定义右侧矩阵重新表述 nls 线性算法的问题，其列分别乘以每个线性参数 a1、a2、a3 和 d。 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 个参数，但以下方差分析表明它与上面的 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))
关于r - 如何在 R 中执行具有共享参数的非线性最小二乘法？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/63530860/