不规则时间序列的滚动回归答案

【问题标题】：Rolling regression on irregular time series不规则时间序列的滚动回归
【发布时间】：2018-04-02 06:55:06
【问题描述】：

摘要 (tldr)

我需要对 不规则 时间序列（即间隔可能甚至不是周期性的，从 0, 1, 2, 3... 到 ...7, 20, 24, 28...）执行滚动回归，这是简单的数字，不一定需要日期/时间，但滚动窗口需要按时间。因此，如果我有一个不规则采样 600 秒且窗口为 30 的时间序列，则每 30 秒执行一次回归，而不是每 30 秒执行一次回归。

我已经阅读了示例，虽然我可以按时间复制滚动总和和中位数，但我似乎无法为回归计算出来。

问题

首先，我阅读了其他一些关于对不规则时间序列数据执行滚动函数的问题，例如：optimized rolling functions on irregular time series with time-based window，以及：Rolling window over irregular time series。

问题在于，到目前为止，所提供的示例对于sum 或median 等方程很简单，但我还没有弄清楚如何执行简单的滚动回归，即使用lm，即仍然基于相同的警告，即窗口基于不规则的时间序列。而且，我的时间序列要简单得多。不需要日期，只是“经过”的时间。

无论如何，正确处理对我来说很重要，因为在不规则时间（例如，时间间隔中的跳过）可能会高估或低估滚动回归中的系数，因为样本窗口将包括 额外时间。

所以我想知道是否有人可以帮助我创建一个以最简单的方式执行此操作的函数？该数据集基于随时间测量一个变量，即 2 个变量：time 和 response。时间以每 x 个经过的时间单位（秒、分钟，因此不是 日期/时间 格式）测量一次，但偶尔会变得不规则。

对于函数中的每一行，它应该基于 n 个时间单位的宽度执行线性回归。宽度不得超过 n 个单位，但可以下限（即减小）以适应不规则的时间采样。例如，如果宽度指定为 20 秒，但时间每 6 秒采样一次，则窗口将四舍五入为 18，而不是 24 秒。

我在这里查看了问题：How to calculate the average slope within a moving window in R，并且我在不规则的时间序列上测试了该代码，但它看起来像是基于规则的时间序列。

样本数据：

sample <- 
structure(list(x = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 
29, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 47, 48, 
49), y = c(50, 49, 48, 47, 46, 47, 46, 45, 44, 43, 44, 43, 42, 
41, 40, 41, 40, 39, 38, 37, 38, 37, 36, 35, 34, 35, 34, 33, 32, 
31, 30, 29, 28, 29, 28, 27, 26, 25, 26, 25, 24, 23, 22, 21, 20, 
19)), .Names = c("x", "y"), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -46L))

我当前的代码（基于我之前提到的问题）。我知道它不是按时间划分的：

library(zoo)
clm <- function(z) coef(lm(y ~ x, as.data.frame(z)))
rollme <- rollapplyr(zoo(sample), 10, clm, by.column = F, fill = NA)

预期输出（手动计算）如下。输出不同于常规的滚动回归——只要时间间隔跳过 29（秒），数字就会不同：

    NA
    NA
    NA
    NA
    NA
    NA
    NA
    NA
    NA
    -0.696969697
    -0.6
    -0.551515152
    -0.551515152
    -0.6
    -0.696969697
    -0.6
    -0.551515152
    -0.551515152
    -0.6
    -0.696969697
    -0.6
    -0.551515152
    -0.551515152
    -0.6
    -0.696969697
    -0.6
    -0.551515152
    -0.551515152
    -0.6
    -0.696969697
    -0.605042017
    -0.638888889
    -0.716981132
    -0.597560976
    -0.528301887
    -0.5
    -0.521008403
    -0.642857143
    -0.566666667
    -0.551515152
    -0.551515152
    -0.6
    -0.696969697
    -0.605042017
    -0.638888889
    -0.716981132

我希望我提供了足够的信息，但让我知道（或给我一个指导某个地方的好例子）让我试试这个？

我尝试过的其他事情：我尝试将时间转换为 POSIXct 格式，但我不知道如何执行 lm：

require(lubridate)    
x <- as.POSIXct(strptime(sample$x, format = "%S"))

更新：添加了 tldr 部分。

【问题讨论】：

