结构不报告中断日期
这是我第一次使用结构，所以请耐心等待。我遇到的问题似乎是 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)
时间序列 [1:21] 从 1996 年到 2016 年：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）。
有什么建议？
出现问题是因为breakpoints()目前只能 (a)NA通过省略它们来处理s，以及 (b) 通过ts类处理时间/日期。这会产生冲突，因为当您NA从 a 中省略 internal 时，ts它会失去其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将需要。