机器算法验证 - 测试时间序列中干预效果的显着性 - 吾爱随笔录

测试时间序列中干预效果的显着性

机器算法验证 r 时间序列干预分析动态回归分段线性

2022-04-13 02:29:22

我正在寻找最好的方法来测试在已知时间发生的干预对时间序列数据的影响的重要性。

以玩具数据集为例，我提出了两种方法。

数据

y <- c(rnorm(10, 10, 0.12), 9.6, 9.4, 9.3, 9.2, 9.15)
x <- seq(1:15)
df <- data.frame(y = y, x = x)

ggplot(df, aes(x,y)) + geom_point() +
  geom_vline(xintercept = 10.5) +
  scale_x_continuous(breaks=df$x)

1. 分段回归

将两个线性回归模型拟合到干预前后的数据子集。

df1 <- subset(df, x <= 10)
m1 <- lm(y ~ x, data = df1) 
summary(m1) #Obviously non-significant

df2 <- subset(df, x > 10)
m2 <- lm(y ~ x, data = df2)
summary(m2) #Obviously significant

使用公式来比较这个答案的斜率。

b1 <- coef(summary(m1))[2, 1]
b2 <- coef(summary(m2))[2, 1]
SEb1 <- coef(summary(m1))[2, 2]
SEb2 <- coef(summary(m2))[2, 2]

Z <- (b1-b2)/sqrt(SEb1^2+SEb2^2)

并计算相应的P值。

2*pnorm(-abs(Z))
[1] 1.395998e-08

（顺便问一下，有没有更优雅的功能可以做到以上几点？）

该 P 值非常显着，是要报告干预效果的值。

结果通过绘制前后两条回归线以图形方式显示。（由于lm表明处的关系斜率与x=1:100 没有显着差异，因此线位于处y=mean(1:10)）

ggplot(df, aes(x,y)) + geom_point() +
  geom_vline(xintercept = 10.5) +
  scale_x_continuous(breaks=df$x) +
  stat_smooth(method="lm", data=df2, se=F, colour="royalblue1", size = 0.75) +
  geom_segment(x = 1, xend = 10, y = mean(df1$y), yend = mean(df1$y),colour="royalblue1", size=0.75)

2. 使用 ARIMA 进行动态回归

拟合两个 ARIMA 模型，一个没有，一个带有回归量，用于编码干预。

library(forecast)
y <- ts(y)
intervention <- c(rep(0,10), rep(1,5))

a1 <- auto.arima(y)
summary(a1)

总结表明auto.arima选择 ARIMA(0,1,0) 作为最佳模型。因此，使用该函数将 ARIMA(0,1,0) 与回归量拟合Arima。

a2 <- Arima(y, order=c(0,1,0), xreg=intervention)

然后使用 LRT 检验比较两个模型，以获得与干预效果相关的 P 值。

library(lmtest)
lrtest(a1, a2)

显然，P 值非常显着。

动态回归的一个优点是它可以用来获得预测。

intf <- c(rep(1,5))
autoplot(forecast(a2, h=5, xreg=intf))

问题

这两种方法是否充分且充分执行？
还有其他方法吗？
哪种方法是首选？

1个回答

您所指的称为结构更改/中断测试或更改点模型。

由于您有一个已知的更改日期，您可以简单地在模型中添加一个交互，并对该系数使用标准 t 检验。如果您测试更多系数，请使用Chow 测试公式（例如，参见这篇文章）。这适用于线性回归或 ARIMA。所以这两个模型之间的选择应该基于一般的考虑，并且不受你想要做测试的影响。

所以你只需建模：

$y_i = \alpha + \alpha^+ 1(x>10) +\beta x_i + \beta^+x_i1(x>10) + \epsilon_i$

和的（联合）显着性进行测试。快速示例，没有 Chow 测试，仅查看单个 coefs： $\beta^+$ $\alpha^+$

y <- c(rnorm(10, 10, 0.12), 9.6, 9.4, 9.3, 9.2, 9.15)
x <- seq(1:15)
df <- data.frame(y = y, x = x)
df$D <- ifelse(x<10, 0,1)

reg <- lm(y~1+D+ x*(1-D)+x*(D), data=df)

然后该summary方法将显示您的系数是否显着：

summary(reg)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.82048    0.08167 120.241  < 2e-16 ***
D            1.62439    0.34875   4.658 0.000696 ***
x            0.03889    0.01451   2.680 0.021417 *  
D:x         -0.19901    0.03054  -6.516 4.33e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

DD:x是和\系数。因此，您可以检查三个测试： $\alpha^+$ $\beta^+$

$H_0: \alpha^+=0$ ，即中断拦截。在这里，的 p 值非常低D，因此被拒绝，因此截距中断。
$H_0: \beta^+=0$ ，即斜率中断。也非常低的 p 值D:x。
$H_0: \alpha^+ \beta^+=0$ 所有系数都没有中断的空值（联合测试）。

为此，请linearHypothesis如下使用。也被拒绝了。

library(car)
linearHypothesis(reg, c("D=0", "D:x=0"))

最后，绘图非常简单：

plot(y~x, data=df)
lines(predict(reg))

总的来说，看strucchangeR中的包，确实不错，也更通用（允许自己搜索日期/更改值）。例如：

library(strucchange)

breakpoints(y~1+x, data=df, h=0.2, breaks=1)
sctest(Fstats(y~1+x, data=df, from=0.2, to=0.2))

第一个估计断点（10，我相信它对应于你的 11），第二个测试恒常性。

其它你可能感兴趣的问题

上一篇围绕 Medoids (PAM) 进行分区是确定性的吗？下一篇马尔可夫随机场与隐马尔可夫模型