结构不报告中断日期

Question

这是我第一次使用结构，所以请多多包涵。我遇到的问题似乎是 strucchange 无法正确识别我的时间序列，但我不知道为什么，也没有在处理这个问题的板上找到答案。这是一个可重现的示例：

require(strucchange)
# time series
nmreprosuccess <- c(0,0.50,NA,0.,NA,0.5,NA,0.50,0.375,0.53,0.846,0.44,1.0,0.285, 
                    0.75,1,0.4,0.916,1,0.769,0.357)
dat.ts <- ts(nmreprosuccess, frequency=1, start=c(1996,1))
str(dat.ts)

从 1996 年到 2016 年的时间序列 [1:21]：0 0.5 NA 0 NA 0.5 NA 0.5 0.375 0.53 ...

对我来说，这意味着时间序列看起来可以使用。

# obtain breakpoints
bp.NMSuccess <- breakpoints(dat.ts~1)
summary(bp.NMSuccess)

这给出了：

Optimal (m+1)-segment partition: 

 Call:
 breakpoints.formula(formula = dat.ts ~ 1)

 Breakpoints at observation number:

 m = 1     6              
 m = 2   3   7            
 m = 3   3           14 16
 m = 4   3   7       14 16
 m = 5   3   7 10    14 16
 m = 6   3   7 10 12 14 16
 m = 7   3 5 7 10 12 14 16

 Corresponding to breakdates:

 m = 1                     0.333333333333333                                                      
 m = 2   0.166666666666667                   0.388888888888889                                    
 m = 3   0.166666666666667                                                                        
 m = 4   0.166666666666667                   0.388888888888889                                    
 m = 5   0.166666666666667                   0.388888888888889 0.555555555555556                  
 m = 6   0.166666666666667                   0.388888888888889 0.555555555555556 0.666666666666667
 m = 7   0.166666666666667 0.277777777777778 0.388888888888889 0.555555555555556 0.666666666666667

 m = 1                                      
 m = 2                                      
 m = 3   0.777777777777778 0.888888888888889
 m = 4   0.777777777777778 0.888888888888889
 m = 5   0.777777777777778 0.888888888888889
 m = 6   0.777777777777778 0.888888888888889
 m = 7   0.777777777777778 0.888888888888889

 Fit:

 m   0       1       2       3       4       5       6       7      
 RSS  1.6986  1.1253  0.9733  0.8984  0.7984  0.7581  0.7248  0.7226
 BIC 14.3728 12.7421 15.9099 20.2490 23.9062 28.7555 33.7276 39.4522

这就是我开始遇到问题的地方。它报告的不是实际的中断日期，而是数字，因此无法将中断线绘制到图表上，因为它们不在中断日期（2002 年），而是在 0.333。

plot.ts(dat.ts, main="Natural Mating")
lines(fitted(bp.NMSuccess, breaks = 1), col = 4, lwd = 1.5)

这张图表中没有任何东西显示出来（我想是因为它对于图表的比例来说太小了）。

此外，当我尝试可能解决此问题的修复程序时，

fm1 <- lm(dat.ts ~ breakfactor(bp.NMSuccess, breaks = 1))

我明白了：

Error in model.frame.default(formula = dat.ts ~ breakfactor(bp.NMSuccess,  : 
  variable lengths differ (found for 'breakfactor(bp.NMSuccess, breaks = 1)')

由于数据中的 NA 值，我收到错误，因此 dat.ts 的长度为 21，breakfactor(bp.NMSuccess, breaks = 1) 的长度为 18（缺少 3 个 NA）。

有什么建议吗？

Answer 1

出现问题是因为breakpoints() 目前只能 (a) 通过省略 NAs 来处理它们，以及 (b) 通过 ts 类处理时间/日期。这会产生冲突，因为当您从 ts 中省略内部 NAs 时，它会丢失其 ts 属性，因此 breakpoints() 无法推断出正确的时间。

解决此问题的“明显”方法是使用可以处理此问题的时间序列类，即zoo。但是，我从来没有将zoo 支持完全集成到breakpoints() 中，因为它可能会破坏当前的一些行为。

长话短说：您目前最好的选择是自己记账，而不是指望breakpoints() 为您做这件事。额外的工作并没有那么大。首先，我们创建一个包含响应和时间向量的时间序列，并省略 NAs：

d <- na.omit(data.frame(success = nmreprosuccess, time = 1996:2016))
d
##    success time
## 1    0.000 1996
## 2    0.500 1997
## 4    0.000 1999
## 6    0.500 2001
## 8    0.500 2003
## 9    0.375 2004
## 10   0.530 2005
## 11   0.846 2006
## 12   0.440 2007
## 13   1.000 2008
## 14   0.285 2009
## 15   0.750 2010
## 16   1.000 2011
## 17   0.400 2012
## 18   0.916 2013
## 19   1.000 2014
## 20   0.769 2015
## 21   0.357 2016

然后我们可以估计断点，然后将观察的“数量”转换回时间尺度。请注意，我在这里明确设置了最小段大小h，因为对于这个简短的系列来说，默认值 15% 可能有点小。 4 仍然很小，但可能足以估计一个恒定的平均值。

bp <- breakpoints(success ~ 1, data = d, h = 4)
bp
##   Optimal 2-segment partition: 
## 
## Call:
## breakpoints.formula(formula = success ~ 1, h = 4, data = d)
## 
## Breakpoints at observation number:
## 6 
## 
## Corresponding to breakdates:
## 0.3333333 

我们忽略了 1/3 观测值的中断“日期”，而是简单地映射回原始时间尺度：

d$time[bp$breakpoints]
## [1] 2004

要使用格式良好的因子水平重新估计模型，我们可以这样做：

lab <- c(
  paste(d$time[c(1, bp$breakpoints)], collapse = "-"),
  paste(d$time[c(bp$breakpoints + 1, nrow(d))], collapse = "-")
)
d$seg <- breakfactor(bp, labels = lab)
lm(success ~ 0 + seg, data = d)
## Call:
## lm(formula = success ~ 0 + seg, data = d)
## 
## Coefficients:
## seg1996-2004  seg2005-2016  
##       0.3125        0.6911  

或用于可视化：

plot(success ~ time, data = d, type = "b")
lines(fitted(bp) ~ time, data = d, col = 4, lwd = 2)
abline(v = d$time[bp$breakpoints], lty = 2)

最后一点：对于如此短的时间序列，只需要简单地改变均值，还可以考虑条件推理（也称为置换测试），而不是 strucchange 中使用的渐近推理。 coin 包正好为此目的提供了 maxstat_test() 函数（= 测试均值单次偏移的短系列）。

library("coin")
maxstat_test(success ~ time, data = d, dist = approximate(99999))
##  Approximative Generalized Maximally Selected Statistics
## 
## data:  success by time
## maxT = 2.3953, p-value = 0.09382
## alternative hypothesis: two.sided
## sample estimates:
##   "best" cutpoint: <= 2004

这会找到相同的断点并提供置换测试 p 值。但是，如果有更多数据并且需要多个断点和/或更多回归系数，则需要strucchange。