如何对分段回归系数进行参数化以表示以下间隔的斜率（而不是斜率的变化）

Question

考虑以下数据集

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

[绘制数据，观察值的分布显然是非线性的，但是在数量= 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的系数表示从前一个间隔开始的斜率变化。

如何设置相关系数以代替第二个间隔的斜率？

我进行了一些研究，但除了在统计状态中实现一些自动化之外，我无法找到所需的规格（请参见此处的“边际”语音https://www.stata.com/manuals13/rmkspline.pdf）。

非常感谢您的帮助。谢谢！

Answer 1

如果您知道断点，那么您几乎有了模型，应该是：

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

Answer 2

这里的关键是使用逻辑变量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 ~ I(Xbar * diff) + I((1 - Xbar) * diff), df)

与nls仔细检查

我们可以使用以下公式用nls仔细检查这些内容，它利用了以下事实：如果我们将两条线都延伸，则在任何数量上使用的线都将是两者中的较低者。

st <- list(a = 0, b1 = 1, b2 = -1)
nls(Sales ~ a + pmin(b1 * (Quantity - 89), b2 * (Quantity - 89)), start = st)
## 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

情节

这里是情节：

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