Кусочная регрессия с прямой линией и горизонтальной линией, соединяющей в точке останова

Question

Я хочу сделать кусочно-линейную регрессию с одной точкой останова, где вторая половина линии регрессии имеет slope = 0. Существуют примеры того, как делать кусочно-линейную регрессию, такую ​​как here. Проблема у меня в том, что я не ясно, как зафиксировать угол наклона половины модели равным 0.

Я попытался

lhs <- function(x) ifelse(x < k, k-x, 0) 
rhs <- function(x) ifelse(x < k, 0, x-k) 
fit <- lm(y ~ lhs(x) + rhs(x)) 

где k это точка разрыва, но сегмент на право не является плоским/горизонтальным.

Я хочу, чтобы ограничить наклон второго отрезка в 0. Я попробовал:

fit <- lm(y ~ x * (x < k) + x * (x > k)) 

но опять же, я не уверен, как получить вторую половину, чтобы иметь нулевой наклон.

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

---

Мое собственное решение

У меня есть решение благодаря комментарий ниже. Вот код, который я использую, чтобы оптимизировать и затем построить припадок:

x <- c(1, 2, 3, 1, 2, 1, 6, 1, 2, 3, 2, 1, 4, 3, 1) 
y <- c(0.041754212, 0.083491254, 0.193129615, 0.104249201, 0.17280516, 
0.154342335, 0.303370501, 0.025503008, 0.123934121, 0.191486527, 
0.183958737, 0.156707866, 0.31019215, 0.281890206, 0.25414608) 

range_x <- max(x) - min(x) 
intervals <- 1000 
coef1 <- c() 
coef2 <- c() 
r2 <- c() 

for (i in 1:intervals) { 
    k <- min(x) + (i-1) * (range_x/intervals)  
    x2 = (x - k) * (x < k) 
    fit <- lm(y ~ x2) 
    coef1[i] <- summary(fit)$coef[1] 
    coef2[i] <- summary(fit)$coef[2] 
    r2[i] <- summary(fit)$r.squared 
    } 

best_r2 <- max(r2) # get best r squared 
pos <- which.max(r2)           
best_k <- min(x) + (pos - 1) * (range_x/intervals) 

plot(x, y) 
curve(coef1[pos] - best_k * coef2[pos] + coef2[pos] * x, 
     from=min(x), to=best_k, add = TRUE) 
segments(best_k, coef1[pos], max(x), coef1[pos]) 

Answer 1

Попробуйте сделать переменные вне выражения.

x2 = (x-k)*(x>k) 
lm(y ~ x2) 

В качестве альтернативы, вы можете использовать I()

lm(y~ I((x-k)*(x>k))) 

I() принимает все, что находится внутри буквально и игнорирует другие возможные (MIS) интерпретации с тем, что функция это внутри.

Если у вас нет четко определенного k, вам придется оптимизировать что-то вроде отклонения от разных значений k.

Answer 2

Существует очень похожий поток на переполнение стека: Piecewise regression with a quadratic polynomial and a straight line joining smoothly at a break point.Единственное отличие состоит в том, что мы теперь рассмотрим:

Оказывается, что функции est, choose.c и pred определены в my answer изменять не требуется вообще; нам нужно только изменить getX вернуть дизайн матрицы для кусочно регрессии:

getX <- function (x, c) cbind("beta0" = 1, "beta1" = pmin(x - c, 0)) 

Теперь мы следуем код в toy example, чтобы подобрать модель для ваших данных:

x <- c(1, 2, 3, 1, 2, 1, 6, 1, 2, 3, 2, 1, 4, 3, 1) 
y <- c(0.041754212, 0.083491254, 0.193129615, 0.104249201, 0.17280516, 
0.154342335, 0.303370501, 0.025503008, 0.123934121, 0.191486527, 
0.183958737, 0.156707866, 0.31019215, 0.281890206, 0.25414608) 

x колеблется от 1 6, поэтому мы рассмотрим

c.grid <- seq(1.1, 5.9, 0.05) 
fit <- choose.c(x, y, c.grid) 
fit$c 
# 4.5 

Наконец мы делаем предсказания сюжет:

x.new <- seq(1, 6, by = 0.1) 
p <- pred(fit, x.new) 
plot(x, y, ylim = c(0, 0.4)) 
matlines(x.new, p[,-2], col = c(1,2,2), lty = c(1,2,2), lwd = 2) 

Мы имеем богатую информацию в подогнанной модели:

str(fit) 
#List of 12 
# $ coefficients : num [1:2] 0.304 0.055 
# $ residuals : num [1:15] -0.06981 -0.08307 -0.02844 -0.00731 0.00624 ... 
# $ fitted.values: num [1:15] 0.112 0.167 0.222 0.112 0.167 ... 
# $ R   : num [1:2, 1:2] -3.873 0.258 9.295 -4.37 
# $ sig2   : num 0.00401 
# $ coef.table : num [1:2, 1:4] 0.3041 0.055 0.0384 0.0145 7.917 ... 
# ..- attr(*, "dimnames")=List of 2 
# .. ..$ : chr [1:2] "beta0" "beta1" 
# .. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)" 
# $ aic   : num -34.2 
# $ bic   : num -39.5 
# $ c   : num 4.5 
# $ RSS   : num 0.0521 
# $ r.squared : num 0.526 
# $ adj.r.squared: num 0.49 

Например, мы можем проверять коэффициенты сводную таблицу:

fit$coef.table 
#  Estimate Std. Error t value  Pr(>|t|) 
#beta0 0.30406634 0.03840657 7.917039 2.506043e-06 
#beta1 0.05500095 0.01448188 3.797915 2.216095e-03