带有断点的数据的简单回归模型

机器算法验证 变化点 回归不连续
2022-03-17 14:51:55

我目前正在研究一个分段回归模型,其中 ) 该模型应采用以下形式:(xi,yi)i=1..N

yi=β0+β1xi forxi<xcrit

yi=α0+β0+β1xi forxixcrit

其中是必须从数据以及系数确定的断点值。xcritα0,β0β1

我在统计/计量经济学论文(即分段回归模型/回归不连续性设计)中看到了类似的想法,但我想获得一些关于我应该使用哪种模型的反馈。

理想情况下,我想使用一个有据可查的框架,我可以在断点上获得某种置信区间。xcrit

3个回答

这将是一个以 R 为中心的答案。一种方法是将调用包装lm在一个传递断点并以该断点为条件构建回归的函数中,然后通过仅迭代断点的可能值来最小化以断点为条件的拟合模型的偏差。这最大化了断点的配置文件对数似然性,并且通常(即,不仅仅是对于这个问题)如果断点迭代内部的函数(在这种情况下为 lm)找到以传递给它的参数为条件的最大似然估计,则整个过程找到所有参数的联合最大似然估计。

例如:

# True model: y = a + b*(obs. no >= shift) + c*x
# a = 0, b = 1, c = 1, shift at observation 31

# Construct sample data
x <- rnorm(100)
shift <- c(rep(0,30),rep(1,70))
y <- shift + x + rnorm(100)

# Find deviance conditional upon breakpoint
lm.shift <- function(y, x, shift.obs) {
  shift.var <- c(rep(0, (shift.obs-1)), rep(1, length(y)-shift.obs+1))
  deviance(lm(y~x+shift.var))
}

# Find deviance of all breakpoint values 
dev.value <- rep(0, length(y))
for (i in 1:length(y)) {
  dev.value[i] <- lm.shift(y, x, i)
}

# Calculate profile-ll based confidence interval
estimate <- which.min(dev.value)
profile.95.dev <- min(dev.value) + qchisq(0.95,1)
est.lb.95 <- max(which(dev.value[1:estimate] > profile.95.dev))
est.ub.95 <- est -1 + min(which(dev.value[estimate:length(y)] > profile.95.dev))

> estimate
[1] 30
> est.lb.95
[1] 28
> est.ub.95
[1] 33

因此,我们的估计值为 30,95% 的置信区间为 28 - 33。非常紧凑,但相对于误差项的标准偏差,这也是一个相当大的转变。

请注意,在计算基于对数似然的置信区间时会出现一些混乱,但基本思想是找到小于估计值的最大指数,其偏差大于下限的截止水平,最小指数大于估计值偏差大于上限的截止水平。

确实应该绘制出偏差曲线,以确保您没有多个彼此接近的局部最小值,这可能会告诉您一些关于假设模型(或数据)的有趣信息:

在此处输入图像描述

这是检测截距变化的示例(在您的符号中为 B0),有时也称为电平或阶跃。这通常发生在时间序列数据中,其中模型中的变量受到 0,0,0,0,0,0,1,1,1,1,1,1 ......未知的任意点。它被称为干预检测,因为断点(干预)是通过反复试验(即搜索过程)找到(检测到)的。如果您的数据不是时间序列,则可以使用时间序列包来识别干预,同时指定频率“1”并禁用任何 ARMA 结构,从而产生您需要的模型。我的数据是可能预期的时间序列那么您需要在存在 ARIMA 结构和 PDL(ADL' s) 在用户建议的输入/因果序列中。如果您希望发布您的数据,我会向您和列表展示这一点。此外,您可能会查看通用时间序列和/或 www.forecastingsolutions.com/publications/Introducing_cart.pdf的异常值检测

听起来你想要一个带有单个结的样条回归。在 SAS 中,请参见 PROC TRANSREG。在 R 中,请参见(例如)Splines 包。