Как параметризировать кусочно-регрессионный коэффициент для представления наклона для следующего интервала (вместо изменения наклона)

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).

Любая помощь очень ценится. Спасибо!

Благодарность: рабочий пример получен из https://towardsdatascience.com/unraveling-spline-regression-in-r-937626bc3d96

Answer 1

Ключом здесь является использование логической переменной is.right, которая ИСТИНА для точек справа от 89 и ЛОЖЬ в противном случае.

Из показанных выходных данных 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

Answer 2

Если вы знаете точки останова, то у вас почти есть модель, она должна быть:

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 ***

А коэффициент 2-го наклона равен 59,572 + (- 80,133) = -20,561