凹时在 R 中使用 NLS 的三段式分段回归的语法
Syntax for three-piece segmented regression using NLS in R when concave
我的目标是拟合一个三段式（即两个断点）回归模型以使用传播的 predictNLS 函数进行预测，确保将节点定义为参数，但我的模型公式似乎不正确.
我已经使用 segmented 包来估计断点位置（在 NLS 中用作起始值），但我想将我的模型保留为 NLS 格式，特别是 nlsLM {minipack.lm} 因为我正在使用 NLS 将其他类型的曲线拟合到我的数据中，希望允许 NLS 优化节点值，有时使用可变权重，并且需要能够从 propagate 轻松计算 Monte Carlo 置信区间].尽管我非常接近 具有正确的公式语法，但我在断点附近没有得到 expected/required 行为。这些段应该直接在断点处相遇（没有任何跳跃），但至少在这个数据上，我在断点处得到了一个奇怪的局部最小值（见下图）。
以下是我的数据和一般过程的示例。我认为我的问题出在 NLS 公式中。
library(minpack.lm)
library(segmented)
y <- c(-3.99448113, -3.82447011, -3.65447803, -3.48447030, -3.31447855, -3.14448753, -2.97447972, -2.80448401, -2.63448380, -2.46448069, -2.29448796, -2.12448912, -1.95448783, -1.78448797, -1.61448563, -1.44448719, -1.27448469, -1.10448651, -0.93448525, -0.76448637, -0.59448626, -0.42448586, -0.25448588, -0.08448548, 0.08551417, 0.25551393, 0.42551411, 0.59551395, 0.76551389, 0.93551398)
x <- c(61586.1711, 60330.5550, 54219.9925, 50927.5381, 48402.8700, 45661.9175, 37375.6023, 33249.1248, 30808.6131, 28378.6508, 22533.3782, 13901.0882, 11716.5669, 11004.7305, 10340.3429, 9587.7994, 8736.3200, 8372.1482, 8074.3709, 7788.1847, 7499.6721, 7204.3168, 6870.8192, 6413.0828, 5523.8097, 3961.6114, 3460.0913, 2907.8614, 2016.1158, 452.8841)
df<- data.frame(x,y)
#Use Segmented to get estimates for parameters with 2 breakpoints
my.seg2 <- segmented(lm(y ~ x, data = df), seg.Z = ~ x, npsi = 2)
#extract knot, intercept, and coefficient values to use as NLS start points
my.knot1 <- my.seg2$psi[1,2]
my.knot2 <- my.seg2$psi[2,2]
my.m_2 <- slope(my.seg2)$x[1,1]
my.b1 <- my.seg2$coefficients[[1]]
my.b2 <- my.seg2$coefficients[[2]]
my.b3 <- my.seg2$coefficients[[3]]
#Fit a NLS model to ~replicate segmented model. Presumably my model formula is where the problem lies
my.model <- nlsLM(y~m*x+b+(b2*(ifelse(x>=knot1&x<=knot2,1,0)*(x-knot1))+(b3*ifelse(x>knot2,1,0)*(x-knot2-knot1))),data=df, start = c(m = my.m_2, b = my.b1, b2 = my.b2, b3 = my.b3, knot1 = my.knot1, knot2 = my.knot2))
应该看起来如何
plot(my.seg2)
看起来如何
plot(x, y)
lines(x=x, y=predict(my.model), col='black', lty = 1, lwd = 1)
我很确定我是“正确的”，但是当用直线绘制 95% 的置信区间并且预测分辨率（例如，x 点的密度）增加时，事情似乎 非常不正确.
谢谢大家的帮助。
它可能部分反映了 segmented 中的限制。 segmented returns 没有量化相关不确定性的单个变化点值。使用 mcp 重做分析，其中 returns 贝叶斯后验，我们看到第二个变化点是双峰分布的：
library(mcp)
model = list(
y ~ 1 + x, # Intercept + slope in first segment
~ 0 + x, # Only slope changes in the next segments
~ 0 + x
)
# Fit it with a large number of samples and plot the change point posteriors
fit = mcp(model, data = data.frame(x, y), iter = 50000, adapt = 10000)
plot_pars(fit, regex_pars = "^cp*", type = "dens_overlay")
仅供参考，mcp 也可以绘制可信区间（红色虚线）：
plot(fit, q_fit = TRUE)
将 g 定义为与 x 具有相同长度的分组向量，对于 X 轴的 3 个部分采用值 1、2、3，并从中创建 nls 模型。结果图看起来不错。
my.knots <- c(my.knot1, my.knot2)
g <- cut(x, c(-Inf, my.knots, Inf), label = FALSE)
fm <- nls(y ~ a[g] + b[g] * x, df, start = list(a = c(1, 1, 1), b = c(1, 1, 1)))
plot(y ~ x, df)
lines(fitted(fm) ~ x, df, col = "red")
（图后续）
约束条件
虽然上面看起来不错并且可能已经足够了，但它并不能保证线段在节点处相交。为此，我们必须施加约束，即双方在节点处相等：
a[2] + b[2] * my.knots[1] = a[1] + b[1] * my.knots[1]
a[3] + b[3] * my.knots[2] = a[2] + b[2] * my.knots[2]
所以
a[2] = a[1] + (b[1] - b[2]) * my.knots[1]
a[3] = a[2] + (b[2] - b[3]) * my.knots[2]
= a[1] + (b[1] - b[2]) * my.knots[1] + (b[2] - b[3]) * my.knots[2]
给予：
# returns a vector of the three a values
avals <- function(a1, b) unname(cumsum(c(a1, -diff(b) * my.knots)))
fm2 <- nls(y ~ avals(a1, b)[g] + b[g] * x, df, start = list(a1 = 1, b = c(1, 1, 1)))
要获得我们可以使用的三个 a 值：
co <- coef(fm2)
avals(co[1], co[-1])
求残差平方和：
deviance(fm2)
## [1] 0.193077
多项式
虽然涉及大量参数，但可以使用多项式拟合代替分段线性回归。 12 次多项式涉及 13 个参数，但残差平方和低于分段线性回归。随着残差平方和的相应增加，可以使用较低的次数。 7次多项式涉及8个参数，虽然残差平方和较高，但视觉上看起来还不错。
fm12 <- nls(y ~ cbind(1, poly(x, 12)) %*% b, df, start = list(b = rep(1, 13)))
deviance(fm12)
## [1] 0.1899218
我的目标是拟合一个三段式（即两个断点）回归模型以使用传播的 predictNLS 函数进行预测，确保将节点定义为参数，但我的模型公式似乎不正确.
我已经使用 segmented 包来估计断点位置（在 NLS 中用作起始值），但我想将我的模型保留为 NLS 格式，特别是 nlsLM {minipack.lm} 因为我正在使用 NLS 将其他类型的曲线拟合到我的数据中，希望允许 NLS 优化节点值，有时使用可变权重，并且需要能够从 propagate 轻松计算 Monte Carlo 置信区间].尽管我非常接近 具有正确的公式语法，但我在断点附近没有得到 expected/required 行为。这些段应该直接在断点处相遇（没有任何跳跃），但至少在这个数据上，我在断点处得到了一个奇怪的局部最小值（见下图）。
以下是我的数据和一般过程的示例。我认为我的问题出在 NLS 公式中。
library(minpack.lm)
library(segmented)
y <- c(-3.99448113, -3.82447011, -3.65447803, -3.48447030, -3.31447855, -3.14448753, -2.97447972, -2.80448401, -2.63448380, -2.46448069, -2.29448796, -2.12448912, -1.95448783, -1.78448797, -1.61448563, -1.44448719, -1.27448469, -1.10448651, -0.93448525, -0.76448637, -0.59448626, -0.42448586, -0.25448588, -0.08448548, 0.08551417, 0.25551393, 0.42551411, 0.59551395, 0.76551389, 0.93551398)
x <- c(61586.1711, 60330.5550, 54219.9925, 50927.5381, 48402.8700, 45661.9175, 37375.6023, 33249.1248, 30808.6131, 28378.6508, 22533.3782, 13901.0882, 11716.5669, 11004.7305, 10340.3429, 9587.7994, 8736.3200, 8372.1482, 8074.3709, 7788.1847, 7499.6721, 7204.3168, 6870.8192, 6413.0828, 5523.8097, 3961.6114, 3460.0913, 2907.8614, 2016.1158, 452.8841)
df<- data.frame(x,y)
#Use Segmented to get estimates for parameters with 2 breakpoints
my.seg2 <- segmented(lm(y ~ x, data = df), seg.Z = ~ x, npsi = 2)
#extract knot, intercept, and coefficient values to use as NLS start points
my.knot1 <- my.seg2$psi[1,2]
my.knot2 <- my.seg2$psi[2,2]
my.m_2 <- slope(my.seg2)$x[1,1]
my.b1 <- my.seg2$coefficients[[1]]
my.b2 <- my.seg2$coefficients[[2]]
my.b3 <- my.seg2$coefficients[[3]]
#Fit a NLS model to ~replicate segmented model. Presumably my model formula is where the problem lies
my.model <- nlsLM(y~m*x+b+(b2*(ifelse(x>=knot1&x<=knot2,1,0)*(x-knot1))+(b3*ifelse(x>knot2,1,0)*(x-knot2-knot1))),data=df, start = c(m = my.m_2, b = my.b1, b2 = my.b2, b3 = my.b3, knot1 = my.knot1, knot2 = my.knot2))
应该看起来如何
plot(my.seg2)
看起来如何
plot(x, y)
lines(x=x, y=predict(my.model), col='black', lty = 1, lwd = 1)
我很确定我是“正确的”，但是当用直线绘制 95% 的置信区间并且预测分辨率（例如，x 点的密度）增加时，事情似乎 非常不正确.
谢谢大家的帮助。
它可能部分反映了 segmented 中的限制。 segmented returns 没有量化相关不确定性的单个变化点值。使用 mcp 重做分析，其中 returns 贝叶斯后验，我们看到第二个变化点是双峰分布的：
library(mcp)
model = list(
y ~ 1 + x, # Intercept + slope in first segment
~ 0 + x, # Only slope changes in the next segments
~ 0 + x
)
# Fit it with a large number of samples and plot the change point posteriors
fit = mcp(model, data = data.frame(x, y), iter = 50000, adapt = 10000)
plot_pars(fit, regex_pars = "^cp*", type = "dens_overlay")
仅供参考，mcp 也可以绘制可信区间（红色虚线）：
plot(fit, q_fit = TRUE)
将 g 定义为与 x 具有相同长度的分组向量，对于 X 轴的 3 个部分采用值 1、2、3，并从中创建 nls 模型。结果图看起来不错。
my.knots <- c(my.knot1, my.knot2)
g <- cut(x, c(-Inf, my.knots, Inf), label = FALSE)
fm <- nls(y ~ a[g] + b[g] * x, df, start = list(a = c(1, 1, 1), b = c(1, 1, 1)))
plot(y ~ x, df)
lines(fitted(fm) ~ x, df, col = "red")
（图后续）
约束条件
虽然上面看起来不错并且可能已经足够了，但它并不能保证线段在节点处相交。为此，我们必须施加约束，即双方在节点处相等：
a[2] + b[2] * my.knots[1] = a[1] + b[1] * my.knots[1]
a[3] + b[3] * my.knots[2] = a[2] + b[2] * my.knots[2]
所以
a[2] = a[1] + (b[1] - b[2]) * my.knots[1]
a[3] = a[2] + (b[2] - b[3]) * my.knots[2]
= a[1] + (b[1] - b[2]) * my.knots[1] + (b[2] - b[3]) * my.knots[2]
给予：
# returns a vector of the three a values
avals <- function(a1, b) unname(cumsum(c(a1, -diff(b) * my.knots)))
fm2 <- nls(y ~ avals(a1, b)[g] + b[g] * x, df, start = list(a1 = 1, b = c(1, 1, 1)))
要获得我们可以使用的三个 a 值：
co <- coef(fm2)
avals(co[1], co[-1])
求残差平方和：
deviance(fm2)
## [1] 0.193077
多项式
虽然涉及大量参数，但可以使用多项式拟合代替分段线性回归。 12 次多项式涉及 13 个参数，但残差平方和低于分段线性回归。随着残差平方和的相应增加，可以使用较低的次数。 7次多项式涉及8个参数，虽然残差平方和较高，但视觉上看起来还不错。
fm12 <- nls(y ~ cbind(1, poly(x, 12)) %*% b, df, start = list(b = rep(1, 13)))
deviance(fm12)
## [1] 0.1899218