明确地说，任务是在协变量 x 上回归时间 y，超过一个滚动时间窗口，比如 20 个单位，时间差不相等。
检查您发布的代码的前几行输出，它给出的斜率系数与您按预期列出的一样。请准确说明问题所在。
对不起，我以为我说得够清楚了（但也许不是）。我会尽快编辑并澄清问题。

标签： r iteration linear-regression zoo

【解决方案1】：

试试这个：

# time interval is 1    
sz=10
    pl2=list()
    for ( i in 1:nrow(sample)){
      if (i<sz) period=sz else
      period=length(sample$x[sample$x>(sample$x[i]-sz) & sample$x<=sample$x[i]])-1
      pl2[[i]]=seq(-period,0)
    }

#update for time interval > 1
sz=10
tint=1
pl2=list()
for ( i in 1:nrow(sample)){
  if (i<sz) period=sz else
  period=length(sample$x[sample$x>(sample$x[i]-sz*tint) & sample$x<=sample$x[i]])-1
  pl2[[i]]=seq(-period,0)
}

rollme3 <- rollapplyr(zoo(sample), pl2, clm, by.column = F, fill = NA)

> tail(rollme3)
   (Intercept)          x
41    47.38182 -0.5515152
42    49.20000 -0.6000000
43    53.03030 -0.6969697
44    49.26050 -0.6050420
45    50.72222 -0.6388889
46    54.22642 -0.7169811

【讨论】：

嗨罗伯特，谢谢！当时间间隔为 1 时效果很好，但大于该时间间隔时（例如每 5 秒）会出现错误：Error in eval(predvars, data, env) : object 'y' not found。我想知道这是否是代码中的一个简单修复......我正在阅读它并试图理解它。
查看更新，定义tint =5，看看它是否有效。
嗨罗伯特，谢谢！我将在更激进的不规则时间序列上对其进行测试，并尽快回复您。

【解决方案2】：

为了完整起见，这里有一个答案，它使用data.table在非等值连接中聚合。

尽管有许多类似的问题，例如 r calculating rolling average with window based on value (not number of rows or date/time variable)，但这个问题值得自己回答，因为 OP 正在寻找 滚动回归的系数。

library(data.table)
ws <- 10   # size of sliding window in time units
setDT(sample)[.(start = x - ws, end = x), on = .(x > start, x <= end),
              as.list(coef(lm(y ~ x.x))), by = .EACHI]

      x  x (Intercept)        x.x
 1: -10  0    50.00000         NA
 2:  -9  1    50.00000 -1.0000000
 3:  -8  2    50.00000 -1.0000000
 4:  -7  3    50.00000 -1.0000000
 5:  -6  4    50.00000 -1.0000000
 6:  -5  5    49.61905 -0.7142857
 7:  -4  6    49.50000 -0.6428571
 8:  -3  7    49.50000 -0.6428571
 9:  -2  8    49.55556 -0.6666667
10:  -1  9    49.63636 -0.6969697
11:   0 10    49.20000 -0.6000000
12:   1 11    48.88485 -0.5515152
13:   2 12    48.83636 -0.5515152
14:   3 13    49.20000 -0.6000000
15:   4 14    50.12121 -0.6969697
16:   5 15    49.20000 -0.6000000
17:   6 16    48.64242 -0.5515152
18:   7 17    48.59394 -0.5515152
19:   8 18    49.20000 -0.6000000
20:   9 19    50.60606 -0.6969697
21:  10 20    49.20000 -0.6000000
22:  11 21    48.40000 -0.5515152
23:  12 22    48.35152 -0.5515152
24:  13 23    49.20000 -0.6000000
25:  14 24    51.09091 -0.6969697
26:  15 25    49.20000 -0.6000000
27:  16 26    48.15758 -0.5515152
28:  17 27    48.10909 -0.5515152
29:  18 28    49.20000 -0.6000000
30:  19 29    51.57576 -0.6969697
31:  22 32    49.18487 -0.6050420
32:  23 33    50.13889 -0.6388889
33:  24 34    52.47170 -0.7169811
34:  25 35    48.97561 -0.5975610
35:  26 36    46.77358 -0.5283019
36:  27 37    45.75000 -0.5000000
37:  28 38    46.34454 -0.5210084
38:  29 39    50.57143 -0.6428571
39:  30 40    47.95556 -0.5666667
40:  31 41    47.43030 -0.5515152
41:  32 42    47.38182 -0.5515152
42:  33 43    49.20000 -0.6000000
43:  34 44    53.03030 -0.6969697
44:  37 47    49.26050 -0.6050420
45:  38 48    50.72222 -0.6388889
46:  39 49    54.22642 -0.7169811
      x  x (Intercept)        x.x

