如何将线性模型与阶跃函数相结合答案

【问题标题】：How to combine linear model with step function如何将线性模型与阶跃函数相结合
【发布时间】：2019-07-22 10:25:13
【问题描述】：

假设我们有这个数据：

library(tidyverse)
library(modelr)

set.seed(42)
d1 <- tibble(x =  0:49, y = 5*x  + rnorm(n = 50))
d2 <- tibble(x = 50:99, y = 10*x + rnorm(n = 50))
data <- rbind(d1, d2)   
ggplot(data, aes(x, y)) +
  geom_point()

如何拟合这些数据？

我尝试了什么：

线性模型

m1 <- lm(y ~ x, data = data)
data %>%
  add_predictions(m1) %>%
  gather(key = cat, value = y, -x) %>%
  ggplot(aes(x, y, color = cat)) +
  geom_point()

阶梯函数



# step model
m2 <- lm(y ~ cut(x, 2), data = data)
data %>%
  add_predictions(m2) %>%
  gather(key = cat, value = y, -x) %>%
  ggplot(aes(x, y, color = cat)) +
  geom_point()

如何将两者结合起来？

【问题讨论】：

标签： r glm lm

【解决方案1】：

在数学上，您的模型采用以下形式

    { a_0 + a_1 x  when x < 50
y = {
    { b_0 + b_1 x when x >= 50

您可以将其与指标函数结合起来，得出单线方程的形式：

y = a_0 + (b_0 - a_0) * 1[x >= 50] + a_1 * x + (b_1 - a_1) * x * 1[x >= 50] + error

简单来说，我们可以这样写：

y = c_0 + c_1 * x + c_2 * z + c_3 * x * z + error

我写z = 1[x >= 50] 是为了强调这个指标函数只是另一个回归量

在 R 中，我们可以这样拟合

lm(y ~ x * I(x >= 50), data = data)

* 将根据需要与x 和1[x >= 50] 完全交互。

with(data, {
  plot(x, y)
  reg = lm(y ~ x * I(x >= 50))

  lines(x, predict(reg, data.frame(x)))
})

如果您不知道跳跃发生在 50 岁，那么道路是敞开的，但您可以例如比较均方误差：

x_range = 1:100
errs = sapply(x_range, function(BREAK) {
  mean(lm(y ~ x * I(x >= BREAK), data = data)$residuals^2)
})
plot(x_range, errs)
x_min = x_range[which.min(errs)]
axis(side = 1L, at = x_min)
abline(v = x_min, col = 'red')

【讨论】：