具有二次多项式和在断点处平滑连接的直线的分段回归
Piecewise regression with a quadratic polynomial and a straight line joining smoothly at a break point
我想用一个断点 xt 拟合分段线性回归，这样对于 x < xt 我们有一个二次多项式，对于 x >= xt 我们有一条直线。两块应该平滑连接，连续性达到 xt 处的一阶导数。这是它可能看起来像的图片：
我将分段回归函数参数化为：
其中 a、b、c 和 xt 是要估计的参数。
我想根据调整后的 R 平方将此模型与整个范围内的二次多项式回归进行比较。
这是我的数据：
y <- c(1, 0.59, 0.15, 0.078, 0.02, 0.0047, 0.0019, 1, 0.56, 0.13,
0.025, 0.0051, 0.0016, 0.00091, 1, 0.61, 0.12, 0.026, 0.0067,
0.00085, 4e-04)
x <- c(0, 5.53, 12.92, 16.61, 20.3, 23.07, 24.92, 0, 5.53, 12.92,
16.61, 20.3, 23.07, 24.92, 0, 5.53, 12.92, 16.61, 20.3, 23.07,
24.92)
我的尝试如下，对于已知的 xt:
z <- pmax(0, x - xt)
x1 <- pmin(x, xt)
fit <- lm(y ~ x1 + I(x1 ^ 2) + z - 1)
但是这条直线在xt处似乎并不与二次多项式相切。我哪里做错了？
类似问题：
- Piecewise regression with a straight line and a horizontal line joining at a break point
- Fitting a V-shape curve to my data（交叉验证）
这是一个很好的练习（可能很难）来消化线性模型的理论和实现。我的回答将包含两部分：
- 第 1 部分（这一部分）介绍了我使用的参数化以及这种分段回归如何简化为普通的最小二乘问题。提供模型估计、断点选择、预测等R函数
- 第 2 部分（另一个）演示了一个可复制的玩具示例，说明如何使用我定义的那些函数。
我必须使用不同的参数化，因为你在问题中给出的参数是错误的！你的参数化只保证了函数值的连续性，而不是一阶导数！这就是为什么您的拟合线与 xt.
处的拟合二次多项式不相切的原因
## generate design matrix
getX <- function (x, c) {
x <- x - c
cbind("beta0" = 1, "beta1" = x, "beta2" = pmin(x, 0) ^ 2)
}
下面的函数 est 在给定的 c:
下总结了模型的估计和推理的 .lm.fit（为了实现最大效率）
## `x`, `y` give data points; `c` is known break point
est <- function (x, y, c) {
## model matrix
X <- getX(x, c)
p <- dim(X)[2L]
## solve least squares with QR factorization
fit <- .lm.fit(X, y)
## compute Pearson estimate of `sigma ^ 2`
r <- c(fit$residuals)
n <- length(r)
RSS <- c(crossprod(r))
sig2 <- RSS / (n - p)
## coefficients summary table
beta <- fit$coefficients
R <- "dimnames<-"(fit$qr[1:p, ], NULL)
Rinv <- backsolve(R, diag(p))
se <- sqrt(rowSums(Rinv ^ 2) * sig2)
tstat <- beta / se
pval <- 2 * pt(abs(tstat), n - p, lower.tail = FALSE)
tab <- matrix(c(beta, se, tstat, pval), nrow = p, ncol = 4L,
dimnames = list(dimnames(X)[[2L]],
c("Estimate", "Std. Error", "t value", "Pr(>|t|)")))
## 2 * negative log-likelihood
nega2logLik <- n * log(2 * pi * sig2) + (n - p)
## AIC / BIC
aic <- nega2logLik + 2 * (p + 1)
bic <- nega2logLik + log(n) * (p + 1)
## multiple R-squared and adjusted R-squared
TSS <- c(crossprod(y - sum(y) / n))
r.squared <- 1 - RSS / TSS
adj.r.squared <- 1 - sig2 * (n - 1) / TSS
## return
list(coefficients = beta, residuals = r, fitted.values = c(X %*% beta),
R = R, sig2 = sig2, coef.table = tab, aic = aic, bic = bic, c = c,
RSS = RSS, r.squared = r.squared, adj.r.squared = adj.r.squared)
}
如您所见，它还 return 各种摘要，就好像 summary.lm 已被调用。现在让我们写另一个包装函数choose.c。它根据 c.grid 和 return 绘制了 RSS 的最佳模型，并选择了 c。
choose.c <- function (x, y, c.grid) {
if (is.unsorted(c.grid)) stop("'c.grid' in not increasing")
## model list
lst <- lapply(c.grid, est, x = x, y = y)
## RSS trace
RSS <- sapply(lst, "[[", "RSS")
## verbose
plot(c.grid, RSS, type = "b", pch = 19)
## find `c` / the model minimizing `RSS`
lst[[which.min(RSS)]]
}
到目前为止一切顺利。为了完成这个故事，我们还想要一个predict例程。
pred <- function (model, x.new) {
## prediction matrix
X <- getX(x.new, model$c)
p <- dim(X)[2L]
## predicted mean
fit <- X %*% model$coefficients
## prediction standard error
Qt <- forwardsolve(t(model$R), t(X))
se <- sqrt(colSums(Qt ^ 2) * model$sig2)
## 95%-confidence interval
alpha <- qt(0.025, length(model$residuals) - p)
lwr <- fit + alpha * se
upr <- fit - alpha * se
## return
matrix(c(fit, se, lwr, upr), ncol = 4L,
dimnames = list(NULL, c("fit", "se", "lwr", "upr")))
}
在本节中，我将演示一个可重现的示例。请确保您在其他答案中定义了源函数。
## we first generate a true model
set.seed(0)
x <- runif(100) ## sample points on [0, 1]
beta <- c(0.1, -0.2, 2) ## true coefficients
X <- getX(x, 0.6) ## model matrix with true break point at 0.6
y <- X %*% beta + rnorm(100, 0, 0.08) ## observations with Gaussian noise
plot(x, y)
现在，假设我们不知道 c，我们想在均匀分布的网格上搜索：
c.grid <- seq(0.1, 0.9, 0.05)
fit <- choose.c(x, y, c.grid)
fit$c
RSS 已选择 0.55。这与真实值 0.6 略有不同，但从图中我们看到 RSS 曲线在 [0.5, 0.6] 之间变化不大，所以我对此很满意。
生成的模型 fit 包含丰富的信息：
#List of 12
# $ coefficients : num [1:3] 0.114 -0.246 2.366
# $ residuals : num [1:100] 0.03279 -0.01515 0.21188 -0.06542 0.00763 ...
# $ fitted.values: num [1:100] 0.0292 0.3757 0.2329 0.1087 0.0263 ...
# $ R : num [1:3, 1:3] -10 0.1 0.1 0.292 2.688 ...
# $ sig2 : num 0.00507
# $ coef.table : num [1:3, 1:4] 0.1143 -0.2456 2.3661 0.0096 0.0454 ...
# ..- attr(*, "dimnames")=List of 2
# .. ..$ : chr [1:3] "beta0" "beta1" "beta2"
# .. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)"
# $ aic : num -240
# $ bic : num -243
# $ c : num 0.55
# $ RSS : num 0.492
# $ r.squared : num 0.913
# $ adj.r.squared: num 0.911
我们可以提取系数的摘要table：
fit$coef.table
# Estimate Std. Error t value Pr(>|t|)
#beta0 0.1143132 0.009602697 11.904286 1.120059e-20
#beta1 -0.2455986 0.045409356 -5.408546 4.568506e-07
#beta2 2.3661097 0.169308226 13.975161 5.730682e-25
最后，我们想看一些预测图。
x.new <- seq(0, 1, 0.05)
p <- pred(fit, x.new)
head(p)
# fit se.fit lwr upr
#[1,] 0.9651406 0.02903484 0.9075145 1.0227668
#[2,] 0.8286400 0.02263111 0.7837235 0.8735564
#[3,] 0.7039698 0.01739193 0.6694516 0.7384880
#[4,] 0.5911302 0.01350837 0.5643199 0.6179406
#[5,] 0.4901212 0.01117924 0.4679335 0.5123089
#[6,] 0.4009427 0.01034868 0.3804034 0.4214819
我们可以画个图：
plot(x, y, cex = 0.5)
matlines(x.new, p[,-2], col = c(1,2,2), lty = c(1,2,2), lwd = 2)
李哲源 is a genius but I would like to suggest another solution, using the Heaviside (unit step) 函数，如果 x>0，则 H(x) = 1； H = 0 如果 x ≤ 0
H <- function(x) as.numeric(x>0)
那么，要拟合的函数就是f(x,c) = b0 + b1 (x-c) + b2 (x-c)^2 H(c-x)，可以和nls一起使用：
fit <- nls(y ~ b0+b1*(x-c)+b2*(x-c)^2*H(c-x),
start = list(b0=0,b1=0,b2=1,c=0.5))
用 李哲源 的玩具示例对其进行测试，得到
summary(fit)$parameters
Estimate Std. Error t value Pr(>|t|)
b0 0.1199124 0.03177064 3.774315 2.777969e-04
b1 -0.2578121 0.07856856 -3.281365 1.440945e-03
b2 2.4316379 0.40105205 6.063148 2.624975e-08
c 0.5400831 0.05287111 10.215089 5.136550e-17
我想用一个断点 xt 拟合分段线性回归，这样对于 x < xt 我们有一个二次多项式，对于 x >= xt 我们有一条直线。两块应该平滑连接，连续性达到 xt 处的一阶导数。这是它可能看起来像的图片：
我将分段回归函数参数化为：
其中 a、b、c 和 xt 是要估计的参数。
我想根据调整后的 R 平方将此模型与整个范围内的二次多项式回归进行比较。
这是我的数据：
y <- c(1, 0.59, 0.15, 0.078, 0.02, 0.0047, 0.0019, 1, 0.56, 0.13,
0.025, 0.0051, 0.0016, 0.00091, 1, 0.61, 0.12, 0.026, 0.0067,
0.00085, 4e-04)
x <- c(0, 5.53, 12.92, 16.61, 20.3, 23.07, 24.92, 0, 5.53, 12.92,
16.61, 20.3, 23.07, 24.92, 0, 5.53, 12.92, 16.61, 20.3, 23.07,
24.92)
我的尝试如下，对于已知的 xt:
z <- pmax(0, x - xt)
x1 <- pmin(x, xt)
fit <- lm(y ~ x1 + I(x1 ^ 2) + z - 1)
但是这条直线在xt处似乎并不与二次多项式相切。我哪里做错了？
类似问题：
- Piecewise regression with a straight line and a horizontal line joining at a break point
- Fitting a V-shape curve to my data（交叉验证）
这是一个很好的练习（可能很难）来消化线性模型的理论和实现。我的回答将包含两部分：
- 第 1 部分（这一部分）介绍了我使用的参数化以及这种分段回归如何简化为普通的最小二乘问题。提供模型估计、断点选择、预测等R函数
- 第 2 部分（另一个）演示了一个可复制的玩具示例，说明如何使用我定义的那些函数。
我必须使用不同的参数化，因为你在问题中给出的参数是错误的！你的参数化只保证了函数值的连续性，而不是一阶导数！这就是为什么您的拟合线与 xt.
## generate design matrix
getX <- function (x, c) {
x <- x - c
cbind("beta0" = 1, "beta1" = x, "beta2" = pmin(x, 0) ^ 2)
}
下面的函数 est 在给定的 c:
.lm.fit（为了实现最大效率）
## `x`, `y` give data points; `c` is known break point
est <- function (x, y, c) {
## model matrix
X <- getX(x, c)
p <- dim(X)[2L]
## solve least squares with QR factorization
fit <- .lm.fit(X, y)
## compute Pearson estimate of `sigma ^ 2`
r <- c(fit$residuals)
n <- length(r)
RSS <- c(crossprod(r))
sig2 <- RSS / (n - p)
## coefficients summary table
beta <- fit$coefficients
R <- "dimnames<-"(fit$qr[1:p, ], NULL)
Rinv <- backsolve(R, diag(p))
se <- sqrt(rowSums(Rinv ^ 2) * sig2)
tstat <- beta / se
pval <- 2 * pt(abs(tstat), n - p, lower.tail = FALSE)
tab <- matrix(c(beta, se, tstat, pval), nrow = p, ncol = 4L,
dimnames = list(dimnames(X)[[2L]],
c("Estimate", "Std. Error", "t value", "Pr(>|t|)")))
## 2 * negative log-likelihood
nega2logLik <- n * log(2 * pi * sig2) + (n - p)
## AIC / BIC
aic <- nega2logLik + 2 * (p + 1)
bic <- nega2logLik + log(n) * (p + 1)
## multiple R-squared and adjusted R-squared
TSS <- c(crossprod(y - sum(y) / n))
r.squared <- 1 - RSS / TSS
adj.r.squared <- 1 - sig2 * (n - 1) / TSS
## return
list(coefficients = beta, residuals = r, fitted.values = c(X %*% beta),
R = R, sig2 = sig2, coef.table = tab, aic = aic, bic = bic, c = c,
RSS = RSS, r.squared = r.squared, adj.r.squared = adj.r.squared)
}
如您所见，它还 return 各种摘要，就好像 summary.lm 已被调用。现在让我们写另一个包装函数choose.c。它根据 c.grid 和 return 绘制了 RSS 的最佳模型，并选择了 c。
choose.c <- function (x, y, c.grid) {
if (is.unsorted(c.grid)) stop("'c.grid' in not increasing")
## model list
lst <- lapply(c.grid, est, x = x, y = y)
## RSS trace
RSS <- sapply(lst, "[[", "RSS")
## verbose
plot(c.grid, RSS, type = "b", pch = 19)
## find `c` / the model minimizing `RSS`
lst[[which.min(RSS)]]
}
到目前为止一切顺利。为了完成这个故事，我们还想要一个predict例程。
pred <- function (model, x.new) {
## prediction matrix
X <- getX(x.new, model$c)
p <- dim(X)[2L]
## predicted mean
fit <- X %*% model$coefficients
## prediction standard error
Qt <- forwardsolve(t(model$R), t(X))
se <- sqrt(colSums(Qt ^ 2) * model$sig2)
## 95%-confidence interval
alpha <- qt(0.025, length(model$residuals) - p)
lwr <- fit + alpha * se
upr <- fit - alpha * se
## return
matrix(c(fit, se, lwr, upr), ncol = 4L,
dimnames = list(NULL, c("fit", "se", "lwr", "upr")))
}
在本节中，我将演示一个可重现的示例。请确保您在其他答案中定义了源函数。
## we first generate a true model
set.seed(0)
x <- runif(100) ## sample points on [0, 1]
beta <- c(0.1, -0.2, 2) ## true coefficients
X <- getX(x, 0.6) ## model matrix with true break point at 0.6
y <- X %*% beta + rnorm(100, 0, 0.08) ## observations with Gaussian noise
plot(x, y)
现在，假设我们不知道 c，我们想在均匀分布的网格上搜索：
c.grid <- seq(0.1, 0.9, 0.05)
fit <- choose.c(x, y, c.grid)
fit$c
RSS 已选择 0.55。这与真实值 0.6 略有不同，但从图中我们看到 RSS 曲线在 [0.5, 0.6] 之间变化不大，所以我对此很满意。
生成的模型 fit 包含丰富的信息：
#List of 12
# $ coefficients : num [1:3] 0.114 -0.246 2.366
# $ residuals : num [1:100] 0.03279 -0.01515 0.21188 -0.06542 0.00763 ...
# $ fitted.values: num [1:100] 0.0292 0.3757 0.2329 0.1087 0.0263 ...
# $ R : num [1:3, 1:3] -10 0.1 0.1 0.292 2.688 ...
# $ sig2 : num 0.00507
# $ coef.table : num [1:3, 1:4] 0.1143 -0.2456 2.3661 0.0096 0.0454 ...
# ..- attr(*, "dimnames")=List of 2
# .. ..$ : chr [1:3] "beta0" "beta1" "beta2"
# .. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)"
# $ aic : num -240
# $ bic : num -243
# $ c : num 0.55
# $ RSS : num 0.492
# $ r.squared : num 0.913
# $ adj.r.squared: num 0.911
我们可以提取系数的摘要table：
fit$coef.table
# Estimate Std. Error t value Pr(>|t|)
#beta0 0.1143132 0.009602697 11.904286 1.120059e-20
#beta1 -0.2455986 0.045409356 -5.408546 4.568506e-07
#beta2 2.3661097 0.169308226 13.975161 5.730682e-25
最后，我们想看一些预测图。
x.new <- seq(0, 1, 0.05)
p <- pred(fit, x.new)
head(p)
# fit se.fit lwr upr
#[1,] 0.9651406 0.02903484 0.9075145 1.0227668
#[2,] 0.8286400 0.02263111 0.7837235 0.8735564
#[3,] 0.7039698 0.01739193 0.6694516 0.7384880
#[4,] 0.5911302 0.01350837 0.5643199 0.6179406
#[5,] 0.4901212 0.01117924 0.4679335 0.5123089
#[6,] 0.4009427 0.01034868 0.3804034 0.4214819
我们可以画个图：
plot(x, y, cex = 0.5)
matlines(x.new, p[,-2], col = c(1,2,2), lty = c(1,2,2), lwd = 2)
李哲源 is a genius but I would like to suggest another solution, using the Heaviside (unit step) 函数，如果 x>0，则 H(x) = 1； H = 0 如果 x ≤ 0
H <- function(x) as.numeric(x>0)
那么，要拟合的函数就是f(x,c) = b0 + b1 (x-c) + b2 (x-c)^2 H(c-x)，可以和nls一起使用：
fit <- nls(y ~ b0+b1*(x-c)+b2*(x-c)^2*H(c-x),
start = list(b0=0,b1=0,b2=1,c=0.5))
用 李哲源 的玩具示例对其进行测试，得到
summary(fit)$parameters
Estimate Std. Error t value Pr(>|t|)
b0 0.1199124 0.03177064 3.774315 2.777969e-04
b1 -0.2578121 0.07856856 -3.281365 1.440945e-03
b2 2.4316379 0.40105205 6.063148 2.624975e-08
c 0.5400831 0.05287111 10.215089 5.136550e-17