请注意，时间序列定期间隔的第 10 到 30 行与 OP 的 rollme 相同。

对as.list() 的调用会强制coef(lm(...)) 的结果出现在不同的列中。

上面的代码使用了一个右对齐的滚动窗口。但是，该代码也可以很容易地调整为支持左对齐窗口：

# left aligned window
setDT(sample)[.(start = x, end = x + ws), on = .(x >= start, x < end),
              as.list(coef(lm(y ~ x.x))), by = .EACHI]

【讨论】：

有没有办法让它与 GLS 回归一起工作？ setDT(sample)[.(start = x, end = x + ws), on = .(x >= start, x < end), as.list(coef(gls(y ~ x.x))), by = .EACHI] 抛出关于无法找到的错误 y Error in eval(predvars, data, env) : object 'y' not found

【解决方案3】：

runner 可以在不规则时间序列中应用任何 R 函数。用户必须将数据指定为x 参数，并将日期向量指定为idx 参数（以使窗口与时间相关）。窗口宽度 k 可以是整数 k = 30 或类似 seq.POSIXt k = "30 secs" 中的字符。

第一个示例展示了如何从 lm 函数获取两个参数 - 输出将是一个矩阵

library(runner)

runner(
  x = sample,
  k = "30 secs",
  idx = sample$datetime,
  function(x) {
    coefficients(lm(y ~ x, data = x))
  }
)

或者可以为每个参数单独执行runner

library(runner)

sample$intercept <- runner(
  sample,
  k = "30 secs",
  idx = sample$datetime,
  function(x) {
    coefficients(lm(y ~ x, data = x))[1]
  }
)

sample$slope <- runner(
  sample,
  k = "30 secs",
  idx = sample$datetime,
  function(x) {
    coefficients(lm(y ~ x, data = x))[2]
  }
)

head(sample, 15)

#               datetime  x  y intercept      slope
# 1  2020-04-13 09:27:20  0 50  50.00000         NA
# 2  2020-04-13 09:27:21  1 49  50.00000 -1.0000000
# 3  2020-04-13 09:27:25  2 48  50.00000 -1.0000000
# 4  2020-04-13 09:27:29  3 47  50.00000 -1.0000000
# 5  2020-04-13 09:27:29  4 46  50.00000 -1.0000000
# 6  2020-04-13 09:27:32  5 47  49.61905 -0.7142857
# 7  2020-04-13 09:27:34  6 46  49.50000 -0.6428571
# 8  2020-04-13 09:27:38  7 45  49.50000 -0.6428571
# 9  2020-04-13 09:27:38  8 44  49.55556 -0.6666667
# 10 2020-04-13 09:27:41  9 43  49.63636 -0.6969697
# 11 2020-04-13 09:27:44 10 44  49.45455 -0.6363636
# 12 2020-04-13 09:27:47 11 43  49.38462 -0.6153846
# 13 2020-04-13 09:27:48 12 42  49.38462 -0.6153846
# 14 2020-04-13 09:27:49 13 41  49.42857 -0.6263736
# 15 2020-04-13 09:27:50 14 40  49.34066 -0.6263736

日期时间列的数据

sample <- structure(
  list(
    datetime = c(3, 1, 4, 4, 0, 3, 2, 4, 0, 3, 3, 3, 1, 1, 1, 3, 0, 2, 4, 2, 2, 
                 3, 0, 1, 2, 4, 0, 1, 4, 4, 1, 2, 1, 3, 0, 4, 4, 1, 3, 0, 0, 2, 
                 1, 0, 2, 0) + Sys.time(),
    x = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 
          20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 32, 33, 34, 35, 36, 37, 38, 
          39, 40, 41, 42, 43, 44, 47, 48, 49), 
    y = c(50, 49, 48, 47, 46, 47, 46, 45, 44, 43, 44, 43, 42, 41, 40, 41, 40, 39,
          38, 37, 38, 37, 36, 35, 34, 35, 34, 33, 32, 31, 30, 29, 28, 29, 28, 27, 
          26, 25, 26, 25, 24, 23, 22, 21, 20,19)
  ), 
  .Names = c("x", "y"), 
  class = c("tbl_df", "tbl", "data.frame"), 
  row.names = c(NA, -46L)
)

【讨论】：