如何参数化分段回归系数以表示以下间隔的斜率（而不是斜率的变化）
How to parametrize piecewise regression coefficient to represent the slope for the following interval (instead of the change in the slope)
考虑以下数据集
Quantity <- c(25,39,45,57,70,85,89,100,110,124,137,150,177)
Sales <- c(1000,1250,2600,3000,3500,4500,5000,4700,4405,4000,3730,3400,3300)
df <- data.frame(Quantity,Sales)
df
绘制数据，观测值的分布显然是非线性的，但在 Quantity = 89 附近呈现出可能的断点（我在此处跳过该图）。因此，我建立了一个联合分段线性模型如下
df$Xbar <- ifelse(df$Quantity>89,1,0)
df$diff <- df$Quantity - 89
reg <- lm(Sales ~ Quantity + I(Xbar * (Quantity - 89)), data = df)
summary(reg)
或者干脆
df$X <- df$diff*df$Xbar
reg <- lm(Sales ~ Quantity + X, data = df)
summary(reg)
然而，根据这个参数化，X 的系数表示斜率相对于前一个区间的变化。
如何将相关系数参数化以表示第二个区间的斜率？
我做了一些研究，但除了 stata 中的一些自动化之外，我无法找到所需的规范（请参阅此处的语音 'marginal' https://www.stata.com/manuals13/rmkspline.pdf）。
非常感谢任何帮助。谢谢！
致谢：
可行的例子是从
https://towardsdatascience.com/unraveling-spline-regression-in-r-937626bc3d96
如果你知道断点，那么模型就差不多了，应该是：
fit=lm(Sales ~ Quantity + Xbar + Quantity:Xbar,data=df)
因为如果你不引入新的截距（Xbar），它会从模型中已有的截距开始，这是行不通的。我们可以绘制它：
plot(df$Quantity,df$Sales)
newdata = data.frame(Quantity=seq(40,200,by=5))
newdata$Xbar= ifelse(newdata$Quantity>89,1,0)
lines(newdata$Quantity,predict(fit,newdata))
系数为：
summary(fit)
Call:
lm(formula = Sales ~ Quantity * Xbar, data = df)
Residuals:
Min 1Q Median 3Q Max
-527.9 -132.2 -15.1 148.1 464.7
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -545.435 327.977 -1.663 0.131
Quantity 59.572 5.746 10.367 2.65e-06 ***
Xbar 7227.288 585.933 12.335 6.09e-07 ***
Quantity:Xbar -80.133 6.856 -11.688 9.64e-07 ***
第二个斜率的系数是59.572+(-80.133) = -20.561
这里的关键是使用一个逻辑变量is.right，它对于 89 右边的点为 TRUE，否则为 FALSE。
从显示的输出来看，60.88 是 89 左侧的斜率，-19.97 是右侧的斜率。这些线相交于数量 = 89，销售额 = 4817.30。
is.right <- df$Quantity > 89
fm <- lm(Sales ~ diff : is.right, df)
fm
## Call:
## lm(formula = Sales ~ diff:is.right, data = df)
##
## Coefficients:
## (Intercept) diff:is.rightFALSE diff:is.rightTRUE
## 4817.30 60.88 -19.97
备选方案
或者，如果您想使用问题中的 Xbar，请按此方式进行。它给出与 fm.
相同的系数
fm2 <- lm(Sales ~ diff : factor(Xbar), df)
或
fm3 <- lm(Sales ~ I(Xbar * diff) + I((1 - Xbar) * diff), df)
用 nls 仔细检查
我们可以使用 nls 和以下公式仔细检查这些，它利用了这样一个事实，即如果我们扩展两条线，则在任何数量处使用的那条线是两条线中较小的一条。
st <- list(a = 0, b1 = 1, b2 = -1)
fm4 <- nls(Sales ~ a + pmin(b1 * (Quantity - 89), b2 * (Quantity - 89)), start = st)
fm4
## Nonlinear regression model
## model: Sales ~ a + pmin(b1 * (Quantity - 89), b2 * (Quantity - 89))
## data: parent.frame()
## a b1 b2
## 4817.30 60.88 -19.97
## residual sum-of-squares: 713120
##
## Number of iterations to convergence: 1
## Achieved convergence tolerance: 2.285e-09
这也行得通：
fm5 <- nls(Sales ~ a + ifelse(Quantity > 89, b2, b1) * diff, df, start = st)
情节
这是一个情节：
plot(Sales ~ Quantity, df)
lines(fitted(fm) ~ Quantity, df)
模型矩阵
这里是线性回归的模型矩阵：
> model.matrix(fm)
(Intercept) diff:is.rightFALSE diff:is.rightTRUE
1 1 -64 0
2 1 -50 0
3 1 -44 0
4 1 -32 0
5 1 -19 0
6 1 -4 0
7 1 0 0
8 1 0 11
9 1 0 21
10 1 0 35
11 1 0 48
12 1 0 61
13 1 0 88
考虑以下数据集
Quantity <- c(25,39,45,57,70,85,89,100,110,124,137,150,177)
Sales <- c(1000,1250,2600,3000,3500,4500,5000,4700,4405,4000,3730,3400,3300)
df <- data.frame(Quantity,Sales)
df
绘制数据，观测值的分布显然是非线性的，但在 Quantity = 89 附近呈现出可能的断点（我在此处跳过该图）。因此，我建立了一个联合分段线性模型如下
df$Xbar <- ifelse(df$Quantity>89,1,0)
df$diff <- df$Quantity - 89
reg <- lm(Sales ~ Quantity + I(Xbar * (Quantity - 89)), data = df)
summary(reg)
或者干脆
df$X <- df$diff*df$Xbar
reg <- lm(Sales ~ Quantity + X, data = df)
summary(reg)
然而，根据这个参数化，X 的系数表示斜率相对于前一个区间的变化。
如何将相关系数参数化以表示第二个区间的斜率？
我做了一些研究，但除了 stata 中的一些自动化之外，我无法找到所需的规范（请参阅此处的语音 'marginal' https://www.stata.com/manuals13/rmkspline.pdf）。
非常感谢任何帮助。谢谢！
致谢： 可行的例子是从 https://towardsdatascience.com/unraveling-spline-regression-in-r-937626bc3d96
如果你知道断点，那么模型就差不多了，应该是：
fit=lm(Sales ~ Quantity + Xbar + Quantity:Xbar,data=df)
因为如果你不引入新的截距（Xbar），它会从模型中已有的截距开始，这是行不通的。我们可以绘制它：
plot(df$Quantity,df$Sales)
newdata = data.frame(Quantity=seq(40,200,by=5))
newdata$Xbar= ifelse(newdata$Quantity>89,1,0)
lines(newdata$Quantity,predict(fit,newdata))
系数为：
summary(fit)
Call:
lm(formula = Sales ~ Quantity * Xbar, data = df)
Residuals:
Min 1Q Median 3Q Max
-527.9 -132.2 -15.1 148.1 464.7
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -545.435 327.977 -1.663 0.131
Quantity 59.572 5.746 10.367 2.65e-06 ***
Xbar 7227.288 585.933 12.335 6.09e-07 ***
Quantity:Xbar -80.133 6.856 -11.688 9.64e-07 ***
第二个斜率的系数是59.572+(-80.133) = -20.561
这里的关键是使用一个逻辑变量is.right，它对于 89 右边的点为 TRUE，否则为 FALSE。
从显示的输出来看，60.88 是 89 左侧的斜率，-19.97 是右侧的斜率。这些线相交于数量 = 89，销售额 = 4817.30。
is.right <- df$Quantity > 89
fm <- lm(Sales ~ diff : is.right, df)
fm
## Call:
## lm(formula = Sales ~ diff:is.right, data = df)
##
## Coefficients:
## (Intercept) diff:is.rightFALSE diff:is.rightTRUE
## 4817.30 60.88 -19.97
备选方案
或者，如果您想使用问题中的 Xbar，请按此方式进行。它给出与 fm.
fm2 <- lm(Sales ~ diff : factor(Xbar), df)
或
fm3 <- lm(Sales ~ I(Xbar * diff) + I((1 - Xbar) * diff), df)
用 nls 仔细检查
我们可以使用 nls 和以下公式仔细检查这些，它利用了这样一个事实，即如果我们扩展两条线，则在任何数量处使用的那条线是两条线中较小的一条。
st <- list(a = 0, b1 = 1, b2 = -1)
fm4 <- nls(Sales ~ a + pmin(b1 * (Quantity - 89), b2 * (Quantity - 89)), start = st)
fm4
## Nonlinear regression model
## model: Sales ~ a + pmin(b1 * (Quantity - 89), b2 * (Quantity - 89))
## data: parent.frame()
## a b1 b2
## 4817.30 60.88 -19.97
## residual sum-of-squares: 713120
##
## Number of iterations to convergence: 1
## Achieved convergence tolerance: 2.285e-09
这也行得通：
fm5 <- nls(Sales ~ a + ifelse(Quantity > 89, b2, b1) * diff, df, start = st)
情节
这是一个情节：
plot(Sales ~ Quantity, df)
lines(fitted(fm) ~ Quantity, df)
模型矩阵
这里是线性回归的模型矩阵：
> model.matrix(fm)
(Intercept) diff:is.rightFALSE diff:is.rightTRUE
1 1 -64 0
2 1 -50 0
3 1 -44 0
4 1 -32 0
5 1 -19 0
6 1 -4 0
7 1 0 0
8 1 0 11
9 1 0 21
10 1 0 35
11 1 0 48
12 1 0 61
13 1 0 88