R中中断时间序列分析的资源

机器算法验证 r 时间序列
2022-03-15 23:01:27

我对 R 相当陌生。我试图阅读时间序列分析并且已经完成

  1. Shumway 和 Stoffer 的时间序列分析及其应用第 3 版，
  2. 海德曼的优秀预测：原则与实践
  3. Avril Coghlan使用 R 进行时间序列分析
  4. A. Ian McLeod 等人使用 R 进行时间序列分析
  5. Marcel Dettling 博士的应用时间序列分析

编辑：我不确定如何处理这个问题，但我在交叉验证之外找到了一个有用的资源。我想把它包括在这里，以防有人偶然发现这个问题。

药物使用研究中中断时间序列研究的分段回归分析

我有一个连续 7 年每天测量的消耗物品数量（计数数据）的单变量时间序列。大约在时间序列的中间对研究人群进行了干预。预计这种干预不会产生立竿见影的效果，而且效果开始的时间基本上是不可知的。

使用 Hyndman 的forecast软件包，我使用auto.arima(). 但我不确定如何使用这种拟合来回答趋势是否发生了统计上的显着变化并量化数量。

# for simplification I will aggregate to monthly counts
# I can later generalize any teachings the community supplies
count <- c(2464, 2683, 2426, 2258, 1950, 1548, 1108,  991, 1616, 1809, 1688, 2168, 2226, 2379, 2211, 1925, 1998, 1740, 1305,  924, 1487, 1792, 1485, 1701, 1962, 2896, 2862, 2051, 1776, 1358, 1110,  939, 1446, 1550, 1809, 2370, 2401, 2641, 2301, 1902, 2056, 1798, 1198,  994, 1507, 1604, 1761, 2080, 2069, 2279, 2290, 1758, 1850, 1598, 1032,  916, 1428, 1708, 2067, 2626, 2194, 2046, 1905, 1712, 1672, 1473, 1052,  874, 1358, 1694, 1875, 2220, 2141, 2129, 1920, 1595, 1445, 1308, 1039,  828, 1724, 2045, 1715, 1840)
# for explanatory purposes
# month <- rep(month.name, 7)
# year <- 1999:2005
ts <- ts(count, start(1999, 1))
train_month <- window(ts, start=c(1999,1), end = c(2001,1))
require(forecast)
arima_train <- auto.arima(train_month)
fit_month <- Arima(train_month, order = c(2,0,0), seasonal = c(1,1,0), lambda = 0)
plot(forecast(fit_month, 36)); lines(ts, col="red")

是否有任何资源专门处理 R 中的中断时间序列分析？我发现这在 SPSS 中处理 ITS，但我无法将其翻译成 R。

3个回答

这称为变化点分析。R 包changepoint可以为您执行此操作：请参阅此处的文档（包括参考文献）：http ://www.lancs.ac.uk/~killick/Pub/KillickEckley2011.pdf

我建议使用重复测量层次模型。这种方法应该提供稳健的结果，因为每个人都将作为他/她自己的控制。尝试从 UCLA查看 此链接。

对于贝叶斯方法，您可以使用mcp对残差应用自回归（在对数空间中）来拟合泊松或二项式模型（因为您有来自固定间隔周期的计数）。然后使用交叉验证将两段模型与单段模型进行比较。

在我们开始之前，请注意，对于这个数据集，这个模型不太适合，交叉验证看起来很不稳定。因此，我会避免在高风险场景中使用以下内容，但它说明了一种通用方法：

# Fit the change point model
library(mcp)
model_full = list(
  count ~ 1 + ar(1),  # intercept and AR(1)
  ~ 1  # New intercept
)
fit_full = mcp(model_full, data = df, family = poisson(), par_x = "year")


# Fit the null model
model_null = list(
  count ~ 1 + ar(1)  # just a stable AR(1)
)
fit_null = mcp(model_null, data = df, family = poisson(), par_x = "year")

# Compare predictive performance using LOO cross-validation
fit_full$loo = loo(fit_full)
fit_null$loo = loo(fit_null)
loo::loo_compare(fit_full$loo, fit_null$loo)

对于当前数据集，这导致

       elpd_diff se_diff
model2    0.0       0.0 
model1 -459.1      64.3

即，有elpd_diff/se_diff利于空模型的比率约为 7（无变化）。可能的改进包括：

  • sin()使用或对周期性趋势建模cos()。
  • 添加有关更改的可能位置的先验信息，例如prior = list(cp_1 = dnorm(1999.8, 0.5)。

在网站上阅读有关建模自回归、进行模型比较和设置先验的mcp更多信息。披露：我是mcp.

其它你可能感兴趣的问题
上一篇什么是普通的，在普通的最小二乘法中？ 下一篇如何在校准概率模型时选择最佳 bin 宽度？